text:
\begin{document}
\maketitle
\begin{abstract}
  We propose a new method for the numerical computation of the cut locus of a compact submanifold of $\R^3$ without boundary. This method is based on a convex variational problem with conic constraints, with proven convergence. We illustrate the versatility of our approach by the approximation of Voronoi cells on embedded surfaces of $\R^3$.
\end{abstract}
\section{Introduction}

Let $S$ be a real analytic surface without boundary embedded in $\R^3$, and $b\in S$ any point of $S$ (that can be thought of as a base point).
\begin{definition}
  The \emph{cut locus} of $b$ in $S$ can be defined as the closure of the set of points $p\in S$ such that there exists at least two minimizing geodesics between $p$ and $b$. We will denote it by $Cut_b(S)$. Equivalently, it is also the set of points around which the distance function to the point $b$ - denoted by $d_b$ - is not smooth.
\end{definition}
\noindent The cut locus is a fundamental object in Riemannian geometry, and it is a natural problem to try and find ways to compute it numerically. In this paper, we propose a numerical approximation of $Cut_b(S)$, based on a convex variational problem on $S$, with proven convergence. It is not trivial to compute $Cut_b(S)$ because it is not stable with respect to $C^1$-small variations of $S$.
See for instance \cite[Example 2]{albano2016stability}. For instance, one can't approximate the cut locus of $S$ with the cut locus of a piecewise linear approximation of $S$.

\textbf{Related works.} Let us review the techniques used in the past by different authors to approximate the cut locus. We may divide them in two categories.

\emph{Geodesic approximation on parametrized surfaces.} This approach was used in \cite{sinclair2002loki} and \cite{misztal2011cut}. In \cite{sinclair2002loki}, on genus $1$ parametrized surfaces, the authors computed a degree $4$ polynomial approximation of the exponential map using the geodesic equation, and deduced an approximation of the cut locus from there. In \cite{misztal2011cut}, the authors used the deformable simplicial complexes (DSC) method and finite differences techniques for geodesic computations, to compute geodesic circles of increasing radius and their self intersection, \textit{i.e.} the cut locus.
They apply the method to genus $1$ surfaces. These papers contain no proof of convergence of the computed cut locus.

\emph{Exact geodesic computation on discretized surfaces.} This approach was used in \cite{itoh2004thaw} and \cite{dey2009cut}. In \cite{itoh2004thaw}, the authors computed the geodesics on a convex triangulated surface. They deduced an approximation of the cut locus of the triangulated surface, and filtered it according to the angle formed by the geodesics meeting at a point of the approximated cut locus, to make their approximation stable. They applied the method to ellipsoids. There is no proof of convergence. In \cite{dey2009cut}, the authors computed shortest curves on a graph obtained from a sufficiently dense sample of points of a surface. From there they deduced an approximation of the cut locus, and filtered it according to the maximal distance (called \emph{spread}) between the geodesics meeting at a point of the approximated cut locus. They proved that the set they compute converges to the cut locus (see \cite[Theorem 4.1]{dey2009cut}).

We may also mention \cite{bonnard2014geometric_cutlocus}, where the authors use some more geometric tools to compute (numerically) the cut locus of an ellipsoid, or a sphere with some particular metric with singularities.

\textbf{Our method}. The strategy we use is quite different. Given $m>0$ a constant, let $u_m$ be the minimizer of the following variational problem
\begin{equation}\label{eq:gradient constraint}
  \min_{\substack{u\in H^1(S)\\ \abs{\nabla_{_S} u} \leq 1 \\u(b)=0}} \int_{S} \left( \abs{\nabla_{_S} u}^2 - mu\right),
\end{equation}
where $\nabla_{_S}$ denotes the gradient operator on the surface $S$.
For $\lambda>0$ to be chosen small, we will use the set $E_{m,\lambda} := \left\{ \abs{\nabla_{_S} u_m}^2 \leq 1 - \frac{\lambda^2}{u_m^2}\right\}$ as an approximation of $Cut_b(S)$.
This is justified by some theoretical results regarding problem \eqref{eq:gradient constraint} obtained in \cite{generauCutLocusCompact2020}, which will be summarized in section \ref{sec:variational problem}. Now the set $E_{m,\lambda}$ can be well approximated using finite elements on a triangulation of the surface $S$.

The rest of the paper is organized as follows. In section \ref{sec:lambda medial axis}, we recall the notion of $\lambda$-medial axis that was introduced in \cite{chazal_lieutier_lambda}, and summarize some of its properties.
In section \ref{sec:lambda cut locus}, following the strategy of the $\lambda$-medial axis, we define a "$\lambda$-cut locus" $Cut_b(S)_{\lambda}$ and show that it can be used as an approximation of the complete cut locus.
In section \ref{sec:variational problem}, we recall the result from \cite{generauCutLocusCompact2020} which states that the set $E_{m,\lambda}$ defined above is a good approximation of $Cut_b(S)_{\lambda}$ if $m$ is big enough.
In section \ref{sec:discretization}, we discretize problem \eqref{eq:gradient constraint} using finite elements, to find a discrete minimizer $v_h$, where $h>0$ is the step of the dicretization. We show that the set
\begin{equation*}
  E_{m,\lambda,h} := \left\{x\in S \setminus\{b\} : \abs{\nabla_{_S} u_{m,h}^l(x)}^2 \leq 1 - \frac{\lambda^2}{(u_{m,h}^l)^2(x)}\right\},
\end{equation*}
is a good approximation of $E_{m,\lambda}$ as $h\to 0$.
In section \ref{sec:numerical}, we present the results of some numerical experiments.

\section{$\lambda$-Medial axis}\label{sec:lambda medial axis}
In this section, we recall briefly the notion of $\lambda$-medial axis introduced by Chazal and Lieutier in \cite{chazal_lieutier_lambda}.
Given an open subset $\Omega$ of $\R^2$, its medial axis $\mathcal{M}(\Omega)$ is defined as the set of points of $\Omega$ that have at least two closest points on the boundary $\partial \Omega$ of $\Omega$:
\[ \mathcal{M}(\Omega) := \big\{ x \in \Omega : \exists y,z \in \partial \Omega, \; y \neq z \; \text{and} \; \dom(x) = \abs{x-y} = \abs{x-z} \big\}, \]
where for any $x\in\Omega$, $\dom(x)$ is the distance from $x$ to the boundary $\partial\Omega$,
\[d_{\partial\Omega}(x)=\min\big\{|x-y|\ :\ y\in\partial\Omega\big\}.\]
The medial axis $\mathcal{M}(\Omega)$ is unstable with respect to small non-smooth perturbations of the boundary of $\Omega$. To deal with this issue, Chazal and Lieutier defined the so called {\it$\lambda$-medial axis} of $\Omega$ by setting, for any $\lambda>0$,
\begin{equation}
  \mathcal{M}_\lambda (\Omega):= \{x\in\Omega : r(x)\geq\lambda\},
\end{equation}
where $r(x)$ is the radius of the smallest ball containing the set of all closest points to $x$ on $\partial \Omega$, \textit{i.e.} the set $\{z\in \partial \Omega : \abs{x-z} = \dom(x) \}$. The map $\lambda \mapsto \mathcal{M}_\lambda (\Omega)$ is non increasing, and
\begin{equation*}
  \mathcal{M}(\Omega) = \bigcup_{\lambda>0}\mathcal{M}_\lambda (\Omega).
\end{equation*}
It is further proved in \cite[section 3, theorem 2]{chazal_lieutier_lambda} that $\calM_\lambda (\Omega)$ has the same homotopy type as $\calM (\Omega)$, for $\lambda$ small enough. These facts justify that $\calM_\lambda (\Omega)$ is a good approximation of $\calM(\Omega)$, for $\lambda$ small enough.
The crucial  difference though is that $\calM_\lambda (\Omega)$ is stable with respect to small variations of $\Omega$, whereas $\calM(\Omega)$ is not. We refer the reader to \cite[section 4]{chazal_lieutier_lambda} for precise statements and proofs.

\section{$\lambda$-Cut locus}\label{sec:lambda cut locus}
We want to define a set similar to the $\lambda$-medial axis, in the case of the cut locus $Cut_b(S)$. To this end, we note that, as it can be seen from \cite[section 2.1]{chazal_lieutier_lambda}, we have
\begin{equation}\label{eq:equivalentdef}
  \mathcal{M}_\lambda(\Omega) = \left\{x\in \Omega : \abs{\nabla \dom(x)}^2 \leq 1 - \frac{\lambda^2}{\dom^2(x)}\right\},
\end{equation}
where $\nabla \dom$ denotes the generalized gradient wherever $\dom$ is not differentiable. Analogously, for $\lambda>0$, we define the $\lambda$-cut locus as
\begin{equation*}
  Cut_b(S)_\lambda := \left\{x\in S \setminus \{b\} : \abs{\nabla_{_S} d_b(x)}^2 \leq 1 - \frac{\lambda^2}{d_b^2(x)}\right\}.
\end{equation*}
Note that, according to \cite[Proposition 3.4]{mantegazza_mennucci_2003}, the function $d_b$ is locally semiconcave on $S\setminus \{b\}$, so it has a generalized gradient everywhere on $S\setminus \{b\}$, whose norm is given by the following formula:
\begin{equation}
  \abs{\nabla_{_S} d_b}(x)=\max(0, \sup\limits_{v\in T_xS, \abs{v}=1} \partial^+_v d_b(x)).
\end{equation}
We have the following proposition from \cite[Proposition 2.9]{generauCutLocusCompact2020}.
\begin{proposition}\label{prop:increasing union}
  The map $\lambda \mapsto Cut_b(S)_\lambda$ is non increasing, and
  \begin{equation*}
    Cut_b(S) = \overline{\bigcup_{\lambda>0}Cut_b(S)_\lambda}.
  \end{equation*}
\end{proposition}

\noindent In addition, the following proposition holds.
\begin{proposition}\label{prop:homotopy type}
  If $S$ is a real analytic surface, then for $\lambda>0$ small enough, one of the connected component of $Cut_b(S)_\lambda$ has the same homotopy type as $Cut_b(S)$, while the other connected components, if any, are contractible.
\end{proposition}

These two propositions justify that $Cut_b(S)_\lambda$ is a good approximation of $Cut_b(S)$, for $\lambda>0$ small enough. Before proving proposition \ref{prop:homotopy type}, we prove the following lemma.
\begin{lemma}\label{lemma:gradient angle}
  Let $x\in Cut_b(S)$ be such that there exists two unit speed minimizing geodesics $\gamma_1, \gamma_2 : [0,d_b(x)] \to S$ such that $\gamma_i (0) = b$ and $\gamma_i (d_b(x)) = x$. Let $\theta \in(0,\pi)$ be the angle between $\gamma_1$ and $\gamma_2$ at $x$. Then, we have
  \begin{equation*}
    \abs{\nabla_{_S} d_b}(x)\leq \cos(\theta/2).
  \end{equation*}
\end{lemma}
\begin{proof}
  For $i = 1,2$, let us set $v_i = -\dot{\gamma_i}(d_b(x))$. Let $v \in T_xS$. Let us denote by $\exp_x$ the Riemannian exponential map at the point $x$. Let $t_0\in(0,d_b(x))$ and $x_i = \exp_{v_i t_0}$. Note that we have $x\notin Cut_{x_i}(S)$, so the function $d_{x_i}$ is smooth at $x$, and its gradient is $-v_i$.
  Given $v\in T_xS$ such that $\abs{v} = 1$, we have
  \begin{align}
    \partial^+_v d_b(x)
    &= \lim\limits_{t\to 0^+} \frac{d_b(\exp_x(vt))-d_b(x)}{t} \nonumber
    \\
    &\leq \lim\limits_{t\to 0^+} \frac{d_{x_i}(\exp_x(vt)) + d_b(x_i) - (d_{x_i}(x) + d_b(x_i))}{t} \nonumber
    \\
    & = \lim\limits_{t\to 0^+} \frac{d_{x_i}(\exp_x(vt)) - d_{x_i}(x)}{t} \nonumber
    \\
    &= -v \cdot v_i \nonumber
  \end{align}
  Given that the angle between $v_1$ and $v_2$ is $\theta$, there exists $i\in\{1,2\}$, such that the angle between $v$ and $v_i$ is at most $\pi - \theta/2$. Thus the last inequality gives $\partial^+_v d_b(x) \leq \cos(\theta/2)$. This concludes the proof.
\end{proof}

Using lemma \ref{lemma:gradient angle}, Proposition \ref{prop:homotopy type} will mainly be a consequence of \cite[Proposition 3.4]{dey2009cut} and the proof of \cite[Proposition 3.5]{dey2009cut}. Following \cite{dey2009cut}, we will use the following terminology.
A point $x$ of a finite graph $G$ is called a tree point if $G\setminus\{x\}$ has a connected component whose closure is a tree. Otherwise, $x$ is called a cycle point. It is a consequence of the proof of \cite[Proposition 3.5]{dey2009cut} that any closed connected subset of $G$ that contains all cycle points is a deformation retract of $G$.
\begin{proof}[Proof of proposition \ref{prop:homotopy type}]
  As $S$ is real analytic, the cut locus $Cut_b(S)$ is a finite graph (see \cite{myers1936connections} in dimension $2$, and \cite{buchner1977simplicial} for the generalization to arbitrary dimensions).
  According to the lemma \ref{lemma:gradient angle}, given any $\theta>0$, if $\lambda$ has been taken small enough, then for any point $x\in Cut_b(S)\setminus Cut_b(S)_\lambda$, the angle between the minimizing geodesics from $b$ to $x$ is smaller than $\theta$.
  Given two unit speed minimizing geodesics $\gamma_1$ and $\gamma_2$, following \cite{dey2009cut}, the \emph{spread} between $\gamma_1$ and $\gamma_2$ is defined as
  \begin{equation*}
    spd(\gamma_1,\gamma_2) = \sup\limits_t d(\gamma_1(t),\gamma_2(t)).
  \end{equation*}
  As geodesics verify a second order differential equation, if their angle at their common starting point is small, then their spread is also small. Therefore, applying \cite[Proposition 3.4]{dey2009cut}, we deduce that if $\lambda$ has been taken small enough, then any point $x\in Cut_b(S)\setminus Cut_b(S)_\lambda$ is a \emph{tree point} of $Cut_b(S)$. It remains to show that $Cut_b(S)_\lambda$ is closed to conclude that it is a deformation retract of $Cut_b(S)$ and conclude the proof. But this is a consequence of the fact that $d_b$ is semiconcave, and the upper semicontinuity of the generalized gradient of convex functions.
\end{proof}
Therefore, we will use $Cut_b(S)_\lambda$ as an approximation of $Cut_b(S)$ for $\lambda$ small enough.

\section{Approximation with a variational problem}\label{sec:variational problem}
For $m>0$, recall that $u_m$ is the minimizer in \eqref{eq:gradient constraint}. For $\lambda>0$, let us define the set $E_{m,\lambda}$ by
\begin{equation*}
  E_{m,\lambda} := \left\{x\in S \setminus\{b\} : \abs{\nabla_{_S} u_m(x)}^2 \leq 1 - \frac{\lambda^2}{u_m^2(x)}\right\}.
\end{equation*}
We have the following theorem (see \cite[Theorem 1.1 and Theorem 1.3]{generauCutLocusCompact2020}):
\begin{theorem}\label{thm:big theorem}
  For any $m>0$, the function $u_m$ is locally $C^{1,1}$ on $S\setminus \{b\}$. For any $m>m'>m_0$,
  \begin{equation}\label{eq:double inclusion}
    Cut_b(S) \subset \{\abs{\nabla_{_S} u_m} < 1 \} \subset \{\abs{\nabla_{_S} u_{m'}} < 1 \}.
  \end{equation}
  Moreover,
  \begin{equation}\label{eq:convergence elastic set}
     \{\abs{\nabla_{_S} u_m} < 1 \}\tend{m}{+\infty}Cut_b(S) \quad \text{in the Hausdorff sense.}
  \end{equation}
  Finally, for any $\varepsilon>0$,
  \begin{equation}
    \sup\limits_{x \in E_{m,\lambda}} d(x,Cut_b(S)_{\lambda}) \tend{m}{+\infty}0, \label{eq:convergence lambda cut locus}
    \quad \text{and} \quad  \sup\limits_{x \in Cut_b(S)_{\lambda+\varepsilon}} d(x,E_{m,\lambda}) \tend{m}{+\infty}0.
  \end{equation}
\end{theorem}

Therefore, we can use $E_{m,\lambda}$ as an approximation of $Cut_b(S)_\lambda$. All in all, we will use $E_{m,\lambda}$ as an approximation of $Cut_b(S)$.

\section{Discretization}\label{sec:discretization}

\subsection{Finite elements of order $r$ on a surface approximation of order $k$}

In this section we introduce a discretization framework adapted to variational problem
\eqref{eq:gradient constraint} based on finite elements. We follow the notations
 of \cite{demlow2009higher, dziuk2013finite}.

Let $S$ be a compact oriented smooth two-dimensional surface embedded in $\R^3$.
For $x \in S$, we denote by $\nu (x)$ the oriented normal vector field on $S$.
Let $d:\R^3 \rightarrow \R$ be the signed distance associated to $S$ and
$U_\eta = \{x\in \R^3, \, |d(x)| < \eta \} $ the tubular
neighborhood of $S$ of width $\eta >0$. It is well known that if $\eta$ is
small enough (for  instance
$0< \eta < \min_{i = 1,2} {\frac{1}{|\kappa_i|}}_{L_{\infty}(S)}$ where the $(\kappa_i)$
stand for the extremal sectional curvatures of S), then for every $x \in U_\eta$
it exists a unique $a(x) \in S$ such that

\begin{equation}
  \label{eq:proj}
  x = a(x) + d(x) \nu(a(x)) = a(x) + d(x) \nabla d(x).
\end{equation}
We consider $S_h^1$ a triangular approximation of $S$ whose vertices lie on $S$
and whose faces are quasi-uniform and shape regular of  diameter at most $h>0$.
Moreover, we will assume that $\calTT_h$, the set of triangular faces of $S_h$,
are contained in some tubular neighborhood $U_\eta$ such that  the map $a$
defined by  \eqref{eq:proj} is unique.

For $k \geq 2$ and for a triangle $T \in \calTT_h$, we consider
the $n_k$ Lagrange basis functions $\Phi_1^k,\dots \Phi_{n_k}^k $ of degree $k$
and define the discrete projection on $S_h$ by:

\begin{equation}
  \label{eq:projk}
  a_k(x) = \sum_{j=1}^{n_k} a(x_j)  \Phi_{j}^k(x)
\end{equation}
where $x_1,\dots,x_{n_k}$ are the nodal points associated to the basis functions.
Now we can define $S_h^k$ a polynomial approximation of order $k$ of $S$ associated
 to $\calTT_h$

\begin{equation}
  \label{eq:Spoly}
  S_h^k = \{ a_k(x),\, x\in S_h \}.
\end{equation}
Observe that by definition the image by $a$ of the nodal points are both on $S$
and on  $S_h^k$. Let us now introduce the finite element spaces on $S_h = S_h^1$
and $S_h^k$ for $k \geq 2$. For every integer $r\geq 1$, let

\begin{equation}

  L^r_h  = \{ \chi  \in \calC^0(S_h),\, \chi|_{T} \in \PP_r,
  \forall  T \in \calTT_h\}
\end{equation}
where $\PP_r$ is the family of polynomials of degree at most $r$. Analogously,
for $k \geq 2$ let

\begin{equation}

  L^{r,k}_h  = \{ \hat{\chi}  \in \calC^0(S_h^k),\, \hat{\chi} =
  \chi \circ  a_k^{-1}, \, \text{for some } \,  \chi \in L^r_h \}.
\end{equation}
Analogously to \eqref{eq:gradient constraint}, we define
\begin{equation} \label{eq:gradient constraint_h}
  \min_{\substack{u\in L_h^{r,k} \\ \abs{\nabla_{_{S_h^k}} u} \leq 1 \\u(b)=0}}
                   F_h^k(u)
\end{equation}
where $ F_h^k(u) = \int_{S_h^k} \left( \abs{\nabla_{_{S_h^k}} u}^2 - mu\right)$ and $b$
 some fixed nodal points of the mesh $\calTT_h$.

\subsection{Convergence of the lifted minimizers}
In order to prove the convergence of our numerical approach, let us first establish
that our discrete problem converges in values in the sense of proposition \ref{prop:values}.
For a function $u$ defined on $S_h^k$, we introduce its lifted  function $u^l$
onto $S$ defined by the relation $u^l(b) = u(x)$ for $b \in S$ where $x$
is the unique point  of $S_h^k$ which satisfies $a(x) = b$.

Below, we focus our analysis in the piecewise linear case $r = k = 1$ which contains
 all the  main ingredients of a proof for the general $(r,k)$ case.
 For every $h>0$, the convex  optimization problem \eqref{eq:gradient constraint_h}
 has a unique solution.
\begin{lemma}\label{lemma:differential of a}
  The differential of the projection $a$ onto $S$, when restricted to the tangent space of $S_h$, is the identity, up to order $2$ in $h$:
  \begin{equation*}
    Da_{|_{TS_h}} = Id + \mathcal{O}(h^2).
  \end{equation*}
  The second differential of $a$, when restricted to the tangent space of $S_h$, is null, up to order $1$ in $h$:
  \begin{equation*}
    D^2a_{|_{TS_h}} =  \mathcal{O}(h).
  \end{equation*}
\end{lemma}
\begin{proof}
  The identity estimate on $Da$ is a direct consequence of \cite[equations (4.12), (4.13) and (4.11)]{dziuk2013finite}, and the fact that, following the notations of \cite[lemma 4.1]{dziuk2013finite}, we have $\nu_{n+1}^2 = 1 - \sum\limits_{j\leq n }\nu_j^2$.
  The estimate on $D^2a$ follows from the same equations, plus the identity $D^2a(x) = -2\nabla d(x)D^2d(x)$.

\end{proof}

  Defining $ F(u) = \int_{S} \left( \abs{\nabla_{_S} u}^2 - mu\right)$, we have

\begin{proposition}
  \label{prop:values}
  Let $u_{m,h}$ be the
  solution of problem \eqref{eq:gradient constraint_h} for $k=r=1$. Let $Lu_{m,h}^l := \frac{u_{m,h}^l}{|\nabla_{_{S}} u_{m,h}^l|_{L_\infty(S)}}$ be the $1$-Lipschitz normalization of $u_{m,h}^l$. Then,  $Lu_{m,h}^l \in H^1(S)$ and
  $$F(Lu_{m,h}^l)  = \min_{\substack{u\in H^1(S)\\ \abs{\nabla_{_{S}} u} \leq 1 \\u(b)=0}}
  F(u)  + \mathcal{O}(h^{\frac{1}{2}}).$$
\end{proposition}

\begin{proof}
    \emph{step 1.} Let $u_m$ be the solution of problem \eqref{eq:gradient constraint}. For $\varepsilon>0$, let $u_{m,\varepsilon}:S\to\R$ be defined by:
    \begin{equation*}
      u_{m,\varepsilon} =
      \begin{cases}
        \frac{d_b(x)^2}{2\varepsilon} \quad &\text{if} \quad d_b(x) \leq \varepsilon
        \\
        u_m(x) - \frac{\varepsilon}{2} \quad &\text{if} \quad d_b(x) \geq \varepsilon.
      \end{cases}
    \end{equation*}
    According to \cite[Lemma 3.3]{generauCutLocusCompact2020}, we have $u_m = d_b$ in a neighborhood of $b$. Therefore, for $\varepsilon>0$ small enough, we have $u_m = d_b$ on $B(b,2\varepsilon)$. In particular, we deduce that $u_{m,\varepsilon}$ is $C^1$ on $S$.
    As $d_b^2$ is smooth in a neighborhood of $b$, the gradient of $d_b^2/2\varepsilon$ is $\mathcal{O}(\varepsilon^{-1})$-Lipschitz on $B(b,\varepsilon)$. Moreover, as $u_m = d_b$ on $B(b,2\varepsilon)$, the gradient of $u_m$ is $\mathcal{O}(\varepsilon^{-1})$-Lipschitz on $B(b,2\varepsilon)\setminus B(b,\varepsilon)$.
    According to lemma \cite[Proposition 3.4]{generauCutLocusCompact2020}, $u_m$ is also locally $C^{1,1}$ on $S\setminus \{b\}$. Therefore its gradient is $\mathcal{O}(\varepsilon^{-1})$-Lipschitz on $S\setminus B(b,\varepsilon)$.
    All in all, we obtain that $u_{m,\varepsilon}$ is $C^{1,1}$ on $S$, and the Lipschitz constant of its gradient is $\mathcal{O}(\varepsilon^{-1})$. Furthermore, as $d_b$ and $u_m$ are both $1$-Lipschitz, we have $\abs{\nabla u_{m,\varepsilon}}\leq1$. Now for $\varepsilon>0$, consider
    \begin{equation*}
      v_{h,\varepsilon} := \frac{I_hu_{m,\varepsilon}}{|\nabla_{_{S_h}} I_hu_{m,\varepsilon}|_{L_\infty(S_h)}},
    \end{equation*}
    where $I_hu_{m,\varepsilon}$ is the $\PP^1$ Lagrange interpolation of $u_{m,\varepsilon}$ on $S_h$. For $x\in S_h$,
    observe that we have the relation $I_h u_{m,\varepsilon}(x) = I_h (u_{m,\varepsilon} \circ a)(x)$ which
    says that $I_h u_{m,\varepsilon}$ is the standard (flat) interpolation of the composed function
    $u_{m,\varepsilon} \circ a$. From lemma \ref{lemma:differential of a}, we know that on every triangle, the differential of $a$ is $\mathcal{O}(h)$-Lipschitz, and $a$ is $\mathcal{O}(1)$-Lipschitz.
    As the gradient of $u_{m,\varepsilon}$ is $\mathcal{O}(\varepsilon^{-1})$-Lipschitz, we deduce that on every triangle, the gradient of $u_m\circ a$ is $\mathcal{O}(\varepsilon^{-1})$-Lipschitz.
    By the quasi uniformity of the mesh, we obtain the uniform interpolation estimates on $S_{h}$:
    \begin{equation}\label{eq:interpolation estimate}
      I_h u_{m,\varepsilon}(x) = (u_{m,\varepsilon} \circ a)(x) + \mathcal{O}(\varepsilon^{-1}h^2)
    \end{equation}
    and
    \begin{equation*}
      \nabla_{_{S_h}} I_h u_{m,\varepsilon}(x) = \nabla_{_{S_h}} (u_{m,\varepsilon} \circ a)(x) + \mathcal{O}(\varepsilon^{-1}h).
    \end{equation*}

    With lemma \ref{lemma:differential of a}, we deduce for all $x\in S_{h}$,
    \begin{equation}\label{eq:interpolation gradient estimate}
      \nabla_{_{S_h}}  I_h u_{m,\varepsilon}(x) = \nabla_{_{S}} u_{m,\varepsilon}(a(x)) + \mathcal{O}(\varepsilon^{-1}h).
    \end{equation}
    Recall that we have $|\nabla_{_{S}} u_{m,\varepsilon}|_{L_\infty(S)} = 1$. Therefore the last identity yields
    $$|\nabla_{_{S_h}}  I_h u_{m,\varepsilon}|_{L_\infty(S_{h,\varepsilon})} = 1 +  \mathcal{O}(\varepsilon^{-1}h).$$
    Thus, $v_{h,\varepsilon} = I_h u_{m,\varepsilon}(1 + \mathcal{O}(\varepsilon^{-1}h))$, and so
    \begin{equation}\label{eq:F_h intermediate 1}
      F_h(v_{h,\varepsilon}) = F_h(I_h u_{m,\varepsilon}) +  \mathcal{O}(\varepsilon^{-1}h).
    \end{equation}
    Applying lemma \ref{lemma:differential of a} again, with a simple change of variable, we find that for any function $f:S_h \to \R$,
    \begin{equation*}
      \int_{S_h}f \circ a = \int_{S} f + \mathcal{O}(h^2).
    \end{equation*}
    Recalling \eqref{eq:interpolation estimate} and \eqref{eq:interpolation gradient estimate}, we obtain
    \begin{equation}\label{eq:F_h intermediate 2}
      F_h(I_h u_{m,\varepsilon}) = F(u_{m,\varepsilon}) +  \mathcal{O}(\varepsilon^{-1}h).
    \end{equation}
    Furthermore, we have
    \begin{equation*}
      \int_S\abs{u_{m,\varepsilon}-u_m} \leq \mathcal{O}(\varepsilon)\quad \text{and} \int_S \abs{\nabla u_{m,\varepsilon} - \nabla u_m}^2 \leq \mathcal{O}(\varepsilon^2),
    \end{equation*}
    so
    \begin{equation*}
      F(u_{m,\varepsilon}) = F(u_m) + \mathcal{O}(\varepsilon).
    \end{equation*}
    Combining this with \eqref{eq:F_h intermediate 1} and \eqref{eq:F_h intermediate 2}, we find
    \begin{equation*}
      F_h(v_{h,\varepsilon}) = F(u_m) + \mathcal{O}(\varepsilon^{-1}h) + \mathcal{O}(\varepsilon).
    \end{equation*}
    Choosing $\varepsilon = h^{\frac{1}{2}}$, this yields
    \begin{equation}\label{eq:min comparison}
      \min_{\substack{u\in H^1(S_h)\\ \abs{\nabla_{_{S_h}} u} \leq 1 \\u(b)=0}} F_h \leq \min_{\substack{u\in H^1(S)\\ \abs{\nabla_{_{S}} u} \leq 1 \\u(b)=0}} F +  \mathcal{O}(h^{\frac{1}{2}}).
    \end{equation}

    \emph{step 2.} Symmetrically, let $u_{m,h}$ the solution of the discrete
    problem  \eqref{eq:gradient constraint_h}, $u^l_h := u_{m,h} \circ (a_{|_{S_h}})^{-1}$ its lifted version on $S$, and $Lu_{m,h}^l :=  \frac{u^l_h}{{|\nabla_{_{S_h}} u^l_h|}_{L_\infty (S_h)}} $. We show as before, using the equation $u_{m,h} = u^l_h \circ a$, that $F(Lu_{m,h}^l) = F_h(u_{m,h}) +  \mathcal{O}(h)$.
    With \eqref{eq:min comparison}, this implies
    \begin{equation*}
      \min_{\substack{u\in H^1(S)\\ \abs{\nabla_{_{S}} u} \leq 1 \\u(b)=0}} F \leq F(Lu_{m,h}^l) \leq \min_{\substack{u\in H^1(S)\\ \abs{\nabla_{_{S}} u} \leq 1 \\u(b)=0}} F +  \mathcal{O}(h^{\frac{1}{2}}),
    \end{equation*}
    which concludes the proof of the proposition.
\end{proof}

We can now establish the convergence of the minimizers:
\begin{proposition}\label{prop:convergence minimizers}
  \begin{equation*}
    \abs{\nabla u_{m,h}^l - \nabla u_m}^2_{L^2(S)} = \mathcal{O}(h^{\frac{1}{2}}) \quad \text{and} \quad \abs{u_{m,h}^l - u_m}_{L^1(S)} = \mathcal{O}(h^{\frac{1}{2}}).
  \end{equation*}
\end{proposition}

\begin{proof} Consider $v = \frac12 (Lu_{m,h}^l + u_m)$. Then, $v$ is admissible for
   problem \eqref{eq:gradient constraint}, so $F(v)\geq F(u_m)$.
    Moreover, the following algebraic identity holds
$$
F(v) = \frac12 F(Lu_{m,h}^l) + \frac12 F(u_m) - \frac14 \int_S |\nabla_{_{S}} u_m - \nabla_{_{S}}  Lu_{m,h}^l|^2.
$$
Therefore, we have
$$
\frac12 F(Lu_{m,h}^l) - \frac12 F(u_m) \geq \frac14 \int_S |\nabla_{_{S}} u_m - \nabla_{_{S}} Lu_{m,h}^l|^2,
$$
which proves, with proposition \ref{prop:values}, that
\begin{equation}\label{eq:l2 gradient}
  \abs{\nabla Lu_{m,h}^l - \nabla u_m}^2_{L^2(S)}  = \mathcal{O}(h^{\frac{1}{2}}).
\end{equation}
Moreover, we have
\begin{equation*}
   F(Lu_{m,h}^l) -  F(u_m) = \int_S \left(\abs{\nabla Lu_{m,h}^l}^2 - \abs{\nabla u_{m}}^2\right) - m \int_S \left(Lu_{m,h}^l - u_m\right).
\end{equation*}
The last two equations imply
\begin{equation}\label{eq:l1 function}
  \abs{ Lu_{m,h}^l -  u_m}_{L^1(S)}  = \mathcal{O}(h^{\frac{1}{2}}).
\end{equation}
As in the proof of proposition \ref{prop:values}, using the relation $u_{m,h} = u_{m,h}^l \circ h$, we show that $Lu_{m,h}^l = u_{m,h}^l(1 + \mathcal{O}(h^2))$. Together with \eqref{eq:l2 gradient} and \eqref{eq:l1 function}, this concludes the proof.
\end{proof}
We just proved that the sequence of the lifted minimizers converges with an order
 at least $1/2$ to the minimizer of problem  \eqref{eq:gradient constraint}. By analogy with the more standard variational context \cite{demlow2009higher, dziuk2013finite},
  we expect a convergence of order $\mathcal{O}((h^r + h^{k+1})^{\frac{1}{2}})$
 using an approximation of orders $(r,k)$.

\subsection{Convergence in measure to the elastic set}
Let us recall that the set $E_{m,\lambda}$ is defined by
\begin{equation*}
  E_{m,\lambda}=\left\{x\in S \setminus\{b\} : \abs{\nabla_{_S} u_m(x)}^2 \leq 1 - \frac{\lambda^2}{u_m^2(x)}\right\}.
\end{equation*}

\begin{proposition}
For any $\lambda >0$ and $\varepsilon>0$ with $\varepsilon<\lambda/2$, let us define
\begin{equation*}
  E_{m,\lambda,h} := \left\{x\in S \setminus\{b\} : \abs{\nabla_{_S} u_{m,h}^l(x)}^2 \leq 1 - \frac{\lambda^2}{(u_{m,h}^l)^2(x)}\right\}.
\end{equation*}
Then, we have
\begin{equation*}
  \abs{E_{m,\lambda+\varepsilon}\setminus E_{m,\lambda,h}} = \mathcal{O}(h^{\frac{1}{2}}) \quad \text{and} \quad \abs{E_{m,\lambda,h}\setminus E_{m,\lambda-\varepsilon}} = \mathcal{O}(h^{\frac{1}{2}}).
\end{equation*}
\end{proposition}

\begin{proof}
By definition of $E_{m,\lambda}$ and $E_{m,\lambda,h}$, we have
\begin{equation*}
  E_{m,\lambda+\varepsilon}\setminus E_{m,\lambda,h} \subset \left\{\abs{\nabla u_{m,h}^l}^2 - \abs{\nabla u_m}^2 > \frac{(\lambda+\varepsilon)^2}{u_m^2} - \frac{\lambda^2}{(u_{m,h}^l)^2} \right\}.
\end{equation*}
Therefore, on $E_{m,\lambda+\varepsilon}\setminus E_{m,\lambda,h}$, we have
\begin{align*}
  \abs{\nabla u_{m,h}^l}^2 - \abs{\nabla u_m}^2
  & > \frac{(\lambda+\varepsilon)^2-\lambda^2}{u_m^2} + \lambda^2 (\frac{1}{u_m^2} - \frac{1}{(u_{m,h}^l)^2})
  \\
  & \geq \frac{2 \varepsilon \lambda + \varepsilon^2}{u_m^2} - \lambda^2 \frac{2}{\min(u_m,u_{m,h}^l)^3}\abs{u_m - u_{m,h}^l}
  \\
  & = \frac{2 \varepsilon \lambda + \varepsilon^2}{(\diam S)^2} - \lambda^2 \frac{2}{(u_m + \mathcal{O}(h^{\frac{1}{2}}))^3}\abs{u_m - u_{m,h}^l},
\end{align*}
where $\diam S$ is the diameter of $S$.
By definition of $E_{m,\lambda}$, we also have $E_{m,\lambda+\varepsilon} \subset \{u_m\geq (\lambda+\varepsilon)\}$, so on $E_{m,\lambda+\varepsilon}\setminus E_{m,\lambda,h}$,
\begin{align*}
  \abs{\nabla u_{m,h}^l}^2 - \abs{\nabla u_m}^2
  > \frac{2 \varepsilon \lambda + \varepsilon^2}{(\diam S)^2} - \lambda^2 \frac{2}{(\lambda + \varepsilon + \mathcal{O}(h^{\frac{1}{2}}))^3}\abs{u_m - u_{m,h}^l}.
\end{align*}
So for $h$ big enough, we have
\begin{align*}
  &E_{m,\lambda+\varepsilon}\setminus E_{m,\lambda,h}
  \\
  &\quad \subset \left\{\abs{\nabla u_{m,h}^l}^2 - \abs{\nabla u_m}^2 + 2 \lambda \abs{u_m - u_{m,h}^l} > \frac{2 \varepsilon \lambda + \varepsilon^2}{(\diam S)^2} \right\}.
  \\
  &\quad \subset \left\{\abs{\nabla u_{m,h}^l}^2 - \abs{\nabla u_m}^2 > \frac{2 \varepsilon \lambda + \varepsilon^2}{2(\diam S)^2} \right\} \cup \left\{\abs{u_m - u_{m,h}^l} > \frac{2 \varepsilon \lambda + \varepsilon^2}{4 \lambda (\diam S)^2} \right\}.
\end{align*}
Now from proposition \ref{prop:convergence minimizers}, we know that for any $\eta>0$, we have the following estimates
\begin{equation*}
  \eta\abs{\left\{\abs{\abs{\nabla u_{m,h}^l}^2 - \abs{\nabla u_m}^2} > \eta \right\}}\leq \int_S \abs{\abs{\nabla u_{m,h}^l}^2 - \abs{\nabla u_m}^2} = \mathcal{O}(h^{\frac{1}{2}}),
\end{equation*}
and
\begin{equation*}
  \eta\abs{\left\{\abs{u_m - u_{m,h}^l} > \eta \right\}}\leq \int_S \abs{u_m - u_{m,h}^l} = \mathcal{O}(h^{\frac{1}{2}}).
\end{equation*}
This gives the estimate $\abs{E_{m,\lambda+\varepsilon}\setminus E_{m,\lambda,h}} = \mathcal{O}(h^{\frac{1}{2}})$. The other estimate is proved by the same method.
\end{proof}

\begin{remark}
  We expect a convergence of order $\mathcal{O}((h^r + h^{k+1})^{\frac{1}{2}})$
 using an approximation of orders $(r,k)$.
\end{remark}

Sections \ref{sec:lambda cut locus}, \ref{sec:variational problem} and \ref{sec:discretization} together justify that the set $E_{m,\lambda,h}$ is a good approximation of the cut locus of $b$ in $M$, if $m$ is big enough, and $\lambda$ and $h$ are small enough.

\section{Numerical illustrations}\label{sec:numerical}

\subsection{Cut locus approximation}

We established the convergence of the minimizers of solutions of problems
\eqref{eq:gradient constraint} when $h$ tends to $0$. For a fixed $h>0$,
this convex discrete problems is of quadratic type with an infinite
 number of conic pointwise constraints.
By the way, it is important to observe that  for $k=r=1$, the gradient pointwise
bounds for a function of $\PP^1$ is
equivalent to a single discrete conic constraint on every triangle with respect
to the degrees of freedom of  $\PP^1(\calTT_h)$.

Nevertheless, we observed in our experiments that using  $\PP^1$ elements may lead to
approximated cut loci with some tiny artificial connected components. Motivated by
 this lack of precision, we use in all following illustrations elements of order
  $r >1$.

 For the general case $r>1$, the bound constraint on the gradient can not be
 easily reduced to a finite set of discrete constraints. In  our experiments, we
approximated the constraint
$| \nabla_{_{S_h^k}} u|_{L^\infty (S_h^k)} \leq 1$  by forcing this constraint only on a
finite number of points of the mesh. In practice, we imposed these constraints
on the Gauss quadrature points of order $g$ on every triangle of $\calTT_h$.

We illustrate in figures \ref{fig:f1}, \ref{fig:f2}, \ref{fig:f3} and \ref{fig:f4}
the approximation of  the cut locus provided by our approach. These computations have been
carried out on meshes of approximated $10^5$ triangles for $k = 2$ and $r = 3$
 using high precision quadrature formula associated to $17$ Gauss points on every element
  of the mesh. Moreover, for $r=3$, we imposed the conic gradient constraints on
  the $g = 9$ Gauss points of every triangle.  In order to solve the resulting
  linear conic constrained quadratic optimization problem,
   we used the \emph{JuMP} modeling language and the finite elements library \emph{Getfem++  }
   \cite{JuMP.jl-2017, renard2006getfem++} combined with \emph{Mosek} optimization
   solver \cite{andersen2000mosek}. For such a precision,  the optimization
   solver identified a solution in less than one hour on a standard computer.

\begin{figure}[htbp]
  \centering
  \begin{tabular}{c c c}
    \includegraphics[width=\widthfig cm]{movie_1Sgenre1.png}
    & \includegraphics[width=\widthfig cm]{movie_120Sgenre1.png}
    & \includegraphics[width=\widthfig cm]{movie_180Sgenre1.png}
  \end{tabular}
  \caption{Three different views of the approximation of a cut locus on a standard torus}
  \label{fig:f1}
\end{figure}

\begin{figure}[htbp]

  \begin{tabular}{c c c}
    \includegraphics[width=\widthfig cm]{movie_1genre1.png}
    & \includegraphics[width=\widthfig cm]{movie_120genre1.png}
    & \includegraphics[width=\widthfig cm]{movie_180genre1.png}
  \end{tabular}
  \caption{Three different views of the approximation of a cut locus on a standard torus,
           without representing the surface}
  \label{fig:f2}
\end{figure}

\begin{figure}[htbp]
  \centering
  \begin{tabular}{ccc}
    \includegraphics[width=\widthfig cm]{movie_50Sgenre2.png}
    & \includegraphics[width=\widthfig cm]{movie_140Sgenre2.png}
    & \includegraphics[width=\widthfig cm]{movie_180Sgenre2.png} \\
  \end{tabular}
  \caption{Three different views of the approximation of a cut locus on a torus of genus $2$}
  \label{fig:f3}
\end{figure}

\begin{figure}[htbp]

  \begin{tabular}{ccc}
    \includegraphics[width=\widthfig cm]{movie_50genre2.png}
    & \includegraphics[width=\widthfig cm]{movie_140genre2.png}
    & \includegraphics[width=\widthfig cm]{movie_180genre2.png}
  \end{tabular}
  \caption{Three different views of the approximation of a cut locus on a torus of genus $2$,
           without representing the surface}
  \label{fig:f4}
\end{figure}
\FloatBarrier

\subsection{Approximation of the boundary of Voronoi cells}

All previous theoretical results still hold if we replace the source point $b$ by any
compact subset of the surface $S$. For instance, if $b$ is replaced by a set of points,
the singular set of the distance function can be decomposed as the union of the
boundary of voronoi cells and the cut loci of every point intersected with its
voronoi cell. As a consequence, if the distribution of source points is homogeneous enough, that is every voronoi cell is small enough, the singular part
of the distance function will be exactly equal to the boundary of the voronoi cells.
We illustrate this remark in the following experiments. We used exactly the
same framework as in previous sections and just replaced the pointwise condition
 at $b$ with the analogous pointwise Dirichlet conditions at every source point.
Figure \ref{fig:vortorus1} and \ref{fig:vortorus2} represent the voronoi
 diagrams obtained with $10$, $30$ and $100$ points for surfaces of genus $2$ and $3$.
  The expected computational complexity is exactly of the same order as with
   a single source point.

\begin{figure}[htbp]
  \begin{tabular}{ccc}
    \includegraphics[width=\widthvoronoi cm]{movie_50V1.png} &
    \includegraphics[width=\widthvoronoi cm]{movie_50V2.png} &
    \includegraphics[width=\widthvoronoi cm]{movie_50V3.png}\\
    \includegraphics[width=\widthvoronoi cm]{movie_150V1.png} &
    \includegraphics[width=\widthvoronoi cm]{movie_150V2.png} &
    \includegraphics[width=\widthvoronoi cm]{movie_150V3.png}
  \end{tabular}
  \caption{Approximation of the voronoi cells on a torus of genus $2$ of $10$, $30$ and $100$ points. Every column represent two different views}
  \label{fig:vortorus1}
\end{figure}

\begin{figure}[htbp]
  \begin{tabular}{ccc}
    \includegraphics[width=\widthvoronoi cm]{movie_50V4.png} &
    \includegraphics[width=\widthvoronoi cm]{movie_50V5.png} &
    \includegraphics[width=\widthvoronoi cm]{movie_50V6.png}\\
    \includegraphics[width=\widthvoronoi cm]{movie_150V4.png} &
    \includegraphics[width=\widthvoronoi cm]{movie_150V5.png} &
    \includegraphics[width=\widthvoronoi cm]{movie_150V6.png}
  \end{tabular}
  \caption{Approximation of the voronoi cells on a torus of genus $3$ of $10$, $30$
  and $100$ points. Every column represent two different views}
  \label{fig:vortorus2}
\end{figure}

\FloatBarrier

\bibliographystyle{plain}
\bibliography{../../biblio/biblio.bib}

\bigskip\noindent
François Générau:
Laboratoire Jean Kuntzmann (LJK),
Universit\'e Joseph Fourier\\
Bâtiment IMAG, 700 avenue centrale,
38041 Grenoble Cedex 9 - FRANCE\\
{\tt francois.generau@univ-grenoble-alpes.fr}

\bigskip\noindent
\'Edouard Oudet:\\
Laboratoire Jean Kuntzmann (LJK),
Universit\'e Grenoble Alpes\\
Bâtiment IMAG, 700 avenue centrale,
38041 Grenoble Cedex 9 - FRANCE\\
{\tt edouard.oudet@univ-grenoble-alpes.fr}

\bigskip\noindent
Bozhidar Velichkov:\\
Dipartimento di Matematica, Universit\`a di Pisa\\
Largo Bruno Pontecorvo, 5, 56127 Pisa - ITALY\\
{\tt bozhidar.velichkov@gmail.com}

\end{document}


TITLE: Binomial and Normal Distribution Problem - Check solution
QUESTION [0 upvotes]: Whooping cough is a highly contagious bacterial infection...About 80% of unvaccinated children who are exposed to whooping cough will develop the infection, as opposed to only about 5% of vaccinated children.
1) Find the probability that exactly 2 out of 20 exposed children develop whooping cough if 17 children have been vaccinated and 3 has not been.
My solution for this question is:
   17 children = 0.05 ^ 17 and 3 children = 0.8 ^ 3
2 out of 20 = (20C2)(0.05^17 + 0.8^3)^2 (1 - (0.05^17 + 0.8^3))^18
i don't know whether my ans is correct
2) 1400 unvaccinated children were exposed to whooping cough. What is the probability that at least 78% of these children develop infections?
My solution:
   this must be approx. to normal
P(Z >= 1092 - 1120 / square_root(224)) = p(Z < 1.87) = 0.9693
Please check my answers and let me know. Thanks in advance.

REPLY [0 votes]: Whooping cough is a highly contagious bacterial infection...About 80% of unvaccinated children who are exposed to whooping cough will develop the infection, as opposed to only about 5% of vaccinated children.
1) Find the probability that exactly 2 out of 20 exposed children develop whooping cough if 17 children have been vaccinated and 3 has not been.

Let $V$ count the number of vaccinated children who develop whooping cough, and $U$ count the number of unvaccinated children who do so.  These can be assumed to be independent and binomially distributed variables.
We have been given: $V\sim\mathcal{B}(17,0.05), U\sim\mathcal{B}(3,0.80), V\bot U$
We are tasked to find: $\mathrm{P}(V+U=2)$
We use the following:
$$V\sim\mathcal{B}(17,0.05)\iff \mathrm{P}(V\!=\!v)={17\choose v}(0.05)^v(0.95)^{17-v}$$
$$U\sim\mathcal{B}(3,0.80)\iff \mathrm{P}(U\!=\!u)={3\choose u}(0.80)^u(0.20)^{3-u}$$
$$\mathrm{P}(V\!+\!U\!=\!2) = \mathrm{P}(V\!=\!0)\mathrm{P}(U\!=\!2)+\mathrm{P}(V\!=\!1)\mathrm{P}(U\!=\!1)+\mathrm{P}(V\!=\!2)\mathrm{P}(U\!=\!0)$$
Put it together.


TITLE: Clarification over Ahlfors page 116, 2.1 about winding numbers
QUESTION [3 upvotes]: Everything on this question is in complex plane.
As the book describes a property of a winding number, it says that:

Outside of the [line segment from $a$ to $b$] the function $(z-a) / (z-b)$ is never real and $\leq 0$.

Here, the above statement should be interpreted as "never (real and $\leq 0$)".
If anyone could explain why this is true that would be great. I do get why any point on the line segment (other than $b$, in which case the denominator is $0$) has to satisfy the condition that $(z-a) / (z-b)$ is real and $\leq 0$, but I am not sure how to prove why any point not on the line has to satisfy the condition also.
Here, $a$ and $b$ are arbitrary complex number in a region determined by a closed curve in the complex plane; both points lie on the same region.

REPLY [1 votes]: The imaginary part of $u=\frac{z-a}{z-b}$ is $$\frac{u-\overline{u}}{2i}=\frac{1}{2i}\left(\frac{z-a}{z-b}-\overline{\frac{z-a}{z-b}}\right)=\frac{1}{2i}\left(\frac{z-a}{z-b}-\frac{\overline{z}-\overline{a}}{\overline{z}-\overline{b}}\right)$$
So, if we want this to be real, we need
$$
\begin{eqnarray*}
\frac{1}{2i}\left(\frac{z-a}{z-b}-\frac{\overline{z}-\overline{a}}{\overline{z}-\overline{b}}\right)&=&0
\\
(z-a)(\overline{z}-\overline{b})&=&(z-b)(\overline{z}-\overline{a}) \\
z\overline{z}-z\overline{b}-a\overline{z}+a\overline{b}&=&z\overline{z}-z\overline{a}-b\overline{z}+b\overline{a}
\\
0&=&z(\overline{a}-\overline{b})+b(\overline{z}-\overline{a})-a(\overline{z}-\overline{b})
\\
0&=&-\Re(b)\Im(z)+\Im(b)\Re(z)+\Re(a)(\Im(z)-\Im(b))+\Im(a)(\Re(b)-\Re(z)) \\
0&=&\Re(b)\Im{a}-\Re(a)\Im(b)-(\Re(a)-\Re(b))\Im(z)+(\Im(b)-\Im(a))\Re(z) \\
\Im(z)&=&\frac{\Im(b) \Re(a) - \Im(a) \Re(b) + (\Im(a) - \Im(b)) \Re(z)}{\Re(a) - \Re(b)}
\end{eqnarray*}
$$
which is the complex point-slope form of the line connecting $a$ and $b$.
So, what about whether the real part is positive or negative? We'll use this formula for $z$ to compute $u$.
$$
\begin{eqnarray*}
u&=&\frac{\Re(z)+i\Im(z)-\left(\Re(a)-i \Im(a)\right)}{\Re(z)+i\Im(z)-\left(\Re(b)-i \Im(b)\right)}\\
&=&
\frac{\Re(z)+i\left(\frac{\Im(b) \Re(a) - \Im(a) \Re(b) + (\Im(a) - \Im(b)) \Re(z)}{\Re(a) - \Re(b)}\right)-(\Re(a)+i \Im(a))}{\Re(z)+i\left(\frac{\Im(b) \Re(a) - \Im(a) \Re(b) + (\Im(a) - \Im(b)) \Re(z)}{\Re(a) - \Re(b)}\right)-(\Re(b)+i \Im(b))}\\
&=&
\frac{\left(\Re(z)-\Re(a)\right)\left(\frac{\Re(a)-\Re(b)+i(-\Im(a)+\Im(b))}{\Re(a)-\Re(b)}\right)}{\left(\Re(z)-\Re(b)\right)\left(\frac{\Re(a)-\Re(b)+i(-\Im(a)+\Im(b))}{\Re(a)-\Re(b)}\right)}\\
&=&
\frac{\Re(z)-\Re(a)}{\Re(z)-\Re(b)}
\end{eqnarray*}
 $$
From there, it's simple algebra to see that $u=\frac{\Re(z)-\Re(a)}{\Re(z)-\Re(b)}$ is negative only between $\Re(a)$ and $\Re(b)$, which corresponds to the line segment joining $a$ and $b$.


TITLE: For permutations with “copies”, in what sense are the elements of the set distinct?
QUESTION [2 upvotes]: Suppose I have a set of elements $$ S = \lbrace A,B,B,D,E \rbrace$$
Then the number of distinct permutations of size 5 is
$$\frac{5!}{2!} = 60$$
This is because $S$ contains a ‘copy’ of B.
But by definition, in a set of elements, each element is distinct. However, now I’m essentially saying “wait actually B and B aren’t distinct.”
So why isn’t the set $$ S = \lbrace A,B,D,E \rbrace$$ if the B elements are not distinct? Notice the permutations of this set are very different
$$4! = 24$$
Are there two different kinds of being distinct here?

REPLY [2 votes]: You should see the counts as the number of outputs of two different experiments:
First we have a bag of 5 letter tiles (a la Scrabble say) with $A,B,B,D,E$ tiles.
We draw 5 tiles from the bag and put them in a row in order we draw them. We have 5 draws, so $5!$ many permutations of the tiles. But in reality we cannot (in the end) distinguish the two B's (We could have marked them and then we'd have $5!$ different results): so the number of end results (sequences of 5 letters) is halved, because for every permutation $\ldots B_1\ldots B_2\ldots$ have a seocnd equivalent one $\ldots B_2\ldots B_1\ldots$.
The second experiment we have a bag with four tiles $A,B,D,E$ and counts the number of possible permutations from a draw of $4$ letters put in a row, so there we count sequences of $4$ letters, not of $5$. Hence the difference in counts.

REPLY [1 votes]: take a simpler case $\{A,B,B\}$ vs $\{A,B\}$.  First one will have $ABB,BBA,BAB$, however the second one only $AB,BA$.  Note that when $BB$ are next to each other it can be mapped to the single instance case one-to-one, but not when there are separated.


TITLE: A pattern appearing in the powers of $\phi$
QUESTION [13 upvotes]: \begin{align}
\phi^5 &= 11,\underline{0}901699\cdots\\
\phi^6 &= 17,\underline{9}44271\cdots\\
\phi^7 &= 29,\underline{6}34441\cdots\\
\phi^8 &= 46,\underline{9}7871\cdots\\
\phi^9 &= 76,\underline{0}1315 \cdots\\
\phi^{10} &= 122,\underline{99}18\cdots\\
\phi^{11} &= 199,\underline{00}502\cdots\\
\phi^{12} &= 321,\underline{99}6894\cdots\\
\phi^{13} &= 521,\underline{00}191\cdots\\
\phi^{14} &= 842,\underline{99}881\cdots\\
\phi^{15} &= 1364,\underline{000}73\cdots\\
\phi^{16} &= 2206,\underline{999}54\cdots\\
\end{align}
Why there is a $0$ $9$ patterns in the powers of the golden ratio

REPLY [11 votes]: $\phi$ is the larger root of $$x^2-x-1=0$$
It's conjugate root is:  $$\overline {\phi}=\frac {1-\sqrt 5}2\approx -0.61803$$
From the quadratic, we see that the sequence $$\{a_n\}=\{\phi^n +\bar {\phi}^n\}$$  satisfies the Fibonacci recursion:  $$a_{n+1}=a_n+a_{n-1}\quad a_1=1\quad a_2=3$$
Of course $\bar {\phi}^n\to 0$ for large $n$ so we must have that $\phi^n$ is nearly an integer for large $n$.  since $\bar {\phi}^n$ alternates sign we see the pattern you have noticed.

REPLY [1 votes]: I cannot provide an explanation for why the $0$ and $9$ patterns appear, but I can give an explanation of why it seems like the powers converge to integers.
If you visit the following website: http://www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/phigits.html
Based on the special properties of $\phi$ (I won't prove it here),
$$\phi^2 -\phi^{-2} = 3,\\\phi^3 -\phi^{-3} = 4,\\\phi^4 -\phi^{-4} = 7,\\\phi^5 -\phi^{-5} = 11\\\cdots$$
It should be evident that as the exponent increases, the second term in each expression will converge to zero, and the entire expression converges to a number.


TITLE: Is there a name for a binomial expansion without coefficients?
QUESTION [3 upvotes]: I am investigating a problem from George E. Andrews Number Theory (Dover, 1971), discussed previously here:
Induction Proof that $x^n-y^n=(x-y)(x^{n-1}+x^{n-2}y+\ldots+xy^{n-2}+y^{n-1})$
I was led astray for a bit, because I misunderstood the rhs to be the difference of $x$ and $y$ times the binomial expansion $(x+y)^{n-1}$.
But this is wrong... the rhs is actually the difference of $x$ and $y$ times the binomial expansion $(x+y)^{n-1}$ "without the coeffients", i.e.,
$$\sum_{i=0}^{n-1}x^{n-1-i}y^i$$
(also briefly discussed here: Binomial summation without coefficient)
Is there a name for this summation?

REPLY [2 votes]: There is no specific name for the bivariate polynomial
  \begin{align*}
x^{n-1}+x^{n-2}y+x^{n-3}y^2+\cdots+x^2y^{n-3}+xy^{n-2}+y^{n-1}
\end{align*}
but the nice expression
  \begin{align*}
(x-&y)(x^{n-1}+x^{n-2}y+x^{n-3}y^2+\cdots+x^2y^{n-3}+xy^{n-2}+y^{n-1})\\
\end{align*}
is called a telescoping sum since all terms besides the first and the last cancel out.
\begin{align*}
(x-&y)(x^{n-1}+x^{n-2}y+x^{n-3}y^2+\cdots+x^2y^{n-3}+xy^{n-2}+y^{n-1})\\
&=x(x^{n-1}+x^{n-2}y+x^{n-3}y^2+\cdots+x^2y^{n-3}+xy^{n-2}+y^{n-1})\\
&\qquad-y(x^{n-1}+x^{n-2}y+x^{n-3}y^2+\cdots+x^2y^{n-3}+xy^{n-2}+y^{n-1})\\
&=\color{blue}{x^{n}}+x^{n-1}y+x^{n-2}y^2+\cdots+x^3y^{n-3}+x^2y^{n-2}+xy^{n-1}\\
&\qquad\ \;-x^{n-1}y-x^{n-2}y^2-\cdots-x^3y^{n-3}-x^2y^{n-2}-xy^{n-1}\color{blue}{-y^{n}}\\
&=\color{blue}{x^n-y^n}\\
\end{align*}


TITLE: Perimeters of random-walk polygons
QUESTION [13 upvotes]: I have a random walk on $\mathbb{Z}^2$ that takes a step
with equal probability in the three directions that avoid
retracing the previous step.
The walk proceeds until it returns to a lattice point
previously visited, at which time it pinches off a simple,
closed loop or polygon:

I'd like to know the distribution of the perimeters of
these polygons.
(In the example above, the perimeter is 10.)
Simulations show that the average perimeter is about 5.6,
with perimeter 4 the overwhelming favorite,
as one would expect:

I feel this distribution must be known to the experts and
not difficult to explicitly detail, but
after looking at hitting times, first-passage times, self-avoidance times,
and various other frequently studied random walk quanities,
I am not finding a close-enough analog to help.
Thanks for any pointers you might provide!
Addendum.
Here is log-plot of the probability of a perimeter of length $L$, based on a simulation of $10^6$ walks. The first point represents 642,225 perimeters for $L=4$, the second point 176,043 perimeters for $L=6$, etc. The last point plotted is 135 instances of $L=38$. (There is one polygon of length $L=74$ in these million trials.) The average perimeter length is 5.62, which occurs after an average of 8.46 steps.

REPLY [3 votes]: I believe the paper below is on this exact question:
http://pre.aps.org/abstract/PRE/v55/i3/pR2093_1


TITLE: Prove that $f:M(m \times n,F) \to Hom(F^n,F^m)$ is an isomorphism of vector spaces
QUESTION [1 upvotes]: So, I'm trying to prove the following result:
The map $f:M(m \times n,F) \to Hom(F^n,F^m)$, which associates a linear map from $F^n$ to $F^m$ to a matrix from the set of $m \times n$ matrices with entries from a field F is an isomorphism of vector spaces.

Proof Attempt:
Now, we just have to prove that it is linear and that it is bijective. Let $A_1,A_2 \in M(m \times n,F)$. Then:
$f(A_1+A_2) = L_0:F^n \to F^m$
We have used the fact that every matrix has a linear map associated with it. In fact, the above means that $f$ is additively linear. So:
$f(A_1+A_2) = f(A_1)+f(A_2)$
Similarly, let $A \in M(m \times n, F)$ and $\lambda \in F$. Then:
$f(\lambda A) = L_1 : F^n \to F^m$
Since $L_1$ is a linear map, it follows that:
$f(\lambda A) = \lambda \cdot f(A)$
We, now, have to prove bijectivity. First, we prove that it is injective. Let $A_1,A_2 \in M(m \times n, F)$ such that:
$f(A_1) = f(A_2)$.
Now, this is going to be a linear map between $F^n$ and $F^m$. From a previously proven result, it is clear that every linear map has a unique matrix associated with it such that $x \to Ax$ for every $x \in F^n$. Hence:
$A_1 = A_2$
To prove surjectivity, we need to show that $f(M(m \times n, F)) = Hom(F^n,F^m)$. Let $L \in Hom(F^n,F^m)$. Since every linear map has a unique matrix associated with it, it suffices to say that there exists an $A \in M(m \times n, F)$ such that $f(A) = L$. However, that means that $L \in f(M(m \times n, F))$. This proves that $Hom(F^n,F^m) \subset f(M(m \times n,F))$ and, thus, establishes equality of the two sets.
That proves the desired result.
^Does my argument above work or no? What could I fix?

REPLY [1 votes]: You haven't defined $f$. Hints:
$(1).\ $ Let $A=(a_{ij})_{1\le i\le n;\ 1\le j\le m}$ and choose bases $\{e_i\}_{1\le i\le n}$ and $\{e'_i\}_{1\le i\le m}$ for $F^n$ and $F^m$ respectively.
$(2).\ f$ will send $A$ to a $\textit{function}\  \phi: F^n\to F^m.$ So, to define $f$, you need to say how $\phi$ acts.
$(3).\ $But the action of $\phi$ is specified by defining it on the basis $\{e_i\}_{1\le i\le n}$ and extending by linearity (this way, you force $\phi$ to be linear), which means that $\phi$ is determined by $(\alpha_{ij})_{1\le i\le n;\ 1\le j\le m}$ such that $\phi(e_i)=\sum^m_{j=1} \alpha_{ij}e_j'.$
$(4).\ $ Make the obvious choice for the $\alpha_{ij}$, and $\textit{now}$ check that $f$ is linear.


text:
\begin{document}

\title{Sparsity in Reservoir Computing Neural Networks
\thanks{\textbf{This paper is currently under review.}
This work has been partially supported by the European Union’s Horizon 2020 Research and Innovation program, under project TEACHING (Grant agreement ID: 871385).
https://www.teaching-h2020.eu}
}

\author{\IEEEauthorblockN{Claudio Gallicchio}
\IEEEauthorblockA{\textit{Department of Computer Science} \\
\textit{University of Pisa}\\
Pisa, Italy \\
gallicch@di.unipi.it}
}

\maketitle

\begin{abstract}
Reservoir Computing (RC) is a well-known strategy for designing Recurrent Neural Networks featured by striking efficiency of training. The crucial aspect of RC is to properly instantiate the hidden recurrent layer that serves as dynamical memory to the system. In this respect, the common recipe is to create a pool of randomly and \emph{sparsely} connected recurrent neurons. While the aspect of sparsity in the design of RC systems has been debated in the literature, it is nowadays understood mainly as a way to enhance the efficiency of computation, exploiting sparse matrix operations.\\
In this paper, we empirically investigate the role of sparsity in RC network design under the perspective of the richness of the developed temporal representations. We analyze both sparsity in the recurrent connections, and in the connections from the input to the reservoir. Our results point out that sparsity, in particular in input-reservoir connections, has a major role in developing internal temporal representations that have a longer short-term memory of past inputs and a higher dimension.
\end{abstract}

\begin{IEEEkeywords}
Reservoir Computing, Echo State Networks, Short-term Memory, Sparse Recurrent Neural Networks
\end{IEEEkeywords}

\section{Introduction}
Recurrent Neural Networks
(RNNs) \cite{kolen2001field} are a fundamental tool for adaptive processing of dynamically evolving information, with excellent performance in fields such as
time-series forecasting \cite{laptev2017time},
machine translation \cite{sutskever2014sequence},
speech and text processing \cite{graves2013speech, nallapati2016abstractive}, just to mention a few.
An increasing number of works are analyzing the role of sparsity in the design of trained (dynamical) neural networks systems, for example through pruning \cite{narang2017exploring} or re-wiring \cite{bellec2017deep} connections. The characterization emerging from these studies is that having sparse connections between neurons
is not only advantageous in computational terms - as it enables fast sparse matrix computations - but can also be beneficial to obtain a better performance in practice.
Moreover, in the context of neurobiologically-inspired information processing systems, a sparse degree of connectivity between neurons has been shown to improve the quality of the developed internal representations \cite{litwin2017optimal}. Interestingly, the optimal amount sparsity in the numerical simulations matched observed properties of cerebellum-like circuits.

Reservoir Computing (RC) neural networks \cite{lukovsevivcius2009reservoir,schrauwen2007overview,jaeger2004harnessing} represent an intriguing development in the field of RNNs.
In RC, the recurrent hidden layer of a RNN is left untrained after initialization subject to asymptotic stability conditions of the corresponding dynamical system. As a result, learning is applied only to a simple readout component with striking advantages in terms of required training times compared to fully trained RNNs. Pushing the involved algorithms towards extreme simplicity and efficiency makes the RC approach very well suited for real-world application scenarios featured by (possibly severe) resource constraints, such as neuromorphic hardware implementations \cite{larger2012photonic} or cyber-physical systems where the learning modules are embedded at the edge \cite{bacciu2014experimental}.

A typical strategy in the design of RC networks is to setup the recurrent layer in a sparse way.
The initial intuition was that sparsity in the recurrent untrained layer could enable a decoupling of state variables and hence richer representations \cite{jaeger2001echo}.
Successively, several authors pointed out empirical evidences contrary to the initial intuition (see, e.g.,
\cite{schrauwen2007overview, xue2007decoupled, gallicchio2011architectural}). Currently, the sparse design of reservoirs is commonly understood mainly as a way to speedup state computations, without a practical effect on the resulting performance.
However, the impact of sparsity on the performance of RC neural networks has been typically studied limited to the recurrent connections only.
In this paper, we intend to shed more light on the role of sparsity in RC by extending the analysis to both recurrent and input connections.
Specifically, we empirically show the effect of recurrent and input sparsity in reservoirs, evaluated by means of short-term memory capacity and effective dimension of the resulting state trajectories.

The rest of this paper is structured as follows.
We introduce the basics of RC methodology in Section~\ref{sec.RC}, discussing initialization and sparsity of reservoirs. Then, in Section~\ref{sec.richness} we present the concepts of short-term memory capacity and effective reservoir dimension. Our experimental analysis is described in Section~\ref{sec.experiments}. Finally, in Section~\ref{sec.conclusions} we draw our conclusions and sketch possible developments.

\section{Reservoir Computing Neural Networks}
\label{sec.RC}
Here we give a brief description of the RC design methodology for RNNs, focusing on the Echo State Network (ESN) \cite{jaeger2004harnessing,jaeger2001echo} model.

An RC network is a neural information processing system that treats data in the form of (temporal) sequences. Architecturally, the neural network is composed by a hidden recurrent layer called \emph{reservoir}, and an output layer called \emph{readout}. Fig.~\ref{fig.esn} illustrates the building blocks of a typical RC network.

\begin{figure}[htbp]
\centering
\includegraphics[width = 0.5\columnwidth]{esn.pdf}
\caption{Architecture of an RC neural network. The dotted arrow indicates trained connections.}
\label{fig.esn}
\end{figure}

In what follows, we denote the number of reservoir neurons, i.e., the reservoir dimension, by $N$, and the state of the reservoir system at time $t$ by $\h(t) \in \R^{N}$.
This state is evolved by following a state update equation:
\begin{equation}
\label{eq.reservoir}
\h(t) = \tanh(\U \x(t) + \W \h(t-1)),
\end{equation}
where $\x(t) \in \R^{M}$ is the M-dimensional input at time $t$, $\U \in \R^{N \times M}$ is the input weight matrix, modulating the influence of the external input on the current state, and $\W \in \R^{N \times N}$ is the recurrent weight matrix, which controls the impact of previous state on the current state. The state is typically set to a zero vector as initial condition, i.e. $\h(0) = \mathbf{0} \in \R^{N}.$
Note that here we dropped from \eqref{eq.reservoir} the reference to bias terms to focus the analysis on the external stimulating input signal alone. Both weight matrices $\U$ and $\W$ remain untrained after initialization (see Section~\ref{sec.initialization}).

The reservoir system is coupled with a linear readout layer that computes an $L$-dimensional output at each time step, i.e. $\y(t) \in \R^{L}$, as an affine transformation of the reservoir state:
\begin{equation}
\label{eq.readout}
\y(t) = \V \h(t) + \mathbf{b},
\end{equation}
where $\V \in \R^{L \times N}$ is a readout weight matrix and $\mathbf{b} \in \R^{L}$ is a bias vector (that assumes a constant unitary input bias for the readout). The readout parameters are the only ones that undergo a training process, typically in closed-form fashion by using pseudo-inversion \cite{lukovsevivcius2009reservoir}.

\subsection{Initialization of Reservoirs}
\label{sec.initialization}
The fundamental characterization of RC neural networks is that all the reservoir parameters remain untrained after initialization. Such initialization is performed in agreement to asymptotic stability conditions expressed by the Echo State Property (ESP) \cite{yildiz2012re,jaeger2001echo,gallicchio2019chasing}, which essentially require to control the magnitude of the weights in $\U$ and $\W$. Usually, both the input weights in $\U$ and the recurrent weights in $\W$ are randomly drawn from a  uniform distribution in $[-1,1]$. After that, the elements in $\U$ are re-scaled by a factor $\omega_{in}$, which takes the role of input scaling. The weights in $\W$ are re-scaled to control the largest absolute eigenvalue, i.e., the spectral radius $\rho$, typically to a value smaller than 1 \cite{jaeger2001echo}.

The design strategy of the reservoir topology (i.e., the way in which the reservoir neurons are connected among each other) has been subject of several studies in literature (see, e.g., \cite{strauss2012design,rodan2010minimum}). While some of the proposed reservoir organizations can be beneficial in specific application circumstances, a random and sparse topological organization of the reservoir is the architecture of choice in general cases. This is the focus of our analysis in this paper.

Making the connections among reservoir neurons sparse has the fundamental practical advantage to reduce the cost of state update operations in \eqref{eq.reservoir}. Actually, for densely connected reservoirs (and assuming $N >> M$) the cost of state updating scales as $\mathcal{O}(N^2)$, i.e. quadratically with the reservoir size.
A first
approach to make the reservoir sparsely connected would be to impose a (small) fixed percentage, say $C$, of non-zero weights in the involved weight matrices matrices. Although reducing the running times in practice for smaller reservoirs, this approach would asymptotically scale as $\mathcal{O}(N^2 \, C/100)$, hence still quadratically with the reservoir size. A more effective approach, which is adopted in this paper, is to fix the number, say $\chi_R$, of incoming recurrent connections for each reservoir unit. This indeed makes the state update cost as small as $\mathcal{O}(N \, \chi_R)$, i.e. scaling only linearly with the number of neurons in the reservoir.
A similar strategy can be adapted for the setup of the input connections. In this case, to ensure that each input dimension is actually forwarded to the reservoir, we fix the number of outgoing connections from each input units, denoted as $\chi_I$. The sparse architectural reservoir setup used in this paper is exemplified in Figure~\ref{fig.sparsity}. Notice that in this case, every row of $\W$ has exactly $\chi_R$ non-zero values, and every column of $\U$ has exactly $\chi_I$ non-zero elements, with both $\chi_R$ and $\chi_I$ being not greater than $N$.

\begin{figure}[htbp]
\centering
\includegraphics[width = 1\columnwidth]{reservoir.pdf}
\caption{Illustration of sparsity in input to reservoir and recurrent reservoir connections. $\chi_R$ indicates the number of incoming recurrent connections for each reservoir unit. $\chi_I$ indicates the number of outgoing connections from each input unit.}
\label{fig.sparsity}
\end{figure}

\section{Short-term Memory and Effective Reservoir Space Dimension}
\label{sec.richness}
The role of the recurrent reservoir system is to embed the  input time-series into an internal ``state'' representation, given by the activation of the reservoir neurons over time. Here we analyze the quality of such internal reservoir representation by quantifying its short-term memory and effective dimension.\\

\noindent
\emph{Short-term Memory Capacity (MC)} \cite{jaeger2001short} tests the ability of a recurrent neural system to reconstruct its driving input time-series from the transient state dynamics. More in detail, the reservoir is driven by a uni-dimensional time-series, $x(t), t = 1, 2, \ldots$, and different readout units are trained to recall progressively delayed versions of the input. I.e., the i-th readout unit $y_i(t)$ should approximate $x(t-i)$. The MC of an RC network is then quantified as follows:
\begin{equation}
    \label{eq.mc}
    MC = \sum_{i = 1}^\infty \frac{cov^2(x(t-i),y_i(t))}{\sigma^2(x(t-i))
    \sigma^2(y_i(t))},
\end{equation}
i.e., as the sum of squared correlation coefficients of the delayed input and reconstructed signals.\\

\noindent
\emph{Effective Dimension ($N_{eff}$)}\cite{abbott2011interactions, litwin2017optimal} is a measure of the number of orthogonal directions in the neuronal system's state trajectory over time. While the evolution of the reservoir system in \eqref{eq.reservoir} is described by an $N$-dimensional state vector $\h(t)$, the actual reservoir trajectory lies into a lower-dimensional manifold whose dimension can be quantified as follows:
\begin{equation}
    \label{eq.dimension}
    N_{eff} = \frac{(\sum_{i = 1}^N \lambda_i)^2}{\sum_{i=1}^N\lambda_i^2},
\end{equation}
where $\lambda_i$, $i = 1, 2, \ldots, N$, denote the eigenvalues of the covariance matrix of the reservoir state activation over time.
When measured for a reservoir under the driving influence of an external time-series, \eqref{eq.dimension} gives an estimate of the number of directions of reservoir state variability that are (linearly) uncorrelated along the observed trajectory.

\section{Experimental Analysis}
\label{sec.experiments}

We measured the short-term memory (MC) and the effective reservoir dimension ($N_{eff}$) introduced in Section~\ref{sec.richness} for RC networks varying the amount of recurrent and input connections.
Our experimental settings are described in Section~\ref{sec.settings}, while the results are reported in Section~\ref{sec.results}.

\subsection{Settings}
\label{sec.settings}
We used a uni-dimensional signal as driving input for the reservoir (i.e., $M = 1$). To maximally test the intrinsic  quality of reservoir representations, we used iid randomly sampled inputs $x(t)$ from a uniform distribution (in $[-0.8, 0.8]$).
The length of the generated input time-series was $6000$, and the number of reservoir neurons was fixed to $N = 100$.
To compute MC, we used the first $5000$ time-steps as training set\footnote{We used
pseudo-inversion to train the readout, discarding the first 1000 time-steps as initial transient.}, using the remaining $1000$ time-steps to assess the MC score. The total number of delays used for the computation of \eqref{eq.mc} was 200, which is in practice sufficient to account for all the non-negligible contributions for 100-dimensional reservoirs. The last $1000$ time-steps of the dataset were also used to compute the effective reservoir dimension $N_{eff}$ (see \eqref{eq.dimension}).

In our experiments, we used RC networks with spectral radius $\rho = 0.9$ and input scaling $\omega_{in} = 1$. While this setup is of common use in RC practice, we also ran preliminary
experiments with
other choices of these hyper-parameters, finding that the outcomes are not qualitatively different.
We varied both the number of recurrent connections ($\chi_R$) and of input connections ($\chi_I$) from 1 to 100 (with step of 1). For each configuration  we averaged the results over 50 reservoir realizations.

\subsection{Results}
\label{sec.results}

\begin{figure*}[htbp]
\centering
\begin{subfigure}{0.48\textwidth}
\centering
\includegraphics[width = \textwidth]{MC.png}
\label{fig.MC}
\end{subfigure}
\begin{subfigure}{0.48\textwidth}
\centering
\includegraphics[width = \textwidth]{DIM.png}
\label{fig.dim}
\end{subfigure}
\caption{Short-term Memory Capacity (MC) and effective dimension ($N_{eff}$) of RC networks. Results corresponds to $N = 100$ reservoir neurons, spectral radius $\rho = 0.9$, and input scaling $\omega_{in} = 1$. Recurrent ($\chi_R$) and input ($\chi_I$) connectivity varied from 1 to 100 with step 1. For each of the 10000 configurations the results are averaged over a number of 50 reservoir realizations.}
\label{fig.results}
\end{figure*}

The achieved values of MC and $N_{eff}$ in correspondence of the possible sparsity settings (values of $\chi_R$ and $\chi_I$) are shown in Fig.~\ref{fig.results}.
We can draw two major observations from the results. First, the number of input connections has a decisive impact on both
the short-term memory and the effective reservoir dimension of the networks.
Indeed, maximally sparse input connections, with $\chi_I = 1$, achieved the highest performances.
Interestingly, simply propagating the input to all the reservoir neurons degrades the performance sensibly. Second, the role of sparsity in recurrent connections seems to be much less important. In fact, the trend in Fig.~\ref{fig.results} indicates that for a given input connectivity, the achieved results are not much sensible to the exact number of recurrent connections (after a minimum number has been exceeded).

The results are further detailed in Fig.~\ref{fig.results2}, which shows
the best result for each choice of input (resp. recurrent) connectivity in Fig.~\ref{fig.results2}(a) (resp. Fig.~\ref{fig.results2}(b)), as well as
the results achieved for maximally sparse input connectivity, i.e. for $\chi_I = 1$,
in Fig.~\ref{fig.results2}(c).
Figs.~\ref{fig.results2}(a)-(b) confirm the already observed trends. On the one hand the performance of the RC networks tends to deteriorate for less sparse input weight matrices. On the other hand, a modest number of recurrent connections is already sufficient to achieve a performance
not far from
the highest possible one. For RC networks with $\chi_I = 1$ (Fig.~\ref{fig.results2}(c)), both MC and $N_{eff}$ saturate for fairly small values of $\chi_R$, without appreciable differences for settings with more than $20$ recurrent connections per reservoir neuron.

\begin{figure*}[htbp]
\centering
\begin{subfigure}{0.32\textwidth}
\centering
\includegraphics[width = \textwidth]{Chi_in.png}
\label{fig.dim}
\end{subfigure}
\begin{subfigure}{0.32\textwidth}
\centering
\includegraphics[width = \textwidth]{Chi_r.png}
\label{fig.dim}
\end{subfigure}
\begin{subfigure}{0.32\textwidth}
\centering
\includegraphics[width = \textwidth]{results_1.png}
\label{fig.MC}
\end{subfigure}
\caption{Short-term Memory Capacity (MC) and effective dimension ($N_{eff}$) of 100 units RC networks, detailed for:
\textbf{(a)} best results  for increasing input connectivity;
\textbf{(b)} best results  for increasing recurrent connectivity;
\textbf{(c)} results  for maximally sparse input connections ($\chi_I = 1$) and increasing recurrent connectivity.
Results are re-scaled to $[0,1]$.}
\label{fig.results2}
\end{figure*}

\section{Conclusions}
\label{sec.conclusions}
We have empirically analyzed the performance of RC neural networks in relation to sparsity of input and recurrent connections. Our results indicate that under commonly used reservoir configurations, the number of non-zero connections can play a decisive role in determining the richness of the developed representations. In particular, while a modest number of recurrent connections is already sufficient to achieve good performance, we found that maximally sparse input to reservoir connections lead to the best results both in terms of short-term memory and in terms of effective dimension of the
state manifold. Overall, our analysis points out a simple rule of thumb for shaping reservoir weight matrices in case of uni-dimensional driving time-series: (i) connect the input to just one reservoir neuron, and (ii)
set a small number of incoming recurrent connections ($\approx 20\%$) for each reservoir neuron.

The study presented in this paper can be seen as preparatory to opening further and deeper lines of research.
First of all, the role of sparsity
can be investigated in synergy with \emph{structured} (rather than random) recurrent reservoir topologies, such as those based on cyclic \cite{rodan2010minimum} or small-world \cite{kawai2019small} connections.
Similarly, the study can be extended towards \emph{deep} RC neural networks \cite{gallicchio2017deep,gallicchio2018design}, where multiple reservoir layers are connected in a pipeline. In this case, the sparsity of input connections for higher layers has the even more intriguing role of modulating the extent of signal propagation between consecutive internal representations. Neuromorphic hardware implementations \cite{moughames2019three,freiberger2019towards,partzsch2011analyzing,larger2012photonic} of deep recurrent neural systems are an important example of a domain where such insights can be capitalized in practice.
Under a broader perspective, and outside the RC world, the analysis presented here pointed out that a sparse setting of RNN connections brings advantages even before learning of the non-zero connections. How these architectural advantages can be further exploited by (supervised or unsupervised) \emph{training} is another exciting open research question.

\bibliographystyle{IEEEtran}
\bibliography{references}

\end{document}

config: arxiv


TITLE: Solve: $4(9^x)+3(4^x)=13(6^x)$
QUESTION [1 upvotes]: I'm new to this site so bear with me if I make some mistakes. Here is the question:
Solve: $4(9^x)+3(4^x)=13(6^x)$
I have absolutely no clue how to attack this problem. I have tried a bit of logs but I come to a dead end. What could I do?

REPLY [2 votes]: Dividing both sides by $4^{x}$:
$4(\frac{9}{4})^{x} + 3 = 13(\frac{3}{2})^{x}$
$4(\frac{3}{2})^{2x} - 13(\frac{3}{2})^{x} + 3 = 0$
Note that this is a quadratic in $(\frac{3}{2})^{x}$, which can be factored as:
$(4(\frac{3}{2})^{x}-1)((\frac{3}{2})^{x}-3) = 0$
Thus, we either have $(\frac{3}{2})^{x} = \frac{1}{4}$ or $(\frac{3}{2})^{x}=3$, which gives a solution set of $\boxed{x = -\log_{\frac{3}{2}}(4), \log_{\frac{3}{2}}(3).}$
If you want, you can use change-of-base and log rules to get everything in normal $\log$s.


TITLE: Optimization problem on graph with weights on nodes and edges
QUESTION [1 upvotes]: I am solving a problem where I have a complete undirected graph with weights on the nodes and on the edges. The weight on the node represents a profit that you obtain if you select that node. The weight on the edge $(i,j)$ represent the distance to go from the node $i$ to node $j$.
The problem is to select a subset of nodes that maximize the profit, such that the distance traveled is less than a threshold $\lambda$.
The mathematical formulation of my idea is:
$$
maximize \; \sum_{i=1}^{n} p_{i}x_{i} \\
subject \;to \; \sum_{i = 1}^{n} \; \sum_{\{j = 1, \; j > i\}}^{n} d_{ij}x_{i}x_{j} \leq \lambda\\
x_{i} \in \{0,1\}, \; \forall i \in \{1,..., n\}
$$
Where $n$ is the number of nodes, $p_{i}$ is the profit of node $i$, and $d_{ij}$ is the distance to go from the node $i$ to node $j$.
The problem is similar to the knapsack 0-1, but having the weights on both nodes and edges, I can not solve it.
Do you have any ideas?

REPLY [1 votes]: You could solve your quadratic problem as a MIP (Mixed Integer Programming) problem using something like:
$$
\begin{array}{ll}
\max&\sum_i p_i x_i \\
    &y_{i,j} \ge x_i+x_j-1  \\
    &\sum_{i,j} d_{i,j} y_{i,j} \le \lambda \\
    & x_i \in \{0,1\} \\
    & y_{i,j} \in [0,1]
\end{array}
$$
where we restrict all $(i,j)\in A$.