# Computing the $L^{1}$ optimal transport density: a FEM approach

Federico Piazzon  Department of Mathematics “Tullio Levi-Civita”, University
of Padua, via Trieste 63, Padova, Italy<EMAIL_ADDRESS>,  Enrico Facca
Univ. Lille, Inria, CNRS, UMR 8524 - Laboratoire Paul Painlevé, F-59000 Lille,
France<EMAIL_ADDRESS>and  Mario Putti  Department of Mathematics
“Tullio Levi-Civita”, University of Padua, via Trieste 63, Padova, Italy
<EMAIL_ADDRESS>
###### Abstract.

The $L^{1}$ optimal transport density $\mu^{*}$ is the unique $L^{\infty}$
solution of the Monge-Kantorovich equations. It has been recently
characterized also as the unique minimizer of the $L^{1}$-transport energy
functional $\mathcal{E}$. In the present work we develop and we prove
convergence of a numerical approximation scheme for $\mu^{*}.$ Our approach
relies upon the combination of a FEM-inspired variational approximation of
$\mathcal{E}$ with a minimization algorithm based on a gradient flow method.

###### Key words and phrases:

Monge-Kantorovich Equations, Optimal Transport, Gradient Flow, Convex
Optimization with Positivity Constraints, Łojasiewicz Inequality

###### 1991 Mathematics Subject Classification:

68Q25 68R10 68U05

##  1\. Introduction

###  1.1. $L^{1}$ optimal transport

Optimal transport is the mathematical problem modeling the quest for an
optimal strategy for moving a mass displacement of a material to a target
distribution. Apart from the original motivation of soil dragging [27],
optimal transport finds a large variety of applications, e.g, metrics in
probability theory, smoothing in image processing, baricentric interpolation
and many others (see for instance [30, 31] and references therein). The modern
formulation of optimal transport is due to Kantorovich, see [23], and can be
stated as follows. Let $\nu^{+},\nu^{-}$ be two Borel measures defined on
$\mathbb{R}^{d}$ with finite equal masses. Let
$c:\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}\cup\\{+\infty\\}$
be a Borel function. Find a non-negative Borel Measure $\gamma$ on
$\mathbb{R}^{d}\times\mathbb{R}^{d}$ (termed _optimal transport plan_)
realizing the following infimum

(1.1)
$\inf\left\\{\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}c(x,y)d\gamma(x,y)\right\\},$

under the constraints

$\displaystyle\gamma(A,\mathbb{R}^{d})=\nu^{+}(A)\ \forall A\mbox{ Borel set
in }\mathbb{R}^{d},$ $\displaystyle\gamma(\mathbb{R}^{d},B)=\nu^{-}(B)\
\forall B\mbox{ Borel set in }\mathbb{R}^{d}.$

Existence of such an optimal transport plan is guaranteed by mild assumptions
on $c$, e.g., lower semicontinuity and boundedness from below.

In this paper we focus on the $L^{1}$ case, i.e., we assume $c(x,y):=\|x-y\|$.
Also, we will work under the (classical) further assumption of
$\nu^{+}=f^{+}dx$, $\nu^{-}=f^{-}dx$ for a given zero-mean function $f\in
L^{\infty}(\mathbb{R}^{d}).$ In this setting not only the problem (1.1) admits
a unique optimal transport plan, also such a plan is defined by an _optimal
transport map_ , i.e., a Borel function
$T:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ realizing the minimum in the Monge
formulation of optimal transport, i.e.

$\inf\left\\{\int_{\mathbb{R}^{d}}c(x,T(x))d\nu^{+},T_{\\#}\nu^{+}=\nu^{-}\right\\},$

where $T_{\\#}\nu^{+}$ is the push-forward measure.

The aim of the present work is to study an approximation algorithm for the
solution of the $L^{1}$-optimal transport. While the numerical solution of
$L^{2}$ optimal transport has been tackled in several previous works (see e.g.
[3, 5, 11] for recent contributions), schemes for the $L^{1}$ case are
sparser. Typical methods for the numerical solution of this latter problem are
based on the entropic regularization of the linear programming problem
associated to the Kantorovich formulation and augmented Lagrangian methods, or
on the reformulation of the problem in Beckmann form [13, 4, 24, 28].

In our work, we take inspiration from the work of [16, 17], where a numerical
solver for the PDE based formulation of this problem was proposed. The PDE
formulation of the $L^{1}$-optimal transport originated in [15] who showed
that the $L^{1}$-optimal transport map $T$ can be constructed starting from
the _optimal transport density_ , i.e., the unique $L^{\infty}$ solution
$\mu^{*}$ of the so-called _Monge-Kantorovich equations_

(1.2) $\begin{cases}-\operatorname{div}{(\mu^{*}\nabla
u^{*})}=f^{+}-f^{-},&\text{ in }\Omega\\\ |\nabla u^{*}|\leq 1,&\text{ in
}\Omega\\\ |\nabla u^{*}|=1\,&\mu^{*}\text{ a.e. in }\Omega\end{cases},$

where $\Omega$ is a bounded convex domain compactly containing
$\operatorname{supp}f$. The approximation of the solution of (1.2) is a very
difficult task. Indeed, the coupling of the _linear_ PDE with the constraints
on the gradient of $u$ give rise to a highly _nonlinear_ problem with
inequality constraints in the variable $\mu$. Also, due to the nonlinearity of
(1.2), it is not _a priori_ clear weather (and in which sense) the solution
$\mu_{n}$ of certain discretized version of (1.2) would converge to the
solution of (1.2).

The fundamental advance introduced in [16, 17] is to address the numerical
solution of the variational reformulation of the Monge-Kantorovich equation
presented by [9] and solve the related optimization problem by a gradient
method. In these works, the authors introduce a Lyapunov-candidate functional
and conjecture the existence of a continuous time gradient flow whose
asymptotic solution coincides with the optimal transport density. Relying on
this conjecture, the authors propose an efficient numerical scheme for the
integration of the flow based on the explicit Euler discretization of the
gradient flow and the finite element approximation of the relevant integrals.
We would like to note that this approach is based on the formulation of a
minimization problem and its subsequent discretization. The resolution of the
conjecture has been hampered by the lack of results related to existence and
uniqueness of the gradient flow and its trajectories, an outcome of the setup
of the minimization problem in non-separable infinite dimensional spaces. As a
consequence, no theoretical convergence of the above numerical scheme have
been obtained and only experimental results are available to justify this
approach.

The aim of the present work is to fill this gap and develop the complete
theoretical convergence analysis of the FEM-based numerical method for the
approximation of the $L^{1}$ optimal transport density. The properties of the
studied discretization naturally suggest that, in contrast with the explicit
nature of the method proposed in [16, 17], implicit Euler and Newton schemes
must combine effectively to provide a super-linearly convergent method for the
solution of the $L^{1}$-optimal transport problem.

The strategy that allows the completion of our goal can be summarized as
follows. Rather than tackling the proof of the conjecture, we reformulate the
problem by first discretizing the relevant functional and then defining a
gradient flow in a finite-dimensional setting. This discretize-minimize
approach, opposite to the minimize-discretize suggested in [17], defines a
sequence of finite-dimensional approximations of the above mentioned Lyapunov-
candidate functional that are designed to provide a sequence of minimizers
that converge towards the optimal transport density. The minimization step of
the proposed strategy exploits the real analyticity of these approximations to
derive a well-posed and smooth gradient flow that converges exponentially
towards the sought solution under mild hypotheses. Within this proposed
framework, we are able to prove the convergence of the proposed algorithm
towards the optimal transport density, to provide a sharp convergence rate
estimate, and to design robust and efficient stopping criteria. The final
scheme includes an adaptive time-stepping strategy that allows the use of a
geometrically increasing sequence of time-step sizes leading to super-
exponential convergence towards the minimizer, thus drastically improving the
computational performance with respect to previous work.

###  1.2. Our study

Our starting point is the variational characterization of the transport
density proposed in [18] and the following definition of a $L^{1}$-_transport
energy_ :

######  Definition 1.1 ($L^{1}$-Transport energy [18]).

Let $\Omega$ be a bounded Lipschitz domain of $\mathbb{R}^{n}$ and let $f\in
L^{\infty}(\Omega)$ be such that
$\operatorname{conv}(\operatorname{supp}f)\subset\subset\Omega$ and
$\int_{\Omega}fdx=0$. We denote by
$\mathcal{E}:\mathcal{M}^{+}(\Omega)\rightarrow[0,+\infty]$ the transport
energy functional defined by

(1.3)
$\mathcal{E}(\mu):=\sup_{u\in\mathscr{C}^{1}(\overline{\Omega}),\int_{\Omega}udx=0}\left(2\int_{\Omega}fu\,dx-\int_{\Omega}|\nabla
u|^{2}d\mu\right)+\int_{\Omega}d\mu.$

and we consider the following variational problem:

######  Problem 1.2 (Minimization of the transport energy).

Given $f,\Omega$ as above, find $\mu_{\mathcal{E}}$ such that

$\mathcal{E}(\mu_{\mathcal{E}})=\inf_{\nu\in\mathcal{M}^{+}(\Omega)}\mathcal{E}(\nu).$

Here and throughout the paper we denote by $\mathcal{M}(\Omega)$ the space of
Borel signed measures on $\Omega),$ by $\mathcal{M}^{+}(\Omega)$ the non-
negative Borel measures, and by $\mathcal{M}^{1}(\Omega)$ the space of Borel
probability measures. In addition, our work relies on the following
fundamental property, see [18, Prop. 2.1].

######  Proposition 1.3.

Problem 1.2 has an unique solution $\mu_{\mathcal{E}}$ which is an absolutely
continuous measure with

$d\mu_{\mathcal{E}}=\mu^{*}dx.$

We start by introducing a double indexed sequence of functionals
$\mathcal{E}_{n,m}$ corresponding to the restriction of $\mathcal{E}$ to
certain finite element spaces $\mathcal{M}_{n}^{+}$ and $\mathcal{W}_{m}$ on
which we approximate $\mu$ and $u$, respectively (Section 2). We consider the
diagonal sequence $\mathcal{E}_{n}:=\mathcal{E}_{n,n}$ and we show that its
$\Gamma$-limit is $\mathcal{E}_{\infty}$, the lower semicontinuous
regularization of the restriction of $\mathcal{E}$ to a dense subspace.
Nevertheless, we prove in Theorem 2.6 that

$\operatorname*{argmin}_{\mathcal{M}^{+}(\Omega)}\mathcal{E}_{\infty}\subseteq\operatorname*{argmin}_{\mathcal{M}^{+}(\Omega)}\mathcal{E},$

moreover _any sequence $\\{\mu_{n}^{*}\\}$ of minimizers of $\mathcal{E}_{n}$
converges to the unique minimizer of $\mathcal{E}$_ in the weak∗ topology of
measures, i.e.,

$\mu_{n}^{*}\rightharpoonup\mu^{*}.$

This result motivates the study of the numerical solution of the problems

$\mu_{n}^{*}\in\operatorname*{argmin}_{\mathcal{M}_{n}^{+}}\mathcal{E}_{n}.$

Next (Section 3), we study conditions that ensure well-posedness and well-
conditioning of these discrete minimum problems. In particular we show in
Proposition 3.2 that under the Hypothesis 3.1 there exists a unique minimizer
of $\mathcal{E}_{n}$ for any $n\in\mathbb{N}$, i.e., the discretized minimum
problems are well-posed. Under a slightly different assumption, namely
Hypothesis 3.3, we prove that such problems are also well-conditioned, see
Proposition 3.4.

The following step (Section 4) is to prove global existence, uniqueness, and
regularity of the following two gradient flows

$\begin{cases}\frac{d}{dt}\boldsymbol{\mu}=-\partial^{\circ}E_{n}(\boldsymbol{\mu}),&\forall
t>0\\\ \boldsymbol{\mu}(0)=\boldsymbol{\mu}_{0}\in\mathbb{R}_{\geq
0}^{N}&\end{cases}\;\;,\;\;\;\begin{cases}\frac{d}{dt}\boldsymbol{\sigma}=-\nabla
F_{n}(\boldsymbol{\sigma}),&\forall t>0\\\
\boldsymbol{\sigma}(0)=\boldsymbol{\sigma}^{0}\end{cases}\;.$

Namely the gradient flow equation (4.1) for the function $E_{n}$ which is
obtained from $\mathcal{E}_{n}:=\mathcal{E}_{n,n}$ with the use of canonical
coordinates on $\mathcal{M}_{n}^{+}\cong\mathbb{R}_{\geq 0}^{N}$, and the
gradient flow equation (4.8) for the function $F_{n}$ which is obtained from
$E_{n}$ by composition with the coordinate square map (i.e.,
$\mathbb{R}^{N}\ni\boldsymbol{\sigma}\mapsto(\sigma_{1}^{2},\sigma_{2}^{2},\dots\sigma_{N}^{2})\in\mathbb{R}_{\geq
0}^{N}$). In particular, in Theorem 4.1, we prove that any trajectory of the
first flow is converging to a global minimizer of $\mathcal{E}_{n}$, while we
prove in Theorem 4.3 that any trajectory of the latter is converging to a
local minimizer of $F_{n}$: applying the coordinate square map we recover a
global minimizer of $\mathcal{E}_{n}.$

The definition and study of the optimization algorithms for $E_{n}$ and
$F_{n}$ derived from the time discretization of the two flows above conclude
our theoretical contributions (Section 5). In the first case one needs to take
into account the non-negativity constraint on $\boldsymbol{\mu}$, so an
explicit time discretization seems to be necessary. Thus the gradient flow
(4.1) of $\mathcal{E}_{n}$ can be integrated by the so-called _projected
forward Euler scheme_ , which can be solved by Algorithm 1. On the contrary,
the gradient flow of $F_{n}$ is defined on $\mathbb{R}^{N}$, and this makes
feasible the use of an implicit scheme: our choice is indeed the _backward
Euler scheme_ coupled with the (zero finding) Newton’s Method with a
prescribed initial guess. We focus on the latter scheme, both because it is
more promising from a numerical point of view, e.g., there are no inequality
constraints, and because -roughly speaking- the higher smoothness of the
continuous time trajectories should provide faster convergence. We prove in
Theorem 5.1 the consistency of this scheme, i.e., the Newton’s Method
converges at each time step and the sequence of the backward Euler iterates
converges to a local minimizer of $F_{n}.$ Then we define Algorithm 2 by
choosing a particular _a posteriori_ stopping criterion that allows us to
prove a stability estimate (see Proposition 5.2). We finally prove in Theorem
5.3 the convergence of Algorithm 2 and give in Proposition 5.5 a sharp
estimate of the rate of such convergence. The combination of Theorem 2.6,
Theorem 5.1, and Theorem 5.3 leads to (see Corollary 5.4)

${\lim_{n}}^{*}\mathcal{I}_{n}\big{(}\lim_{\ell}\boldsymbol{\mu}_{n}^{\ell}\big{)}=\mu^{*}dx,$

where we denoted by $\lim^{*}$ the limit in the weak∗ topology of measures,
and by $\mathcal{I}_{n}$ the embedding $\mathbb{R}_{\geq
0}^{N}\hookrightarrow\mathcal{M}^{+}$ deriving from the FEM coordinates.

Finally (Section 2), we perform some numerical tests of Algorithm 2 in few
cases where the optimal transport density is known. Also we check numerically
the assumptions we made for the well-posedness and well-conditioning of the
discrete variational problems.

##  2\. $\Gamma$-convergence of finite elements discretization

###  2.1. Finite elements discretization of $\mathcal{E}$

We want to introduce and study a sequence of functionals that provides a well-
designed sequence of approximations of $\mathcal{E}$ for minimization
purposes. In order to be termed ”well designed” the sequence needs to fulfill
certain desirable properties as preserving the convexity of $\mathcal{E}$,
having (at each stage) finite dimensional domain, possibly being smooth, and
being relatively easy to be evaluated at any given point. In addition, since
we are interested in the minimization of $\mathcal{E}$, we need to be able to
prove that any sequence of minimizers converges to the minimizer of
$\mathcal{E}$. We introduce a suitable FEM-based sequence of approximations of
$\mathcal{E}$ in the present subsection and we will prove its convergence in
the next subsection.

Let us remark that, as pointed out in [15] and in [18], the role of $\Omega$
is not relevant in our work as long as it is an open bounded convex set
compactly containing the set
$\operatorname{conv}\left(\operatorname{supp}f^{+}\cup\operatorname{supp}f^{-}\right),$
where $\operatorname{conv}$ stands for the convex envelope operator and
$\operatorname{supp}$ for the support. Therefore, possibly up to translations
an dilations, we can assume $\Omega:=(0,1)^{d},$ a very convenient choice for
our purposes. This choice prevents the need of the analysis of the so-called
geometric error that occurs when considering approximations of $\Omega$ by
unions of simplices which are the elements of choice in our approach.

First, we pick a sequence of nested triangulations
$\mathcal{T}^{n}=\\{T_{j}^{n}\\}_{j=1}^{N}$ of $\Omega$ indexed over
$n\in\mathbb{N}$ and obtained by uniform refinements. Let us denote by
$\mathcal{P}_{1}(\mathcal{T}^{n})$ the space of continuous real functions $v$
such that $v|_{T_{j}^{n}}$ is an affine function for any $j=1,2,\dots,N.$
Similarly, $\mathcal{P}_{0}(\mathcal{T}^{n})$ is the space of functions $v\in
L^{\infty}(\Omega)$ such that $v|_{T_{j}^{n}}$ is equivalent to a constant
function for any $j=1,2,\dots,N.$ Then we can consider the finite dimensional
function spaces

$\mathcal{M}_{n}:=\mathcal{P}_{0}(\mathcal{T}^{n}),\;\;\mathcal{W}_{m}:=\left\\{u\in\mathcal{P}_{1}(\mathcal{T}^{m+1}):\int_{\Omega}udx=0\right\\}\;,$

and denote by $\mathcal{M}_{n}^{+}$ the non-negative cone in
$\mathcal{M}_{n}$, i.e.,

$\mathcal{M}_{n}^{+}:=\left\\{\mu\in\mathcal{P}_{0}(\mathcal{T}^{n}):\mu|_{T_{j}^{n}}\geq
0,\,\forall T_{j}^{n}\in\mathcal{T}^{n}\right\\}.$

Second, we pick a monotone sequence $\\{\delta_{m}\\}$ such that
$\delta_{m}\downarrow 0$ as $m\to+\infty$.

We can now introduce a double sequence
$\\{\mathcal{E}_{n,m}\\}_{n,m\in\mathbb{N}}$ of approximations of the
functional
$\mathcal{E}:\mathcal{M}^{+}(\Omega)\rightarrow\mathbb{R}\cup\\{+\infty\\}$
that is a good candidate for the variational approximation of $\mathcal{E}.$
Notice that we can re-write $\mathcal{E}$, for notational convenience, as

$\mathcal{E}(\mu):=\mathcal{L}(\mu)\;+\;\int_{\Omega}d\mu,$

where $\mathcal{L}(\mu):=\sup_{u\in
L^{2}(\Omega),\int_{\Omega}udx=0}\mathcal{D}(\mu,u)$ and

$\mathcal{D}(\mu,u):=\begin{cases}2\int_{\Omega}fu\,dx-\int_{\Omega}|\nabla
u|^{2}d\mu&\text{ if
}u\in\mathscr{C}^{1}(\overline{\Omega}),\;\int_{\Omega}udx=0\\\ +\infty&\text{
otherwise}\end{cases}.$

Then we define

(2.1)
$\mathcal{E}_{n,m}(\mu):=\begin{cases}\sup_{u\in\mathcal{W}_{m}}\mathcal{D}_{m}(\mu,u)+\int_{\Omega}d\mu&\text{
if }\mu\in\mathcal{M}_{n}^{+}\\\ +\infty&\text{ otherwise}\end{cases},$

where, for any $\mu\in L^{\infty}_{+}(\Omega):=\\{\mu\in
L^{\infty}(\Omega):\,\mu\geq 0\text{ a.e.}\\}$, we set

(2.2)
$\mathcal{D}_{m}(\mu,u):=\begin{cases}2\int_{\Omega}fu\,dx-\int_{\Omega}(\mu+\delta_{m})|\nabla
u|^{2}dx&\text{ if }u\in\mathcal{W}_{m}\\\ +\infty&\text{
otherwise}\end{cases}.$

######  Remark 2.1.

We would like to stress here that the use of the two indices $n$ and $m$ is
dictated by notational convenience. Indeed, as we will see in the next
sections, the principal object of our study will be the diagonal sequence
$\mathcal{E}_{n,n}$, and in this case it produces the wanted combination of
meshes and spaces $\mathcal{T}^{n}$ for $\mathcal{M}_{n}$ and
$\mathcal{T}^{n+1}$ for $\mathcal{W}_{n}$.

In our work we will use the functionals

(2.3) $\displaystyle\mathcal{E}_{\infty,m}(\mu)$
$\displaystyle:=\inf_{n}\mathcal{E}_{n,m}$ (2.4)
$\displaystyle\widetilde{\mathcal{E}}_{m}(\mu)$
$\displaystyle:=\operatorname*{\textsc{sc}^{\,-}}\mathcal{E}_{\infty,m}(\mu),$
(2.5) $\displaystyle\mathcal{E}_{\infty}(\mu)$
$\displaystyle:=\sup_{m}\widetilde{\mathcal{E}}_{m}(\mu),$

where $\operatorname*{\textsc{sc}^{\,-}}$ stands for the lower semicontinuous
envelope with respect to the weak∗ topology of measures, i.e.,

$\operatorname*{\textsc{sc}^{\,-}}\mathcal{F}(\mu):=\sup\left\\{\mathcal{G}(\mu),\;\mathcal{G}\leq\mathcal{F},\;\mathcal{G}\text{
is l.s.c. in the weak${}^{*}$ topology}\right\\}.$

Here and throughout the paper, when we consider measures that are absolutely
continuous with respect to the Lebesgue measure, we use (by a slight abuse of
notation) the same symbol both for the density of a measure and the measure
itself. Accordingly, $\mu+\delta$ denotes the measure $\mu+\delta dx.$

Note that the definition of the functionals $\mathcal{E}_{n,m}$ is taylored to
the need of treating the convergence of minimizers. A more explicit
representation is obtained by means of the linear isomorphism
$\mathcal{I}_{n}$ mapping $\mathbb{R}^{N}$ onto $\mathcal{M}_{n}$. Precisely,
for any
$\boldsymbol{\mu}=\\{{\boldsymbol{\mu}}_{\\{1\\}},{\boldsymbol{\mu}}_{\\{2\\}},\dots,{\boldsymbol{\mu}}_{\\{N\\}}\\}\in\mathbb{R}^{N}$
we can define

$\mathcal{I}_{n}(\boldsymbol{\mu}):=\sum_{i=1}^{N}{\boldsymbol{\mu}}_{\\{i\\}}\chi_{T_{i}^{n}}(x),$

and, for any $\mu\in\mathcal{M}_{n}$, we denote by $\boldsymbol{\mu}$ the
vector

(2.6)
$\boldsymbol{\mu}:=\mathcal{I}_{n}^{-1}(\mu)=\left(\frac{\int_{T_{1}^{n}}\mu
dx}{|T_{1}^{n}|},\frac{\int_{T_{2}^{n}}\mu
dx}{|T_{2}^{n}|},\dots,\frac{\int_{T_{N}^{n}}\mu dx}{|T_{N}^{n}|}\right)^{t}.$

Analogously, given a basis $\\{\phi_{m,1},\phi_{m,2},\dots,\phi_{m,M}\\}$ of
$\mathcal{W}_{m}$ (where $M:=\operatorname{dim}\mathcal{W}_{m}$), we can
identify any function $u\in\mathcal{W}_{m}$ with its coordinates
$\boldsymbol{u}\in\mathbb{R}^{M}$, where
$u(x)=\sum_{p=1}^{M}{\boldsymbol{u}}_{\\{p\\}}\phi_{m,p}(x).$

Given $n,m\in\mathbb{N}$ and $i\in\\{1,2,\dots,N\\}$, we us introduce the
matrices

(2.7)
$A^{(n,m,i)}:=\left(\int_{T_{i}^{n}}\nabla\phi_{m,p}\cdot\nabla\phi_{m,q}dx\right)_{p,q=1,\dots,M}$

and, for any $\boldsymbol{\mu}\in\mathbb{R}_{\geq 0}^{N}$, the matrix-valued
function (the stiffness matrix)

(2.8)
$A^{(n,m)}(\boldsymbol{\mu}):=\sum_{i=1}^{N}({\boldsymbol{\mu}}_{\\{i\\}}+\delta_{m})A^{(n,m,i)}.$

We extend such a definition to each $\mu\in L^{\infty}(\Omega)$, $\mu\geq 0$
a.e., by setting

(2.9)
$A^{(n,m)}(\mu):=\left(\int_{\Omega}(\mu+\delta_{m})\nabla\phi_{m,p}\cdot\nabla\phi_{m,q}dx\right)_{p,q=1,\dots,M}.$

We also denote by $\boldsymbol{f}_{m}$ the load vector
$({\boldsymbol{f}_{m}}_{\\{1\\}},\dots,{\boldsymbol{f}_{m}}_{\\{M\\}})^{t}$,
where ${\boldsymbol{f}_{m}}_{\\{p\\}}=\int_{\Omega}f\phi_{m,p}dx$,
$p=1,2,\dots,M.$

######  Remark 2.2.

We remark that, in order to distinguish indeces denoting an element of a
sequence from indeces refering to a component of a vector, we use the
following convention: in the latter instance we add parentheses to the index,
e.g., given the sequence of vectors
$\\{\boldsymbol{\mu}_{n}\\}_{n\in\mathbb{N}}$, we denote by
${\boldsymbol{\mu}_{n}}_{\\{i\\}}$ the $i$-th component of the vector
$\boldsymbol{\mu}_{n}$.

It is worth pointing out that the operation of taking the supremum of
$\mathcal{D}_{m}$ among $u\in\mathcal{W}_{m}$ (used in the definition of
$\mathcal{E}_{n,m}$ given in eq. (2.1)) corresponds to calculating the
$\mathcal{P}_{1}$-finite element approximation of the solution of

(2.10) $\begin{cases}-\operatorname{div}((\mu+\delta_{m})\nabla u)=f&\text{ in
}\Omega\\\ \partial_{n}u=0&\text{ on }\partial\Omega\\\ \int
udx=0\end{cases}.$

More precisely, given $\mu\in\mathcal{M}_{n}^{+}$, we have

(2.11)
$\sup_{u\in\mathcal{W}_{m}}\mathcal{D}_{m}(\mu,u)=\max_{u\in\mathcal{W}_{m}}\mathcal{D}_{m}(\mu,u)=\int
fu[\mu]dx=\int(\mu+\delta_{m})|\nabla u[\mu]|^{2}dx,$

where $u[\mu]$ is the unique element of $\mathcal{W}_{m}$ such that

(2.12) $\int_{\Omega}\nabla u[\mu]\cdot\nabla
v(\mu+\delta_{m})dx=\int_{\Omega}fvdx,\;\forall v\in\mathcal{W}_{m},$

in other words,

(2.13)
$A^{(n,m)}(\boldsymbol{\mu})\boldsymbol{u}[\boldsymbol{\mu}]=\boldsymbol{f}_{m}.$

In view of this we set $E_{n,m}(\boldsymbol{\mu}):=\mathcal{E}_{n,m}(\mu)$,
that is

(2.14)
$E_{n,m}(\boldsymbol{\mu})=\begin{cases}(\boldsymbol{f}_{m})^{t}[A^{(n,m)}(\boldsymbol{\mu})]^{-1}\boldsymbol{f}_{m}+\langle\boldsymbol{\mu};\boldsymbol{a}_{n}\rangle&\text{
if }\ {\boldsymbol{\mu}}_{\\{i\\}}\geq 0,\,i=1,2,\dots,N\\\ +\infty&\text{
otherwise}\end{cases},$

where we denoted by $\boldsymbol{a}_{n}$ the area vector
$\left(\int_{T_{1}^{n}}dx,\int_{T_{2}^{n}}dx,\dots,\int_{T_{N}^{N}}dx\right)^{t}.$

Let us consider the coordinate square map
$\cdot^{\boldsymbol{2}}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ defined by
setting

(2.15)
${\boldsymbol{\sigma}}^{\boldsymbol{2}}:=({\boldsymbol{\sigma}}_{\\{1\\}}^{2},{\boldsymbol{\sigma}}_{\\{2\\}}^{2},\dots,{\boldsymbol{\sigma}}_{\\{N\\}}^{2}).$

We can introduce the functional

$F_{n,m}(\boldsymbol{\sigma}):=E_{n,m}(\boldsymbol{\sigma^{2}}).$

Notice that the image under the map $\cdot^{\boldsymbol{2}}$ of any local
minimizer $\boldsymbol{\hat{\sigma}}$ of $F_{n,m}$ on $\mathbb{R}^{N}$ is a
local minimizer of $E_{n,m}$ on $\mathbb{R}_{\geq 0}^{N}.$ Moreover
${\boldsymbol{\hat{\sigma}}}^{\boldsymbol{2}}$ is a global minimizer of
$E_{n}$ on $\mathbb{R}_{\geq 0}^{N}$ due to the convexity of the objective.

The next proposition collects some properties of $E_{n,m}$ and $F_{n,m}$ that
will be useful later on.

######  Proposition 2.3 (Differential properties of $E_{n,m}$ and $F_{n,m}$).

For any $n,m\in\mathbb{N}$ the functional $E_{n,m}$ is convex on
$\mathbb{R}^{N}$ and real analytic on

$\mathbb{R}_{\delta_{m}}^{N}:=\\{\boldsymbol{\mu}\in\mathbb{R}^{N}:\,{\boldsymbol{\mu}}_{\\{i\\}}>-\delta_{m},\,i=1,2,\dots,N\\}\supset\mathbb{R}_{\geq
0}^{N}.$

The functional $F_{n,m}$ is real analytic on $\mathbb{R}^{N}.$ In particular,
$\forall\boldsymbol{\mu}\in\mathbb{R}_{\geq 0}^{N}$ and any
$\boldsymbol{\sigma}\in\mathbb{R}^{N}$ such that
${\boldsymbol{\sigma}}^{\boldsymbol{2}}=\boldsymbol{\mu}$, we have

(2.16) $\displaystyle\frac{\partial
E_{n,m}}{\partial{\boldsymbol{\mu}}_{\\{i\\}}}(\boldsymbol{\mu})=1-(\boldsymbol{u}[\boldsymbol{\mu}])^{t}A^{(n,m,i)}\boldsymbol{u}[\boldsymbol{\mu}]$
(2.17) $\displaystyle\frac{\partial
F_{n,m}}{\partial{\boldsymbol{\sigma}}_{\\{i\\}}}(\boldsymbol{\sigma})=2{\boldsymbol{\sigma}}_{\\{i\\}}\frac{\partial
E_{n,m}}{\partial{\boldsymbol{\mu}}_{\\{i\\}}}({\boldsymbol{\sigma}}^{\boldsymbol{2}})$
(2.18)
$\displaystyle\frac{\partial^{2}E_{n,m}}{\partial{\boldsymbol{\mu}}_{\\{i\\}}\partial{\boldsymbol{\mu}}_{\\{j\\}}}(\boldsymbol{\mu})=(\boldsymbol{u}[\boldsymbol{\mu}])^{t}A^{(n,m,i)}[A^{(n,m)}(\boldsymbol{\mu})]^{-1}A^{(n,m,j)}\boldsymbol{u}[\boldsymbol{\mu}]$
(2.19)
$\displaystyle\frac{\partial^{2}F_{n,m}}{\partial{\boldsymbol{\mu}}_{\\{i\\}}\partial{\boldsymbol{\mu}}_{\\{j\\}}}(\boldsymbol{\sigma})=4{\boldsymbol{\sigma}}_{\\{i\\}}{\boldsymbol{\sigma}}_{\\{j\\}}\frac{\partial^{2}E_{n,m}}{\partial{\boldsymbol{\mu}}_{\\{i\\}}\partial{\boldsymbol{\mu}}_{\\{j\\}}}({\boldsymbol{\sigma}}^{\boldsymbol{2}})+2\delta_{i,j}\frac{\partial
E_{n,m}}{\partial{\boldsymbol{\mu}}_{\\{i\\}}}({\boldsymbol{\sigma}}^{\boldsymbol{2}}),$

where $\boldsymbol{u}[\boldsymbol{\mu}]$ is defined in (2.13). The convex
subdifferential $\partial E_{n,m}(\boldsymbol{\mu})$ of $E_{n,m}$ can be
characterized on the closure of $\mathbb{R}_{\delta_{m}}^{N}$ by

(2.20) $\partial
E_{n,m}(\boldsymbol{\mu})=\left\\{\boldsymbol{\xi}\in\mathbb{R}^{N}:{\boldsymbol{\xi}}_{\\{i\\}}\leq\frac{\partial
E_{n,m}}{\partial{\boldsymbol{\mu}}_{\\{i\\}}}(\boldsymbol{\mu})\;\forall
i,\;{\boldsymbol{\xi}}_{\\{i\\}}=\frac{\partial
E_{n,m}}{\partial{\boldsymbol{\mu}}_{\\{i\\}}}(\boldsymbol{\mu})\,\forall
i:{\boldsymbol{\mu}}_{\\{i\\}}>0\right\\}.$

The minimal norm subdifferential
$\partial^{\circ}E_{n,m}(\boldsymbol{\mu}):=\operatorname*{argmin}_{\xi\in\partial
E_{n,m}(\boldsymbol{\mu})}\|\xi\|_{2}$ satisfies

(2.21)
$\big{(}\partial^{\circ}E_{n,m}(\boldsymbol{\mu})\big{)}_{i}=\begin{cases}\frac{\partial}{\partial{\boldsymbol{\mu}}_{\\{i\\}}}E_{n,m}(\boldsymbol{\mu})&\text{
if }{\boldsymbol{\mu}}_{\\{i\\}}>0\\\
\left(\frac{\partial}{\partial{\boldsymbol{\mu}}_{\\{i\\}}}E_{n,m}(\boldsymbol{\mu})\right)^{-}&\text{
if }{\boldsymbol{\mu}}_{\\{i\\}}=0\end{cases}$

###### Proof.

The convexity of $\mathcal{E}_{n,m}$ is elementary because it is defined as
the supremum of affine functionals and the composition with $\mathcal{I}_{n}$
clearly preserves such a property.

Notice that $A(\boldsymbol{\mu}+\delta_{n})$ is a positive definite (and thus
invertible) symmetric matrix for any $\boldsymbol{\mu}\in\mathbb{R}_{\geq
0}^{N}$ and linearly depending on $\boldsymbol{\mu}$. Since the matrix
inversion is a real analytic operation the function $E_{n,m}$ is real analytic
on $\mathbb{R}^{N}_{\delta_{m}}.$

Using $\partial
A^{-1}(\boldsymbol{\mu}+\delta_{h})/\partial{\boldsymbol{\mu}}_{\\{i\\}}=-A^{-1}(\boldsymbol{\mu}+\delta_{h})\partial
A^{-1}(\boldsymbol{\mu}+\delta_{h})/\partial{\boldsymbol{\mu}}_{\\{i\\}}A^{-1}(\boldsymbol{\mu}+\delta_{h})$
and (2.13), we obtain eq. 2.16 and similarly eq. 2.18. By the chain rule we
obtain eq. 2.17 and eq. 2.19. Notice that $E_{n,m}$ is the restriction to the
non-negative cone of a (similarly defined) real analytic function defined on
$\mathbb{R}_{\delta_{m}}^{N}.$ For $\epsilon>0$ small enough, using the
definition of convex subdifferential, we can write

$\displaystyle
E_{n,m}(\boldsymbol{\mu}+\epsilon{\boldsymbol{e}}_{\\{i\\}})-E_{n,m}(\boldsymbol{\mu})\geq\epsilon\langle\boldsymbol{\xi};{\boldsymbol{e}}_{\\{i\\}}\rangle=\epsilon{\boldsymbol{\xi}}_{\\{i\\}},\;\;\forall\boldsymbol{\xi}\in\partial
E_{n,m}(\boldsymbol{\mu}),\,\forall i$ $\displaystyle
E_{n,m}(\boldsymbol{\mu}-\epsilon{\boldsymbol{e}}_{\\{i\\}})-E_{n,m}(\boldsymbol{\mu})\geq\epsilon\langle\boldsymbol{\xi};-{\boldsymbol{e}}_{\\{i\\}}\rangle=-\epsilon{\boldsymbol{\xi}}_{\\{i\\}},\;\;\forall\boldsymbol{\xi}\in\partial
E_{n,m}(\boldsymbol{\mu}),\,\forall i:{\boldsymbol{\mu}}_{\\{i\\}}>0.$

Dividing by $\epsilon$ and passing to the limit as $\epsilon\to 0^{+}$ we get

$\displaystyle\frac{\partial
E_{n,m}}{\partial{\boldsymbol{\mu}}_{\\{i\\}}}(\boldsymbol{\mu})$
$\displaystyle\geq{\boldsymbol{\xi}}_{\\{i\\}},\;\forall\boldsymbol{\xi}\in\partial
E_{n,m}(\boldsymbol{\mu}),\,\forall i$ $\displaystyle\frac{\partial
E_{n,m}}{\partial{\boldsymbol{\mu}}_{\\{i\\}}}(\boldsymbol{\mu})$
$\displaystyle\leq{\boldsymbol{\xi}}_{\\{i\\}},\;\forall\boldsymbol{\xi}\in\partial
E_{n,m}(\boldsymbol{\mu}),\,\forall i:{\boldsymbol{\mu}}_{\\{i\\}}>0$

Thus (2.20) follows. Equation (2.21) is immediately obtained by minimizing the
$L^{2}$ norm over this set. ∎

######  Corollary 2.4.

Any critical point $\boldsymbol{\mu}_{n,m}^{*}\in\mathbb{R}_{\geq 0}^{N}$ for
$E_{n,m}$ (i.e., $\partial^{\circ}E_{n,m}(\boldsymbol{\mu}_{n,m}^{*})=0$) is a
global minimizer for $E_{n,m}$. If
$\boldsymbol{\sigma}_{n,m}^{*}\in\mathbb{R}^{N}$ is a local minimizer for
$F_{n,m}$, then
$\boldsymbol{\mu}_{n,m}^{*}:={(\boldsymbol{\sigma}_{n,m}^{*})}^{\boldsymbol{2}}$
is a global minimizer for $E_{n,m}$.

Also, using Proposition 2.3, we can derive a discrete version of Monge
Kantorovich equations (1.2).

######  Corollary 2.5 (Discrete Monge-Kantorovich equations).

For any $n,m\in\mathbb{N}$ there exists at least one minimizer $\mu_{n,m}^{*}$
of $\mathcal{E}_{n,m}$ on $\mathcal{\mathcal{missing}}\mathcal{M}_{n}^{+}.$ We
have $\mu\in\operatorname*{argmin}\mathcal{E}_{n,m}$ if and only if there
exists a (unique) $u\in\mathcal{W}_{m}$ such that the following equations hold
true

(2.22)
$\begin{cases}A^{(n,m)}(\boldsymbol{\mu})\boldsymbol{u}=\boldsymbol{f}_{m}&\\\
\boldsymbol{u}^{t}A^{(n,m,i)}\boldsymbol{u}=1&\,\forall
i:\,{\boldsymbol{\mu}}_{\\{i\\}}>0\\\
\boldsymbol{u}^{t}A^{(n,m,i)}\boldsymbol{u}\leq 1&\,\forall
i:\,{\boldsymbol{\mu}}_{\\{i\\}}=0\end{cases}$

###### Proof.

The function $E_{n,m}$ is convex and coercive (i.e.,
$\liminf_{\|\boldsymbol{\mu}\|_{1}\to+\infty}E_{n,m}(\boldsymbol{\mu})=+\infty$)
and thus the existence of a minimizer follows by the direct method.

Again by convexity the condition $0\in\partial E_{n,m}(\boldsymbol{\mu})$ is
_equivalent_ to $\boldsymbol{\mu}\in\operatorname*{argmin}E_{n,m}$ and thus
equivalent to $0=\partial^{\circ}E_{n,m}(\boldsymbol{\mu}).$ The latter
equations can be written precisely as (2.22) using eq. 2.21. ∎

###  2.2. Convergence of minimizers of $\mathcal{E}_{n}:=\mathcal{E}_{n,n}$ to
$\mu^{*}$

The aim of the present section is to prove the following result.

######  Theorem 2.6.

Let $\mu_{n}^{*}\in\operatorname*{argmin}\mathcal{E}_{n}.$ Then

(2.23) $\mu_{n}^{*}\rightharpoonup\mu^{*},$

the optimal transport density of the Monge-Kantorovich equations.

######  Remark 2.7.

We choose to focus our study on the ”diagonal sequence” of functionals
$\mathcal{E}_{n}:=\mathcal{E}_{n,n}$ to simplify our notation and avoid few
technicalities in our proofs. However, the reader can easily check that our
results still hold for any subsequence
$\\{\mathcal{E}_{n_{j},m_{j}}\\}_{j\in\mathbb{N}}$ such that $\\{n_{j}\\}$ and
$\\{m_{j}\\}$ are non-decreasing diverging sequences.

$\Gamma$-convergence is the main tool that we are to going to use in the proof
of Theorem 2.6, and, to fix our notation, we recall in the following
paragraphs the needed definitions and properties. Let $(X,\tau)$ be a
topological space and, for any $x\in X$, let us denote by $\mathcal{N}(x)$ the
filter of the neighborhoods of $x.$ Let
$f_{n}:X\rightarrow\overline{\mathbb{R}},$ $n\in\mathbb{N}.$ We define

$\displaystyle\left(\operatorname*{\Gamma-\liminf}_{n\to+\infty}f_{n}\right)(x)$
$\displaystyle:=\sup_{U\in\mathcal{N}(x)}\liminf_{n\to\infty}\inf_{y\in
U}f_{n}(y),$
$\displaystyle\left(\operatorname*{\Gamma-\limsup}_{n\to+\infty}f_{n}\right)(x)$
$\displaystyle:=\sup_{U\in\mathcal{N}(x)}\limsup_{n\to\infty}\inf_{y\in
U}f_{n}(y).$

If there exists a function $f:X\rightarrow\mathbb{R}\cup\\{-\infty,+\infty\\}$
such that

(2.24) $\left(\operatorname*{\Gamma-\limsup}_{n\to+\infty}f_{n}\right)(x)\leq
f(x)\leq\left(\operatorname*{\Gamma-\liminf}_{n\to+\infty}f_{n}\right)(x),\;\;\forall
x\in X,$

then we say that $f_{n}$ $\Gamma$-converges to $f$ with respect to the
topology $\tau$ and we write $f_{n}\xrightarrow[]{\Gamma}f$ or
$\operatorname*{\Gamma-\lim}f_{n}=f.$ Our main interest on this notion of
convergence is given by the following property (cfr. for instance [14, Cor.
7.20]). Assume that $f_{n}\xrightarrow[]{\Gamma}f$ and $x_{n}$ is a minimizer
of $f_{n}$. Then any cluster point $x$ of $\\{x_{n}\\}$ is a minimizer of $f$
and $f(x)=\limsup_{n}f_{n}(x_{n}).$ If moreover $x_{n}$ converges to $x$ in
the topology $\tau,$ then $f(x)=\lim_{n}f_{n}(x_{n}).$

A useful property of $\Gamma$-convergence is that it is well-behaving under
monotone limits of lower semicontinuous functionals, [14, Prop. 5.4, Rem. 5.5,
Prop. 5.7]. Indeed we can prove a preliminary result exploiting the two
monotonicities of the double sequence $\mathcal{E}_{n,m}.$

######  Proposition 2.8.

Let the functionals $\mathcal{E}$, $\mathcal{E}_{n,m}$
$\widetilde{\mathcal{E}}_{m}$ and $\mathcal{E}_{\infty}$ be defined as in
(2.1), (2.4), and (2.5). Then

(2.25)
$\operatorname*{\Gamma-\lim}_{m}\widetilde{\mathcal{E}}_{m}=\mathcal{E}_{\infty}$

and

(2.26)
$\operatorname*{\Gamma-\lim}_{n}\mathcal{E}_{n}=\operatorname*{\Gamma-\lim}_{m}\operatorname*{\Gamma-\lim}_{n}\mathcal{E}_{n,m}=\operatorname*{\Gamma-\lim}_{n}\operatorname*{\Gamma-\lim}_{m}\mathcal{E}_{n,m}=\mathcal{E}_{\infty}$

###### Proof.

Equation 2.25 follows easily observing that
$\\{\widetilde{\mathcal{E}}_{m}\\}$ is an increasing sequence of lower
semicontinuous functionals and hence its $\Gamma$-limit coincides with its
point-wise limit.

To prove the second statement, we first notice that

(2.27)
$\sup_{m}\mathcal{E}_{\infty,m}=\sup_{m}\inf_{n}\mathcal{E}_{n,m}=\inf_{n}\sup_{m}\mathcal{E}_{n,m}.$

The first equality follows by definition. For the second equality we proceed
as follows. We first notice that, if
$\mu\in\cup_{n\in\mathbb{N}}\mathcal{M}_{n}^{+},$ then there exists $\bar{n}$
such that, for any $n\geq\bar{n}$, $\mu\in\mathcal{M}_{n}^{+}$ and
$\mathcal{E}_{n,m}(\mu)=\mathcal{E}_{\bar{n},m}(\mu)$ for all
$m\in\mathbb{N}.$ Then we can write

$\displaystyle\sup_{m}\inf_{n}\mathcal{E}_{n,m}(\mu)=$
$\displaystyle\sup_{m}\begin{cases}\mathcal{E}_{\bar{n},m}(\mu)&\text{ if
}\mu\in\cup_{n\in\mathbb{N}}\mathcal{M}_{n}^{+}\\\ +\infty&\text{
otherwise}\end{cases}$ $\displaystyle=$
$\displaystyle\begin{cases}\sup_{m}\mathcal{E}_{\bar{n},m}(\mu)&\text{ if
}\mu\in\cup_{n\in\mathbb{N}}\mathcal{M}_{n}^{+}\\\ +\infty&\text{
otherwise}\end{cases}=\inf_{n}\sup_{m}\mathcal{E}_{n,m}(\mu)\;.$

If we pick the lower semicontinuous regularization of eq. 2.27 we get

(2.28)
$\operatorname*{\textsc{sc}^{\,-}}\sup_{m}\mathcal{E}_{\infty,m}=\operatorname*{\textsc{sc}^{\,-}}\sup_{m}\inf_{n}\mathcal{E}_{n,m}=\operatorname*{\textsc{sc}^{\,-}}\inf_{n}\sup_{m}\mathcal{E}_{n,m}.$

Note that, since iterated $\sup$ operators commute, using eq. 2.25 we have

(2.29)
$\operatorname*{\textsc{sc}^{\,-}}\sup_{m}\mathcal{E}_{\infty,m}=\sup_{m}\operatorname*{\textsc{sc}^{\,-}}\mathcal{E}_{\infty,m}=\operatorname*{\Gamma-\lim}_{m}\widetilde{\mathcal{E}}_{m}=\mathcal{E}_{\infty}.$

On the other hand, using again the monotonicity we have

$\displaystyle\operatorname*{\textsc{sc}^{\,-}}\sup_{m}\inf_{n}\mathcal{E}_{n,m}=$
$\displaystyle\sup_{m}\operatorname*{\textsc{sc}^{\,-}}\inf_{n}\mathcal{E}_{n,m}=\sup_{m}\operatorname*{\Gamma-\lim}_{n}\mathcal{E}_{n,m}$
(2.30) $\displaystyle=$
$\displaystyle\operatorname*{\Gamma-\lim}_{m}\operatorname*{\Gamma-\lim}_{n}\mathcal{E}_{n,m},$
(2.31)
$\displaystyle\operatorname*{\textsc{sc}^{\,-}}\inf_{n}\sup_{m}\mathcal{E}_{n,m}=$
$\displaystyle\operatorname*{\textsc{sc}^{\,-}}\inf_{n}\operatorname*{\Gamma-\lim}_{m}\mathcal{E}_{n,m}=\operatorname*{\Gamma-\lim}_{n}\operatorname*{\Gamma-\lim}_{m}\mathcal{E}_{n,m}$

The combination of equations (2.28), (2.29), (2.30), and (2.31) leads to

$\operatorname*{\Gamma-\lim}_{n}\operatorname*{\Gamma-\lim}_{m}\mathcal{E}_{n,m}=\operatorname*{\Gamma-\lim}_{m}\operatorname*{\Gamma-\lim}_{n}\mathcal{E}_{n,m}=\mathcal{E}_{\infty}.$

Since the iterated Gamma-limits exist and coincide, also the diagonal Gamma-
limit exists and it is equal to the iterated limits, i.e., the first equality
of (2.26) holds true. ∎

The above proven $\Gamma$-convergence ensures that any cluster point of a
sequence of minimizers of $\mathcal{E}_{n}$ is a minimizer of
$\mathcal{E}_{\infty},$ but it does not imply by itself the existence of the
cluster point. However, we are able to show that indeed this cluster point
does exist.

######  Proposition 2.9.

For any $n\in\mathbb{N}$, let
$\mu_{n}^{*}\in\operatorname*{argmin}\mathcal{E}_{n}.$ Then
$\\{\mu_{n}^{*}\\}$ is pre-compact in the weak∗ topology of measures. Any
cluster point $\mu_{\infty}^{*}$ of $\\{\mu_{n}^{*}\\}$ lies in
$\operatorname*{argmin}\mathcal{E}_{\infty}$ and we can extract a subsequence
$l\mapsto n_{l}$ such that

(2.32) $\mu_{n_{l}}^{*}\rightharpoonup\mu_{\infty}^{*}\;\;\text{ and
}\;\;\mathcal{E}_{n}(\mu_{n}^{*})\to\mathcal{E}_{\infty}(\mu_{\infty}^{*})\;\text{
as }l\to+\infty.$

###### Proof.

Notice that, denoting by $dx$ the $d$-dimensional Lebesgue measure,

(2.33)
$\sup_{n}\int\mu_{n}^{*}\,dx\leq\sup_{n}\mathcal{E}_{n}(\mu_{n}^{*})\leq\sup_{h}\mathcal{E}_{n}(dx)=\sup_{m}\mathcal{E}_{1,m}(dx).$

Let $u_{m}\in\mathcal{W}_{m}$ be the solution of (2.12) with $\mu=dx.$ Notice
that, denoting by $C$ the Poincaré constant of $\Omega$ and using (2.11), we
have

$\displaystyle\sup_{m}\mathcal{E}_{1,m}(dx)=$
$\displaystyle\sup_{m}2\int_{\Omega}fu_{m}dx-\int_{\Omega}(1+\delta_{m})|\nabla
u_{m}|^{2}dx+\int_{\Omega}dx$ $\displaystyle=$
$\displaystyle\sup_{m}\int_{\Omega}fu_{m}dx+\int_{\Omega}dx\leq\sup_{m}\|f\|_{2}\|u_{m}\|_{2}+m(\Omega)$
$\displaystyle\leq$
$\displaystyle\sup_{m}\left(\frac{C}{1+\delta_{m}}\|f\|_{\infty}^{2}+1\right)dx(\Omega)=(C\|f\|_{\infty}^{2}+1)dx(\Omega)<+\infty.$

Here we used the Poincaré Inequality and the weak formulation (2.10). It
follows that the sequence $\mu_{n}$ has bounded mass:

$\sup_{n}\int\mu_{n}^{*}<+\infty.$

Thus we can extract a weak∗ converging subsequence, and the limit is a
minimizer of the $\Gamma$-limit functional (i.e., of $\mathcal{E}_{\infty}$ by
2.25) due to basic properties of $\Gamma$-limits, cfr. e.g., [14, Cor. 7.20].
∎

The proof of Theorem 2.6 essentially relies on the uniqueness and the
$L^{\infty}$ regularity of the optimal transport density $\mu^{*}.$ Indeed,
under our assumptions, the optimal transport density is uniquely determined
and it is an absolutely continuous measure having $L^{\infty}$ density with
respect to the Lebesgue measure, cfr. [15, 1, 19]. This is a key element of
our construction because the functional $\mathcal{E}$ and
$\mathcal{E}_{\infty}$ coincide on $L^{\infty}_{+}$, as we state in
Proposition 2.11 below. We need to prove a continuity property first.

######  Lemma 2.10.

Given $\mu\in L^{\infty}_{+}(\Omega)$,then, for any $\delta>0$, there exists a
sequence $\\{\mu_{n}\\}$, with $\mu_{n}\in\mathcal{M}_{n}^{+}$ such that

  1. i)

$\mu_{n}\to\mu$ almost everywhere in $\Omega$,

  2. ii)

$\max\\{\|\mu\|_{\infty},\sup_{n}\|\mu_{n}\|_{\infty}\\}<+\infty$,

  3. iii)

for any $m\in\mathbb{N}$,

(2.34) $\lim_{n}\
\mathcal{D}_{m}(\mu_{n},u)+\int\mu_{n}dx=\mathcal{D}_{m}(\mu,u)+\int\mu.$

###### Proof.

Let us define

(2.35)
$\mu_{n}:=\sum_{i=1}^{N}{\boldsymbol{\mu}_{n}}_{\\{i\\}}\chi_{T^{n}_{i}}(x):=\sum_{i=1}^{N}\frac{\int_{T^{n}_{i}}\mu
dx}{|T^{n}_{i}|}\chi_{T^{n}_{i}}(x).$

By the Lebesgue Differentiation Theorem it follows that $\mu_{n}\to\mu$ almost
everywhere in $\Omega$. Clearly

$\|\mu_{n}\|_{\infty}\leq\max_{i=1,\dots,N}|{\boldsymbol{\mu}_{n}}_{\\{i\\}}|\leq\|\mu\|_{\infty}.$

By the Lebesgue Dominated Convergence Theorem we have that

$\lim_{n}\int_{\Omega}\mu_{n}\psi dx=\int_{\Omega}\mu\psi dx,\;\forall\psi\in
L^{1}(\Omega).$

Note that in particular $\lim_{n}\int_{\Omega}\mu_{n}dx=\int_{\Omega}\mu dx$
follows. Thus, noticing that $\nabla\phi_{m,p}\cdot\nabla\phi_{m,q}\in
L^{1}(\Omega)$, and using the notation
$A^{(\infty,m)}_{p,q}(\mu):=\int_{\Omega}\nabla\phi_{m,p}\cdot\nabla\phi_{m,q}(\mu+\delta_{m})dx,$
for any $\mu\in L^{\infty}_{+}(\Omega)$, we have

$\lim_{n}A^{(n,m)}(\mu_{n})=\lim_{n}A^{(\infty,m)}(\mu_{n})=A^{(\infty,m)}(\mu)$

for any $m\in\mathbb{N}.$ Also notice that, denoting by
$A^{(\infty,m)}(1-\delta_{m})$ the (strictly positive definite) stiffness
matrix of the FEM $\mathcal{P}^{1}$ discretization of the Neumann Laplacian on
zero-mean functions, we have

$\displaystyle
A^{(\infty,m)}(\mu)\succeq\delta_{m}A^{(\infty,m)}(1-\delta_{m})\succ 0\,,$
$\displaystyle A^{(n,m)}(\mu_{n})\succeq\delta_{m}A^{(n,m)}(1-\delta_{m})\succ
0\text{ uniformly in }n,$

where $\succeq$ denotes the canonical ordering of symmetric semidefinite
matrices. Hence all the considered matrices are invertible and we have

$\lim_{n}[A^{(n,m)}(\mu_{n})]^{-1}=\lim_{n}[A^{(\infty,m)}(\mu_{n})]^{-1}=[A^{(\infty,m)}(\mu)]^{-1}.$

Finally we have

$\displaystyle\lim_{n}\mathcal{D}_{m}(\mu_{n},u)=\lim_{n}\boldsymbol{f}_{m}^{t}[A^{(n,m)}(\mu_{n})]^{-1}A^{(n,m)}(\mu_{n})[A^{(n,m)}(\mu_{n})]^{-1}\boldsymbol{f}_{m}$
$\displaystyle=$
$\displaystyle\boldsymbol{f}_{m}^{t}[A^{(\infty,m)}(\mu)]^{-1}A^{(\infty,m)}(\mu)[A^{(\infty,m)}(\mu)]^{-1}\boldsymbol{f}_{m}=\mathcal{D}_{m}(\mu,u).$

∎

We wold like to remark that the above result can be interpreted as a
simplified application of $\Gamma$-convergence of quadratic forms on
$\mathcal{W}_{m}$.

######  Proposition 2.11.

Under our assumptions we have

(2.36)
$\displaystyle\mathcal{E}(\mu)\leq\mathcal{E}_{\infty}(\mu),\;\forall\mu\in\mathcal{M}^{+}(\Omega):\mathcal{L}(\mu)<+\infty,$
(2.37)
$\displaystyle\mathcal{E}(\mu)\geq\mathcal{E}_{\infty}(\mu),\;\forall\mu\in
L^{\infty}_{+}(\Omega).$

###### Proof.

For notational convenience we define the sets

$\mathcal{M}_{\infty}:=\cup_{n\in\mathbb{N}}\mathcal{M}_{n}^{+}\;\;\;\mathcal{W}_{\infty}:=\cup_{m\in\mathbb{N}}\mathcal{W}_{m}\,.$

We begin proving (2.36). Let us pick $\mu\in\mathcal{M}^{+}$. Then, using the
notation $\llcorner$ for the restriction of a function to a set, i.e.,

$f\llcorner_{A}(x):=\begin{cases}f(x)&\text{ if }x\in A\\\
+\infty&\text{otherwise}\end{cases}\;,$

we can write

$\displaystyle\mathcal{E}_{\infty}(\mu)=\sup_{m}\operatorname*{\textsc{sc}^{\,-}}\inf_{n}\mathcal{E}_{n,m}(\mu)$
$\displaystyle=$
$\displaystyle\operatorname*{\textsc{sc}^{\,-}}\left[\sup_{m}\left(\sup_{u\in\mathcal{W}_{m}}2\int
fudx-\int|\nabla u|^{2}(d\nu+\delta_{m}dx)+\int
d\nu\right)\llcorner_{\nu\in\mathcal{M}_{\infty}^{+}}\right](\mu)$
$\displaystyle=$
$\displaystyle\operatorname*{\textsc{sc}^{\,-}}\left[\left(\sup_{m}\sup_{u\in\mathcal{W}_{m}}2\int
fudx-\int|\nabla u|^{2}(d\nu+\delta_{m}dx)+\int
d\nu\right)\llcorner_{\nu\in\mathcal{M}_{\infty}^{+}}\right](\mu)$
$\displaystyle\geq$
$\displaystyle\operatorname*{\textsc{sc}^{\,-}}\left[\sup_{m}\sup_{u\in\mathcal{W}_{m}}2\int
fudx-\int|\nabla u|^{2}(d\nu+\delta_{m_{0}}dx)+\int
d\nu\right]\llcorner_{\nu\in\mathcal{M}_{\infty}^{+}}(\mu)$ $\displaystyle=$
$\displaystyle\operatorname*{\textsc{sc}^{\,-}}\left[\sup_{u\in\mathcal{W}_{\infty}}2\int
fudx-\int|\nabla
u|^{2}(d\nu+\delta_{m_{0}}dx)+\int(d\nu+\delta_{m_{0}}dx)\right]\llcorner_{\nu\in\mathcal{M}_{\infty}^{+}}(\mu)-\delta_{m_{0}}|\Omega|.$

Now, if $\nu\in\mathcal{M}_{\infty}^{+}$, we have

$\displaystyle\sup_{u\in\mathcal{W}_{\infty}}2\int fudx-\int|\nabla
u|^{2}(d\nu+\delta_{m_{0}}dx)+\int(d\nu+\delta_{m_{0}}dx)$ $\displaystyle=$
$\displaystyle\sup_{u\in H^{1}(\Omega),\;\int u=0}2\int fudx-\int|\nabla
u|^{2}(d\nu+\delta_{m_{0}}dx)+\int(d\nu+\delta_{m_{0}}dx)$ $\displaystyle\geq$
$\displaystyle\sup_{u\in\mathscr{C}^{1}(\overline{\Omega}),\;\int u=0}2\int
fudx-\int|\nabla u|^{2}(d\nu+\delta_{m_{0}}dx)+\int(d\nu+\delta_{m_{0}}dx)$
$\displaystyle=$ $\displaystyle\mathcal{E}(\nu+\delta_{m_{0}}),$

if instead $\nu\in\mathcal{M}^{+}\setminus\mathcal{M}_{\infty}^{+}$, by our
definition of restriction, we have

$\left(\sup_{u\in\mathcal{W}_{\infty}}2\int fudx-\int|\nabla
u|^{2}(d\nu+\delta_{m_{0}}dx)+\int(d\nu+\delta_{m_{0}}dx)\right)\llcorner_{\mathcal{M}_{\infty}^{+}}=+\infty.$

In other words, due to the arbitrariness of $m_{0}$, for any $\delta>0,$ we
have

(2.38)
$\mathcal{E}_{\infty}(\mu)\geq\operatorname*{\textsc{sc}^{\,-}}\left(\mathcal{E}\big{|}_{\mathcal{M}_{\infty}^{+}}\right)(\mu+\delta)-\delta|\Omega|.$

Taking the $\liminf$ as $\delta\downarrow 0$ and using the lower
semicontinuity of
$\operatorname*{\textsc{sc}^{\,-}}\left(\mathcal{E}\big{|}_{\mathcal{M}_{\infty}^{+}}\right)$,
we obtain

(2.39)
$\mathcal{E}_{\infty}(\mu)\geq\operatorname*{\textsc{sc}^{\,-}}\left(\mathcal{E}\big{|}_{\mathcal{M}_{\infty}^{+}}\right)(\mu).$

Since $\mathcal{M}^{+}(\Omega)$ is first countable, the relaxed functional
$\operatorname*{\textsc{sc}^{\,-}}\left(\mathcal{E}\big{|}_{\mathcal{M}_{\infty}^{+}}\right)$
has the following equivalent characterization (see [14, Prop. 3.6]).

(2.40)
$\displaystyle\forall\mu\in\mathcal{M}^{+}(\Omega)\;\exists\\{\mu_{j}\\}_{j\in\mathbb{N}}\rightharpoonup^{*}\mu:\;\mathcal{\operatorname*{\textsc{sc}^{\,-}}}\left(\mathcal{E}\big{|}_{\mathcal{M}_{\infty}^{+}}\right)(\mu)\geq\limsup_{j}\left(\mathcal{E}\big{|}_{\mathcal{M}_{\infty}^{+}}\right)(\mu_{j}),$
(2.41)
$\displaystyle\operatorname*{\textsc{sc}^{\,-}}\left(\mathcal{E}\big{|}_{\mathcal{M}_{\infty}^{+}}\right)(\mu)\leq\liminf_{j}\left(\mathcal{E}\big{|}_{\mathcal{M}_{\infty}^{+}}\right)(\mu_{j}),\;\;\forall\\{\mu_{j}\\}_{j\in\mathbb{N}}\rightharpoonup^{*}\mu.$

Let us pick $\mu_{j}$ as in (2.40): it is evident that, if
$\mathcal{E}_{\infty}(\mu)<+\infty$, then
$\mu_{j}\in\mathcal{M}_{\infty}^{+}$. As a consequence and because of the
lower semi-continuity of $\mathcal{E}$, we obtain

(2.42)
$\mathcal{E}_{\infty}(\mu)\geq\limsup_{j}\left(\mathcal{E}\big{|}_{\mathcal{M}_{\infty}^{+}}\right)(\mu_{j})\geq\liminf_{j}\mathcal{E}(\mu_{j})\geq\mathcal{E}(\mu).$

This concludes the proof of (2.36).

Before proving (2.37), we claim that, for any $\mu\in\mathcal{M}^{+}(\Omega),$

(2.43) $\lim_{\delta\to 0^{+}}\mathcal{E}(\mu+\delta)=\mathcal{E}(\mu).$

In order to prove this, we first notice that, again due to lower
semicontinuity, we have

(2.44) $\liminf_{\delta\to 0^{+}}\mathcal{E}(\mu+\delta)\geq\mathcal{E}(\mu).$

On the other hand, we can prove the reverse inequality for the $\limsup$ using
the definition of $\mathcal{E}$. Note that we assume that either the supremum
defining $\mathcal{E}(\mu+\delta)$ is indeed a maximum that is achieved for
some $u_{\delta},$ or there exists a sequence
$u_{s}\in\mathscr{C}^{1}(\overline{\Omega})$ such that

(2.45) $\mathcal{E}(\mu+\delta)=\lim_{s}2\int fu_{s}dx-\int|\nabla
u_{s}|^{2}(d\mu+\delta dx)+\int(d\mu+\delta dx).$

In the latter case we have

$\displaystyle\mathcal{E}(\mu+\delta)=$ $\displaystyle\lim_{s}2\int
fu_{s}dx-\int|\nabla u_{s}|^{2}(d\mu+\delta dx)+\int(d\mu+\delta dx)$
$\displaystyle\leq$ $\displaystyle\lim_{s}2\int fu_{s}dx-\int|\nabla
u_{s}|^{2}d\mu+\int d\mu+\delta|\Omega|$ $\displaystyle\leq$
$\displaystyle\sup_{u\in\mathscr{C}^{1}(\overline{\Omega}),\;\int udx=0}2\int
fudx-\int|\nabla u|^{2}d\mu+\int d\mu+\delta|\Omega|$ $\displaystyle=$
$\displaystyle\mathcal{E}(\mu)+\delta|\Omega|.$

In the former case we have

$\displaystyle\mathcal{E}(\mu+\delta)$ $\displaystyle=2\int
fu_{\delta}dx-\int|\nabla u_{\delta}|^{2}(d\mu+\delta dx)+\int(d\mu+\delta
dx)$ $\displaystyle\leq\sup_{u\in\mathscr{C}^{1}(\overline{\Omega}),\;\int
udx=0}2\int fudx-\int|\nabla u|^{2}d\mu+\int d\mu+\delta|\Omega|$
$\displaystyle=$ $\displaystyle\mathcal{E}(\mu)+\delta|\Omega|.$

Thus in both cases we have

(2.46) $\limsup_{\delta\to 0^{+}}\mathcal{E}(\mu+\delta)\leq\mathcal{E}(\mu).$

Equations (2.44) and (2.46) prove (2.43).

We are now ready to finish the proof of (2.37). To this aim, let us recall
that, if $g$ is a continuous function on a topological space $X$ and $Y$ is a
dense subset of $X$, then $\sup_{x\in Y}g(x)=\sup_{x\in X}g(x)$. Notice also
that, for $\mu\in L^{\infty}_{+}(\Omega)$ and any $\delta>0$ the function
$u\mapsto 2\int fudx-\int(\mu+\delta)|\nabla u|^{2}dx$ is continuous on
$X:=H^{1}(\Omega)\cap\\{u:\int_{\Omega}udx=0\\}$, while
$Y:=\\{u\in\mathscr{C}^{1}(\overline{\Omega}),\;\int_{\Omega}udx=0\\}$ is
dense in $X.$ Hence we have

$\sup_{u\in Y}\mathcal{D}(\mu+\delta,u)=\sup_{u\in X}2\int
fudx-\int(\mu+\delta)|\nabla u|^{2}dx.$

Using this last fact and (2.43), for any $\mu\in L^{\infty}_{+}(\Omega)$, we
have the following.

$\displaystyle\mathcal{E}(\mu)$
$\displaystyle=\lim_{m}\sup_{u\in\mathscr{C}^{1}(\overline{\Omega}),\;\int_{\Omega}udx=0}\mathcal{D}(\mu+\delta_{m},u)+\int_{\Omega}(\mu+\delta_{m})dx$
$\displaystyle=\lim_{m}\sup_{u\in
H^{1}(\Omega),\;\int_{\Omega}udx=0}\mathcal{D}(\mu+\delta_{m},u)+\int_{\Omega}(\mu+\delta_{m})dx$
$\displaystyle\geq\lim_{m}\sup_{u\in\mathcal{W}_{m}}\mathcal{D}(\mu+\delta_{m},u)+\int_{\Omega}\mu
dx$
$\displaystyle=\sup_{m}\sup_{u\in\mathcal{W}_{m}}\mathcal{D}_{m}(\mu,u)+\int_{\Omega}\mu
dx.$

Lemma 2.10 allows us to select $\mu_{n}\in\mathcal{M}_{n}^{+}$ such that
$\mu_{n}\to\mu$ almost everywhere in $\Omega$ and

$\sup_{u\in\mathcal{W}_{m}}\mathcal{D}_{m}(\mu,u)+\int_{\Omega}\mu
dx=\lim_{n}\sup_{u\in\mathcal{W}_{m}}\mathcal{D}_{m}(\mu_{n},u)+\int_{\Omega}\mu_{n}dx=\lim_{n}\mathcal{E}_{n,m}(\mu_{n}).$

Now we observe that, since we are defining the spaces $\mathcal{M}_{n}^{+}$ by
uniform refinements,

$\mathcal{E}_{n,m}(\mu_{n})=\begin{cases}\mathcal{E}_{s,m}(\mu_{n})&\text{ if
}s\geq n\\\ +\infty&\text{ otherwise}\end{cases}\,,$

then $\mathcal{E}_{n,m}(\mu_{n})=\inf_{s}\mathcal{E}_{s,m}(\mu_{n})$ for any
$n\in\mathbb{N}$. Hence

$\displaystyle\mathcal{E}(\mu)$
$\displaystyle\geq\sup_{m}\lim_{n}\inf_{s}\mathcal{E}_{s,m}(\mu_{n})\geq\sup_{m}\lim_{n}\operatorname*{\textsc{sc}^{\,-}}\inf_{s}\mathcal{E}_{s,m}(\mu_{n})$
$\displaystyle\geq\sup_{m}\liminf_{n}(\operatorname*{\textsc{sc}^{\,-}}\inf_{s}\mathcal{E}_{s,m})(\mu_{n})=\sup_{m}\liminf_{n}\widetilde{\mathcal{E}}_{m}(\mu_{n})\geq\sup_{m}\widetilde{\mathcal{E}}_{m}(\mu)$
$\displaystyle=\mathcal{E}_{\infty}(\mu).$

This concludes the proof of (2.37). ∎

We are now ready to conclude the proof of Theorem 2.6.

###### of Theorem 2.6.

Let us pick $\mu_{n}\in\operatorname*{argmin}\mathcal{E}_{n}$ for any
$h\in\mathbb{N}.$ Due to Proposition 2.9 there exists a weak∗ convergent
subsequence and a limit point $\bar{\mu}$ such that
$\bar{\mu}\in\operatorname*{argmin}\mathcal{E}_{\infty}.$ Recall that
$\mu^{*}\in L^{\infty}_{+}(\Omega)$ and notice that

$\displaystyle\mathcal{E}(\mu^{*})$ $\displaystyle\leq\mathcal{E}(\bar{\mu})$
$\displaystyle\text{ by the }\mathcal{E}\text{-optimality of }\mu^{*}$
$\displaystyle\mathcal{E}(\bar{\mu})$
$\displaystyle\leq\mathcal{E}_{\infty}(\bar{\mu})$ by (2.36)
$\displaystyle\mathcal{E}_{\infty}(\bar{\mu})$
$\displaystyle\leq\mathcal{E}_{\infty}(\mu^{*})$ $\displaystyle\text{ by the
}\mathcal{E}_{\infty}\text{-optimality of }\bar{\mu}$
$\displaystyle\mathcal{E}_{\infty}(\mu^{*})$
$\displaystyle\leq\mathcal{E}(\mu^{*})$ $\displaystyle\text{ by
\eqref{technicallemmaclaim2}}.$

Thus

$\mathcal{E}(\mu^{*})\leq\mathcal{E}(\bar{\mu})\leq\mathcal{E}_{\infty}(\bar{\mu})\leq\mathcal{E}_{\infty}(\mu^{*})\leq\mathcal{E}(\mu^{*}).$

Hence the above inequalities are indeed equalities and $\bar{\mu}$ is a global
minimizer of $\mathcal{E}.$

Since the unique global minimizer of $\mathcal{E}$ is $\mu^{*}$ then
$\bar{\mu}=\mu^{*}$ and the whole sequence $\mu_{n}^{*}$ converges to
$\mu^{*}$ in the weak∗ topology of measures. ∎

##  3\. Well-posedness and conditioning of finite dimensional minimum problems

We have seen that the sequence of minimum problems
$\boldsymbol{\mu}_{n}^{*}\in\operatorname*{argmin}_{\mathbb{R}_{\geq
0}^{N}}E_{n,n}$ converges to the optimal transport density, it is quite
natural to ask whether each of this problems is well-posed and well-
conditioned. We briefly discuss these questions under two different sets of
assumptions, namely Hypothesis 3.1 and Hypothesis 3.3 below. We will look at
both $E_{n}:=E_{n,n}$ and $F_{n}:=F_{n,n}$ defined in (2.14) and lines below.

######  Hypothesis 3.1.

There exists
$\boldsymbol{\mu}_{n}^{*}\in\operatorname*{argmin}_{\mathbb{R}_{\geq
0}^{N}}E_{n}$ such that

(3.1)
$\left(\boldsymbol{\mu}_{n}^{*}+\operatorname{Ker}\operatorname{Hess}E_{n}(\boldsymbol{\mu}^{*}_{h})\right)\cap\mathbb{R}_{\geq
0}^{N}=\\{\boldsymbol{\mu}_{n}^{*}\\}.$

######  Proposition 3.2 (Well-posedness under Hypothesis 3.1).

If Hypothesis 3.1 holds true, then the functional $E_{n}$ admits the unique
minimizer $\boldsymbol{\mu}_{n}^{*}$.

###### Proof.

We first prove that there exists $\lambda>0$ such that, for any
$\nu\in\mathbb{R}_{\geq 0}^{N}$, we have

(3.2) $\inf_{\boldsymbol{\nu}\in\mathbb{R}_{\geq
0}^{N}}(\boldsymbol{\nu}-\boldsymbol{\mu}_{n}^{*})^{t}\operatorname{Hess}E_{n}(\boldsymbol{\mu}_{n}^{*})(\boldsymbol{\nu}-\boldsymbol{\mu}_{n}^{*})\geq\lambda\|\boldsymbol{\nu}-\boldsymbol{\mu}_{n}^{*}\|^{2}.$

For, let us consider the eigendecomposition of the symmetric positive
semidefinite matrix $\operatorname{Hess}E_{n}(\boldsymbol{\mu}_{n}^{*})$. More
precisely we pick an orthonormal basis
$\\{\boldsymbol{\theta}^{1},\dots,\boldsymbol{\theta}^{N}\\}$ of
$\mathbb{R}_{\geq 0}^{N}$, where
$K:=\operatorname{Ker}\operatorname{Hess}E_{n}(\boldsymbol{\mu}_{n}^{*})=\operatorname{span}\\{\boldsymbol{\theta}^{R+1},\dots\boldsymbol{\theta}^{N}\\}$,
$K^{\perp}=\operatorname{span}\\{\boldsymbol{\theta}^{1},\dots\boldsymbol{\theta}^{R}\\}$,
and, for suitable $\lambda_{1},\dots\lambda_{R}>0$, we have

$({\boldsymbol{\theta}}_{\\{i\\}})^{t}\operatorname{Hess}E_{n}(\boldsymbol{\mu}_{n}^{*})\boldsymbol{\theta}_{j}=\begin{cases}\lambda_{i}\delta_{i,j},&\text{
if }\max\\{i,j\\}\leq R\\\ 0,&\text{ otherwise}\end{cases}.$

Therefore we have

(3.3)
$(\boldsymbol{\nu}-\boldsymbol{\mu}_{n}^{*})^{t}\operatorname{Hess}E_{n}(\boldsymbol{\mu}_{n}^{*})(\boldsymbol{\nu}-\boldsymbol{\mu}_{n}^{*})=\sum_{i=1}^{R}\lambda_{i}|\langle\boldsymbol{\nu}-\boldsymbol{\mu}_{n}^{*};\boldsymbol{\theta}^{i}\rangle|^{2}\geq\min_{i\leq
R}\lambda_{i}\|\pi_{K^{\perp}}(\boldsymbol{\nu}-\boldsymbol{\mu}_{n}^{*})\|^{2},$

where we denoted by $\pi_{K^{\perp}}$ the orthogonal projection onto
$K^{\perp}.$ We claim that there exists $C>0$ such that

(3.4)
$\sqrt{C}\leq\inf_{\boldsymbol{\mu}\geq-\boldsymbol{\mu}_{n}^{*},\;\|\boldsymbol{\mu}\|=1}\|\pi_{K^{\perp}}\boldsymbol{\mu}\|=\inf_{\boldsymbol{\mu}\geq-\boldsymbol{\mu}_{n}^{*}}\frac{\|\pi_{K^{\perp}}\boldsymbol{\mu}\|}{\|\boldsymbol{\mu}\|}=\inf_{\boldsymbol{\nu}\in\mathbb{R}_{\geq
0}^{N}}\frac{\|\pi_{K^{\perp}}(\boldsymbol{\nu}-\boldsymbol{\mu}_{n}^{*})\|}{\|\boldsymbol{\nu}-\boldsymbol{\mu}_{n}^{*}\|}.$

The two equalities are quite obvious, while the inequality can be proven by
contradiction using the hypothesis (3.1) and the fundamental fact that we are
working in a finite dimensional space. Using (3.3) and (3.4) we can conclude
that (3.2) holds true with $\lambda:=C\min_{i\leq R}\lambda_{i}.$

Now we want to conclude that $\boldsymbol{\mu}_{n}^{*}$ is the unique
minimizer of $E_{n}$ in $\mathbb{R}_{\geq 0}^{N}$. Assume that we can pick
$\boldsymbol{\nu}_{h}^{*}\in\mathbb{R}_{\geq 0}^{N}$ with
$E_{n}(\boldsymbol{\nu}_{n}^{*})=E_{n}(\boldsymbol{\mu}_{n}^{*}),$ and define
$g:[0,1]\rightarrow\mathbb{R}$ with
$g(t):=E_{n}((1-t)\boldsymbol{\mu}_{n}^{*}+t\boldsymbol{\nu}_{n}^{*}).$ On one
hand, since $E_{n}$ is convex, so it is $g$. On the other hand, since
$E_{n}(\boldsymbol{\nu}_{n}^{*})=E_{n}(\boldsymbol{\mu}_{n}^{*})\geq g(t)$,
$g$ needs to be constant. However we also have

$g^{\prime\prime}(0)=(\boldsymbol{\nu}_{n}^{*}-\boldsymbol{\mu}_{n}^{*})^{t}\operatorname{Hess}E_{n}(\boldsymbol{\mu}_{n}^{*})(\boldsymbol{\nu}_{n}^{*}-\boldsymbol{\mu}^{*}_{n})\geq\lambda\|\boldsymbol{\nu}_{n}^{*}-\boldsymbol{\mu}_{n}^{*}\|^{2}>0.$

This is a contraddiction, so $\boldsymbol{\mu}_{n}^{*}$ is the unique
minimizer of $E_{n}$ in $\mathbb{R}_{\geq 0}^{N}$. ∎

######  Hypothesis 3.3.

There exists
$\boldsymbol{\mu}_{n}^{*}\in\operatorname*{argmin}_{\mathbb{R}_{\geq
0}^{N}}E_{n}$ and $\lambda>0$ such that,

(3.5)
$\displaystyle[\operatorname{Hess}E_{n}(\boldsymbol{\mu}_{n}^{*})]_{\llcorner\\{{\boldsymbol{e}}_{\\{i\\}}:{\boldsymbol{\mu}_{n}^{*}}_{\\{i\\}}>0\\}}\succeq\lambda\mathbb{I},$
(3.6) $\displaystyle\partial_{i}E_{n}(\boldsymbol{\mu}_{n}^{*})\neq
0,\;\;\forall i:{\boldsymbol{\mu}_{n}^{*}}_{\\{i\\}}=0.$

######  Proposition 3.4 (Well-posedness and well-conditioning under Hypothesis
3.3).

If Hypothesis 3.3 holds true, then the functional $E_{n}$ admits a unique
minimizer $\boldsymbol{\mu}_{n}^{*}$. Moreover, for any

(3.7)
$0<\Lambda<\min\left\\{4\lambda\min_{i:{\boldsymbol{\mu}_{n}^{*}}_{\\{i\\}}\neq
0}{\boldsymbol{\mu}_{n}^{*}}_{\\{i\\}}\;,\;2\min_{i:{\boldsymbol{\mu}_{n}}_{\\{i\\}}^{*}=0}\partial_{i}E_{n}(\boldsymbol{\mu}_{n}^{*})\right\\},$

there exists $R>0$ such that, for any $\boldsymbol{\sigma}_{n}^{*}$ such that
${(\boldsymbol{\sigma}_{n}^{*})}^{\boldsymbol{2}}=\boldsymbol{\mu}_{n}^{*}$
and for any $\boldsymbol{\mu}:={\boldsymbol{\sigma}}^{\boldsymbol{2}}$ with
$\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{n}^{*}\|<R,$ we have

(3.8)
$\frac{\|\boldsymbol{\mu}-\boldsymbol{\mu}_{n}^{*}\|}{\|\boldsymbol{\mu}_{n}^{*}\|}\leq\frac{16}{\Lambda}\max_{i:{\boldsymbol{\sigma}^{*}_{h}}_{\\{i\\}}\neq
0}|\partial_{i}E_{n}(\boldsymbol{\mu})|.$

###### Proof.

Using equations (2.19), (3.5), and (3.6), we obtain

(3.9)
$\operatorname{Hess}F_{n}(\boldsymbol{\sigma}_{n}^{*})\succeq\min\left\\{4\lambda\min_{i:{\boldsymbol{\mu}_{n}^{*}}_{\\{i\\}}\neq
0}{\boldsymbol{\mu}_{n}^{*}}_{\\{i\\}}\;,\;2\min_{i:{\boldsymbol{\mu}_{n}^{*}}_{\\{i\\}}=0}\partial_{i}E_{n}(\boldsymbol{\mu}_{n}^{*})\right\\}\mathbb{I}.$

Thus we can pick $0<R<\|\boldsymbol{\sigma}_{n}^{*}\|_{\infty}$ such that, for
any $\boldsymbol{\sigma}\in\mathbb{R}^{N}$, with
$|\boldsymbol{\sigma}_{n}^{*}-\boldsymbol{\sigma}|<R$, we have

(3.10)
$\operatorname{Hess}F_{n}(\boldsymbol{\sigma})\succeq\Lambda\mathbb{I}.$

Using the first order Taylor expansion of $\nabla F_{n}(\boldsymbol{\sigma})$
centered at $\boldsymbol{\sigma}_{n}^{*}$, we can write, for any
$\boldsymbol{\sigma}$ as above,

$\displaystyle\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{n}^{*}\|=\left\|[\operatorname{Hess}F_{n}(\boldsymbol{\eta})]^{-1}\nabla
F_{n}(\boldsymbol{\sigma})\right\|\leq\frac{1}{\Lambda}\|\nabla
F_{n}(\boldsymbol{\sigma})\|$ (3.11) $\displaystyle\leq$
$\displaystyle\frac{2\|\boldsymbol{\sigma}\|_{\infty}}{\Lambda}\max_{i:{\boldsymbol{\sigma}}_{\\{i\\}}\neq
0}|\partial_{i}E_{n}(\boldsymbol{\mu})|,$

where $\boldsymbol{\eta}\in[\boldsymbol{\sigma},\boldsymbol{\sigma}_{n}^{*}].$
Notice that

$\|\boldsymbol{\mu}-\boldsymbol{\mu}_{n}^{*}\|^{2}\leq\|\boldsymbol{\sigma}+\boldsymbol{\sigma}_{n}^{*}\|_{\infty}^{2}\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{n}^{*}\|^{2}\leq
4\max\left\\{\|\boldsymbol{\sigma}\|_{\infty}^{2}\;,\;\|\boldsymbol{\sigma}_{n}^{*}\|_{\infty}^{2}\right\\}\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{n}^{*}\|^{2},$

thus

(3.12) $\|\boldsymbol{\mu}-\boldsymbol{\mu}_{n}^{*}\|\leq
2(R+\|\boldsymbol{\sigma}_{n}^{*}\|_{\infty})\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{n}^{*}\|\leq
4\|\boldsymbol{\sigma}_{n}^{*}\|_{\infty}\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{n}^{*}\|.$

Therefore

$\displaystyle\|\boldsymbol{\mu}-\boldsymbol{\mu}_{n}^{*}\|\leq\frac{8\|\boldsymbol{\sigma}_{n}^{*}\|_{\infty}\|\boldsymbol{\sigma}\|_{\infty}}{\Lambda}\max_{i:{\boldsymbol{\sigma}}_{\\{i\\}}\neq
0}|\partial_{i}E_{n}(\boldsymbol{\mu})|$ $\displaystyle\leq$
$\displaystyle\frac{16\|\boldsymbol{\sigma}_{n}^{*}\|_{\infty}^{2}}{\Lambda}\max_{i:{\boldsymbol{\sigma}}_{\\{i\\}}\neq
0}|\partial_{i}E_{n}(\boldsymbol{\mu})|\leq\frac{16\|\boldsymbol{\mu}_{n}^{*}\|}{\Lambda}\max_{i:{\boldsymbol{\sigma}}_{\\{i\\}}\neq
0}|\partial_{i}E_{n}(\boldsymbol{\mu})|.$

∎

Note that, as a byproduct of this proof, we can state the following.

######  Corollary 3.5.

Hypothesis (3.3) holds true at a given
$\boldsymbol{\mu}_{n}^{*}\in\mathbb{R}_{\geq 0}^{N}$ for some $\lambda>0$ if
and only if $\operatorname{Hess}F_{n}(\boldsymbol{\sigma}_{n}^{*})$ is
strictly positive definite at any $\boldsymbol{\sigma}_{n}^{*}$ such that
${\boldsymbol{\sigma}_{n}^{*}}^{\boldsymbol{2}}=\boldsymbol{\mu}_{n}^{*}.$

##  4\. Existence, uniqueness and long time behavior of minimizing flows for
$E_{n}$

In order to minimize the functional $E_{n}$ we consider two different gradient
flows. The more natural choice is to design a gradient flow for $E_{n}$
including explicitly the non-negativity constraint
$\boldsymbol{\mu}\in\mathcal{M}_{n}^{+}$. It turns out that the latter
constraint makes the flow trajectories possibly non-smooth. The second
gradient flow we propose is applied directly to $F_{n}$ yelidng smooth
trajectories and ideal characteristics for the design of a more efficient
numerical scheme. For this reason this will be our method of choice described
in the next section.

###  4.1. The Euclidean gradient flow of $E_{n}$

Consider the following:

(4.1)
$\begin{cases}\frac{d}{dt}\boldsymbol{\mu}=-\partial^{\circ}E_{n}(\boldsymbol{\mu})&t\in[0,+\infty[\\\
\boldsymbol{\mu}(0)=\boldsymbol{\mu}^{0}\in\mathbb{R}_{\geq
0}^{N}&\end{cases}\;,$

where $\partial^{\circ}E_{n}(\boldsymbol{\mu})$ is the (unique) element of
minimal norm of the convex subdifferential of $E_{n}$ at $\boldsymbol{\mu}.$
We show next that every trajectory of this flow converges to a minimizer of
$E_{n}$, as $t\to+\infty$, regardless of the choice of
$\boldsymbol{\mu}^{0}\in\mathbb{R}_{\geq 0}^{N}.$

######  Theorem 4.1 (Existence,uniqueness and long time behavior of the flow).

Given $\boldsymbol{\mu}_{n}^{0}\in\mathbb{R}_{\geq 0}^{N}$ with
$\partial^{\circ}E_{n}(\boldsymbol{\mu}^{0})\neq 0$, there exists a unique
locally absolutely continuous and a.e. diffferentiable curve $[0,+\infty[\ni
t\rightarrow\boldsymbol{\mu}_{n}(t;\boldsymbol{\mu}_{n}^{0})\in\mathbb{R}_{\geq
0}^{N}$, with
$\boldsymbol{\mu}_{n}(0;\boldsymbol{\mu}_{n}^{0})=\boldsymbol{\mu}_{n}^{0}$,
that is an energy solution of (4.1), i.e.,

(4.2)
$E_{n}(\boldsymbol{\mu}_{n}(t;\boldsymbol{\mu}_{n}^{0}))=E_{n}(\boldsymbol{\mu}_{n}^{0})-\int_{0}^{t}\|\partial^{\circ}E_{n}\|^{2}(\boldsymbol{\mu}_{n}(s;\boldsymbol{\mu}_{n}^{0}))\,ds,\;\forall
t>0.$

In particular
$\boldsymbol{\mu}_{n}^{\prime}(t;\boldsymbol{\mu}_{n}^{0})=-\partial^{\circ}E_{n}(\boldsymbol{\mu}_{n}(t;\boldsymbol{\mu}_{n}^{0}))$
for almost any $t>0.$

Moreover, for any $\boldsymbol{\mu}^{0}\in\mathcal{M}_{n}^{+}$ with
$\partial^{\circ}E_{n}(\boldsymbol{\mu}^{0})\neq 0$, we have

(4.3)
$\lim_{t\to+\infty}\boldsymbol{\mu}_{n}(t;\boldsymbol{\mu}_{n}^{0})\in\operatorname*{argmin}_{\mathbb{R}_{\geq
0}^{N}}E_{n}(\boldsymbol{\mu}_{n}^{*}).$

###### Proof.

The proof of existence and uniqueness of the flow (together with the local
absolute continuity of the trajectories) relies on the application of [2, Th.
4.0.4] and [2, Prop. 1.4.1]. We first need to verify that the assumptions of
[2, Th. 4.0.4] are satisfied, i.e., Namely:

  1. (i)

$(\mathbb{R}_{\geq 0}^{N},d)$, with
$d(\boldsymbol{\mu},\boldsymbol{\nu}):=\|\boldsymbol{\mu}-\boldsymbol{\nu}\|_{2}$,
is a complete metric space.

  2. (ii)

$E_{n}$ is proper, lower semicontinuous, and bounded from below.

  3. (iii)

For any
$\boldsymbol{\mu},\boldsymbol{\mu}^{0},\boldsymbol{\mu}^{1}\in\mathbb{R}_{\geq
0}^{N}$ there exists a curve $\gamma:[0,1]\rightarrow\mathbb{R}_{\geq 0}^{N}$,
with $\gamma(0)=\boldsymbol{\mu}^{0},$ $\gamma(1)=\boldsymbol{\mu}^{1}$, such
that, for any $\tau\in(0+\infty)$,

(4.4)
$\Phi(\tau,\boldsymbol{\mu};\gamma(t))\leq(1-t)\Phi(\tau,\boldsymbol{\mu};\boldsymbol{\mu}^{0})+t\Phi(\tau,\boldsymbol{\mu};\boldsymbol{\mu}^{1})-\frac{1}{2\tau}t(1-t)d^{2}(\boldsymbol{\mu},\boldsymbol{\mu}^{1}),\;\forall
t\in(0,1),$

where we use the notation

$\Phi(\boldsymbol{u},\tau;\boldsymbol{v}):=E_{n}(\boldsymbol{v})+\frac{d^{2}(\boldsymbol{u},\boldsymbol{v})}{2\tau}.$

The property (_i_) hods because the positive cone $\mathbb{R}_{\geq 0}^{N}$ is
a closed set of the complete metric space $(\mathbb{R}^{N},d)$ and thus it is
complete. Property (_ii_) follows directly from Proposition 2.3 and the fact
that $E_{n}$ is non-negative. Finally (4.4) can be verified on line
$\gamma(t):=(1-t)\boldsymbol{\mu}^{0}+t\boldsymbol{\mu}^{1}$ employing the
convexity of $E_{n}$ and the $2$-convexity of $d^{2}$ along straight lines,
namely

$d^{2}(\boldsymbol{\mu},\gamma(t))\leq(1-t)d^{2}(\boldsymbol{\mu},\boldsymbol{\mu}^{0})+td^{2}(\boldsymbol{\mu},\boldsymbol{\mu}^{1})-t(1-t)d^{2}(\boldsymbol{\mu}^{0},\boldsymbol{\mu}^{1}).$

By [2, Th. 4.0.4 (_ii_)], there exists a locally Lipschitz curve of maximal
slope for $E_{n}$ with respect to the metric slope of $E_{n}.$ In our rather
simple setting the metric slope corresponds to the norm of the minimal
subdifferential (i.e., $\|\partial^{\circ}E_{n}\|_{2}$) and the definition of
a curve of maximal slope simplifies to

(4.5)
$|\boldsymbol{\mu}_{n}^{\prime}(t;\boldsymbol{\mu}_{n}^{0})|^{2}=\|\partial^{\circ}E_{n}\|_{2}^{2}(\boldsymbol{\mu}_{n}(t;\boldsymbol{\mu}_{n}^{0}))\\\
=-\frac{d}{dt}\left(E_{n}(\boldsymbol{\mu}_{n}(t;\boldsymbol{\mu}_{n}^{0}))\right),\;\text{
a.e. in }[0,+\infty[.$

Since $\mathbb{R}^{N}$ is a Hilbert space, (4.5) implies (see [2, Prop. 1.4.1
and Cor. 1.4.2]) that
$\boldsymbol{\mu}_{n}^{\prime}(t;\boldsymbol{\mu}_{n}^{0})=-\partial^{\circ}E_{n}(\boldsymbol{\mu}_{n}(t;\boldsymbol{\mu}_{n}^{0}))$
for almost any $t>0.$

We are left to prove that
$\lim_{t\to+\infty}\boldsymbol{\mu}_{n}(t;\boldsymbol{\mu}_{n}^{0})$ exists
and is a minimizer of $E_{n}$ regardless of the choice of the initial data. To
this aim we apply [2, Cor. 4.0.6] and, using again the fact that the metric
slope of $E_{n}$ at $\boldsymbol{\mu}_{n}$ is equal to
$\|\partial^{\circ}E_{n}(\boldsymbol{\mu}_{n})\|_{2}$, we obtain

(4.6) $\displaystyle
E_{n}(\boldsymbol{\mu}_{n}(t;\boldsymbol{\mu}_{n}^{0}))-\min_{\boldsymbol{\mu}_{n}\in\mathbb{R}_{\geq
0}^{N}}E_{n}(\boldsymbol{\mu}_{n})\leq\min_{\boldsymbol{\mu}_{n}\in\operatorname*{argmin}E_{n}}\frac{\|\boldsymbol{\mu}_{n}-\boldsymbol{\mu}_{n}^{0}\|_{2}^{2}}{2t}\to
0,$ (4.7)
$\displaystyle\|\boldsymbol{\mu}_{n}^{\prime}(t;\boldsymbol{\mu}_{n}^{0})\|_{2}=\|\partial^{\circ}E_{n}(\boldsymbol{\mu}_{n}(t;\boldsymbol{\mu}_{n}^{0}))\|_{2}\leq\min_{\boldsymbol{\mu}_{n}\in\operatorname*{argmin}E_{n}}\frac{\|\boldsymbol{\mu}_{n}-\boldsymbol{\mu}_{n}^{0}\|_{2}}{t}\to
0.$

Note that (4.6) implies that, if the curve
$\boldsymbol{\mu}_{n}(t;\boldsymbol{\mu}_{n}^{0})$ admits a limit point
$\boldsymbol{\mu}_{n}\in\mathbb{R}_{\geq 0}^{N}$ for $t\to+\infty$, then
$\boldsymbol{\mu}_{n}\in\operatorname*{argmin}E_{n}.$ Such a limit indeed
exists due to eq. 4.7. ∎

###  4.2. The Euclidean gradient flow of $F_{n}$

We recall here the definition of the coordinate square map
$\cdot^{\boldsymbol{2}}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ introduced in
(2.15) of Sec. 2.1:

$\boldsymbol{\sigma}^{\boldsymbol{2}}:=({\boldsymbol{\sigma}}_{\\{1\\}}^{2},{\boldsymbol{\sigma}}_{\\{2\\}}^{2},\dots,{\boldsymbol{\sigma}}_{\\{N\\}}^{2}),$

and the $F_{n}$ functional:

$F_{n}(\boldsymbol{\sigma}):=E_{n}(\boldsymbol{\sigma^{2}}).$

The minimization of $F_{n}$ on $\mathbb{R}^{N}$ has at least two advantages
with respect to the minimization of $E_{n}$ over $\mathbb{R}_{\geq 0}^{N},$
namely it is an unconstrained optimization problem and the objective function
is globally real analytic. On the other hand a number of potetial drawbacks
can be listed. The objective function is not convex anymore and, by the
symmetry of 2, uniqueness of the minimizer is lost.

We attack the problem of minimizing $F_{n}$ considering the following gradient
flow:

(4.8) $\begin{cases}\frac{d}{dt}\boldsymbol{\sigma}_{n}=-\nabla
F_{n}(\boldsymbol{\sigma}_{n}),&\forall t>0\\\
\boldsymbol{\sigma}_{n}(0)=\boldsymbol{\sigma}_{n}^{0}\end{cases}\;.$

We first prove existence, uniqueness, and regularity of the solutions of
(4.8), then we prove that all the trajectories converge to a global minimizer
of $F_{n}$ in Theorem 4.3.

######  Theorem 4.2 (Global existence, uniqueness, and regularity for (4.8)).

Let $\boldsymbol{\sigma}_{n}^{0}\in\mathbb{R}^{N}.$ There exists a unique
real-analytic classical solution $]0,+\infty[\ni
t\mapsto\boldsymbol{\sigma}_{n}(t;\boldsymbol{\sigma}_{n}^{0})\in\mathbb{R}^{N}$
of (4.8).

###### Proof.

For any $c\in\mathbb{R}$ and any $\epsilon>0$ we consider the compact sets

$S_{c}:=\\{\boldsymbol{\sigma}\in\mathbb{R}^{N}:F_{n}(\boldsymbol{\sigma})\leq
c^{2}\\},\;\;\;S_{c}^{\epsilon}:=\cup_{\boldsymbol{\sigma}\in
S_{c}}B(\boldsymbol{\sigma},\epsilon].$

Note that $\max_{\boldsymbol{\sigma}\in S_{c}}\|\boldsymbol{\sigma}\|\leq|c|$
and $\max_{\boldsymbol{\sigma}\in
S_{c}^{\epsilon}}\|\boldsymbol{\sigma}\|\leq|c|+\epsilon$, since
$F_{n}\geq\|\cdot\|^{2}$ by definition. Let us pick $c\in\mathbb{R}$ such that
$\boldsymbol{\sigma}_{n}^{0}\in S_{c}.$ It is clear that

$\displaystyle\max_{\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{n}^{0}\|\leq\epsilon}\|\nabla
F_{n}(\boldsymbol{\sigma})\|\leq\max_{\boldsymbol{\sigma}\in
S_{c}^{\epsilon}}\|\nabla F_{n}(\boldsymbol{\sigma})\|$ $\displaystyle\leq$
$\displaystyle 2\max_{\boldsymbol{\sigma}\in
S_{c}^{\epsilon}}\left(\|\boldsymbol{\sigma}\|_{2}\max_{i}\left||T_{i}^{n}|-(\boldsymbol{f}_{n})^{t}(A^{(n,n)}({\boldsymbol{\sigma}}^{\boldsymbol{2}}))^{-1}A^{(i,n,n)}(A^{(n,n)}({\boldsymbol{\sigma}}^{\boldsymbol{2}}))^{-1}\boldsymbol{f}_{n}\right|\right)$
$\displaystyle\leq$ $\displaystyle 2\max_{\boldsymbol{\sigma}\in
S_{c}^{\epsilon}}\|\boldsymbol{\sigma}\|_{2}\max_{i}\left(|T_{i}^{n}|+\frac{C^{2}\lambda_{max}(A^{(i,n,n)})}{\delta_{n}^{2}}\|\boldsymbol{f}_{n}\|^{2}\right)$
$\displaystyle\leq$ $\displaystyle 2\max_{\boldsymbol{\sigma}\in
S_{c}^{\epsilon}}\|\boldsymbol{\sigma}\|_{2}\max_{i}\left(|T_{i}^{n}|+\frac{C^{2}\lambda_{max}(A^{(i,n,n)})}{\delta_{n}^{2}}\|\boldsymbol{f}_{n}\|_{2}^{2}\right)$
$\displaystyle\leq$
$\displaystyle(|c|+\epsilon)\max_{i}\left(|T_{i}^{h}|+\frac{C^{2}\lambda_{max}(A^{(i,n,n)})}{\delta_{n}^{2}}\|\boldsymbol{f}_{n}\|_{2}^{2}\right)$
$\displaystyle=:$ $\displaystyle(|c|+\epsilon)K_{n},$

Where we denoted by $C$ the Poincare’ constant of $\Omega$ and we used that
$A^{(n,n)}({\boldsymbol{\sigma}}^{\boldsymbol{2}})$ is symmetric positive
definite and its minimal eigenvalue is at least $\delta_{n}/C$ for any
$\boldsymbol{\mu}\in\mathbb{R}_{+}^{N}.$ Moreover $A^{(i,n,n)}$ is symmetric
and positive semi-definite for any $i=\,2,\dots,N.$

Since $\nabla F_{n}$ is a real analytic function on $\mathbb{R}^{N}$, it is in
particular a uniformly Lipschitz function on $S_{c}^{\epsilon}\supseteq
B]\boldsymbol{\sigma}_{n}^{0},\epsilon].$ By the classical Picard-Lindelöf
Theorem, there exists a unique solution
$\boldsymbol{\sigma}_{n}^{(0)}(t;\boldsymbol{\sigma}_{n}^{0})$ of (4.8) in the
interval $[-\tau,\tau]$ with

$\tau:=\frac{\epsilon}{(|c|+\epsilon)K_{n}}.$

Let us introduce the sequence

$t_{j}:=j\frac{2}{3}\tau,\;\;\;j=0,1,\dots$

and the notation
$\boldsymbol{\sigma}_{n}^{j+1}:=\boldsymbol{\sigma}_{n}^{(j)}(t_{j+1};\boldsymbol{\sigma}_{n}^{s})$
Now notice that

$\displaystyle F_{n}(\boldsymbol{\sigma}_{n}^{1})=$ $\displaystyle
F_{n}(\boldsymbol{\sigma}_{n}^{0})+\int_{0}^{t_{1}}\langle\nabla
F_{n}(\boldsymbol{\sigma}_{n}^{(0)}(s;\boldsymbol{\sigma}_{n}^{0})),\frac{d}{dt}\boldsymbol{\sigma}_{n}^{(0)}(s;\boldsymbol{\sigma}_{n}^{0}))\rangle
ds$ $\displaystyle=$ $\displaystyle
F_{n}(\boldsymbol{\sigma}_{n}^{0})-\int_{0}^{t_{1}}\|\nabla
F_{n}(\boldsymbol{\sigma}_{n}^{(0)}(s;\boldsymbol{\sigma}_{n}^{0}))\|^{2}ds\leq
c^{2}.$

Therefore $\boldsymbol{\sigma}_{n}^{1}\in S_{c}.$ Thus we can repeat the
argument to construct a curve
$\boldsymbol{\sigma}_{n}^{1}(\cdot;\boldsymbol{\sigma}_{n}^{1}):[t^{1}-\tau,t^{1}+\tau]\rightarrow\mathbb{R}^{N}$
which is the uique solution of the Cauchy problem

$\begin{cases}\dfrac{d}{dt}\boldsymbol{\sigma}=-\nabla
F_{n}(\boldsymbol{\sigma})&t>t^{1}\\\
\boldsymbol{\sigma}(t^{1})=\boldsymbol{\sigma}_{n}^{1}\end{cases}.$

Uniqueness implies that
$\boldsymbol{\sigma}_{n}^{(0)}(t;\boldsymbol{\sigma}_{n}^{0})=\boldsymbol{\sigma}_{n}^{(1)}(t-t^{1};\boldsymbol{\sigma}_{n}^{1})$
for any $t\in[-1/3\tau,\tau]$. This provides existence and uniqueness of the
solution of (4.8) in $[-\tau,\frac{5}{3}\tau].$ Iterating the same argument we
can construct a solution
$\boldsymbol{\sigma}_{n}(t;\boldsymbol{\sigma}_{n}^{0})$ of (4.8) in
$[0,+\infty[$ that turns out to be unique. Analyticity of the flow follows
from the Cauchy Kovalevskaya Theorem and Proposition 2.3. ∎

######  Theorem 4.3 (Long-time asymptotics of (4.8)).

Let $\boldsymbol{\sigma}_{n}^{0}\in\mathbb{R}^{N}$ be such that $\nabla
F_{n}(\boldsymbol{\sigma}_{n}^{0})\neq 0$. Then there exists
$\boldsymbol{\mu}_{n}^{*}\in\operatorname*{argmin}_{\mathbb{R}_{\geq
0}^{N}}E_{n}$ such that

(4.9)
$\lim_{t\to+\infty}(\boldsymbol{\sigma}_{n}(t;\boldsymbol{\sigma}_{n}^{0}))^{\boldsymbol{2}}=\boldsymbol{\mu}_{n}^{*}.$

###### Proof.

For a given $\boldsymbol{\sigma}_{n}^{0}$ we denote by
$\boldsymbol{\sigma}_{n}(t):=\boldsymbol{\sigma}_{n}(t;\boldsymbol{\sigma}_{n}^{0})$
the curve that solves (4.8). By the Fundamental Theorem of Integral Calculus
we can write, for any $t>0$,

$\min F_{n}-c^{2}\leq
F_{n}(\boldsymbol{\sigma}_{n}(t))-F_{n}(\boldsymbol{\sigma}_{n}^{0})\\\
=\int_{0}^{t}\langle\nabla
F_{n}(\boldsymbol{\sigma}_{n}(s));\boldsymbol{\sigma}_{n}^{\prime}(s)\rangle
ds=-\int_{0}^{t}\|\boldsymbol{\sigma}_{n}^{\prime}(s)\|^{2}ds.$

Thus we have

$\int_{0}^{+\infty}\|\boldsymbol{\sigma}_{n}^{\prime}(t)\|^{2}dt\leq
c^{2}-\min F_{n}<+\infty.$

In particular this implies that for any $\epsilon>0$ there exists
$T_{\epsilon}>0$ such that
$\int_{T_{\epsilon}}^{+\infty}\|\boldsymbol{\sigma}_{n}^{\prime}(t)\|^{2}dt\leq\epsilon$,
and that any sequence $\\{\boldsymbol{\sigma}_{n}(t_{k})\\}_{k\in\mathbb{N}}$,
with $t_{k}\to+\infty$, is a Cauchy sequence. Hence we can define

$\hat{\boldsymbol{\sigma}}_{n}:=\lim_{t\to+\infty}\boldsymbol{\sigma}_{n}(t)=\boldsymbol{\sigma}_{n}^{0}+\int_{0}^{+\infty}\boldsymbol{\sigma}_{n}^{\prime}(t)dt.$

Now we want to show that
$\lim_{t\to+\infty}\boldsymbol{\sigma}_{n}^{\prime}(t)=0$. Indeed thanks to
the computations above (see the Proof of Theorem 4.2) we can write

$\displaystyle\sup_{t>0}\|\boldsymbol{\sigma}_{n}^{\prime}\|$
$\displaystyle=\max_{\boldsymbol{\sigma}\in S_{c}}\|\nabla
F_{n}(\boldsymbol{\sigma})\|\leq|c|K_{n}$
$\displaystyle\sup_{t>0}\|\boldsymbol{\sigma}_{n}^{\prime\prime}\|^{2}$
$\displaystyle=\sup_{t>0}\boldsymbol{\sigma}_{n}^{\prime}(t)^{t}\operatorname{Hess}F_{n}(\boldsymbol{\sigma}_{n}(t))\boldsymbol{\sigma}_{n}^{\prime}(t)\leq
R^{2}\sup_{t>0}\|\boldsymbol{\sigma}_{n}^{\prime}(t)\|^{2},$

where

$R^{2}=\sup_{\boldsymbol{\sigma}\in
S_{c}}\rho(\operatorname{Hess}F_{n}(\boldsymbol{\sigma}))<\infty,$

and $\rho(\cdot)$ denotes the spectral radius. It follows that
$\sup_{t>0}\|\boldsymbol{\sigma}_{n}^{\prime\prime}(t)\|\leq R|c|K_{n}.$ Let
us assume by contradiction that
$\limsup_{t\to+\infty}\|\boldsymbol{\sigma}_{n}^{\prime}(t)\|^{2}=l>0.$ Then
we can pick $\\{t_{k}\\}$ such that $\inf_{k}|t_{k}-t_{k+1}|=:\alpha>0$ and
$s_{k}:=\|\boldsymbol{\sigma}_{n}^{\prime}(t_{k})\|^{2}\to l,\;\;\text{ as
}k\to+\infty.$ Notice that for any $t>0$ such that

$|t-t_{k}|\leq\inf_{j}\\{\sqrt{s_{j}}/(4RcK_{n}),\alpha\\},\;k=1,2,\dots,$

then

$\displaystyle\|\boldsymbol{\sigma}_{n}^{\prime}(t)\|^{2}=$
$\displaystyle\|\boldsymbol{\sigma}_{n}^{\prime}(t_{k})+(t-t_{k})\boldsymbol{\sigma}_{n}^{\prime\prime}(\xi)\|^{2}$
$\displaystyle=$
$\displaystyle\|\boldsymbol{\sigma}_{n}^{\prime}(t_{k})\|^{2}+2(t-t_{k})\langle\boldsymbol{\sigma}_{n}^{\prime\prime}(\xi),\boldsymbol{\sigma}_{n}^{\prime}(t_{k})\rangle+(t-t_{k})^{2}\|\boldsymbol{\sigma}_{n}^{\prime\prime}(\xi)\|^{2}$
$\displaystyle\geq$ $\displaystyle
s_{k}-2|t-t_{k}|\sup_{s>0}\|\boldsymbol{\sigma}_{n}^{\prime\prime}(s)\|\sqrt{s_{k}}$
$\displaystyle\geq$ $\displaystyle
s_{k}-2|t-t_{k}|R|c|K_{n}\sqrt{s_{k}}\geq\frac{s_{k}}{2}\;.$

Therefore

$\displaystyle+\infty$
$\displaystyle>\int_{0}^{+\infty}\|\boldsymbol{\sigma}_{n}^{\prime}(s)\|^{2}ds\geq\sum_{k=1}^{+\infty}\int_{t_{k}-\alpha}^{t_{k}+\alpha}\|\boldsymbol{\sigma}_{n}^{\prime}(s)\|^{2}ds\geq\frac{\alpha}{2}\sum_{k=1}^{+\infty}s_{k}\;,$

contradicting the hypothesis (i.e., $\lim_{k}s_{k}=l>0$). Thus
$\limsup_{t\to+\infty}\|\boldsymbol{\sigma}_{n}^{\prime}(t)\|=0$ and, since
$\|\boldsymbol{\sigma}_{n}^{\prime}(t)\|\geq 0$ $\forall t>0$,
$\lim_{t\to+\infty}\|\boldsymbol{\sigma}_{n}^{\prime}(t)\|=0.$ In particular,
$\hat{\boldsymbol{\sigma}}_{n}$ is a stationary point for $F_{n}.$

We observe that, by the Karush-Kuhn-Tucker (KKT) Theorem and the convexity of
$E_{n}$ the following conditions are necessary and sufficient for
$\boldsymbol{\mu}_{n}^{*}\in\mathbb{R}_{\geq 0}^{N}$ to be a minimizer of
$E_{n}$ on $\mathbb{R}_{\geq 0}^{N}$

(4.10) $\begin{cases}\partial_{i}E_{n}(\boldsymbol{\mu}_{n}^{*})=0&\forall
i:{\boldsymbol{\mu}_{n}}_{\\{i\\}}^{*}>0\\\
\partial_{i}E_{n}(\boldsymbol{\mu}_{n}^{*})\geq 0&\forall
i:{\boldsymbol{\mu}_{n}}_{\\{i\\}}^{*}=0\end{cases}\;.$

Since

(4.11) $0=\nabla
F_{n}(\hat{\boldsymbol{\sigma}}_{n})={2(\hat{\boldsymbol{\sigma}}_{n})}_{\\{1\\}}\partial_{1}E_{n}({\hat{\boldsymbol{\sigma}}_{n}}^{\boldsymbol{2}}),\dots,2(\hat{\boldsymbol{\sigma}}_{n})_{N}\partial_{N}E_{n}({\hat{\boldsymbol{\sigma}}_{n}}^{\boldsymbol{2}})),$

the first equation in (4.10) is always satisfied by
$\boldsymbol{\mu}_{n}^{*}:={\hat{\boldsymbol{\sigma}}_{n}}^{\boldsymbol{2}}.$

It is not difficult to prove (using an argument similar to the one used above
while proving $\lim_{t\to+\infty}\|\boldsymbol{\sigma}_{n}^{\prime}(t)\|=0$)
that we can pick a sequence $0<t_{k}\to+\infty$ such that
${\boldsymbol{\sigma}_{n}}_{\\{i\\}}^{\prime}(t_{k})\cdot{\boldsymbol{\sigma}_{n}}_{\\{i\\}}(t_{k})\leq
0$ for any such $i$. Thus, using the definition of the gardient flow,

$\begin{cases}{\boldsymbol{\sigma}_{n}}_{\\{i\\}}^{\prime}(t_{k})\leq 0&\text{
if }{\boldsymbol{\sigma}_{n}}_{\\{i\\}}(t_{k})\geq 0\\\
{\boldsymbol{\sigma}_{n}}_{\\{i\\}}^{\prime}(t_{k})\geq 0&\text{ if
}{\boldsymbol{\sigma}_{n}}_{\\{i\\}}(t_{k})\leq
0\end{cases}\;\Longrightarrow\;\\\
\begin{cases}-2{\boldsymbol{\sigma}_{n}}_{\\{i\\}}(t_{k})\partial_{i}E_{n}(\boldsymbol{\sigma}_{n}(t_{k}))\leq
0&\text{ if }{\boldsymbol{\sigma}_{n}}_{\\{i\\}}(t_{k})\geq 0\\\
-2{\boldsymbol{\sigma}_{n}}_{\\{i\\}}(t_{k})\partial_{i}E_{n}(\boldsymbol{\sigma}_{n}(t_{k}))\geq
0&\text{ if }{\boldsymbol{\sigma}_{n}}_{\\{i\\}}(t_{k})\leq 0\end{cases},$

$\forall k$ and for any $i$ such that
$\hat{\boldsymbol{\sigma}}_{n,\\{i\\}}=0.$ Hence
$\partial_{i}E_{n}(\boldsymbol{\sigma}_{n}(t_{k}))\geq 0$ $\forall k$ and for
any such $i$, i.e., also the second equation of (4.10) holds true. ∎

Theorem 4.3 is a qualitative result, but we can give a sharp quantitative
estimate of the rate of convergence exploiting the real analyticity of
$F_{n}$. The key element here is the _Łojasiewicz Theorem_ (see [25]).
Łojasiewicz’s. Precisely, if $D\subseteq\mathbb{R}^{N}$ is an open set and
$f:D\rightarrow\mathbb{R}$ is real analytic, then, for any $x\in D$, there
exists $R>0$, a real number $\vartheta\in(0,1/2]$ (called the _Łojasiewicz
exponent_ of $f$ at $x$), and $L>0$ such that

(4.12) $|f(x)-f(y)|^{1-\vartheta}\leq L\|\nabla f(y)\|,\;\;\forall y\in
B(x,R).$

The above is the so-called Łojasiewicz Inequality, which has already been used
whithin similar contexts for proving the convergence of finite and infinite
dimensional gradient flows, see e.g. [8], and for estimating the rate of
convergence of the flow to a critical point of the objective, see [12, 22] and
references therein. Inequality (4.12) has been successfully used also for
proving the convergence of the time discretization of gradient flows by e.g.
backward Euler schemes, see [26]. We can specialize this quite classical tool
to our setting using the results of [21] as follows.

######  Theorem 4.4 (Rate of convergence to equilibria of (4.8), Thm. 2.2 of
[21]).

Let
$\boldsymbol{\sigma}_{n}(t):=\boldsymbol{\sigma}_{n}(t;\boldsymbol{\sigma}_{n}^{0})$
be a trajectory of (4.8) and
$\boldsymbol{\sigma}_{n}^{*}:=\lim_{t\to+\infty}\boldsymbol{\sigma}_{n}(t)\in\operatorname*{argmin}F_{n}$.
Then:

  1. i)

If the Łojasiewicz Inequality (4.12) holds for $F_{n}$ at
$\boldsymbol{\sigma}_{n}^{*}$ with $\vartheta<1/2$, then there exists
$t^{0}\geq 0$ such that, for any $t>t^{0}$, we have

(4.13)
$\|\boldsymbol{\sigma}_{n}(t)-\boldsymbol{\sigma}_{n}^{*}\|\leq\left(\|\boldsymbol{\sigma}_{n}(t^{0})-\boldsymbol{\sigma}_{n}^{*}\|^{\frac{2\vartheta-1}{\vartheta}}+\vartheta^{\frac{1-2\vartheta}{\vartheta}}\frac{1-2\vartheta}{L^{\frac{1}{\vartheta}}}(t-t^{0})\right)^{\frac{\vartheta}{2\vartheta-1}}.$

  2. ii)

If the Łojasiewicz Inequality (4.12) holds for $F_{n}$ at
$\boldsymbol{\sigma}_{n}^{*}$ with $\vartheta=1/2$ (and in particular if
Hypothesis 3.3 holds), then there exists $t^{0}\geq 0$ such that, for any
$t>t^{0}$, we have

(4.14)
$\|\boldsymbol{\sigma}_{n}(t)-\boldsymbol{\sigma}_{n}^{*}\|\leq\|\boldsymbol{\sigma}_{n}(t^{0})-\boldsymbol{\sigma}_{n}^{*}\|\exp\left(-\frac{1}{2L^{2}}(t-t^{0})\right).$

###  4.3. Refined conditioning estimates by Łojasiewicz Inequality

In this subsection we exploit the existence and regularity of the gradient
flow of $F_{n}$ together with the Łojasiewicz Inequality to derive a
conditioning estimate similar to (3.8) (obtained in Proposition 3.4) valid
also when Hypothesis 3.3 does not hold. This result will be used later in the
design of the stopping criterion for our main minimization algorithm 2.

In the next proposition we use the following notation for the Euclidean
distance of the point $\boldsymbol{\sigma}$ from the set minimizers of
$F_{n}$:

$d(\boldsymbol{\sigma},\operatorname*{argmin}F_{n})=\min\\{\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{n}^{*}\|,\;\boldsymbol{\sigma}_{n}^{*}\in\operatorname*{argmin}F_{n}\\}.$

Note that the minimum is well defined because $\operatorname*{argmin}F_{n}$ is
compact due to the coercivity and lower semicontinuity of $F_{n}$.

######  Proposition 4.5 (Refined conditioning estimates).

For any $h\in\mathbb{N}$ there exists $R>0$, $L>0$, and $\vartheta\in(0,1/2]$,
such that, for any $\boldsymbol{\sigma}$ such that
$d(\boldsymbol{\sigma},\operatorname*{argmin}F_{n})<R$, there exists
$\boldsymbol{\sigma}_{n}^{*}\in\operatorname*{argmin}F_{n}$ with

(4.15)
$\|\boldsymbol{\sigma}_{n}-\boldsymbol{\sigma}_{n}^{*}\|<\frac{L^{\frac{1}{1-\vartheta}}}{\vartheta}\|\nabla
F_{n}(\boldsymbol{\sigma}_{n}(t))\|^{\frac{\vartheta}{1-\vartheta}}\,.$

If Hypothesis 3.3 holds, we can take $\vartheta=1/2$, i.e. we have

$\|\boldsymbol{\sigma}_{n}-\boldsymbol{\sigma}_{n}^{*}\|<2L^{2}\|\nabla
F_{n}(\boldsymbol{\sigma}_{n}(t))\|\,.$

###### Proof.

We notice that the constants in the Łojasiewicz Inequality can be chosen
uniformly with respect to $\boldsymbol{\sigma}_{n}^{*}$ because of the
aforementioned compactness of $\operatorname*{argmin}F_{n}.$ We construct the
trajectory $t\mapsto\boldsymbol{\sigma}_{n}(t)$ by solving the gradient flow
of $F_{n}$ starting at $\boldsymbol{\sigma}_{n}$. This is a well-defined
unique analytic curve due to Theorem 4.2 that has a limit point
$\boldsymbol{\sigma}_{n}^{*}\in\operatorname*{argmin}F_{n}$ due to Theorem
4.3. Then, using eq. (2.10) found in the proof of Thm. 2.2 of [21], we can
write

$\|\boldsymbol{\sigma}_{n}(t)-\boldsymbol{\sigma}_{n}^{*}\|^{2}\leq\frac{L^{\frac{1}{1-\vartheta}}}{\vartheta}\|\nabla
F_{n}(\boldsymbol{\sigma}_{n}(t))\|^{\frac{\vartheta}{1-\vartheta}}$

If Hypothesis 3.3 holds true, then the Hessian of $F_{n}$ at
$\boldsymbol{\sigma}_{n}^{*}$ is not degenerate. Due to e.g., [26, Prop. 2.2],
the Łojasiewicz Inequality holds for $F_{n}$ in a suitable neighbourhood of
$\boldsymbol{\sigma}_{n}^{*}$ with $\vartheta=1/2.$ ∎

######  Remark 4.6.

It is worth pointing out that Proposition 4.5 is a slightly stronger result
than the condition estimate (3.8) that we proved in Section 3, since the
latter depends on Hypothesis 3.3. On the contrary, Proposition 4.5 still leads
to the (weaker) estimate (4.15) even if Hypothesis 3.3 does not hold true.

##  5\. Minimization algorithms from numerical integration of the gradient
flows

Since the gradient flow of $\mathcal{E}_{n}$ converges to
$\boldsymbol{\mu}_{n}^{*}\in\operatorname*{argmin}\mathcal{E}_{n}$ as
$t\to+\infty$, we may try to numerically integrate the flow, in order to
design an approximation algorithm for $\boldsymbol{\mu}_{n}^{*}$. The most
straighforward idea is to adapt the forward Euler scheme to our constrained
setting by enforcing the positivity constraint at each iteration. When the
feasible set is a convex subset of an Hilbert space, one can modify such a
scheme including an orthogonal projection of the Euler update on the feasible
set. This algorithm (and its many variants) can be casted within the framework
of _proximal algorithms_ (see [29]) and is termed _projected gradient descent_
or _projected forward Euler scheme_ , and its application to our case is
described in Algorithm 1.

Algorithm 1 Projected forward Euler Scheme

Input $\boldsymbol{\mu}\in\mathbb{R}_{+}^{N}$, $\tau>0$,
$n_{step}\in\mathbb{N}$, $toll>0$

while $s<n_{step}$ and $err>toll$ do

$s=s+1$

Compute $\boldsymbol{u}=A^{-1}(\boldsymbol{\mu}+\delta_{n})\boldsymbol{f}$

Compute $\boldsymbol{v}=(\boldsymbol{u}^{t}A^{(j)}\boldsymbol{u}^{t})_{j}-1$

Set $v^{(j)}=(v^{(j)})^{+}$ if $\mu^{(j)}=0$

Set $\boldsymbol{\mu}:=(\boldsymbol{\mu}+\tau\boldsymbol{v})^{+}$

$err=\|\boldsymbol{v}\|_{2}$

end while

return $\boldsymbol{\mu}$

The convergence analysis of Algorithm 1 is rather standard. It can be proven
that, for sufficiently small and constant time steps, any cluster point of the
sequence computed by Algorithm 1 is a stationary point of the objective
functional, see for instance [7]. Unfortunately, the use of the orthogonal
projection in the definition of the iterates implies a non-monotone
convergence and may cause few stability issues, expecially if the objective
tends to be non-strongly convex.

For these reasons we pursue a different strategy by working with $F_{n}$, thus
avoiding the explicit implementation of the non-negativity constraint. This
entails the derivation of a numerical algorithm that mimics De Giorgi’s
_minimizing movements_ which have been introduced for the construction of
gradient flows in very general contexts, such as Banach or metric spaces, [2].
This approach applied to the construction of the gradient flow of the
functional $F_{n}$ consists in picking an initial point
$\boldsymbol{\sigma}_{n}^{0}$ and a time step $\tau$ (or a time step sequence
$\\{\tau^{\ell}\\}$), then computing iteratively
$\boldsymbol{\sigma}_{n}^{\ell+1}$ such that

(5.1)
$\boldsymbol{\sigma}_{n}^{\ell+1}\in\operatorname*{argmin}_{\boldsymbol{\sigma}\in\mathbb{R}^{N}}G_{n}(\boldsymbol{\sigma};\boldsymbol{\sigma}_{n}^{\ell},\tau):=\operatorname*{argmin}_{\boldsymbol{\sigma}\in\mathbb{R}^{N}}F_{n}(\boldsymbol{\sigma})+\frac{\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{n}^{\ell}\|^{2}}{2\tau}.$

Then, one shrinks the size of the time step (or the size of the partition in
the case of variable time stepping) to zero and tries to prove that the linear
interpolation of the discrete time trajectories is converging to a solution of
the gradient flow equation.

Here we are interested in approximating a minimizer of $F_{n}$ rather than
computing the trajectory of the flow. Therefore we are not going to shrink the
time step to zero, instead we keep iterating eq. (5.1) and consider
$\lim_{\ell}\boldsymbol{\sigma}_{n}^{\ell+1}.$ The derived optimization
technique is called _variational backward Euler method_ [26] and is still
casted by some authors (see e.g., [20]) in the framework of proximal
algorithms.

Note that equation (5.1) does not define yet a numerical algorithm, since one
needs to choose a method for solving the optimization problem and a stopping
criterion for $\ell$, possibly with error bounds. Starting from (5.1), the
proposed numerical algorithm is based on the following choices:

  * •

the variational problem (5.1) is replaced by the first-order optimality
condition

(5.2) $\nabla
G_{n}(\boldsymbol{\sigma}_{n}^{\ell+1};\boldsymbol{\sigma}_{n}^{\ell},\tau)=\nabla
F(\boldsymbol{\sigma}_{n}^{\ell+1})+\frac{\boldsymbol{\sigma}_{n}^{\ell+1}-\boldsymbol{\sigma}_{n}^{\ell}}{\tau}=0,$

which in principle may be not equivalent to eq. (5.1).

  * •

the solution of the non-linear problem (5.2) is approximated by Newton method
starting at $\boldsymbol{\sigma}_{n}^{\ell}$,

  * •

Newton method is stopped with a specific criterion based on the relative
residual and on the sign of the new iterate, quantities that are both easily
computed _a posteriori_ ,

  * •

the algorithm stops when $\|\nabla F(\boldsymbol{\sigma}_{n}^{\ell})\|$ is
smaller than a prescribed tollerance, as suggested by Proposition 4.5.

The final scheme is reported in Algorithm 2.

Algorithm 2 Backward Euler Scheme with Newton solver

Input $\boldsymbol{\sigma}_{n}^{0}\in\mathbb{R}^{N}$, $\tau>0$,
$n_{step}\in\mathbb{N}$, $toll>0$, $\epsilon>0$

Set $\ell:=0$

Compute $res=\|\nabla F_{n}(\boldsymbol{\sigma}_{n}^{0})\|$

if $res=0$ then

Exit with error.

end if

while $\ell<n_{step}$ and $res>toll$ do

Set $\ell=\ell+1$, $\boldsymbol{\sigma}^{new}:=\boldsymbol{\sigma}^{old}$

Compute $res_{Newton}:=\nabla
G_{n}(\boldsymbol{\sigma}^{new};\boldsymbol{\sigma}^{old},\tau)$

while
$|{res_{Newton}}_{\\{i\\}}|>\epsilon|{\boldsymbol{\sigma}^{new}-\boldsymbol{\sigma}^{old}}_{\\{i\\}}|$
for some $i$ or
$\operatorname{sign}\boldsymbol{\sigma}^{new}\neq\operatorname{sign}\boldsymbol{\sigma}^{old}$
do

Compute
$\boldsymbol{\sigma}^{new}=\boldsymbol{\sigma}^{new}-[\operatorname{Hess}G_{n}(\boldsymbol{\sigma}^{new};\boldsymbol{\sigma}^{old},\tau)]^{-1}\nabla
G_{n}(\boldsymbol{\sigma}^{new};\boldsymbol{\sigma}^{old},\tau)$

Compute $res_{Newton}:=\nabla
G_{n}(\boldsymbol{\sigma}^{new};\boldsymbol{\sigma}^{old},\tau)$

end while

Compute $res=\|\nabla F_{n}(\boldsymbol{\sigma}^{new})\|$

end while

return $\boldsymbol{\sigma}^{new}$

Passing from the study of the iterated solution of (5.1) to the study of
Algorithm 2 some issues come into play. First, since the function
$G_{n}(\cdot;\boldsymbol{\sigma}_{n}^{\ell},\tau)$ may be not globally convex,
the first order optimality condition (5.2), though necessary, is not in
general sufficient for (5.1) holding true. Moreover, equations (5.1) and (5.2)
may have more than one solution. Lastly, convergence of Newton method is
guaranteed only for a suitable choice of the initial guess. In Theorem 5.1 we
will prove that, for a sufficiently small $\tau$ depending only on
$\boldsymbol{\sigma}^{0}$, we can iterativelly compute
$\boldsymbol{\sigma}_{n}^{\ell+1}$ by solving (5.2) by Newton method (possibly
with an infinite number of steps) with $\boldsymbol{\sigma}_{n}^{\ell}$ as
initial guess. The resulting sequence
$\\{\boldsymbol{\sigma}_{n}^{\ell}\\}_{\ell\in\mathbb{N}}$ satisfies (5.1) at
each step and the sequence
$\\{{(\boldsymbol{\sigma}_{n}^{\ell})}^{\boldsymbol{2}}\\}_{\ell\in\mathbb{N}}$
converges to a global minimizer of $E_{n}.$ Note that the use of the Newton
method and the prescribed choice of the initial data provide a natural well-
defined sequence $\\{\boldsymbol{\sigma}_{n}^{\ell}\\}_{\ell\in\mathbb{N}}$
even when the solution of (5.1) is not unique.

Theorem 5.1 refers to an exact solution of (5.2), typically obtained by
infinite iterations. On the other hand, the propagation of the error
introduced when Newton method is stopped erlier must be taken into account in
the numerical analysis of the scheme. In Algorithm 2 we choose a specific
stopping criterion that allows us to provide a _stability estimate_ in terms
of the decrease of the functional $F_{n}$ along the computed discrete
trajectory, see Proposition 5.2. Finally, using these last results, we are
able to prove in Theorem 5.3 the _convergence_ of Algorithm 2 as
$n_{step}\to+\infty$ towards the solution of the discrete Monge-Kantorovich
equations (2.22).

######  Theorem 5.1 (Consistency of Algorithm 2).

Let $\min F_{n}<c<+\infty$. Then, there exists $\tau^{*}(c)>0$ such that, for
any $\boldsymbol{\sigma}_{n}^{0}$ with $F_{n}(\boldsymbol{\sigma}_{n}^{0})<c$
and $\nabla F_{n}(\boldsymbol{\sigma}_{n}^{0})\neq 0$, and for any
$0<\tau<\tau^{*}(c),$ the following sequences are well defined

$\displaystyle\boldsymbol{\sigma}_{n}^{\ell+1,0}:=\boldsymbol{\sigma}_{n}^{\ell}$
$\displaystyle\boldsymbol{\sigma}_{n}^{\ell+1,r+1}:=\boldsymbol{\sigma}_{n}^{\ell+1,r}-\left[\operatorname{Hess}G_{n}(\boldsymbol{\sigma}^{\ell+1,r};\boldsymbol{\sigma}_{n}^{\ell},\tau)\right]^{-1}\nabla
G_{n}(\boldsymbol{\sigma}^{\ell+1,r};\boldsymbol{\sigma}_{n}^{\ell},\tau)$
$\displaystyle\boldsymbol{\sigma}_{n}^{\ell+1}:=\lim_{r}\boldsymbol{\sigma}_{n}^{\ell+1,r}.$

For any $\ell\in\mathbb{N}$ we have

(5.3) $\nabla
F_{n}(\boldsymbol{\sigma}_{n}^{\ell+1})=-\frac{\boldsymbol{\sigma}_{n}^{\ell+1}-\boldsymbol{\sigma}_{n}^{\ell}}{\tau}.$

Moreover there exists $\boldsymbol{\mu}_{n}^{*}\in\operatorname*{argmin}E_{n}$
such that

(5.4)
$\lim_{\ell}{(\boldsymbol{\sigma}_{n}^{\ell})}^{\boldsymbol{2}}=\boldsymbol{\mu}_{n}^{*}.$

###### Proof.

The first part of the proof is a sharp version of the classical proof of the
local convergence of Newton method. Let $\Omega_{\ell}$ be the connected
component of the interior of the set $\\{F_{n}\leq
F_{n}(\boldsymbol{\sigma}_{n}^{\ell})\\}$ containing
$\boldsymbol{\sigma}_{n}^{\ell}$. Let us set

(5.5) $\displaystyle d_{\ell}:=\operatorname{diam}(\Omega_{\ell}),$ (5.6)
$\displaystyle
U_{\ell}:=\cup_{\boldsymbol{\sigma}\in\Omega_{\ell}}B(\boldsymbol{\sigma},d_{\ell}),$
(5.7)
$\displaystyle\gamma_{\ell}:=\max_{\boldsymbol{\varsigma}\in\operatorname{conv}U_{\ell}}\max_{\boldsymbol{\sigma}\in\Omega_{\ell}}\|\nabla
G_{n}(\boldsymbol{\varsigma};\boldsymbol{\sigma},\tau)\|,$ (5.8)
$\displaystyle\lambda_{\ell}:=\min_{\boldsymbol{\varsigma}\in\operatorname{conv}U_{\ell}}\lambda_{min}(\operatorname{Hess}F_{n}(\boldsymbol{\varsigma})),$
(5.9)
$\displaystyle\Lambda:=\max_{\boldsymbol{\varsigma}\in\overline{\Omega}}\lambda_{max}(\operatorname{Hess}F_{n}(\boldsymbol{\varsigma})),$
(5.10) $\displaystyle
P_{i}(\boldsymbol{\varsigma}):=\operatorname{Hess}\partial_{i}F_{n}(\boldsymbol{\varsigma})),\;i=1,2,\dots,N,$
(5.11) $\displaystyle
R_{\ell}:=\left(\max_{\boldsymbol{\varsigma}\in\operatorname{conv}U_{\ell}}\rho(P_{1}(\boldsymbol{\varsigma})),\dots,\max_{\boldsymbol{\varsigma}\in\operatorname{conv}U_{\ell}}\rho(P_{N}(\boldsymbol{\varsigma}))\right).$

Notice that the function $G_{n}(\cdot;\boldsymbol{\sigma},\tau)$ has Hessian
matrix independent by $\boldsymbol{\sigma}$:

$\operatorname{Hess}G_{n}(\boldsymbol{\xi};\boldsymbol{\sigma},\tau)=\operatorname{Hess}F_{n}(\boldsymbol{\xi})+\frac{1}{\tau}\mathbb{I}$

and it is strongly convex on $U_{\ell}$ for any
$\boldsymbol{\sigma}\in\Omega_{\ell}$, provided

(5.12) $\tau<\frac{1}{\lambda_{\ell}^{-}}=:\tau_{1,\ell}.$

Assuming (5.12), we can denote by $\hat{\boldsymbol{\sigma}}_{n}^{\ell+1}$ the
unique minimizer of $G_{n}(\cdot;\boldsymbol{\sigma}_{n}^{\ell},\tau)$ in
$U_{\ell}.$ Note that, since

$F_{n}(\hat{\boldsymbol{\sigma}}_{n}^{\ell+1})+\frac{\|\hat{\boldsymbol{\sigma}}_{n}^{\ell+1}-\boldsymbol{\sigma}_{n}^{\ell}\|^{2}}{2\tau}=G_{n}(\hat{\boldsymbol{\sigma}}_{n}^{\ell+1};\boldsymbol{\sigma}_{n}^{\ell},\tau)\leq
G_{n}(\boldsymbol{\sigma}_{n}^{\ell};\boldsymbol{\sigma}_{n}^{\ell},\tau)=F_{n}(\boldsymbol{\sigma}_{n}^{\ell}),$

then

(5.13)
$F_{n}(\boldsymbol{\sigma}_{n}^{\ell})-F_{n}(\hat{\boldsymbol{\sigma}}_{n}^{\ell+1})\geq\frac{\|\hat{\boldsymbol{\sigma}}_{n}^{\ell+1}-\boldsymbol{\sigma}_{n}^{\ell}\|^{2}}{2\tau}.$

This in particular implies that
$F_{n}(\hat{\boldsymbol{\sigma}}_{n}^{\ell+1})\leq
F_{n}(\boldsymbol{\sigma}_{n}^{\ell})$, but
$\hat{\boldsymbol{\sigma}}_{n}^{\ell+1}$ may lie in a component of
$\\{F_{n}<F(\boldsymbol{\sigma}_{n}^{\ell})\\}$ different from
$\Omega_{\ell}.$ This cannot occour because
$G_{n}(\cdot):=G_{n}(\cdot\,;\boldsymbol{\sigma}_{n}^{\ell},\tau)$ is
$(\lambda_{\ell}+1/\tau)$-convex in $U_{\ell}$, for any
$\boldsymbol{\sigma},\boldsymbol{\varsigma}\in U_{\ell}$, i.e.,

$G_{n}((1-t)\boldsymbol{\sigma}+t\boldsymbol{\varsigma})\leq(1-t)G_{n}(\boldsymbol{\sigma})+tG_{n}(\boldsymbol{\varsigma})-t(1-t)\frac{\lambda_{\ell}\tau+1}{2\tau}|\boldsymbol{\sigma}-\boldsymbol{\varsigma}|^{2},\;\;\forall
t\in(0,1)\;.$

Indeed, if we assume that $\hat{\boldsymbol{\sigma}}_{n}^{\ell+1}$ lies in the
disjoint component $\tilde{\Omega}_{\ell}$, take
$\boldsymbol{\xi}\in[\hat{\boldsymbol{\sigma}}_{n}^{\ell+1},\boldsymbol{\sigma}_{n}^{\ell}]\cap\partial\tilde{\Omega}_{\ell}$
so that $F_{n}(\boldsymbol{\xi})=F_{n}(\boldsymbol{\sigma}_{n}^{\ell})$, and
set $\boldsymbol{\sigma}=\boldsymbol{\sigma}_{n}^{\ell}$ and
$\boldsymbol{\varsigma}=\hat{\boldsymbol{\sigma}}_{n}^{\ell+1}$ to yeld

$\displaystyle G(\boldsymbol{\xi})$
$\displaystyle=F_{n}(\boldsymbol{\xi})+\frac{\|\boldsymbol{\xi}-\boldsymbol{\sigma}_{n}^{\ell}\|^{2}}{2\tau}=F_{n}(\boldsymbol{\sigma}_{n}^{\ell})+\frac{\|\boldsymbol{\xi}-\boldsymbol{\sigma}_{n}^{\ell}\|^{2}}{2\tau}$
$\displaystyle\leq(1-t)F_{n}(\boldsymbol{\sigma}_{n}^{\ell})+tF_{n}(\hat{\boldsymbol{\sigma}}_{n}^{\ell+1})$
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;+t\frac{\|\hat{\boldsymbol{\sigma}}_{n}^{\ell+1}-\boldsymbol{\sigma}_{n}^{\ell}\|^{2}}{2\tau}-t(1-t)\frac{\lambda_{\ell}\tau+1}{2\tau}\|\hat{\boldsymbol{\sigma}}_{n}^{\ell+1}-\boldsymbol{\sigma}_{n}^{\ell}\|^{2}\,,$

where
$t=\|\boldsymbol{\xi}-\boldsymbol{\sigma}_{n}^{\ell}\|/\|\hat{\boldsymbol{\sigma}}_{n}^{\ell+1}-\boldsymbol{\sigma}_{n}^{\ell}\|$.
From the previous equation we can write:

$\displaystyle
t(F_{n}(\boldsymbol{\sigma}_{n}^{\ell})-F_{n}(\hat{\boldsymbol{\sigma}}_{n}^{\ell+1}))$
$\displaystyle\leq$
$\displaystyle-\frac{\|\boldsymbol{\xi}-\boldsymbol{\sigma}_{n}^{\ell}\|}{2\tau}+t\frac{\|{\hat{\boldsymbol{\sigma}}_{n}}^{\ell+1}-\boldsymbol{\sigma}_{n}^{\ell}\|}{2\tau}-t(1-t)(\lambda_{\ell}\tau+1)\frac{\|{\hat{\boldsymbol{\sigma}}_{n}}^{\ell+1}-\boldsymbol{\sigma}_{n}^{\ell}\|}{2\tau}$
$\displaystyle\leq$
$\displaystyle-t(1-t)\lambda_{\ell}\tau\frac{\|\hat{\boldsymbol{\sigma}}_{n}^{\ell+1}-\boldsymbol{\sigma}_{n}^{\ell}\|^{2}}{2\tau}\leq
0$

which is in contrast with (5.13). As a consequence
$\hat{\boldsymbol{\sigma}}_{n}^{\ell+1}\in\Omega_{\ell}.$

We want to show that using Newton method with $\boldsymbol{\sigma}_{n}^{\ell}$
as initial guess we compute $\hat{\boldsymbol{\sigma}}_{n}^{\ell+1}$ (possibly
with an infinite number of steps). It is convenient to introduce the following
notations in order to simplify the computations:

$\displaystyle
e^{\ell+1,r}:=\boldsymbol{\sigma}_{n}^{\ell+1,r}-\hat{\boldsymbol{\sigma}}_{n}^{\ell+1}$
$\displaystyle
s^{\ell+1,r}:=\boldsymbol{\sigma}_{n}^{\ell+1,r}-\hat{\boldsymbol{\sigma}}_{n}^{\ell+1,r-1}\,.$

Then, we can rewrite the Newton step as

$s^{\ell+1,r+1}=-[\operatorname{Hess}G_{n}(\boldsymbol{\sigma}_{n}^{\ell+1,r})]^{-1}\nabla
G_{n}(\boldsymbol{\sigma}_{n}^{\ell+1,r}).$

We now assume that, for a given $r\in\mathbb{N}$ and for $s=0,1,\ldots,r$,
$\boldsymbol{\sigma}_{n}^{\ell+1,s}\in U_{\ell}.$ Writing the second order
Taylor expansion for $0=\nabla G_{n}(\hat{\boldsymbol{\sigma}}_{n}^{\ell+1})$
centered at $\boldsymbol{\sigma}_{n}^{\ell+1,r}$, multipling both sides of the
equation by
$[\operatorname{Hess}G_{n}(\boldsymbol{\sigma}_{n}^{\ell+1,r})]^{-1}$ and
using Newton equation, we obtain

$\displaystyle\|e^{\ell+1,r+1}\|$ $\displaystyle\leq$
$\displaystyle\frac{1}{2}\left\|[\operatorname{Hess}G_{n}(\boldsymbol{\sigma}_{n}^{\ell+1,r})]^{-1}((e^{\ell+1,r})^{t}P_{1}(\boldsymbol{\xi})e^{\ell+1,r},\dots,(e^{\ell+1,r})^{t}P_{N}(\boldsymbol{\xi})e^{\ell+1,r})^{t}\right\|$
$\displaystyle\leq$
$\displaystyle\frac{\tau}{2(\lambda_{\ell}\tau+1)}\|R\|\|e^{\ell+1,r}\|^{2}.$

Thus

(5.14)
$\|e^{\ell+1,r+1}\|\leq\frac{\|R\|\tau}{2(\lambda_{\ell}\tau+1)}\|e^{\ell+1,r}\|^{2}.$

Using the first order Taylor expansion of $0=\nabla
G_{n}(\hat{\boldsymbol{\sigma}}_{n}^{\ell+1})$ centered at
$\boldsymbol{\sigma}_{n}^{\ell}=\boldsymbol{\sigma}_{n}^{\ell+1,0}$ we obtain

$\|e^{\ell+1,0}\|\leq\frac{\tau}{\lambda_{\ell}\tau+1}\|\nabla
G_{n}(\boldsymbol{\sigma}_{n}^{\ell})\|\leq\frac{\tau\gamma}{\lambda_{\ell}\tau+1}.$

It is not difficult to see that, if

(5.15)
$\tau<\frac{1}{(\sqrt{\|R_{\ell}\|\gamma_{\ell}}-\lambda_{\ell})^{+}}=:\tau_{2,\ell},$

then $\frac{\|R\|\tau}{2(\lambda_{\ell}\tau+1)}\|e^{\ell+1,0}\|<1/2$. Thus,
due to (5.14), we have

(5.16) $\|e^{\ell+1,1}\|\leq\frac{1}{2}\|e^{\ell+1,0}\|\,.$

In particular,

$\frac{\|R\|\tau}{2(\lambda_{\ell}\tau+1)}\|e^{\ell+1,1}\|\leq\frac{1}{2}\frac{\|R\|\tau}{2(\lambda_{\ell}\tau+1)}\|e^{\ell+1,0}\|$

and, by finite induction,
$\|e^{\ell+1,r+1}\|\leq\frac{1}{2^{r+1}}\|e^{\ell+1,0}\|$. Therefore, we can
conclude that $\boldsymbol{\sigma}_{n}^{\ell+1,r+1}\in
B(\hat{\boldsymbol{\sigma}}_{n}^{\ell+1},d_{\ell})\subseteq U_{\ell}$ and that
we can iterate the above reasoning to get

$\|e^{\ell+1,r+1}\|\leq\frac{1}{2^{r+1}}d_{\ell}\qquad\mbox{for all
}r\in\mathbb{N}\,.$

Using in addition (5.14), we can conclude that
$\boldsymbol{\sigma}_{n}^{\ell+1,r}$ is converging quadratically to
$\boldsymbol{\sigma}_{n}^{\ell+1}:=\hat{\boldsymbol{\sigma}}_{n}^{\ell+1}$ and
(5.3) holds true.

Note that the bounds (5.12) and (5.15) on the time-step size can be uniformly
satisfied by setting

(5.17)
$\tau^{*}(c):=\min\\{\tau^{*}_{1}(c),\tau^{*}_{2}(c)\\}\leq\min\\{\inf_{\ell}\tau_{1,\ell},\inf_{\ell}\tau_{2,\ell}\\}\,,$

where

$\tau^{*}_{1}(c):=\frac{1}{\lambda^{-}_{c}}\qquad\tau^{*}_{2}(c):=\frac{1}{(\sqrt{\|R_{c}\|\gamma_{c}}-\lambda_{c})^{+}}\,,$

and

$\displaystyle d:=\operatorname{diam}(\Omega_{0})\,,$ $\displaystyle
U:=\cup_{\boldsymbol{\sigma}\in\Omega_{0}}B(\boldsymbol{\sigma},d)\,,$
$\displaystyle\lambda_{c}:=\min_{\boldsymbol{\varsigma}\in\operatorname{conv}U}\lambda_{min}(\operatorname{Hess}F_{n}(\boldsymbol{\varsigma}))\,,$
$\displaystyle
R_{c}:=\left(\max_{\boldsymbol{\varsigma}\in\operatorname{conv}U}\rho(P_{1}(\boldsymbol{\varsigma})),\dots,\max_{\boldsymbol{\varsigma}\in\operatorname{conv}U}\rho(P_{N}(\boldsymbol{\varsigma}))\right)\,,$
$\displaystyle\gamma_{c}:=\max_{\boldsymbol{\varsigma}\in\operatorname{conv}U}\max_{\boldsymbol{\sigma}\in\Omega_{0}}\|\nabla
G_{n}(\boldsymbol{\varsigma};\boldsymbol{\sigma},\tau)\|\,.$

Now we show that the sequence
$\\{\boldsymbol{\sigma}_{n}^{\ell}\\}_{\ell\in\mathbb{N}}$ is convergent to a
point $\boldsymbol{\sigma}_{n}^{*}$, that needs to be a critical point for
$F_{n}.$ To this aim, since $0=\tau\nabla
G_{n}(\boldsymbol{\sigma}_{n}^{\ell};\boldsymbol{\sigma}^{\ell-1},\tau)=\tau\nabla
F_{n}(\boldsymbol{\sigma}_{n}^{\ell})+\boldsymbol{\sigma}_{n}^{\ell}-\boldsymbol{\sigma}^{\ell-1}$,
if $\nabla F_{n}(\boldsymbol{\sigma}_{n}^{\ell})=0$ for some $\ell$, then
$\boldsymbol{\sigma}_{n}^{\ell}\equiv\boldsymbol{\sigma}^{\ell-1}\equiv\dots\equiv\boldsymbol{\sigma}^{0}$.
In such a case $\nabla F_{n}(\boldsymbol{\sigma}^{0})=\nabla
F_{n}(\boldsymbol{\sigma}_{n}^{\ell})=0$ but this violates the hypothesis
Theorem 5.1. Thus we can assume $\nabla
F_{n}(\boldsymbol{\sigma}_{n}^{\ell})\neq 0$ $\forall\ell\in\mathbb{N}.$

Since $\boldsymbol{\sigma}_{n}^{\ell}\in\Omega_{\ell}$
$\forall\ell\in\mathbb{N}$, using (5.3) and the definintion of $\lambda_{c}$,
we can write the second order Taylor expansion of
$F_{n}(\boldsymbol{\sigma}_{n}^{j})$ centered at
$\boldsymbol{\sigma}_{n}^{j+1}$ and get the estimate

$F_{n}(\boldsymbol{\sigma}_{n}^{j})-F_{n}(\boldsymbol{\sigma}_{n}^{j+1})\geq\frac{\lambda_{c}\tau+2}{2\tau}\|\boldsymbol{\sigma}_{n}^{j+1}-\boldsymbol{\sigma}_{n}^{j}\|^{2}=\tau\frac{\lambda_{c}\tau+2}{2}\|\nabla
F_{n}(\boldsymbol{\sigma}_{n}^{j+1})\|^{2}.$

Summing up over $j$ ranging from $0$ to $\ell-1$ we obtain

(5.18)
$F_{n}(\boldsymbol{\sigma}^{0})-F_{n}(\boldsymbol{\sigma}_{n}^{\ell})\geq\frac{\lambda_{c}\tau+2}{2\tau}\sum_{j=0}^{\ell}\|\boldsymbol{\sigma}^{j+1}-\boldsymbol{\sigma}^{j}\|^{2}=\tau\frac{\lambda_{c}\tau+2}{2}\sum_{j=0}^{\ell}\|\nabla
F_{n}(\boldsymbol{\sigma}^{j+1})\|^{2}.$

Since $F_{n}$ is bounded from below, then the sequences
$\\{\|\boldsymbol{\sigma}_{n}^{\ell+1}-\boldsymbol{\sigma}_{n}^{\ell}\|\\}_{\ell\in\mathbb{N}}$
and $\\{\|\nabla
F_{n}(\boldsymbol{\sigma}_{n}^{\ell})\|\\}_{\ell\in\mathbb{N}}$ are square-
summable, and we can conclude that
$\\{\boldsymbol{\sigma}_{n}^{\ell}\\}_{\ell\in\mathbb{N}}$ is Cauchy. Denoting
by $\boldsymbol{\sigma}_{n}^{*}$ the limit of
$\\{\boldsymbol{\sigma}_{n}^{\ell}\\}_{\ell\in\mathbb{N}},$ the continuity of
$\nabla F_{n}$ allows us to state that $\boldsymbol{\sigma}_{n}^{*}$ is
critical for $F_{n}.$

We are still left to prove that
$\boldsymbol{\mu}_{n}^{*}:={\boldsymbol{\sigma}_{n}^{*}}^{\boldsymbol{2}}$ is
a minimizer of $E_{n}$. To this aim we introduce the set
$I_{0}:=\\{i\in\\{1,2,\dots,N\\}:\,{\boldsymbol{\mu}_{n}}_{\\{i\\}}^{*}=0\\}$
and we distinguish the two cases $I_{0}=\emptyset$ and $I_{n}\neq\emptyset.$
In the first case, if for any $i=1,2,\dots,N$
${\boldsymbol{\mu}_{n}}_{\\{i\\}}^{*}\neq 0$ and
$0=\partial_{i}F_{n}(\boldsymbol{\sigma}_{n}^{*})=2{\boldsymbol{\sigma}_{n}}_{\\{i\\}}^{*}\partial_{i}E_{n}(\boldsymbol{\sigma}_{n}^{*})$,
then we have $\nabla E_{n}(\boldsymbol{\mu}_{n}^{*})=0$. The convexity of
$E_{n}$ allows to conclude.

In the second case we first notice that, for any $i=1,2,\dots,N$, the sequence
$\\{{\boldsymbol{\sigma}_{n}}_{\\{i\\}}^{\ell}\\}_{\ell\in\mathbb{N}}$ has
constant sign, being either constantly non-negative or constantly non-
positive. FMoreover we recall that, $\boldsymbol{\sigma}_{n}^{\ell+1}$ is
characterized as the unique minimizer of
$G_{n}(\cdot;\boldsymbol{\sigma}_{n}^{\ell},\tau)$ in $\Omega_{\ell}$. If for
one $i\in\\{1,2,\dots N\\}$ we have
$\operatorname{sign}{\boldsymbol{\sigma}_{n}}_{\\{i\\}}^{\ell}\neq\operatorname{sign}({\boldsymbol{\sigma}_{n}}_{\\{i\\}}^{\ell+1})$,
then the point

$\tilde{\boldsymbol{\sigma}}_{n}^{\ell+1}:=(\operatorname{sign}{\boldsymbol{\sigma}_{n}^{\ell}}_{\\{1\\}}|{\boldsymbol{\sigma}_{n}^{\ell+1}}_{\\{1\\}}|,\dots,\operatorname{sign}{(\boldsymbol{\sigma}_{n}^{\ell}}_{\\{N\\}}|{\boldsymbol{\sigma}_{n}^{\ell+1}}_{\\{N\\}}|)$

is such that
$F_{n}(\tilde{\boldsymbol{\sigma}}_{n}^{\ell+1})=F_{n}(\boldsymbol{\sigma}_{n}^{\ell+1})$
and
$\|\tilde{\boldsymbol{\sigma}}_{n}^{\ell+1}-\boldsymbol{\sigma}_{n}^{\ell}\|^{2}\leq\|\boldsymbol{\sigma}_{n}^{\ell+1}-\boldsymbol{\sigma}_{n}^{\ell}\|^{2}$.
Thus $\tilde{\boldsymbol{\sigma}}_{n}^{\ell+1}$ minimizes
$G_{n}(\cdot;\boldsymbol{\sigma}_{n}^{\ell},\tau)$ in $\Omega_{\ell}$ as well,
leading to a contradiction. Thus
$\\{{\boldsymbol{\sigma}_{n}}_{\\{i\\}}^{\ell}\\}_{\ell\in\mathbb{N}}$ needs
to be either non-negative for any $\ell$ or non-positive for any $\ell$.

Then, we can pick a subsequence $j\mapsto\ell_{j}$ such that
$\operatorname{sign}{\boldsymbol{\sigma}_{n}}_{\\{i\\}}^{\ell_{j}}=\operatorname{sign}\partial_{i}F_{n}(\boldsymbol{\sigma}^{\ell_{j}}),$
for any $i\in I_{0}.$ This last claim can be proven again by contradiction.
Indeed, writing
${\boldsymbol{\sigma}_{n}}_{\\{i\\}}^{\ell}={\boldsymbol{\sigma}_{n}}_{\\{i\\}}^{\ell-1}-\tau\partial_{i}F_{n}(\boldsymbol{\sigma}_{n}^{\ell})$,
and using the fact that the sign of
${\boldsymbol{\sigma}_{n}}_{\\{i\\}}^{\ell}$ is constant, shows that
$\\{{\boldsymbol{\sigma}_{n}}_{\\{i\\}}^{\ell}\\}_{\ell\in\mathbb{N}}$ must be
a non-increasing sequence for $(\boldsymbol{\sigma}_{n}^{0})_{i}<0$ or a non-
decreasing sequence for $(\boldsymbol{\sigma}_{n}^{0})_{i}>0$, and cannot
converge to $0$.

From the fact that

$\operatorname{sign}{(\boldsymbol{\sigma}_{n}^{\ell_{j}})}_{\\{i\\}}=\operatorname{sign}\partial_{i}F_{n}(\boldsymbol{\sigma}_{n}^{\ell_{j}})=\operatorname{sign}\left(2{(\boldsymbol{\sigma}_{n}^{\ell_{j}})}_{\\{i\\}}\partial_{i}E_{n}({(\boldsymbol{\sigma}_{n}^{\ell_{j}})}^{\boldsymbol{2}})\right),$

we easily obtain

$\partial_{i}E_{n}(\boldsymbol{\mu}_{n}^{*})=\lim_{j}\partial_{i}E_{n}({(\boldsymbol{\sigma}_{n}^{\ell_{j}})}^{\boldsymbol{2}})\geq
0,\;\;\forall i\in I.$

On the other hand, $\partial_{i}E_{n}(\boldsymbol{\mu}_{n}^{*})=0$ for any
$i\notin I_{0}$ due to (5.18) and the definition of $I_{0}$. Thus we have

(5.19) $\begin{cases}\partial_{i}E_{n}(\boldsymbol{\mu}_{n}^{*})\geq
0,&\forall i\\\ \partial_{i}E_{n}(\boldsymbol{\mu}_{n}^{*})=0,&\forall
i:{\boldsymbol{\mu}_{n}}_{\\{i\\}}^{*}>0\end{cases}\;\;.$

These are precisely the KKT necessary optimality conditions for $E_{n}$ on
$\mathbb{R}_{\geq 0}^{N}.$ Since the function $E_{n}$ is convex, such
conditions are also sufficient, i.e.,
$\boldsymbol{\mu}_{n}^{*}\in\operatorname*{argmin}E_{n}.$ ∎

In order to continue our study of Algorithm 2, it is convenient to introduce
some notations. Let us denote by

$\Phi(\cdot;\tau):\mathbb{R}^{N}_{\geq}\rightarrow\mathbb{R}^{N}_{\geq}$

the map that, for any $\boldsymbol{\sigma}_{n}^{\ell}\in\mathbb{R}^{N}$,
returns the exact solution of (5.2) provided by the Newton method with
$\boldsymbol{\sigma}_{n}^{\ell}$ as initial guess. As a biproduct of Theorem
5.1 this map is well defined, provided $\tau>0$ is sufficiently small. We also
define the map
$\Phi_{\epsilon}(\cdot;\tau):\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ as

$\Phi_{\epsilon}(\boldsymbol{\sigma}^{\ell};\tau):=\boldsymbol{\sigma}^{\ell+1,\hat{r}},$

where

(5.20)
$\hat{r}:=\min\left\\{r:|\partial_{i}G_{n}(\boldsymbol{\sigma}_{n}^{\ell+1,r})|\leq\epsilon|{\boldsymbol{\sigma}_{n}^{\ell+1,r}}_{\\{i\\}}-{\boldsymbol{\sigma}_{n}^{\ell}}_{\\{i\\}}|,\;\operatorname{sign}{\boldsymbol{\sigma}_{n}^{\ell+1,r}}_{\\{i\\}}=\operatorname{sign}{\boldsymbol{\sigma}_{n}^{\ell}}_{\\{i\\}}\,\forall
i\right\\}.$

We remark that, given an intial point $\boldsymbol{\sigma}_{n}^{0}$ and
suitable $\epsilon,\tau>0$, the Algorithm 2 computes the finite sequence

$(\boldsymbol{\sigma}_{n}^{0},\boldsymbol{\sigma}_{n}^{1},\boldsymbol{\sigma}_{n}^{2},\boldsymbol{\sigma}_{n}^{3},\dots):=\left(\boldsymbol{\sigma}_{n}^{0},\Phi_{\epsilon}(\boldsymbol{\sigma}_{n}^{0};\tau),\Phi^{2}_{\epsilon}(\boldsymbol{\sigma}_{n}^{0};\tau),\Phi^{3}_{\epsilon}(\boldsymbol{\sigma}_{n}^{0};\tau),\dots\right)$

of length at most $n_{step}+1$.

######  Proposition 5.2 (Stability estimate for Algorithm 2).

Let $\min F_{n}<c<+\infty$ and let $\tau<\tau^{*}(c)$ be as above. There
exists $\epsilon^{*}>0$, depending only on $c$ and $\tau$, such that

(5.21)
$F_{n}(\boldsymbol{\sigma})-F_{n}(\Phi_{\epsilon}(\boldsymbol{\sigma};\tau))\geq\frac{F_{n}(\boldsymbol{\sigma})-F_{n}(\Phi(\boldsymbol{\sigma};\tau))}{2},\;\;\forall\boldsymbol{\sigma}:F_{n}(\boldsymbol{\sigma})<c,$

for any $0<\epsilon\leq\epsilon^{*}$.

###### Proof.

It is convenient to introduce the notation

$\displaystyle\Delta$
$\displaystyle:=F_{n}(\boldsymbol{\sigma})-F_{n}(\Phi(\boldsymbol{\sigma};\tau))\;,$
$\displaystyle\Delta_{\epsilon}$
$\displaystyle:=F_{n}(\boldsymbol{\sigma})-F_{n}(\Phi_{\epsilon}(\boldsymbol{\sigma};\tau))\;,$
$\displaystyle\boldsymbol{s}$
$\displaystyle:=\Phi(\boldsymbol{\sigma};\tau)-\boldsymbol{\sigma}\;,$
$\displaystyle\boldsymbol{s}_{\epsilon}$
$\displaystyle:=\Phi_{\epsilon}(\boldsymbol{\sigma};\tau)-\boldsymbol{\sigma}\;.$

Let us pick $0<\epsilon<\lambda_{c}+1/\tau$. Using the standard error bound
for the Newton method we can write

$\|\boldsymbol{e}_{\epsilon}\|:=\|\Phi_{\epsilon}(\boldsymbol{\sigma};\tau)-\Phi(\boldsymbol{\sigma};\tau)\|\leq\max_{\boldsymbol{\xi}\in\Omega_{0}}\|[\operatorname{Hess}G_{n}(\boldsymbol{\xi})]^{-1}\|\|\nabla
G_{n}(\Phi_{\epsilon}(\boldsymbol{\sigma};\tau))\|.$

Then, by the $(\lambda_{c}+1/\tau)$-convexity of $G_{n}$ and the stopping
criterion (5.20), we get

(5.22)
$\|\boldsymbol{e}_{\epsilon}\|\leq\frac{\tau\epsilon}{\lambda_{c}\tau+1}.$

By the triangular inequality we have

$\|\boldsymbol{s}_{\epsilon}\|\leq\frac{\lambda_{c}\tau+1}{(\lambda_{c}-\epsilon)\tau+1}\|\boldsymbol{s}\|:=C_{\epsilon}\|\boldsymbol{s}\|\,,$

which yields:

(5.23)
$\|\boldsymbol{e}_{\epsilon}\|\leq\frac{C_{\epsilon}\epsilon\tau}{\lambda_{c}\tau+1}\|\boldsymbol{s}\|\,.$

Using the second order Taylor expansion of $G_{n}(\boldsymbol{\sigma})$
centered at $\Phi(\boldsymbol{\sigma};\tau)$ we can obtain

(5.24) $\Delta\geq\frac{\lambda_{c}\tau+2}{2\tau}\|\boldsymbol{s}\|^{2}.$

On the other hand, due to (5.23) and (5.24), we can write

$\displaystyle\Delta_{\epsilon}$
$\displaystyle=\Delta-(F_{n}(\Phi(\boldsymbol{\sigma};\tau))-F_{n}(\Phi_{\epsilon}(\boldsymbol{\sigma};\tau)))=\Delta+\langle\boldsymbol{e}_{\epsilon},\frac{\boldsymbol{s}}{\tau}\rangle-\frac{1}{2}\boldsymbol{e}_{\epsilon}^{T}\operatorname{Hess}F_{n}(\boldsymbol{\xi})\boldsymbol{e}_{\epsilon}$
$\displaystyle\geq\Delta-\frac{\|\boldsymbol{s}\|\|\boldsymbol{e}_{\epsilon}\|}{\tau}-\frac{\Lambda}{2}\|\boldsymbol{e}_{\epsilon}\|^{2}\geq\Delta-\|\boldsymbol{s}\|^{2}\left(C_{\epsilon}\frac{\epsilon}{\lambda_{c}\tau+1}+\frac{\Lambda
C_{\epsilon}^{2}}{2}\frac{(\epsilon\tau)^{2}}{(\lambda_{c}\tau+1)^{2}}\right)$
$\displaystyle\geq\Delta\left[1-\frac{2\tau
C_{\epsilon}\epsilon}{(\lambda_{c}\tau+2)(\lambda_{c}\tau+1)}\left(1+\frac{\Lambda\tau^{2}C_{\epsilon}\epsilon}{\lambda_{c}\tau+1}\right)\right]\;.$

In order to conclude the proof, we are left to verify that for small
$\epsilon>0$ we have

$1-\frac{2\tau
C_{\epsilon}\epsilon}{(\lambda_{c}\tau+2)(\lambda_{c}\tau+1)}\left(1+\frac{\Lambda\tau^{2}C_{\epsilon}\epsilon}{\lambda_{c}\tau+1}\right)\geq\frac{1}{2}.$

Thus

(5.25)
$C_{\epsilon}\epsilon\leq\left(\sqrt{1+\Lambda\tau(\lambda_{c}\tau+2)}-1\right)\frac{\lambda_{c}\tau+1}{2\Lambda\tau^{2}}\,.$

Note that $C_{\epsilon}\epsilon=\epsilon+o(\epsilon)$ as $\epsilon\to 0^{+}$,
so (5.25) is certainly satisfied for $\epsilon>0$ small enough. ∎

######  Theorem 5.3 (Convergence of Algorithm 2).

Let $\min F_{n}<c<+\infty$ and $\boldsymbol{\sigma}_{n}^{0}$ such that
$F_{n}(\boldsymbol{\sigma}_{n}^{0})<c$ and $\nabla
F_{n}(\boldsymbol{\sigma}_{n}^{0})\neq 0$. Let $\tau>0$ and $\epsilon$ satisfy
the hypothesis of Proposition 5.2. Then the sequence
$\\{\boldsymbol{\sigma}_{n}^{\ell}\\}_{\ell\in\mathbb{N}}:=\\{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)\\}_{\ell\in\mathbb{N}}$
admits a limit $\boldsymbol{\sigma}_{n}^{*}$. Moreover

(5.26)
$\boldsymbol{\mu}_{n}^{*}:={(\boldsymbol{\sigma}_{n}^{*})}^{\boldsymbol{2}}\in\operatorname*{argmin}E_{n}.$

###### Proof.

If $\nabla
F_{n}(\Phi_{\epsilon}^{(\ell)}(\boldsymbol{\sigma}_{n}^{0};\tau))=0$ for
certain $\ell$, then the use of
$\Phi_{\epsilon}^{(\ell)}(\boldsymbol{\sigma}_{n}^{0};\tau)$ as initial guess
in the Newton method for the computation of
$[\Phi_{\epsilon}]^{(\ell+1)}(\boldsymbol{\sigma}_{n}^{0};\tau)$ guarantees
that
$\Phi_{\epsilon}^{(\ell+j)}(\boldsymbol{\sigma}_{n}^{0};\tau)=\Phi_{\epsilon}^{(\ell)}(\boldsymbol{\sigma}_{n}^{0};\tau)$
for any $j\in\mathbb{N}.$ More importantly, using (5.20), we can repeat the
argument of the proof of Theorem 5.1, to show that
$\Phi_{\epsilon}^{(\ell)}(\boldsymbol{\sigma}_{n}^{0};\tau)\equiv\boldsymbol{\sigma}_{n}^{0}$,
which is not possible since we are assuming $\nabla
F(\boldsymbol{\sigma}_{n}^{0})\neq 0.$ Therefore we can assume without lost of
generality that $\nabla
F_{n}(\Phi_{\epsilon}^{(\ell)}(\boldsymbol{\sigma}_{n}^{0};\tau))\neq 0$
$\forall\ell\in\mathbb{N}.$

From Proposition 5.2 and the optimality of the exact step, i.e.,
$\Phi^{(\ell)}(\boldsymbol{\sigma}_{n}^{0};\tau)-\Phi^{(\ell-1)}(\boldsymbol{\sigma}_{n}^{0};\tau)=-\tau\nabla
F_{n}(\Phi^{(\ell)}(\boldsymbol{\sigma}_{n}^{0};\tau))$, to get

$\displaystyle
F_{n}(\Phi_{\epsilon}^{(\ell)}(\boldsymbol{\sigma}_{n}^{0};\tau))-F_{n}(\Phi_{\epsilon}^{(\ell+1)}(\boldsymbol{\sigma}_{n}^{0};\tau))$
$\displaystyle\geq$
$\displaystyle\frac{1}{2}F_{n}(\Phi_{\epsilon}^{(\ell)}(\boldsymbol{\sigma}_{n}^{0};\tau))-F_{n}(\Phi(\Phi_{\epsilon}^{(\ell)}(\boldsymbol{\sigma}_{n}^{0};\tau));\tau))$
$\displaystyle\geq$
$\displaystyle\frac{\lambda_{c}\tau+1}{4\tau}\left\|\Phi(\Phi_{\epsilon}^{(\ell)}(\boldsymbol{\sigma}_{n}^{0};\tau));\tau)-\Phi_{\epsilon}^{(\ell)}(\boldsymbol{\sigma}_{n}^{0};\tau)\right\|^{2}.$

Notice that, using the notation introduced in the proof of Proposition 5.2 and
in particular eq. (5.24), this last inequality can be written in the compact
form as

(5.27)
$\Delta_{\epsilon}\geq\frac{1}{2}\Delta\geq\frac{\lambda_{c}\tau+2}{2\tau}\|\boldsymbol{s}\|^{2}>\frac{\lambda_{c}\tau+1}{4\tau}\|\boldsymbol{s}\|^{2}.$

On the other hand, using (5.22), we get

(5.28)
$\|\boldsymbol{s}_{\epsilon}\|=\|\boldsymbol{s}+e_{\epsilon}\|\leq\frac{\lambda_{c}\tau+1}{(\lambda_{c}-\epsilon)\tau+1}\|\boldsymbol{s}\|.$

Therefore we have

$\displaystyle
F_{n}(\Phi_{\epsilon}^{(\ell)}(\boldsymbol{\sigma}_{n}^{0};\tau))-F_{n}(\Phi_{\epsilon}^{(\ell+1)}(\boldsymbol{\sigma}_{n}^{0};\tau))\geq\frac{[(\lambda_{c}-\epsilon)\tau+1]^{2}}{4\tau(\lambda_{c}\tau+1)}\|\boldsymbol{s}_{\epsilon}\|^{2}$
$\displaystyle=$
$\displaystyle\frac{[(\lambda_{c}-\epsilon)\tau+1]^{2}}{4\tau(\lambda_{c}\tau+1)}\|\Phi_{\epsilon}^{(\ell+1)}(\boldsymbol{\sigma}_{n}^{0};\tau)-\Phi_{\epsilon}^{(\ell)}(\boldsymbol{\sigma}_{n}^{0};\tau)\|^{2}.$

It is clear that, for any $\ell\in\mathbb{N}$ we have
$F_{n}(\boldsymbol{\sigma}_{n}^{0})-F_{n}(\Phi_{\epsilon}^{(\ell)}(\boldsymbol{\sigma}_{n}^{0};\tau))\leq
c-\min F_{n}<+\infty$. Thus we can write

$\displaystyle+\infty>c-\min
F_{n}\geq\lim_{\ell}\sum_{j=0}^{\ell-1}\left(F_{n}(\Phi_{\epsilon}^{(j)}(\boldsymbol{\sigma}_{n}^{0};\tau))-F_{n}(\Phi_{\epsilon}^{(j+1)}(\boldsymbol{\sigma}_{n}^{0};\tau))\right)$
$\displaystyle\geq$
$\displaystyle\frac{[(\lambda_{c}-\epsilon)\tau+1]^{2}}{2\tau(\lambda_{c}\tau+1)}\sum_{j=0}^{+\infty}\|\Phi_{\epsilon}^{(j+1)}(\boldsymbol{\sigma}_{n}^{0};\tau)-\Phi_{\epsilon}^{(j)}(\boldsymbol{\sigma}_{n}^{0};\tau)\|^{2}$

Thus in particular
$F_{n}(\Phi_{\epsilon}^{(\ell)}(\boldsymbol{\sigma}_{n}^{0};\tau))-F_{n}(\Phi_{\epsilon}^{(\ell+1)}(\boldsymbol{\sigma}_{n}^{0};\tau))\to
0$ as $\ell\to+\infty$ and
$\Phi_{\epsilon}^{(\ell)}(\boldsymbol{\sigma}_{n}^{0};\tau)$ is a Cauchy
sequence, whose limit we denote by $\boldsymbol{\sigma}_{n}^{*}$.

By a similar reasoning, starting from (5.27) we can show that

$\|\Phi(\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau);\tau)-\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)\|\to
0,$

thus
$\Phi(\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau);\tau)\to\boldsymbol{\sigma}_{n}^{*}.$
It follows, by the definition of the map $\Phi$, that we have

$-\tau\nabla
F_{n}(\Phi(\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau);\tau))=\Phi(\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau);\tau)-\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau).$

Therefore we have $\|\nabla
F_{n}(\boldsymbol{\sigma}_{n}^{*})\|=\lim_{\ell}\|\nabla
F_{n}(\Phi(\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau);\tau))\|=0,$
i.e., $\boldsymbol{\sigma}_{n}^{*}$ is critical for $F_{n}.$

We are left to show that $\boldsymbol{\sigma}_{n}^{*}$ is actually a local
minimizer of $F_{n}$ and thus $\boldsymbol{\mu}_{n}^{*}$ is a global minimizer
for $E_{n}.$ We reason as in the final step of the proof of Theorem 5.1. If
none of the components of
$\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0})$ tends to zero, then we
have $\nabla E_{n}({(\boldsymbol{\sigma}_{n}^{*})}^{\boldsymbol{2}})=0$ and we
can conclude using the convexity of $E_{n}$.

Let us assume that for a given $i\in\\{1,2,\dots,N\\}$ we have
$[\Phi_{\epsilon}^{\ell}{\boldsymbol{\sigma}_{n}}_{\\{i\\}}^{0};\tau)\to 0.$
Now we show that we can pick a subsequence $j\mapsto\ell_{j}$ such that

(5.29)
$\operatorname{sign}\partial_{i}F\left(\Phi_{\epsilon}^{\ell_{j}}(\boldsymbol{\sigma}_{n}^{0};\tau)\right)=\operatorname{sign}{\Phi_{\epsilon}^{\ell_{j}}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}},\;\;\forall
j\in\mathbb{N}.$

We can prove it by contradiction. Assume that

(5.30)
$\operatorname{sign}\partial_{i}F\left(\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)\right)=-\operatorname{sign}{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}},\;\;\forall\ell\in\mathbb{N}.$

Then, from the stopping criterion we obtain:

$\left|{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}-{\Phi_{\epsilon}^{\ell-1}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}+\tau\partial_{i}F\left(\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)\right)\right|\\\
\leq\tau\epsilon\left|{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}-{\Phi_{\epsilon}^{\ell-1}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}\right|\,.$

This inequality holds only if one of the following cases occurs:

$\displaystyle{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}$
$\displaystyle\leq{\Phi_{\epsilon}^{\ell-1}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}-\frac{\tau}{1-\epsilon\tau}\partial_{i}F\left(\Phi_{\epsilon}^{\ell_{j}}(\boldsymbol{\sigma}_{n}^{0};\tau)\right)$
(A)
$\displaystyle{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}$
$\displaystyle\geq{\Phi_{\epsilon}^{\ell-1}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}$
$\displaystyle{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}$
$\displaystyle\geq{\Phi_{\epsilon}^{\ell-1}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}-\tau\partial_{i}F(\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)))$
$\displaystyle{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}$
$\displaystyle\leq{\Phi_{\epsilon}^{\ell-1}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}-\frac{\tau}{1-\epsilon\tau}\partial_{i}F\left(\Phi_{\epsilon}^{\ell_{j}}(\boldsymbol{\sigma}_{n}^{0};\tau)\right)$
(B)
$\displaystyle{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}$
$\displaystyle<{\Phi_{\epsilon}^{\ell-1}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}$
$\displaystyle{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}$
$\displaystyle\geq{\Phi_{\epsilon}^{\ell-1}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}-\tau\partial_{i}F(\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)))$
$\displaystyle{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}$
$\displaystyle\geq{\Phi_{\epsilon}^{\ell-1}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}-\frac{\tau}{1+\epsilon\tau}\partial_{i}F\left(\Phi_{\epsilon}^{\ell_{j}}(\boldsymbol{\sigma}_{n}^{0};\tau)\right)$
(C)
$\displaystyle{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}$
$\displaystyle\geq{\Phi_{\epsilon}^{\ell-1}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}$
$\displaystyle{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}$
$\displaystyle<{\Phi_{\epsilon}^{\ell-1}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}-\tau\partial_{i}F(\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)))$
$\displaystyle{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}$
$\displaystyle\geq{\Phi_{\epsilon}^{\ell-1}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}-\frac{\tau}{1+\epsilon\tau}\partial_{i}F\left(\Phi_{\epsilon}^{\ell_{j}}(\boldsymbol{\sigma}_{n}^{0};\tau)\right)$
(D)
$\displaystyle{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}$
$\displaystyle<{\Phi_{\epsilon}^{\ell-1}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}$
$\displaystyle{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}$
$\displaystyle<{\Phi_{\epsilon}^{\ell-1}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}-\tau\partial_{i}F(\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)))$

Notice that, due to (5.30), we have the following implications:

Cases (A) or (C) $\displaystyle\Rightarrow$
$\displaystyle{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}\geq{\Phi_{\epsilon}^{\ell-1}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}\text{
and }{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}\geq
0\,,$ Cases (B) or (D) $\displaystyle\Rightarrow$
$\displaystyle{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}<{\Phi_{\epsilon}^{\ell-1}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}\text{
and
}{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}<0\,.$

Since the (component-wise) sign of the sequence
$\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)$ is constant by
construction, then we obtain:

$\begin{cases}{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}\text{
is non-decreasing }&\text{ if }{\boldsymbol{\sigma}_{n}}_{\\{i\\}}^{0}\geq
0\\\
{\Phi_{\epsilon}^{\ell}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}\text{ is
decreasing }&\text{ if
}{\boldsymbol{\sigma}_{n}}_{\\{i\\}}^{0}<0\end{cases}\,.$

This implies, in particular, that the sequence cannot approach $0$ as
$\ell\to+\infty$, leading to a contradiction. Therefore we can pick a
subsequence as in (5.29) and relabel it. Notice that

$\operatorname{sign}\partial_{i}F\left(\Phi_{\epsilon}^{\ell_{j}}(\boldsymbol{\sigma}_{n}^{0};\tau)\right)=\operatorname{sign}{\Phi_{\epsilon}^{\ell_{j}}(\boldsymbol{\sigma}_{n}^{0};\tau)}_{\\{i\\}}\partial_{i}E_{n}\left({\Phi_{\epsilon}^{\ell_{j}}(\boldsymbol{\sigma}_{n}^{0};\tau)}^{\boldsymbol{2}}\right),$

so we have
$\partial_{i}E_{n}\left({\Phi_{\epsilon}^{\ell_{j}}(\boldsymbol{\sigma}_{n}^{0};\tau)}^{\boldsymbol{2}}\right)\geq
0$, $\forall j\in\mathbb{N}.$ Hence

$\partial_{i}E_{n}(\boldsymbol{\mu}_{n}^{*})=\lim_{j}\partial_{i}E_{n}\left({\Phi_{\epsilon}^{\ell_{j}}(\boldsymbol{\sigma}_{n}^{0};\tau)}^{\boldsymbol{2}}\right)\geq
0.$

Finally, since we proved that the KKT conditions hold for the convex function
$E_{n}$ at $\boldsymbol{\mu}_{n}^{*}$, i.e.,

$\begin{cases}\partial_{i}E_{n}(\boldsymbol{\mu}_{n}^{*})=0,&\forall
i:{\boldsymbol{\mu}_{n}}_{\\{i\\}}^{*}\neq 0\\\
\partial_{i}E_{n}(\boldsymbol{\mu}_{n}^{*})\geq 0&\forall
i\in\\{1,2,\dots,N\\}\end{cases},$

then $\boldsymbol{\mu}_{n}^{*}$ is a global minimum for $E_{n}$ on
$\mathbb{R}^{N}_{\geq}.$ ∎

The combination of Theorem 5.3 with Theorem 2.6 leads to the following result,
which summarize the outcome of our construction.

######  Corollary 5.4.

If we assume that $\boldsymbol{\sigma}_{n}^{0}\equiv\boldsymbol{\sigma}^{0}$
with $\boldsymbol{\sigma}^{0}\in\mathcal{M}_{n}^{+}$ for any $n>\bar{n}$,
then, under the hypothesis of Theorem 5.3 and Theorem 2.6, we have

${\lim_{n}}^{*}\lim_{\ell}\mathcal{I}_{n}({\boldsymbol{\sigma}_{n}^{\ell}}^{\boldsymbol{2}})=\mu^{*},$