# A regularized Interior Point Method for sparse Optimal Transport on Graphs

Stefano Cipolla111School of Mathematics, University of Edinburgh, Edinburgh,
UK<EMAIL_ADDRESS>Jacek Gondzio222School of Mathematics, University
of Edinburgh, Edinburgh, UK<EMAIL_ADDRESS>Filippo Zanetti333School of
Mathematics, University of Edinburgh, Edinburgh, UK<EMAIL_ADDRESS>
###### Abstract

In this work, the authors address the Optimal Transport (OT) problem on graphs
using a proximal stabilized Interior Point Method (IPM). In particular,
strongly leveraging on the induced primal-dual regularization, the authors
propose to solve large scale OT problems on sparse graphs using a bespoke IPM
algorithm able to suitably exploit primal-dual regularization in order to
enforce scalability. Indeed, the authors prove that the introduction of the
regularization allows to use sparsified versions of the normal Newton
equations to inexpensively generate IPM search directions. A detailed
theoretical analysis is carried out showing the polynomial convergence of the
inner algorithm in the proposed computational framework. Moreover, the
presented numerical results showcase the efficiency and robustness of the
proposed approach when compared to network simplex solvers.

Keywords: Convex programming, primal-dual regularized interior point methods,
optimal transport on graphs, polynomial complexity, inexact interior point
methods.

##  1 Introduction

The Optimal Transport (OT) problem requires to move a certain distribution of
mass from one configuration into another, minimizing the total cost required
for the operation. It has been studied extensively, from the early work of
Kantorovich [22], to the development of ever faster algorithms for various OT
formulations, e.g. [8, 19, 27, 31, 36, 41]. Recently, there has been a growing
interest in using Interior Point Methods (IPMs) [17] in applications that
involve optimal transport, in particular for very large scale instances of
such problems, see e.g. [33, 46, 48].

A particularly interesting problem is the optimal transport over sparse
graphs: in this case, the transport of mass is only possible along a specific
subset of connections, which is noticeably smaller than the full list of edges
of a fully connected bipartite graph, as it would happen in a standard
discrete OT formulation. The use of OT and the Wasserstein distance (i.e. the
optimal objective function of the OT problem) is becoming more and more common
in many practical applications, e.g. neural networks [20], image processing
[21], inverse problems [32] and in the analysis of large complex networks
[34].

The specific formulation of the problem is the following: suppose that
$G=(V,E)$ is a connected graph with directed edges $E\subset V\times V$ and
weights $\mathbf{c}\in\mathbb{R}_{+}^{|E|}$. Define the incidence matrix
$A\in\\{-1,0,1\\}^{|V|\times|E|}$ as

$A_{ve}:=\begin{cases}-1,&\hbox{ if }\ e=(v,w)\hbox{ for some }w\in V\\\
1,&\hbox{ if }\ e=(w,v)\hbox{ for some }w\in V\\\ 0,&\hbox{ otherwise.
}\end{cases}$

We consider the optimal transport problem in the Beckmann form [12]:

$\mathcal{W}_{1}(\boldsymbol{\rho_{0}},\boldsymbol{\rho_{1}}):=\begin{cases}\min_{\mathbf{x}\in\mathbb{R}^{|E|}}&\sum_{e\in
E}c_{e}\mathbf{x}_{e}\\\
s.t.&A\mathbf{x}=\boldsymbol{\rho_{1}}-\boldsymbol{\rho_{0}}\\\
&\mathbf{x}\geq 0\par\end{cases},$ (1)

where
$\boldsymbol{\rho_{0}},\boldsymbol{\rho_{1}}\in\\{\boldsymbol{\rho}\in\mathbb{R}^{|V|}\,:\,\mathbf{1}^{T}\boldsymbol{\rho}=1\hbox{
and }\boldsymbol{\rho}\geq 0\\}=:Prob(V)$. In the following we will define
$|E|:=n$ and $|V|:=m$. OT on graphs has been recently studied in [12, 13] and,
in this formulation, it is similar to the more general minimum cost flow
problem on networks [1], which has also seen extensive use of IPMs, e.g. [5,
14, 37, 38].

Sparse graphs have on average very few edges per node, which can lead to
nearly disconnected regions and seriously limit the possible paths where mass
can be moved. As a result, finding a solution to the optimal transport problem
on a sparse graph requires more sophisticated algorithms and may be more
computationally challenging compared to solving the same problem on a denser
graph. In particular, first order methods like the network simplex may
struggle and move slowly towards optimality, due to the limited number of
edges available, while an interior point method manages to identify quickly
the subset of basic variables (i.e. the subset of edges with non-zero flow)
and converges faster.

In this work, the authors address the efficient solution of the optimal
transport problem (1) considering the Proximal-Stabilized Interior Point
framework (PS-IPM), recently introduced and analysed in [6].

As originally observed in [37], when IPMs are used to solve the minimum cost
flow problem on networks, the normal form of the related Newton systems is
structured as a Laplacian matrix of the graph (defined as the difference of
the diagonal matrix of the vertex degrees minus the adjacency matrix) and the
iterates of IPM determine the associate weights of this matrix, see also eq.
(38) . In [9], this observation was exploited to solve such Laplacian linear
systems (which are, in turn, particular instances of symmetric M-matrices)
through the fast specialized solution of $O(\ln m)$ linear systems involving
symmetric diagonally dominant matrices [42]. We refer the interested reader to
[45] for a survey on fast Laplacian solvers and to [15] for information
concerning the distribution of Laplacian’s singular values.

###  1.1 Contribution and organization

This work focuses on the efficient solution of large scale OT problems on
sparse graphs using a bespoke IPM algorithm able to suitably exploit primal-
dual regularization in order to enforce scalability. The organization of the
work and its main contributions can be summarized as follows:

  * •

In Section 2, the authors briefly recall the proximal stabilized framework
responsible for the primal-dual regularization of the IPMs here considered.

  * •

In Section 3, the authors provide a detailed convergence analysis of the
inexact infeasible primal-dual regularized IPM, when a proximal stabilization
procedure is used. Moreover, they prove its polynomial complexity.

  * •

In Section 4, the authors prove that the normal form of the related Newton
system is naturally structured as a shifted Laplacian matrix characterized by
a strict diagonal dominance. Such feature consistently simplifies the
factorization of the normal equations and allows the use of standard libraries
for the solution of the corresponding linear systems. On the other hand, such
factorizations could incur a significant fill-in even when the original graph
is sparse, hence limiting the applicability of the proposed approach for the
solution of large scale problems.

  * •

In Section 5, to overcome potential scalability issues related to the fill-in
mentioned above, the authors propose to generate IPM search directions using
sparsified versions of the IPM normal equations. In particular, the original
normal matrix takes the form $A(\Theta^{-1}+\rho I)^{-1}A^{T}+\delta I$, where
$\rho,\delta$ are regularization parameters and $\Theta$ is a diagonal matrix
related to the IPM barrier parameter $\mu$; the authors propose to use a
perturbed normal matrix, where the entries of $(\Theta^{-1}+\rho I)^{-1}$ that
are sufficiently small (when compared to $\mu$) are set to zero (completely
ignoring the corresponding columns of matrix $A$). This strategy reduces the
time required to assemble and solve the normal equations systems, providing a
fundamental advantage to the algorithm.

The resulting sparsified linear systems are solved either using a Cholesky
factorization (if that displays only negligible fill-in) or using the
conjugate gradient method and employing a simple and inexpensive incomplete
Cholesky preconditioner. In both these cases either the complete or the
incomplete Cholesky factorization remains very sparse, and this translates
into an outstanding efficiency of the proposed method. Moreover, the authors
are able to interpret the Newton directions generated using sparsified Newton
matrices as inexact Newton directions. Relying on the convergence theory
developed in Section 3, the authors are able to prove that, under suitable
choice of the sparsification parameters, the above described approach gives
rise to a polynomially convergent algorithm.

  * •

In Section 6, the authors present experimental results which demonstrate the
efficiency and robustness of the proposed approach. Extensive numerical
experiments, involving very large and sparse graphs coming from public domain
random generators as well as from real world applications, show that, for
sufficiently large problems, the approach presented in this work consistently
outperforms, in terms of computational time, the Lemon network simplex
implementation [26], one of the state-of-the-art solvers available for network
problems.

###  1.2 Notation

In the paper, vectors are indicated with bold letters. $\|\cdot\|$ indicates
the Euclidean norm. $I$ represents the identity matrix and $\mathbf{e}$ the
vector of all ones. Given a vertex $v$ of a graph $G$, we denote as $deg(v)$
its degree, i.e. the number of edges that are incident to $v$. Concerning the
variables inside the algorithm, we use a subscript $k$ to indicate the
external proximal iteration and a superscript $j$ to indicate the internal IPM
iteration. Given a sequence $\\{\mu^{j}\\}_{j\in\mathbb{N}}$ and a continuos
function $f$, the big-O notation $O(\cdot)$ is used as follows:

$\\{u^{j}\\}_{j\in\mathbb{N}}\in O\big{(}f(\mu^{j})\big{)}\hbox{ iff
}\exists\;\;C>0\hbox{ s.t. }u^{j}\leq Cf(\mu^{j})\hbox{ for all
}j\in\mathbb{N}.$

##  2 Computational Framework

###  2.1 Proximal-Stabilized Interior Point Method

Let us consider the following primal-dual formulation of a Linear Program
(LP):

$\displaystyle\begin{aligned}
\min_{\mathbf{x}\in\mathbb{R}^{n}}\;&\mathbf{c}^{T}\mathbf{x}\\\
\hbox{s.t.}\;&A\mathbf{x}=\mathbf{b}\\\ &\mathbf{x}\geq 0\\\
\end{aligned}\;\;\;\;\;\;\;\;\begin{aligned}
\max_{\mathbf{s}\in\mathbb{R}^{n},\;\mathbf{y}\in\mathbb{R}^{m}}\;&\mathbf{b}^{T}\mathbf{y}\\\
\hbox{s.t.}\;&\mathbf{c}-A^{T}\mathbf{y}-\mathbf{s}=0\\\ &\mathbf{s}\geq 0\\\
\end{aligned}$ (2)

where $A\in\mathbb{R}^{m\times n}$ with $m\leq n$ is not required to have full
rank. Notice that problem (1) is indeed formulated in this way.

We solve this problem using PS-IPM [6], which is a proximal-stabilized version
of classic Interior Point Method. Broadly speaking, PS-IPM resorts to the
Proximal Point Method (PPM) [39] to produce primal-dual regularized forms of
problem (2). Indeed, given an approximation $(\mathbf{x}_{k},\mathbf{y}_{k})$
of the solution of such problem, PS-IPM uses interior point methods to produce
the next PPM step $(\mathbf{x}_{k+1},\mathbf{y}_{k+1})$, which, in turn,
represents a better approximation of the solution of problem (2).

In this regard, the problem that needs to be solved at every PPM step takes
the form

$\displaystyle\begin{aligned}
\min_{\begin{subarray}{c}\mathbf{x}\in\mathbb{R}^{n}\\\
\mathbf{y}\in\mathbb{R}^{m}\end{subarray}}\;&\mathbf{c}^{T}\mathbf{x}+\frac{\rho}{2}\|\mathbf{x}-\mathbf{x}_{k}\|^{2}+\frac{\delta}{2}\|\mathbf{y}\|^{2}\\\
\hbox{s.t.}\;&A\mathbf{x}+\delta(\mathbf{y}-\mathbf{y}_{k})=\mathbf{b}\\\
&\mathbf{x}\geq 0,\\\ \end{aligned}\qquad\begin{aligned}
\max_{\begin{subarray}{c}\mathbf{x},\;\mathbf{s}\in\mathbb{R}^{n}\\\
\mathbf{y}\in\mathbb{R}^{m}\end{subarray}}\;&\mathbf{y}^{T}\mathbf{b}{-\frac{\rho}{2}\|\mathbf{x}\|^{2}-\frac{\delta}{2}\|\mathbf{y}-\mathbf{y}_{k}\|^{2}}\\\
\hbox{s.t.}\;&\rho(\mathbf{x}-\mathbf{x}_{k})-A^{T}\mathbf{y}-\mathbf{s}+\mathbf{c}=0\\\
&\mathbf{s}\geq 0\\\ \end{aligned}.$ (PPM$(k)$)

######  Definition 1.

Solution of problem (PPM$(k)$)
Using standard duality theory, we say that
$(\mathbf{x}_{k}^{*},\mathbf{y}_{k}^{*},\mathbf{s}_{k}^{*})$ is a solution of
problem (PPM$(k)$) if the following identities hold

$\displaystyle
A\mathbf{x}_{k}^{*}+\delta(\mathbf{y}_{k}^{*}-\mathbf{y}_{k})-\mathbf{b}=0$
$\displaystyle\rho(\mathbf{x}-\mathbf{x}_{k})-A^{T}\mathbf{y}_{k}^{*}-\mathbf{s}+\mathbf{c}=0$
(3) $\displaystyle(\mathbf{x}_{k}^{*})^{T}\mathbf{s}_{k}^{*}=0\hbox{ and
}(\mathbf{x}_{k}^{*},\mathbf{s}_{k}^{*})\geq 0$

More in particular, the PS-IPM here considered uses two nested cycles to solve
problem (2). The outer loop uses an inexact proximal point method [28], as
shown in Algorithm 1: the current approximate solution
$(\mathbf{x}_{k},\mathbf{y}_{k})$ is used to regularize the LP problem, which
is then solved using an IPM to find the next approximate solution
$(\mathbf{x}_{k+1},\mathbf{y}_{k+1})\approx(\mathbf{x}^{*}_{k},\mathbf{y}^{*}_{k})$.
And indeed, at the inner loop level, an inexact infeasible interior point
method is used to solve the PPM sub-problems, see Algorithm 2. Notice that
both methods are inexact: the outer cycle is inexact because the sub-problems
are solved approximately by an IPM; the IPM is inexact because the Newton
systems are also solved inexactly (see Section 2.2 for more details). Notice
also that the IPM is referred to as infeasible because the intermediate
iterates are not required to be inside the feasible region. We also call the
inner loop regularized, because it is a primal-dual regularized version of the
original LP (2).

Regularization in interior point methods was originally introduced in [40] and
extensively used in [2], as a tool to stabilize and improve the linear algebra
routines needed for their efficient implementation. In this work and in [6],
the regularization is introduced as a result of the application of the PPM at
the outer cycle level. To summarize, in the following we use three acronyms:
PPM refers to the outer cycle; IPM refers to the inner cycle; PS-IPM refers to
the overall procedure, combining PPM and IPM.

Input: $tol>0$, $\sigma_{r}\in(0,1)$, $\tau_{1}>0$.

Initialization: Iteration counter $k=0$; initial point
$(\mathbf{x}_{0},\mathbf{y}_{0})$

1 while  _Stopping Criterion ( 42) False_ do

2 Use Algorithm 2 with starting point
$(\mathbf{x}^{0}_{k},\mathbf{y}^{0}_{k})=(\mathbf{x}_{k},\mathbf{y}_{k})$ to
find $(\mathbf{x}_{k+1},\mathbf{y}_{k+1})$ s.t.
$\|\mathbf{r}_{k}(\mathbf{x}_{k+1},\mathbf{y}_{k+1})\|<\frac{{\sigma_{r}^{k}}}{\tau_{1}}\min\\{1,\|(\mathbf{x}_{k+1},\mathbf{y}_{k+1})-(\mathbf{x}_{k},\mathbf{y}_{k})\|\\}$
(4)

3

4 Update the iteration counter: $k:=k+1$.

5 end while

Algorithm 1 PPM, outer loop of PS-IPM

Concerning the stopping criteria, we finally highlight that Algorithm 1 is
stopped based on the criterion (42). Algorithm 2 instead, is stopped according
to the accuracy that is required for the solution of current sub-problem and
based on the following natural residual, see [28], of problem (PPM$(k)$):

######  Definition 2 (Natural Residual).

$\mathbf{r}_{k}(\textbf{x},\mathbf{y}):=\begin{bmatrix}\mathbf{x}\\\
\mathbf{y}\end{bmatrix}-\Pi_{D}\Big{(}\begin{bmatrix}\mathbf{x}\\\
\mathbf{y}\end{bmatrix}-\begin{bmatrix}\rho(\mathbf{x}-\mathbf{x}_{k})+\mathbf{c}-A^{T}\mathbf{y}\\\
A\mathbf{x}-\mathbf{b}+\delta(\mathbf{y}-\mathbf{y}_{k})\end{bmatrix}\Big{)},$

where

$D:=\mathbb{R}_{\geq 0}^{n}\times\mathbb{R}^{m}$

and where $\Pi_{D}$ is the corresponding projection operator. Moreover, it is
easy to verify that
$(\mathbf{x}_{k}^{*},\mathbf{y}_{k}^{*},\mathbf{s}_{k}^{*})$ is a solution of
problem (PPM$(k)$) if and only if
$\mathbf{r}_{k}(\mathbf{x}_{k}^{*},\mathbf{y}_{k}^{*})=0$, [6, Sec. 2.3].

###  2.2 Interior point method

We now focus on the inner cycle and give a brief description of the IPM used
to solve problem (PPM$(k)$). To this aim, we introduce the following
Lagrangian function which uses a logarithmic barrier to take into account the
inequality constraints

$\begin{split}L_{k}(\mathbf{x},\mathbf{y})=&\frac{1}{2}[\mathbf{x}^{T},\mathbf{y}^{T}]\begin{bmatrix}\rho
I&0\\\ 0&\delta I\end{bmatrix}\begin{bmatrix}\mathbf{x}\\\
\mathbf{y}\end{bmatrix}+[\mathbf{c}^{T}-\rho\mathbf{x}_{k}^{T},0]\begin{bmatrix}\mathbf{x}\\\
\mathbf{y}\end{bmatrix}\\\
&-\mathbf{y}^{T}(A\mathbf{x}+\delta(\mathbf{y}-\mathbf{y}_{k})-\mathbf{b})-\mu\sum_{i=1}^{n}\ln(x_{i}).\end{split}$
(5)

The KKT conditions that arise from the gradients of the Lagrangian (5) are

$\displaystyle\nabla_{\mathbf{x}}L_{k}(\mathbf{x},\mathbf{y})=\rho\mathbf{x}-A^{T}\mathbf{y}+\mathbf{c}-\rho{\mathbf{x}_{k}}-\begin{bmatrix}\frac{\mu}{x_{1}}\\\
\vdots\\\ \frac{\mu}{x_{n}}\end{bmatrix}=0;$
$\displaystyle-\nabla_{\mathbf{y}}L_{k}(\mathbf{x},\mathbf{y})=(A\mathbf{x}+\delta(\mathbf{y}-\mathbf{y}_{k})-\mathbf{b})=0.$

Setting $s_{i}=\frac{\mu}{x_{i}}$ for $i\in\\{1,\dots,n\\}$, we consider the
following function

$F_{k}^{\mu,\sigma}(\mathbf{x},\mathbf{y},\mathbf{s}):=\begin{bmatrix}\rho(\mathbf{x}-{\mathbf{x}_{k}})-A^{T}\mathbf{y}-\mathbf{s}+\mathbf{c}\\\
A\mathbf{x}+\delta(\mathbf{y}-\mathbf{y}_{k})-\mathbf{b}\\\
SX\mathbf{e}-\sigma\mu\mathbf{e}\end{bmatrix},$ (6)

where $\sigma\in(0,1)$ is the barrier reduction parameter,
$S=\text{diag}(\mathbf{s})$ and $X=\text{diag}(\mathbf{x})$. A primal–dual
interior point method applied to problem (PPM$(k)$) relies on the use of
Newton iterations to solve a nonlinear problem of the form

$F_{k}^{\mu,\sigma}(\mathbf{x},\mathbf{y},\mathbf{s})=0,\;\;\mathbf{x},\;\mathbf{s}>0.$

A Newton step for (6) from the current iterate
$(\mathbf{x},\mathbf{y},\mathbf{s})$ is obtained by solving the system

$\begin{bmatrix}\rho I&-A^{T}&-I\\\ A&\delta I&0\\\
S&0&X\end{bmatrix}\begin{bmatrix}\Delta\mathbf{x}\\\ \Delta\mathbf{y}\\\
\Delta\mathbf{s}\end{bmatrix}=-F_{k}^{\mu,\sigma}(\mathbf{x},\mathbf{y},\mathbf{s})=:\begin{bmatrix}\boldsymbol{\xi}_{d}\\\
\boldsymbol{\xi}_{p}\\\ \boldsymbol{\xi}_{\mu,\sigma}\end{bmatrix},$ (7)

i.e., the following relations hold:

$\displaystyle\rho\Delta\mathbf{x}-A^{T}\Delta\mathbf{y}-\Delta\mathbf{s}=\boldsymbol{\xi}_{d}$
(8) $\displaystyle
A\Delta\mathbf{x}+\delta\Delta\mathbf{y}=\boldsymbol{\xi}_{p}$ (9)
$\displaystyle
S\Delta\mathbf{x}+X\Delta\mathbf{s}=\boldsymbol{\xi}_{\mu,\sigma},$ (10)

where $(\Delta\mathbf{x},\Delta\mathbf{y},\Delta\mathbf{s})$ is the Newton
direction to be taken at each iteration (with an appropriate stepsize).

The solution of (7) is delivered by the following computational procedure

$\displaystyle\big{(}A(\Theta^{-1}+\rho I)^{-1}A^{T}+\delta
I\big{)}\Delta\mathbf{y}=\boldsymbol{\xi}_{p}-A(\Theta^{-1}+\rho
I)^{-1}(X^{-1}\boldsymbol{\xi}_{\mu,\sigma}+\boldsymbol{\xi}_{d})$ (11)
$\displaystyle(\Theta^{-1}+\rho
I)\Delta\mathbf{x}=A^{T}\Delta\mathbf{y}+\boldsymbol{\xi}_{d}+X^{-1}\boldsymbol{\xi}_{\mu,\sigma}$
(12) $\displaystyle
X\Delta\mathbf{s}=(\boldsymbol{\xi}_{\mu,\sigma}-S\Delta\mathbf{x}).$ (13)

where $\Theta:=XS^{-1}$. Before continuing let us give basic definitions used
in the remainder of this work.

######  Definition 3.

Normal Matrix:

$S_{\rho,\delta}:=A(\Theta^{-1}+\rho I)^{-1}A^{T}+\delta I.$ (14)

Neighbourhood of the infeasible central path:

$\begin{split}\mathcal{N}_{k}(\bar{\gamma},\underline{\gamma},\gamma_{p},\gamma_{d}):=&\\{(\mathbf{x},\mathbf{y},\textbf{s})\in{\color[rgb]{0,0,0}\mathbb{R}^{n}_{>0}\times\mathbb{R}^{m}\times\mathbb{R}^{n}_{>0}}\;:\;\\\
&\bar{\gamma}\mathbf{x}^{T}\mathbf{s}/n\geq
x_{i}s_{i}\geq\underline{\gamma}\mathbf{x}^{T}\mathbf{s}/n\;\hbox{ for
}\;i=1,\dots,n;\\\
&\mathbf{x}^{T}\mathbf{s}\geq\gamma_{p}\|A\mathbf{x}+\delta(\mathbf{y}-\mathbf{y}_{k})-\mathbf{b}\|;\\\
&\mathbf{x}^{T}\mathbf{s}\geq\gamma_{d}\|\rho(\mathbf{x}-\mathbf{x}_{k})-A^{T}\mathbf{y}-\mathbf{s}+\mathbf{c}\|\\},\end{split}$
(15)

where $\bar{\gamma}>1>\underline{\gamma}>0$ and $(\gamma_{p},\gamma_{d})>0$.

The neighbourhood here considered is standard in the analysis of infeasible
IPMs [47]: it requires the iterates to be close enough to the central path
(according to parameters $\bar{\gamma}$ and $\underline{\gamma}$), and the
primal-dual constraint violations to be reduced at the same rate as the
complementarity product $\mathbf{x}^{T}\mathbf{s}$. Within this neighbourhood,
$\mathbf{x}^{T}\mathbf{s}\to 0$ guarantees convergence to a primal-dual
optimal solution.

Moreover, we consider an inexact solution of the linear system (11):

######  Assumption 1.

$S_{\rho,\delta}\Delta\mathbf{y}=\bar{\boldsymbol{\xi}}_{p}+\boldsymbol{\boldsymbol{\zeta}}\hbox{
where }\|\boldsymbol{\boldsymbol{\zeta}}\|\leq
C_{\text{inexact}}\,{\mathbf{x}}^{T}\mathbf{s},$ (16)

where $C_{\text{inexact}}\in(0,1)$ and we defined

$\bar{\boldsymbol{\xi}}_{p}:=\boldsymbol{\xi}_{p}-A(\Theta^{-1}+\rho
I)^{-1}(X^{-1}\boldsymbol{\xi}_{\mu,\sigma}+\boldsymbol{\xi}_{d}).$

It is important to note that the above Assumption 1 is a non-standard
requirement in inexact Newton methods [23, 11]. Its particular form is
motivated by the use of IPM and the needs of the complexity analysis in
Section 3. It is chosen in agreement with the definition of the infeasible
neighbourhood (15) of the central path of the sub-problem considered. Using
(12) and (16) in (9), we have

$\|A\Delta\mathbf{x}+\delta\Delta\mathbf{y}-\boldsymbol{\xi}_{p}\|=\|S_{\rho,\delta}\Delta\mathbf{y}-\bar{\boldsymbol{\xi}}_{p}\|=\|\boldsymbol{\boldsymbol{\zeta}}\|,$

whereas equations (8) and (10) are satisfied exactly. Therefore the inexact
Newton directions computed according to (16) satisfy:

$\begin{bmatrix}\rho I&-A^{T}&-I\\\ A&\delta I&0\\\
S&0&X\end{bmatrix}\begin{bmatrix}\Delta\mathbf{x}\\\ \Delta\mathbf{y}\\\
\Delta\mathbf{s}\end{bmatrix}=\begin{bmatrix}\boldsymbol{\xi}_{d}\\\
\boldsymbol{\xi}_{p}\\\
\boldsymbol{\xi}_{\mu,\sigma}\end{bmatrix}+\begin{bmatrix}0\\\
\boldsymbol{\boldsymbol{\zeta}}\\\ 0\end{bmatrix}.$ (17)

Define

$\begin{bmatrix}\mathbf{x}^{j}_{k}(\alpha)\\\ \mathbf{y}^{j}_{k}(\alpha)\\\
\mathbf{s}^{j}_{k}(\alpha)\end{bmatrix}:=\begin{bmatrix}\mathbf{x}_{k}^{j}\\\
\mathbf{y}_{k}^{j}\\\
\mathbf{s}_{k}^{j}\end{bmatrix}+\begin{bmatrix}\alpha\Delta\mathbf{x}_{k}^{j}\\\
\alpha\Delta\mathbf{y}_{k}^{j}\\\
\alpha\Delta\mathbf{s}_{k}^{j}\end{bmatrix},$ (18)

i.e. $\mathbf{x}^{j}_{k}(\alpha)$ is the point reached from
$\mathbf{x}^{j}_{k}$ after a step of length $\alpha$ along the Newton
direction. Notice that, after selecting the correct stepsize $\alpha^{j}_{k}$,
we define $\mathbf{x}^{j+1}_{k}:=\mathbf{x}^{j}_{k}(\alpha^{j}_{k})$.

We report in Algorithm 2 a prototype IPM scheme for the solution of problem
(PPM$(k)$). The fundamental steps involved in the algorithm are: computing the
Newton direction by solving (17) with a level of inexactness that satisfies
(16), see Line 2; finding the largest stepsize that guarantees to remain
inside the neighbourhood and to sufficiently reduce the complementarity
products, see Line 2; preparing the quantities to be used in the next
iteration, see Lines 2-2.

We study the convergence of Algorithm 2 in Section 3. Concerning the notation,
recall that the subscript $k$ is related to the iteration count of the outer
Algorithm 1 (PPM) whereas the superscript $j$ is related to the iteration of
the inner Algorithm 2 (IPM). To avoid over-complicating the notation, notice
that in the following we use $\boldsymbol{\xi}_{p,k}^{j}$ and
$\boldsymbol{\xi}_{d,k}^{j}$ instead of $(\boldsymbol{\xi}_{p})_{k}^{j}$ and
$(\boldsymbol{\xi}_{d})_{k}^{j}$.

Input: $0<\sigma,\bar{\sigma}<1$, barrier reduction parameters s.t.
$\sigma<\bar{\sigma}$;

$C_{\text{inexact}}\in(0,1)$ inexactness parameter s.t.
$\gamma_{p}C_{\text{inexact}}<{\sigma}$;

Initialization:

Iteration counter $j=0$;

Primal–dual point
$(\mathbf{x}_{k}^{0},\mathbf{y}_{k}^{0},\mathbf{s}_{k}^{0})\in\mathcal{N}_{k}(\bar{\gamma},\underline{\gamma},\gamma_{p},\gamma_{d})$

Compute $\mu_{k}^{0}:={\mathbf{x}_{k}^{0}}^{T}\mathbf{s}_{k}^{0}/n$,
$\boldsymbol{\xi}^{0}_{d,k}$, and $\boldsymbol{\xi}^{0}_{p,k}$.

1 while  _Stopping Criterion ( 4) False_ do

2 Solve the KKT system (17) using
$[\boldsymbol{\xi}^{j}_{d,k},\boldsymbol{\xi}^{j}_{p,k},\boldsymbol{\xi}^{j}_{\mu_{k}^{j},\sigma}]^{T}$
with $\|\boldsymbol{\boldsymbol{\zeta}}_{k}^{j}\|\leq
C_{\text{inexact}}(\mathbf{x}^{j}_{k})^{T}\mathbf{s}_{k}$ to find
$[\Delta\mathbf{x}_{k}^{j},\;\Delta\mathbf{y}_{k}^{j},\;\Delta\mathbf{s}^{j}_{k}]^{T}$
;

3 Compute
$\alpha_{p}^{*,j}=\sup\\{\alpha\in\mathbb{R}\;:\;\mathbf{x}_{k}^{j}(\alpha)\geq
0\\}$
$\alpha_{d}^{*,j}=\sup\\{\alpha\in\mathbb{R}\;:\;\mathbf{s}_{k}^{j}(\alpha)\geq
0\\}$ and define $\alpha^{*,j}:=\\{\alpha_{p}^{*,j},\alpha_{d}^{*,j}\\}$;

4 if  _
$(\mathbf{x}_{k}(\alpha^{*,j}),\mathbf{y}_{k}(\alpha^{*,j}),\mathbf{s}_{k}(\alpha^{*,j}))$
is a solution of (PPM$(k)$) _ then

5 Stop

6 end if

7

Find $\alpha_{k}^{j}$ as the maximum $\alpha\in[0,1]$ s.t. for all
$\alpha\in[0,\alpha_{k}^{j}]$
$\displaystyle(\mathbf{x}_{k}^{j}(\alpha),\mathbf{y}^{j}_{k}(\alpha),\mathbf{s}^{j}_{k}(\alpha))\in\mathcal{N}_{k}(\bar{\gamma},\underline{\gamma},\gamma_{p},\gamma_{d})\hbox{
and }$ 8(19)
$\displaystyle\mathbf{x}^{j}_{k}(\alpha)^{T}\mathbf{s}_{k}^{j}(\alpha)\leq(1-(1-\bar{\sigma})\alpha){\mathbf{x}_{k}^{j}}^{T}\mathbf{s}_{k}^{j}\hbox{;\lx@prepend@indentation@{\hskip
5.0pt\lx@algo@rule\hskip
10.00002pt}\algocf@printnl{8}\lx@algo@@endline\lx@algo@@startline }$

9

10 Set $\begin{bmatrix}\mathbf{x}_{k}^{j+1}\\\ \mathbf{y}_{k}^{j+1}\\\
\mathbf{s}_{k}^{j+1}\end{bmatrix}=\begin{bmatrix}\mathbf{x}_{k}^{j}\\\
\mathbf{y}_{k}^{j}\\\
\mathbf{s}_{k}^{j}\end{bmatrix}+\begin{bmatrix}\alpha_{k}^{j}\Delta\mathbf{x}_{k}^{j}\\\
\alpha_{k}^{j}\Delta\mathbf{y}_{k}^{j}\\\
\alpha_{k}^{j}\Delta\mathbf{s}_{k}^{j}\end{bmatrix}$ ;

11 Compute the infeasibilities $\boldsymbol{\xi}^{j+1}_{d,k}$,
$\boldsymbol{\xi}^{j+1}_{p,k}$ and barrier parameter
$\mu^{j+1}_{k}:={\mathbf{x}^{j+1}_{k}}^{T}\mathbf{s}^{j+1}_{k}/n$ ;

12 Update the iteration counter: $j:=j+1$.

13 end while

Algorithm 2 IPM: inner loop of PS-IPM

##  3 Convergence and complexity

In this section, we show that the particular inexact IPM in Algorithm 2, used
as inner solver in Algorithm 1, is convergent. Moreover, at the end of the
present section, we show that such IPM converges to an $\varepsilon-$accurate
solution in a polynomial number of iterations. Our implant of the proof is
inspired by the works [3, 4, 7, 16, 24, 25] but consistently differs from the
hypothesis and techniques used there.

The PPM iteration counter $k$ is fixed through this section and, for the sake
of readability, is used only when writing the fixed PPM iteration
$(\mathbf{x}_{k},\mathbf{y}_{k},\mathbf{s}_{k})$ and not in the context of the
IPM iterations $(\mathbf{x}^{j}_{k},\mathbf{y}^{j}_{k},\mathbf{s}^{j}_{k})$.

We start from analysing the progress made in a single Newton iteration. Using
(6), (7), and (8) we obtain

$\begin{split}&\rho(\mathbf{x}^{j}(\alpha)-\mathbf{x}_{k})-A^{T}\mathbf{y}^{j}{(\alpha)}-\mathbf{s}^{j}{(\alpha)}+\mathbf{c}\\\
&=(\rho(\mathbf{x}^{j}-\mathbf{x}_{k})-A^{T}\mathbf{y}^{j}-\mathbf{s}^{j}+\mathbf{c})+\alpha(\rho\Delta\mathbf{x}^{j}-A^{T}\Delta\mathbf{y}^{j}-\Delta\mathbf{s}^{j})\\\
&=(1-\alpha)(\rho(\mathbf{x}^{j}-\mathbf{x}_{k})-A^{T}\mathbf{y}^{j}-\mathbf{s}^{j}+\mathbf{c}),\end{split}\;\;\;$
(20)

whereas, using (6) and (17) we have

$\begin{split}&A\mathbf{x}^{j}(\alpha)+\delta(\mathbf{y}^{j}(\alpha)-\mathbf{y_{k}})-\mathbf{b}\\\
&=(A\mathbf{x}^{j}+\delta(\mathbf{y}^{j}-\mathbf{y_{k}})-\mathbf{b})+\alpha(A\Delta\mathbf{x}^{j}+\delta\Delta\mathbf{y}^{j})\\\
&=(1-\alpha)(A\mathbf{x}^{j}+\delta(\mathbf{y}^{j}-\mathbf{y_{k}})-\mathbf{b})+\alpha\boldsymbol{\boldsymbol{\zeta}}^{j}.\end{split}$
(21)

The last block equation in (17) yields

$(\mathbf{s}^{j})^{T}\Delta\mathbf{x}^{j}+(\mathbf{x}^{j})^{T}\Delta\mathbf{s}^{j}=-(\mathbf{x}^{j})^{T}\mathbf{s}^{j}+\sigma
n\mu^{j}=(\sigma-1)(\mathbf{x}^{j})^{T}\mathbf{s}^{j}$ (22)

and

$s_{i}\Delta x_{i}+x_{i}\Delta
s_{i}=\sigma\frac{\mathbf{x}^{T}\mathbf{s}}{n}-x_{i}s_{i}\;\;.$ (23)

Finally, using (22), we state the following identity

$(\mathbf{x}^{j}+\alpha\Delta\mathbf{x}^{j})^{T}(\mathbf{s}^{j}+\alpha\Delta\mathbf{s}^{j})=(\mathbf{x}^{j})^{T}\mathbf{s}^{j}(1+\alpha(\sigma-1))+\alpha^{2}(\Delta\mathbf{x}^{j})^{T}\Delta\mathbf{s}^{j}.$
(24)

With the next Theorem 1 we prove that Algorithm 2 is well-defined: at each
iteration, there exist a non-empty interval of values for the stepsize
$\alpha$ such that the next iterate still lies in the neighbourhood
$\mathcal{N}_{k}(\bar{\gamma},\underline{\gamma},\gamma_{p},\gamma_{d})$ and
such that the complementarity product
$(\mathbf{x}_{k}^{j})^{T}\mathbf{s}_{k}^{j}$ is reduced by a sufficient
amount, as required in (19).

######  Theorem 1.

Let us suppose that
$(\mathbf{x}^{j},\mathbf{y}^{j},\mathbf{s}^{j})\in\mathcal{N}_{k}(\bar{\gamma},\underline{\gamma},\gamma_{p},\gamma_{d})$
s.t. $(\mathbf{x}^{j})^{T}\mathbf{s}^{j}>0$ is given. If the stopping
conditions at Line 2 of Algorithm 2 are not satisfied, then there exists
$0<\hat{\alpha}^{j}{\color[rgb]{0,0,0}<\alpha^{*,j}}$ such that conditions
(19) are satisfied for all $\alpha\in[0,\hat{\alpha}^{j}]$.

###### Proof.

In this proof we omit also the IPM iterate counter $j$, i.e.
$(\mathbf{x}^{j},\mathbf{y}^{j},\mathbf{s}^{j})\equiv(\mathbf{x},\mathbf{y},\mathbf{s})$.
Let us define the following functions, for all $i=1,\dots,n$

$\displaystyle f_{i}(\alpha):=$ $\displaystyle(x_{i}+\alpha\Delta
x_{i})(s_{i}+\alpha\Delta
s_{i})-\underline{\gamma}(\mathbf{x}+\alpha\Delta\mathbf{x})^{T}(\mathbf{s}+\alpha\Delta\mathbf{s})/n,$
$\displaystyle\bar{f}_{i}(\alpha):=$
$\displaystyle\bar{\gamma}(\mathbf{x}+\alpha\Delta\mathbf{x})^{T}(\mathbf{s}+\alpha\Delta\mathbf{s})/n-(x_{i}+\alpha\Delta
x_{i})(s_{i}+\alpha\Delta s_{i}),$ $\displaystyle h(\alpha):=$
$\displaystyle(1-(1-\bar{\sigma})\alpha)\mathbf{x}^{T}\mathbf{s}-(\mathbf{x}+\alpha\Delta\mathbf{x})^{T}(\mathbf{s}+\alpha\Delta\mathbf{s}),$
$\displaystyle g_{d}(\alpha):=$
$\displaystyle(\mathbf{x}+\alpha\Delta\mathbf{x})^{T}(\mathbf{s}+\alpha\Delta\mathbf{s})$
$\displaystyle-\gamma_{d}\|\rho(\mathbf{x}+\alpha\Delta\mathbf{x}-\mathbf{x}_{k})-A^{T}(\mathbf{y}+\alpha\Delta\mathbf{y})-(\mathbf{s}+\alpha\Delta\mathbf{s})+\mathbf{c}\|,$
$\displaystyle g_{p}(\alpha):=$
$\displaystyle(\mathbf{x}+\alpha\Delta\mathbf{x})^{T}(\mathbf{s}+\alpha\Delta\mathbf{s})$
$\displaystyle-\gamma_{p}\|A(\mathbf{x}+\alpha\Delta\mathbf{x})+\delta(\mathbf{y}+\alpha\Delta\mathbf{y}-\mathbf{y}_{k})-\mathbf{b}\|.$

Using (20) in the expressions of $g_{d}(\alpha)$ we have

$g_{d}(\alpha)=(\mathbf{x}+\alpha\Delta\mathbf{x})^{T}(\mathbf{s}+\alpha\Delta\mathbf{s})-\gamma_{d}(1-\alpha)\|\rho(\mathbf{x}-\mathbf{x}_{k})-A^{T}\mathbf{y}-\mathbf{s}+\mathbf{c}\|,$

whereas using (21) in the expressions of $g_{p}(\alpha)$ we have

$g_{p}(\alpha)\geq(\mathbf{x}+\alpha\Delta\mathbf{x})^{T}(\mathbf{s}+\alpha\Delta\mathbf{s})-\gamma_{p}((1-\alpha)\|A\mathbf{x}+\delta(\mathbf{y}-\mathbf{y}_{k})-\mathbf{b}\|+\alpha\|\boldsymbol{\boldsymbol{\zeta}}\|).$

We start proving that there exists $\hat{\alpha}^{j}>0$ such that

$f_{i}(\alpha)\geq 0,\;\bar{f}_{i}(\alpha)\geq 0,\;h(\alpha)\geq
0,\;g_{p}(\alpha)\geq 0,\;g_{d}(\alpha)\geq 0$

for all $i=1,\dots,n$ and for all $\alpha\in[0,\hat{\alpha}^{j}]$. In the
following we will use exensively the identity (24). We have

$\begin{split}f_{i}(\alpha)=&\underbrace{(1-\alpha)(x_{i}s_{i}-\underline{\gamma}\frac{\mathbf{x}^{T}\mathbf{s}}{n})}_{\geq
0}+\alpha^{2}(\Delta x_{i}\Delta
s_{i}-\underline{\gamma}\frac{(\Delta\mathbf{x})^{T}\Delta\mathbf{s}}{n})+\alpha\sigma(1-\underline{\gamma})\frac{\mathbf{x}^{T}\mathbf{s}}{n}\\\
&\geq\alpha^{2}(\Delta x_{i}\Delta
s_{i}-\underline{\gamma}\frac{(\Delta\mathbf{x})^{T}\Delta\mathbf{s}}{n})+\alpha\sigma(1-\underline{\gamma})\frac{\mathbf{x}^{T}\mathbf{s}}{n}.\end{split}$
(25)

Since $\mathbf{x}^{T}\mathbf{s}>0$, using a simple continuity argument, we can
infer the existence of a small enough $\underline{f}_{i}>0$ s.t.
$f_{i}(\alpha)\geq 0$ for all $\alpha\in[0,\underline{f}_{i}]$. Reasoning
analogously, we have

$\begin{split}&\bar{f}_{i}(\alpha)=\underbrace{(1-\alpha)(\bar{\gamma}\frac{\mathbf{x}^{T}\mathbf{s}}{n}-x_{i}s_{i})}_{\geq
0}+\alpha^{2}(\bar{\gamma}\frac{(\Delta\mathbf{x})^{T}\Delta\mathbf{s}}{n}-\Delta
x_{i}\Delta
s_{i})+\alpha\sigma(\bar{\gamma}-1)\frac{\mathbf{x}^{T}\mathbf{s}}{n}\\\
&\geq\alpha^{2}(\bar{\gamma}\frac{(\Delta\mathbf{x})^{T}\Delta\mathbf{s}}{n}-\Delta
x_{i}\Delta
s_{i})+\alpha\sigma(\bar{\gamma}-1)\frac{\mathbf{x}^{T}\mathbf{s}}{n},\end{split}$
(26)

and hence there exists a small enough $\bar{f}_{i}>0$ s.t.
$\bar{f}_{i}(\alpha)\geq 0$ for all $\alpha\in[0,\bar{f}_{i}]$.

Concerning $h(\alpha)$, we have

$h(\alpha)=\mathbf{x}^{T}\mathbf{s}(\bar{\sigma}-\sigma)\alpha-\alpha^{2}(\Delta\mathbf{x})^{T}\Delta\mathbf{s},$
(27)

and, since $\mathbf{x}^{T}\mathbf{s}(\bar{\sigma}-\sigma)>0$, there exists
$\hat{h}>0$ small enough s.t. $h(\alpha)\geq 0$ for all
$\alpha\in[0,\hat{h}]$.

Concerning $g_{d}(\alpha)$, we have

$\begin{split}g_{d}(\alpha)=&\underbrace{(1-\alpha)(\mathbf{x}^{T}\mathbf{s}-\gamma_{d}\|(\rho(\mathbf{x}-\mathbf{x}_{k})-A^{T}\mathbf{y}-\mathbf{s}+\mathbf{c})\|)}_{\geq
0}+\alpha\sigma\mathbf{x}^{T}\mathbf{s}+\alpha^{2}(\Delta\mathbf{x})^{T}\Delta\mathbf{s}\\\
&\geq\alpha\sigma\mathbf{x}^{T}\mathbf{s}+\alpha^{2}(\Delta\mathbf{x})^{T}\Delta\mathbf{s},\end{split}$
(28)

and hence there exists $\hat{g}_{d}>0$ small enough s.t. $g_{d}(\alpha)\geq 0$
for all $\alpha\in[0,\hat{g}_{d}]$.

Finally, concerning $g_{p}(\alpha)$, we have

$\begin{split}g_{p}(\alpha)\geq&\underbrace{(1-\alpha)(\mathbf{x}^{T}\mathbf{s}-\gamma_{p}\|(A\mathbf{x}+\delta(\mathbf{y}-\mathbf{y}_{k})-\mathbf{b})\|)}_{\geq
0}+\alpha\sigma\mathbf{x}^{T}\mathbf{s}+\\\
&+\alpha^{2}(\Delta\mathbf{x})^{T}\Delta\mathbf{s}-\alpha\gamma_{p}\|\boldsymbol{\boldsymbol{\zeta}}\|\\\
&\geq\alpha(\sigma-\gamma_{p}C_{\text{inexact}})\mathbf{x}^{T}\mathbf{s}+\alpha^{2}(\Delta\mathbf{x})^{T}\Delta\mathbf{s},\end{split}$
(29)

and hence there exists $\hat{g}_{p}>0$ small enough s.t. $g_{p}(\alpha)\geq 0$
for all $\alpha\in[0,\hat{g}_{p}]$. Let us define

$\hat{\alpha}^{j}=\min\\{\min_{i}\underline{f}_{i},\;\;\min_{i}\bar{f}_{i},\;\;\hat{h},\;\;\hat{g}_{d},\;\;\hat{g}_{p},\;1\\}>0.$

To prove the thesis, it remains to show that $\alpha^{*,j}>\hat{\alpha}^{j}$,
i.e. that

$(\mathbf{x}(\alpha),\mathbf{y}(\alpha),\mathbf{s}(\alpha))\in\mathbb{R}^{n}_{>0}\times\mathbb{R}^{m}\times\mathbb{R}^{n}_{>0}\hbox{
for all }\alpha\in[0,\hat{\alpha}^{j}].$

To this aim, let us suppose by contradiction that
$\alpha^{*,j}\leq\hat{\alpha}^{j}$. By definition of $\alpha^{*,j}$, there
exists $\bar{\ell}\in\\{1,\dots,n\\}$ s.t. $(x_{\bar{\ell}}+\alpha^{*,j}\Delta
x_{\bar{\ell}})(s_{\bar{\ell}}+\alpha^{*,j}\Delta s_{\bar{\ell}})=0$. We have
hence

$f_{\bar{\ell}}(\alpha^{*,j})=-\underline{\gamma}(\mathbf{x}\big{(}\alpha^{*,j})\big{)}^{T}\mathbf{s}(\alpha^{*,j})/n\geq
0\Rightarrow(\mathbf{x}\big{(}\alpha^{*,j})\big{)}^{T}\mathbf{s}(\alpha^{*,j})=0.$

From the above implication, using $g_{d}(\alpha^{*,j})$ and
$g_{p}(\alpha^{*,j})$, we obtain that

$\begin{split}&A\mathbf{x}(\alpha^{*,j})+\delta(\mathbf{y}(\alpha^{*,j})-\mathbf{y}_{k})-\mathbf{b}=0\\\
&\rho(\mathbf{x}(\alpha^{*,j})-\mathbf{x}_{k})-A^{T}\mathbf{y}(\alpha^{*,j})-\mathbf{s}(\alpha^{*,j})+\mathbf{c}=0,\end{split}$

i.e.
$(\mathbf{x}(\alpha^{*,j}),\mathbf{y}(\alpha^{*,j}),\mathbf{s}(\alpha^{*,j}))$
is a solution of problem (PPM$(k)$). We have hence obtained a contradiction
since we are supposing that Algorithm 2 did not stop at Line 2.  ∎

Before proving the convergence of Algorithm 2 we would like to emphasize that
the above proof complements and expands [24, Remark 3.1].

The next two results establish that Algorithm 2 converges to a solution of
problem (PPM$(k)$). This is done by establishing that the right-hand sides of
the Newton systems are uniformly bounded and by showing that the
complementarity product $(\mathbf{x}^{j})^{T}\mathbf{s}^{j}$ cannot be bounded
away from zero.

######  Corollary 1.

The right-hand sides of the Newton systems are uniformly bounded.

###### Proof.

As a consequence of Theorem 1, we can suppose the existence of a sequence of
iterates
$\\{(\mathbf{x}^{j},\mathbf{y}^{j},\mathbf{s}^{j})\\}_{j\in\mathbb{N}}$
produced by Algorithm 2 s.t.

$(\mathbf{x}^{j},\mathbf{y}^{j},\mathbf{s}^{j})\in\mathcal{N}_{k}(\bar{\gamma},\underline{\gamma},\gamma_{p},\gamma_{d}).$

Since by construction
$(\mathbf{x}^{j})^{T}\mathbf{s}^{j}\leq(\mathbf{x}^{0})^{T}\mathbf{s}^{0}$, we
have from (15)

$\begin{split}&\|A\mathbf{x}^{j}+\delta(\mathbf{y}^{j}-\mathbf{y}_{k})-\mathbf{b}\|\leq(\mathbf{x}^{0})^{T}\mathbf{s}^{0}/\gamma_{p},\\\
&\|\rho(\mathbf{x}^{j}-\mathbf{x}_{k})-A^{T}\mathbf{y}^{j}-\mathbf{s}^{j}+\mathbf{c}\|\leq(\mathbf{x}^{0})^{T}\mathbf{s}^{0}/\gamma_{d}.\end{split}$

Moreover, we have

$\begin{split}&\|S^{j}X^{j}\mathbf{e}-\sigma\mu^{j}\mathbf{e}\|\leq\|S^{j}X^{j}\mathbf{e}\|+\sigma\mu^{j}\|\mathbf{e}\|\\\
&\leq\frac{\bar{\gamma}}{\sqrt{n}}(\mathbf{x}^{j})^{T}\mathbf{s}^{j}+\frac{\sigma}{\sqrt{n}}(\mathbf{x}^{j})^{T}\mathbf{s}^{j}\leq\frac{\bar{\gamma}+\sigma}{\sqrt{n}}(\mathbf{x}^{0})^{T}\mathbf{s}^{0}.\end{split}$

∎

######  Theorem 2.

Algorithm 2 produces a sequence of iterates in
$\mathcal{N}_{k}(\bar{\gamma},\underline{\gamma},\gamma_{p},\gamma_{d})$ s.t.
$\lim(\mathbf{x}^{j})^{T}\mathbf{s}^{j}=0$, i.e.
$(\mathbf{x}^{j},\mathbf{y}^{j},\mathbf{s}^{j})$ converges to a solution of
problem (PPM$(k)$).

###### Proof.

Let us argue by contradiction supposing that there exists $\varepsilon^{*}>0$
s.t. $(\mathbf{x}^{j})^{T}\mathbf{s}^{j}>\varepsilon^{*}$ for all
$j\in\mathbb{N}$.

Claim 1 There exists a constant $C_{1}$ dependent only on $n$ s.t.

$\|[\Delta\mathbf{x}^{j},\;\Delta\mathbf{y}^{j},\;\Delta\mathbf{s}^{j}]^{T}\|\leq
C_{1}\hbox{ for all }j\in\mathbb{N}.$

The proof of this fact follows observing that the Newton matrices in (17)
satisfy all the hypothesis of [3, Th. 1], i.e. they have a uniformly bounded
inverse, and that the right-hand sides are uniformly bounded, see Corollary 1.
As a consequence, there exists another constant $C_{2}$ s.t.

$\begin{split}&|\Delta x_{i}^{j}\Delta
s_{i}^{j}-\underline{\gamma}\frac{(\Delta\mathbf{x}^{j})^{T}\Delta\mathbf{s}^{j}}{n}|\leq
C_{2},\\\
&|\bar{\gamma}\frac{(\Delta\mathbf{x}^{j})^{T}\Delta\mathbf{s}^{j}}{n}-\Delta
x_{i}^{j}\Delta s_{i}^{j}|\leq C_{2},\\\
&|(\Delta\mathbf{x}^{j})^{T}\Delta\mathbf{s}^{j}|\leq C_{2}.\end{split}$ (30)

Claim 2 There exists $\alpha^{*}>0$ s.t. $\alpha^{j}\geq\alpha^{*}$ for all
$j\in\mathbb{N}$.

Using (30) in equations (25), (26), (27), (28), (29) we have

$\displaystyle f_{i}(\alpha)\geq$ $\displaystyle\;\alpha^{2}(\Delta
x_{i}^{j}\Delta
s_{i}^{j}-\underline{\gamma}\frac{(\Delta\mathbf{x}^{j})^{T}\Delta\mathbf{s}^{j}}{n})+\alpha\sigma(1-\underline{\gamma})\frac{(\mathbf{x}^{j})^{T}\mathbf{s}^{j}}{n}$
$\displaystyle\geq$
$\displaystyle\;-C_{2}\alpha^{2}+\alpha\sigma(1-\underline{\gamma})\varepsilon^{*}/n,$
$\displaystyle\bar{f}_{i}(\alpha)\geq$
$\displaystyle\;(\bar{\gamma}\frac{(\Delta\mathbf{x}^{j})^{T}\Delta\mathbf{s}^{j}}{n}-\Delta
x_{i}^{j}\Delta
s_{i}^{j})+\alpha\sigma(\bar{\gamma}-1)\frac{(\mathbf{x}^{j})^{T}\mathbf{s}^{j}}{n}$
$\displaystyle\geq$
$\displaystyle\;-C_{2}\alpha^{2}+\alpha\sigma(\bar{\gamma}-1)\varepsilon^{*}/n,$
$\displaystyle h(\alpha)=$
$\displaystyle\;(\mathbf{x}^{j})^{T}\mathbf{s}^{j}(\bar{\sigma}-\sigma)\alpha-\alpha^{2}(\Delta\mathbf{x}^{j})^{T}\Delta\mathbf{s}^{j}$
$\displaystyle\geq$
$\displaystyle\;\varepsilon^{*}(\bar{\sigma}-\sigma)\alpha-C_{2}\alpha^{2},$
$\displaystyle g_{d}(\alpha)\geq$
$\displaystyle\;\alpha\sigma(\mathbf{x}^{j})^{T}\mathbf{s}^{j}+\alpha^{2}(\Delta\mathbf{x}^{j})^{T}\Delta\mathbf{s}^{j}\geq\varepsilon^{*}\sigma\alpha-
C_{2}\alpha^{2},$ $\displaystyle g_{p}(\alpha)\geq$
$\displaystyle\;\alpha(\sigma-q\gamma_{p})(\mathbf{x}^{j})^{T}\mathbf{s}^{j}+\alpha^{2}(\Delta\mathbf{x}^{j})^{T}\Delta\mathbf{s}^{j}\geq\varepsilon^{*}\alpha(\sigma-\gamma_{p}C_{\text{inexact}})-C_{2}\alpha^{2}.$

Hence $\alpha^{j}\geq\alpha^{*}$, where

$\alpha^{*}:=\min\bigg{\\{}{\color[rgb]{0,0,0}1},\frac{\sigma(1-\underline{\gamma})\varepsilon^{*}/n}{C_{2}},\frac{\sigma(\bar{\gamma}-1)\varepsilon^{*}/n}{C_{2}},\frac{(\bar{\sigma}-\sigma)\varepsilon^{*}}{C_{2}},\frac{\sigma\varepsilon^{*}}{C_{2}},\frac{(\sigma-\gamma_{p}C_{\text{inexact}})\varepsilon^{*}}{C_{2}}\bigg{\\}}.$
(31)

The convergence claim follows observing that the inequality

$\varepsilon^{*}\leq(\mathbf{x}^{j})^{T}\mathbf{s}^{j}\leq(1-(1-\bar{\sigma})\alpha^{*})^{j}(\mathbf{x}^{0})^{T}\mathbf{s}^{0}$

leads to a contradiction for $j\to\infty$. ∎

###  3.1 Polynomial Complexity

In this section we show that the number of iterations needed to reduce
$\mu^{j}$ below a certain tolerance $\varepsilon$ grows polynomially with the
size of the problem. Moreover, it is important to note that, from the
definition of the central path
$\mathcal{N}_{k}(\bar{\gamma},\underline{\gamma},\gamma_{p},\gamma_{d})$, the
same number of iterations is sufficient to reduce also the primal and dual
infeasibility below the tolerance $\varepsilon$.

In the following we will omit the index $j$ when this does not lead to
ambiguities. It is important to note that the linear system in (17) can be
written in an alternative form as follows

$\underbrace{\begin{bmatrix}\rho I&A^{T}&-I\\\ A&-\delta I&0\\\
S&0&X\end{bmatrix}}_{=:J}\begin{bmatrix}\Delta\mathbf{x}\\\
-\Delta\mathbf{y}\\\
\Delta\mathbf{s}\end{bmatrix}=\begin{bmatrix}\boldsymbol{\xi}_{d}\\\
\boldsymbol{\xi}_{p}+\boldsymbol{\boldsymbol{\zeta}}\\\
\boldsymbol{\xi}_{\mu,\sigma}\end{bmatrix}.$

Using [3, Remark 1], we have that

${J}^{-1}=\begin{bmatrix}H^{-1}&\frac{1}{\delta}H^{-1}A^{T}&H^{-1}X^{-1}\\\
\frac{1}{\delta}AH^{-1}&\frac{1}{\delta^{2}}AH^{-1}A^{T}-\frac{1}{\delta}I&\frac{1}{\delta}AH^{-1}X^{-1}\\\
-\Theta^{-1}H^{-1}&-\frac{1}{\delta}\Theta^{-1}H^{-1}A^{T}&(I-\Theta^{-1}H^{-1})X^{-1}\end{bmatrix},$
(32)

where $H:=\rho I+\Theta^{-1}+\frac{1}{\delta}A^{T}A$.

To prove polynomial complexity, we start by bounding the terms that appear in
the expression (32). The next two technical results are useful in this sense.

######  Lemma 1.

We have that

$\|H^{-1}\|\in O(1).$

###### Proof.

Using the Sherman-Morrison-Woodbury formula, we get

$H^{-1}=(\rho I+\Theta^{-1})^{-1}-\frac{1}{\delta}(\rho
I+\Theta^{-1})^{-1}A^{T}\big{(}I+\frac{1}{\delta}A(\rho
I+\Theta^{-1})^{-1}A^{T}\big{)}^{-1}A(\rho I+\Theta^{-1})^{-1}.$ (33)

We observe

$(\rho
I+\Theta^{-1})^{-1}_{ii}=\frac{\Theta_{ii}}{\rho\Theta_{ii}+1}=\frac{\rho\Theta_{ii}}{\rho\Theta_{ii}+1}\frac{1}{\rho}<\frac{1}{\rho}$
(34)

and hence

$\begin{split}\|H^{-1}\|\;\leq\;&\|(\rho I+\Theta^{-1})^{-1}\|\cdot\\\
&\cdot\Big{(}1+\frac{1}{\delta}\|A^{T}\|\|A\|\|\big{(}I+\frac{1}{\delta}A(\rho
I+\Theta^{-1})^{-1}A^{T})^{-1}\|\|(\rho I+\Theta^{-1}\big{)}^{-1}\|\Big{)}\\\
\leq\;&\frac{1}{\rho}(1+\frac{1}{\delta\rho}\|A^{T}\|\|A\|),\end{split}$

where we used that $\|\big{(}I+\frac{1}{\delta}A(\rho
I+\Theta^{-1})^{-1}A^{T}\big{)}^{-1}\|\leq 1$. ∎

######  Corollary 2.

We have that

$\|\Theta^{-1}H^{-1}\|\in O(1).$

###### Proof.

Using (34), we observe that

$\big{(}\Theta^{-1}(\rho
I+\Theta^{-1})^{-1}\big{)}_{ii}=\frac{1}{\rho\Theta_{ii}+1}$

and hence, using (33), we have

$\begin{split}\|\Theta^{-1}H^{-1}\|\leq\;&\|\Theta^{-1}(\rho
I+\Theta^{-1})^{-1}\|\;\cdot\\\
&\Big{(}1+\frac{1}{\delta}\|A^{T}\|\|A\|\|\big{(}1+\frac{1}{\delta}A(\rho
I+\Theta^{-1})^{-1}A^{T})^{-1}\|\|(\rho I+\Theta^{-1}\big{)}^{-1}\|\Big{)}\\\
\leq\;&1+\frac{1}{\delta\rho}\|A^{T}\|\|A\|,\\\ \end{split}$

where we used (34). ∎

Therefore, we know that

$\|H^{-1}\|\leq C_{3}\hbox{ and }\|\Theta^{-1}H^{-1}\|\leq C_{4},$

for some positive constants $C_{3}$ and $C_{4}$. Let us define
$C_{5}:=\max\\{C_{3},C_{4}\\}$.

Now that we have bounded the terms in (32), we look for an upper bound on the
norms of the Newton directions that depend polynomially on the size of the
problem $n$. This is a crucial step to find a polynomial lower bound on the
minimum stepsize (31) which leads to the polynomial complexity result
mentioned at the beginning of this section.

######  Theorem 3.

There exists a positive constant $C_{6}$ s.t.

$\|[\Delta\mathbf{x}^{j},\;\Delta\mathbf{y}^{j},\;\Delta\mathbf{s}^{j}]^{T}\|\leq
C_{6}n\sqrt{\mu^{j}}\hbox{ for all }j\in\mathbb{N}.$

###### Proof.

Using (32), we have

$\begin{split}&\Delta\mathbf{x}^{j}={(H^{j})^{-1}\boldsymbol{\xi}^{j}_{d}+\frac{1}{\delta}(H^{j})^{-1}A^{T}(\boldsymbol{\xi}^{j}_{p}+\boldsymbol{\boldsymbol{\zeta}}^{j})}+{(H^{j})^{-1}(X^{j})^{-1}\boldsymbol{\xi}^{j}_{\sigma,\mu^{j}}},\\\
&-\Delta\mathbf{y}^{j}={\frac{1}{\delta}A(H^{j})^{-1}\boldsymbol{\xi}^{j}_{d}+\Big{(}\frac{1}{\delta^{2}}A(H^{j})^{-1}A^{T}-\frac{1}{\delta}I\Big{)}(\boldsymbol{\xi}^{j}_{p}+\boldsymbol{\boldsymbol{\zeta}}^{j})}+{\frac{1}{\delta}A(H^{j})^{-1}(X^{j})^{-1}\boldsymbol{\xi}^{j}_{\sigma,\mu^{j}}},\\\
\end{split}$

and hence

$\small\begin{split}\|\Delta\mathbf{x}^{j}\|\leq\;&\|(H^{j})^{-1}\|\|\boldsymbol{\xi}^{j}_{d}\|+\frac{1}{\delta}\|(H^{j})^{-1}\|\|A^{T}\|\|\boldsymbol{\xi}^{j}_{p}+\boldsymbol{\zeta}^{j}\|\\\
+&\|(H^{j})^{-1/2}\|\|(H^{j})^{-1/2}(X^{j})^{-1/2}(S^{j})^{1/2}\|\|(X^{j})^{1/2}(S^{j})^{1/2}\mathbf{e}-\sigma\mu^{j}(X^{j})^{-1/2}(S^{j})^{-1/2}\mathbf{e}\|\\\
\leq\;&\frac{C_{5}n\mu^{j}}{\gamma_{d}}+\frac{C_{5}\|A^{T}\|(1+\gamma_{p}C_{\text{inexact}})n\mu^{j}}{\delta\gamma_{p}}+C_{5}\Big{(}\sqrt{\bar{\gamma}n}+\sigma\frac{\sqrt{n}}{\sqrt{\underline{\gamma}}}\Big{)}\sqrt{\mu^{j}}\\\
\leq&{\color[rgb]{0,0,0}\bigg{(}\frac{C_{5}\sqrt{\mu^{0}}}{\gamma_{d}}+\frac{C_{5}\|A^{T}\|(1+\gamma_{p}C_{\text{inexact}})\sqrt{\mu^{0}}}{\delta\gamma_{p}}+C_{5}\Big{(}\sqrt{\bar{\gamma}}+\sigma\frac{1}{\sqrt{\underline{\gamma}}}\Big{)}\bigg{)}n\sqrt{\mu^{j}}}\\\
\leq\;&{\color[rgb]{0,0,0}C_{\Delta x}n\sqrt{\mu^{j}},}\end{split}$ (35)

where we used

$\begin{split}&\|(X^{j})^{1/2}(S^{j})^{1/2}\mathbf{e}-\sigma\mu^{j}(X^{j})^{-1/2}(S^{j})^{-1/2}\mathbf{e}\|\\\
\leq\;&\|(X^{j})^{1/2}(S^{j})^{1/2}\mathbf{e}\|+\|\sigma\mu^{j}(X^{j})^{-1/2}(S^{j})^{-1/2}\mathbf{e}\|,\\\
\end{split}$

$\mu^{j}=\sqrt{\mu^{j}}\,\sqrt{\mu^{j}}\leq\sqrt{\mu^{0}}\,\sqrt{\mu^{j}}$ and
$\sqrt{n}\leq n$, and where $C_{\Delta x}$ is a positive constant.

Analogously

$\begin{split}\|\Delta\mathbf{y}^{j}\|\leq\;&\frac{1}{\delta}\|A\|\|(H^{j})^{-1}\|\|\boldsymbol{\xi}^{j}_{d}\|+\Big{(}\frac{1}{\delta^{2}}\|(H^{j})^{-1}\|\|A\|\|A^{T}\|+\frac{1}{\delta}\Big{)}\|\boldsymbol{\xi}^{j}_{p}+\boldsymbol{\zeta}^{j}\|\\\
+&\frac{1}{\delta}\|A\|\|(H^{j})^{-1/2}\|\|(H^{j})^{-1/2}(X^{j})^{-1/2}(S^{j})^{1/2}\|\cdot\\\
&\cdot\|(X^{j})^{1/2}(S^{j})^{1/2}\mathbf{e}-\sigma\mu^{j}(X^{j})^{-1/2}(S^{j})^{-1/2}\mathbf{e}\|\\\
\leq\;&\frac{\|A\|C_{5}n\mu^{j}}{\delta\gamma_{d}}+\frac{(C_{5}\|A^{T}\|\|A\|+\delta)(1+\gamma_{p}C_{\text{inexact}})n\mu^{j}}{\delta^{2}\gamma_{p}}+\\\
&+\frac{C_{5}\|A\|}{\delta}\Big{(}\sqrt{\bar{\gamma}n}+\sigma\frac{\sqrt{n}}{\sqrt{\underline{\gamma}}}\Big{)}\sqrt{\mu^{j}}\\\
\leq\;&{\color[rgb]{0,0,0}\bigg{(}\frac{\|A\|C_{5}\sqrt{\mu^{0}}}{\delta\gamma_{d}}+\frac{(C_{5}\|A^{T}\|\|A\|+\delta)(1+\gamma_{p}C_{\text{inexact}})\sqrt{\mu^{0}}}{\delta^{2}\gamma_{p}}+}\\\
&{\color[rgb]{0,0,0}+\frac{C_{5}\|A\|}{\delta}\Big{(}\sqrt{\bar{\gamma}}+\sigma\frac{1}{\sqrt{\underline{\gamma}}}\Big{)}\bigg{)}n\sqrt{\mu^{j}}}\\\
\leq\;&{\color[rgb]{0,0,0}C_{\Delta y}n\sqrt{\mu^{j}},}\end{split}$ (36)

where $C_{\Delta y}$ is a positive constant. Finally, using the fact that
$\Delta\mathbf{s}^{j}=\rho\Delta\mathbf{x}^{j}-A^{T}\Delta\mathbf{y}^{j}-\boldsymbol{\xi}_{d}^{j}$
and using (35), (36), and the definition of
$\mathcal{N}_{k}(\bar{\gamma},\underline{\gamma},\gamma_{p},\gamma_{d})$, we
have

$\begin{split}\|\Delta\mathbf{s}^{j}\|=&\rho\|\Delta\mathbf{x}^{j}\|+\|A^{T}\|\|\Delta\mathbf{y}^{j}\|+\|\boldsymbol{\xi}_{d}^{j}\|\\\
\leq\;&{\color[rgb]{0,0,0}\rho C_{\Delta x}n\sqrt{\mu^{j}}+\|A^{T}\|C_{\Delta
y}n\sqrt{\mu^{j}}+\frac{n\mu^{j}}{\gamma_{d}}}\\\
\leq\;&{\color[rgb]{0,0,0}\Big{(}\rho C_{\Delta x}+\|A^{T}\|C_{\Delta
y}+\frac{\sqrt{\mu^{0}}}{\gamma_{d}}\Big{)}n\sqrt{\mu^{j}}}\\\
\leq\;&{\color[rgb]{0,0,0}C_{\Delta s}n\sqrt{\mu^{j}},}\end{split}$

where $C_{\Delta s}$ is a positive constant.

The thesis follows setting $C_{6}=\max(C_{\Delta x},C_{\Delta y},C_{\Delta
s})$. ∎

The next straightforward technical result specializes the polynomial bound of
Theorem 3 to the terms that appear in equations (25)-(29).

######  Corollary 3.

There exists a positive constant $C_{7}$ such that, for all $i$

$\begin{split}&\Big{|}\Delta x_{i}\Delta
s_{i}-\underline{\gamma}\frac{(\Delta\mathbf{x})^{T}\Delta\mathbf{s}}{n}\Big{|}\leq
C_{7}n^{2}\mu,\\\
&\Big{|}\bar{\gamma}\frac{(\Delta\mathbf{x})^{T}\Delta\mathbf{s}}{n}-\Delta
x_{i}\Delta s_{i}\Big{|}\leq C_{7}n^{2}\mu,\\\
&|(\Delta\mathbf{x})^{T}\Delta\mathbf{s}|\leq C_{7}n^{2}\mu.\end{split}$ (37)

We can now apply the previous results and obtain a bound similar to (31), but
that depends polynomially on the size of the problem $n$. This is the last
fundamental step before the polynomial complexity result can be stated.

######  Theorem 4.

There exists a constant $\tilde{\alpha}$ s.t. $\alpha^{j}\geq\tilde{\alpha}$
for all $j\in\mathbb{N}$ and

$\tilde{\alpha}\geq C_{8}n^{-2},$

where $C_{8}$ is a positive constant.

###### Proof.

Using (37) in equations (25), (26), (27), (28), (29) we have:

$\displaystyle f_{i}(\alpha)\geq$ $\displaystyle\;\alpha^{2}\Big{(}\Delta
x_{i}\Delta
s_{i}-\underline{\gamma}\frac{(\Delta\mathbf{x})^{T}\Delta\mathbf{s}}{n}\Big{)}+\alpha\sigma(1-\underline{\gamma})\frac{\mathbf{x}^{T}\mathbf{s}}{n}$
$\displaystyle\geq$
$\displaystyle\;-C_{7}n^{2}\mu\alpha^{2}+\alpha\sigma(1-\underline{\gamma})\mu,$
$\displaystyle\bar{f}_{i}(\alpha)\geq$
$\displaystyle\;\alpha^{2}\Big{(}\bar{\gamma}\frac{(\Delta\mathbf{x})^{T}\Delta\mathbf{s}}{n}-\Delta
x_{i}\Delta
s_{i}\Big{)}+\alpha\sigma(\bar{\gamma}-1)\frac{\mathbf{x}^{T}\mathbf{s}}{n}$
$\displaystyle\geq$
$\displaystyle\;-C_{7}n^{2}\mu\alpha^{2}+\alpha\sigma(\bar{\gamma}-1)\mu,$
$\displaystyle h(\alpha)=$
$\displaystyle\;\mathbf{x}^{T}\mathbf{s}(\bar{\sigma}-\sigma)\alpha-\alpha^{2}(\Delta\mathbf{x})^{T}\Delta\mathbf{s}\geq
n\mu(\bar{\sigma}-\sigma)\alpha-C_{7}n^{2}\mu\alpha^{2},$ $\displaystyle
g_{d}(\alpha)\geq$
$\displaystyle\;\alpha\sigma\mathbf{x}^{T}\mathbf{s}+\alpha^{2}(\Delta\mathbf{x})^{T}\Delta\mathbf{s}\geq
n\mu\sigma\alpha-C_{7}n^{2}\mu\alpha^{2},$ $\displaystyle g_{p}(\alpha)\geq$
$\displaystyle\;\alpha(\sigma-\gamma_{p}C_{\text{inexact}})\mathbf{x}^{T}\mathbf{s}+\alpha^{2}(\Delta\mathbf{x})^{T}\Delta\mathbf{s}\geq
n\mu\alpha(\sigma-\gamma_{p}C_{\text{inexact}})-C_{7}n^{2}\mu\alpha^{2}.$

Hence, defining

${\color[rgb]{0,0,0}\tilde{\alpha}:=}\min\bigg{\\{}{\color[rgb]{0,0,0}1},\frac{\sigma(1-\underline{\gamma})}{C_{7}n^{2}},\;\frac{\sigma(\bar{\gamma}-1)}{C_{7}n^{2}},\;\frac{(\bar{\sigma}-\sigma)}{C_{7}n},\;\frac{\sigma}{C_{7}n},\;\frac{(\sigma-\gamma_{p}C_{\text{inexact}})}{C_{7}n}\bigg{\\}},$

the thesis follows observing that, by definition,
$\alpha^{j}\geq\tilde{\alpha}$.

∎

Finally, we are ready to show that the number of iterations required to reduce
$\mu$ below a certain tolerance $\varepsilon$ is proportional to $n^{2}$.

######  Theorem 5.

Algorithm 2 has polynomial complexity, i.e. given $\varepsilon>0$ there exists
$K\in O(n^{2}\ln(\frac{1}{\varepsilon}))$ s.t. $\mu^{j}\leq\varepsilon$ for
all $j\geq K$.

###### Proof.

Thesis follows observing that

$(\mathbf{x}^{j})^{T}\mathbf{s}^{j}\leq\big{(}1-(1-\bar{\sigma})\tilde{\alpha}\big{)}^{j}(\mathbf{x}^{0})^{T}\mathbf{s}^{0}\leq\Big{(}1-(1-\bar{\sigma})\frac{C_{8}}{n^{2}}\Big{)}^{j}(\mathbf{x}^{0})^{T}\mathbf{s}^{0}.$

∎

##  4 Properties of the regularized normal equations system

In this section we show some properties of matrix $S_{\rho,\delta}$ that are
useful for the analysis performed in the next section.

For the original graph $G(V,E)$, we define the adjacency matrix
$\mathcal{A}\in\mathbb{R}^{|V|\times|V|}$ such that $\mathcal{A}_{ij}=1$ if
there exists an edge between nodes $i$ and $j$ and $\mathcal{A}_{ij}=0$
otherwise. Notice that for an undirected graph $\mathcal{A}$ is symmetric; we
assume that there are no self-loops, so that the diagonal of $\mathcal{A}$ is
made of zeros. Let us define the degree matrix of a graph
$\mathcal{D}\in\mathbb{R}^{|V|\times|V|}$ such that $\mathcal{D}$ is diagonal
and $\mathcal{D}_{jj}$ is the degree of node $j$. Notice that
$\mathcal{D}_{jj}=(\mathcal{A}\mathbf{e})_{j}$.

Let us define also the Laplacian matrix of a graph as
$\mathcal{L}\in\mathbb{R}^{|V|\times|V|}$ such that
$\mathcal{L}=\mathcal{D}-\mathcal{A}$. An important relation between the
Laplacian $\mathcal{L}$ and the node-arc incidence matrix $A$ is that
$\mathcal{L}=AA^{T}$.

Given a diagonal matrix $\Theta$ and a parameter $\rho$, we define the re-
weighted graph $G_{\Theta}$ as the graph with the same connectivity of $G$, in
which the weight of every edge $j$ is scaled by a factor
$\sqrt{\frac{\Theta_{jj}}{1+\rho\Theta_{jj}}}$. The adjacency matrix of the
new graph $\mathcal{A}_{\Theta}$ has the same sparsity pattern of
$\mathcal{A}$, but takes into account the new weight of the edges. The same
happens for the degree matrix $\mathcal{D}_{\Theta}$. The incidence matrix
therefore becomes $A_{\Theta}=A(\Theta^{-1}+\rho I)^{-1/2}$.

The new Laplacian matrix thus reads
$\mathcal{L}_{\Theta}=\mathcal{D}_{\Theta}-\mathcal{A}_{\Theta}$ and can be
written as

$\mathcal{L}_{\Theta}=A_{\Theta}A_{\Theta}^{T}=A(\Theta^{-1}+\rho
I)^{-1}A^{T}.$

We have just shown that the normal equations matrix can be interpreted as the
Laplacian matrix of the original graph, where the edges are weighted according
to the diagonal entries of matrix $\Theta$. This result is important because
solving linear systems that involve Laplacian matrices is much easier than
solving general linear systems. We summarize this result in the following
Lemma.

######  Lemma 2.

The matrix $A(\Theta^{-1}+\rho I)^{-1}A^{T}$ is the Laplacian of a weighted
undirected graph and hence, for every $\Theta\in\mathbb{R}^{n\times n}_{+}$

$A(\Theta^{-1}+\rho I)^{-1}A^{T}=\mathcal{D}_{\Theta}-\mathcal{A}_{\Theta},$
(38)

where $\mathcal{D}_{\Theta}$ and $\mathcal{A}_{\Theta}$ are the degree and
adjacency matrices of the weighted graph.

The next result shows that the normal matrix is strictly diagonally dominant,
due to the presence of dual regularization. This property is significant
because it assures that the incomplete Cholesky factorization of the normal
equations matrix $S_{\rho,\delta}$ can always be computed without the
algorithm breaking down (in exact arithmetic), see e.g. [29].

######  Lemma 3.

If $\delta>0$ the matrix $S_{\rho,\delta}$ is strictly diagonally dominant.

###### Proof.

From Lemma 2, we have that $A(\Theta^{-1}+\rho
I)^{-1}A^{T}=\mathcal{D}_{\Theta}-\mathcal{A}_{\Theta}$ and hence

$\sum_{j\neq i}|(S_{\rho,\delta})_{ij}|=\sum_{j\neq
i}|-(\mathcal{A}_{\Theta})_{ij}|=\sum_{j\neq
i}(\mathcal{A}_{\Theta})_{ij}<({\mathcal{D}_{\Theta}})_{ii}+\delta=(S_{\rho,\delta})_{ii}.$

∎

The next two technical results are related to the distribution of eigenvalues
of the normal matrix. They are used in the next section to show that the
inexactness introduced when sparsifying the normal matrix remains bounded.

######  Lemma 4.

$\lambda_{max}(AA^{T})\leq 2\max_{v\in V}deg(v)$.

###### Proof.

The proof is straightforward using Gershgorin’s circle Theorems. ∎

######  Lemma 5.

The eigenvalues $\lambda$ of matrix $S_{\rho,\delta}$ satisfy

$\delta\leq\lambda<\delta+\frac{2}{\rho}\max_{v\in V}deg(v).$

###### Proof.

Using the Rayleigh quotient, for some vector $\mathbf{u}$ and
$\mathbf{v}=A^{T}\mathbf{u}$, the eigenvalues can be written as

$\lambda=\frac{\mathbf{v}^{T}(\Theta^{-1}+\rho
I)^{-1}\mathbf{v}}{\mathbf{v}^{T}\mathbf{v}}\frac{\mathbf{u}^{T}AA^{T}\mathbf{u}}{\mathbf{u}^{T}\mathbf{u}}+\delta$

The lower bound $\lambda\geq\delta$ is trivial; the upper bound follows from
Lemma 4 and (34). ∎

##  5 Sparsification of the reduced matrix

We now propose a technique to reduce the number of nonzeros in the normal
equations $S_{\rho,\delta}$, based on the weights of the edges in the re-
weighted graph (according to Lemma 2). We then show that this sparsification
strategy is sound and produces a polynomially convergent interior point
algorithm.

In this section we omit the IPM iteration counter $j$ and we consider all the
IPM-related quantities as a function of $\mu\to 0$. As IPMs progress towards
optimality, we expect the following partition of the diagonal matrix $\Theta$
contributed by the barrier term:

$\begin{split}&\mathcal{B}:=\\{i=1,\dots,n\hbox{ s.t. }x_{i}\to
x_{i}^{*}>0,\;s_{i}\to s_{i}^{*}=0\\}\\\ &\mathcal{N}:=\\{i=1,\dots,n\hbox{
s.t. }x_{i}\to x_{i}^{*}=0,\;s_{i}\to s_{i}^{*}>0\\},\end{split}$

where the optimal solution $(x^{*},y^{*},s^{*})$ was defined in (1). Notice
that $(\mathcal{B},\mathcal{N})$ is the partition corresponding to the optimal
solution.

######  Assumption 2.

We suppose that the following asymptotic estimates hold

$\begin{split}&s_{i}\in O(\mu)\hbox{ and }x_{i}\in O(1)\hbox{ for
}i\in\mathcal{B}\\\ &x_{i}\in O(\mu)\hbox{ and }s_{i}\in O(1)\hbox{ for
}i\in\mathcal{N}\\\ \end{split}$ (39)

and, since $x^{-1}_{i}s_{i}\approx\mu x^{-2}_{i}$ when an IPM iterate is
sufficiently close to the central path, using (39), we suppose

$\Theta_{ii}^{-1}=x^{-1}_{i}s_{i}=O({\mu})\hbox{ for }i\in\mathcal{B}\hbox{
and }\;\;\Theta_{ii}^{-1}=x^{-1}_{i}s_{i}=O(\mu^{-1})\hbox{ for
}i\in\mathcal{N}.$

This assumption makes sense given the neighbourhood that is considered, see
e.g. [18].

Due to Assumption 2, we consider the following asymptotic estimates of
$(\Theta^{-1}+\rho I)^{-1}$

$(\Theta^{-1}+\rho I)^{-1}_{ii}=\begin{cases}O(\frac{1}{\rho+\mu})&\hbox{ if
}i\in\mathcal{B}\\\ O(\frac{\mu}{1+\rho\mu})&\hbox{ if
}i\in\mathcal{N}.\end{cases}$

The diagonal entries of $\Theta$ give a specific weight to each column of
matrix $A$ (or equivalently, give a weight to each edge of the original sparse
graph, as shown in Lemma 2). The columns for which the corresponding
$\Theta_{ii}$ is $O(\mu)$ have a very small impact on the normal matrix, but
still contribute to its sparsity pattern. In order to save time and memory
when forming (complete or incomplete) Cholesky factorizations, we propose the
following sparsification strategy: we introduce a suitable threshold
$C_{t}\in\mathbb{R}_{+}$ and define

$(\widehat{\Theta}_{C_{t}\mu,\rho}^{\dagger})_{ii}:=\begin{cases}(\Theta^{-1}+\rho
I)^{-1}_{ii}&\hbox{ if }(\Theta^{-1}+\rho
I)^{-1}_{ii}\geq\frac{C_{t}\mu}{1+\rho\mu}\\\ 0&\hbox{ if }(\Theta^{-1}+\rho
I)^{-1}_{ii}<\frac{C_{t}\mu}{1+\rho\mu}.\end{cases}$ (40)

We define the $\mu$-sparisified version $S^{C_{t}\mu}_{\rho,\delta}$ of
$S_{\rho,\delta}$ as

$S^{C_{t}\mu}_{\rho,\delta}:=A\widehat{\Theta}_{C_{t}\mu,\rho}^{\dagger}A^{T}+\delta
I.$ (41)

Notice that this matrix completely ignores some of the columns of $A$ (and
some of the edges of the graph). The dual regularization $\delta I$ guarantees
that the resulting matrix is non-singular, irrespective of the level of
sparsification chosen. In this paper, we consider using inexact Newton
directions produced by solving linear systems with matrix
$S^{C_{t}\mu}_{\rho,\delta}$, rather than $S_{\rho,\delta}$.

######  Remark 1.

It is important to note that, in general, the sparsity pattern of the matrix
$S^{C_{t}\mu}_{\rho,\delta}$ depends on the choice of the parameter $C_{t}$
and on the partitioning $(\mathcal{B},\mathcal{N})$. Indeed, when $\mu$ is
sufficiently small, we expect that

$\Big{|}\Big{\\{}i\in\\{1,\dots,n\\}\hbox{ s.t. }(\Theta^{-1}+\rho
I)^{-1}_{ii}\geq\frac{C_{t}\mu}{1+\rho\mu}\Big{\\}}\Big{|}=|\mathcal{B}|.$

Let us now show how the algorithm is affected by the use of the proposed
sparsified normal matrix. Notice that the results presented below depend
strongly on two facts: the optimization problem evolves on a graph and thus
the normal matrix is a Laplacian, with very desirable properties; the interior
point method employs primal-dual regularization.

We start by showing how much the normal matrix deviates from its sparsified
counterpart.

######  Theorem 6.

The sparsification strategy in (40) produces a matrix which is close to the
original $S_{\rho,\delta}$, in the sense that

$\|S_{\rho,\delta}-S^{C_{t}\mu}_{\rho,\delta}\|=O\Big{(}\frac{C_{t}\mu}{1+\rho\mu}\Big{)}.$

###### Proof.

We have that $S_{\rho,\delta}-S^{C_{t}\mu}_{\rho,\delta}=AE^{\mu}A^{T}$ where
$E^{\mu}$ is the diagonal matrix defined as

$E^{\mu}_{ii}:=(\Theta^{-1}+\rho
I)^{-1}_{ii}-(\widehat{\Theta}_{C_{t}\mu,\rho}^{\dagger})_{ii}=\begin{cases}0&\hbox{
if }(\Theta^{-1}+\rho I)^{-1}_{ii}\geq\frac{C_{t}\mu}{1+\rho\mu}\\\
(\Theta^{-1}+\rho I)^{-1}_{ii}&\hbox{ if }(\Theta^{-1}+\rho
I)^{-1}_{ii}<\frac{C_{t}\mu}{1+\rho\mu}\end{cases}.$

Hence we have $\lambda_{max}(E^{\mu})\leq\frac{C_{t}\mu}{1+\rho\mu}$. Thesis
follows using Lemma 4 and observing that

$\frac{\mathbf{v}^{T}(S_{\rho,\delta}-S^{C_{t}\mu}_{\rho,\delta})\mathbf{v}}{\mathbf{v}^{T}\mathbf{v}}=\frac{\mathbf{v}^{T}(AE^{\mu}A^{T})\mathbf{v}}{\mathbf{v}^{T}\mathbf{v}}\leq\lambda_{max}(E^{\mu})\lambda_{max}(AA^{T}).$

∎

We now show that the condition number of both matrices is uniformly bounded.
This is an important property when using iterative Krylov solvers to find the
Newton directions.

######  Lemma 6.

When $\mu$ is sufficiently small, the condition numbers of $S_{\rho,\delta}$
and $S^{C_{t}\mu}_{\rho,\delta}$ satisfy

${k}_{2}(S_{\rho,\delta})\in
O\bigg{(}1+\frac{1}{\delta(\rho+\mu)}\bigg{)}\;\;\hbox{ and
}\;\;{k}_{2}(S^{C_{t}\mu}_{\rho,\delta})\in
O\bigg{(}1+\frac{1}{\delta(\rho+\mu)}\bigg{)}.$

###### Proof.

The thesis follows from Lemma 4 and observing that for
$\mathbf{v}=A^{T}\mathbf{w}$

$\begin{split}&\delta\leq\frac{\mathbf{w}^{T}S_{\rho,\delta}\mathbf{w}}{\mathbf{w}^{T}\mathbf{w}}\leq\delta+\frac{\mathbf{v}^{T}(\Theta^{-1}+\rho
I)^{-1}\mathbf{v}}{\mathbf{v}^{T}\mathbf{v}}\frac{\mathbf{w}^{T}AA^{T}\mathbf{w}}{\mathbf{w}^{T}\mathbf{w}}\\\
&\leq\delta+O\Big{(}\frac{1}{\rho+\mu}\big{(}2\max_{v\in
V}deg(v)\big{)}\Big{)}.\end{split}$

∎

We now show that the solution of the sparsified linear system is “close" to
the solution of the original one, and the bound depends on $\mu$. This result
depends on the spectral distribution that was shown in the previous section.

######  Theorem 7.

For all $\mathbf{v}\in\mathbb{R}^{m}$ we have that

$(S^{C_{t}\mu}_{\rho,\delta})^{-1}\mathbf{v}=S^{-1}_{\rho,\delta}\mathbf{v}+\boldsymbol{\boldsymbol{\psi}}$

where $\|\boldsymbol{\boldsymbol{\psi}}\|\in
O(\frac{C_{t}\mu}{\delta^{2}(1+\rho\mu)}\|\mathbf{v}\|)$.

###### Proof.

Using Lemma 5 and Theorem 6, we have that

$\|I-S_{\rho,\delta}^{-1}S_{\rho,\delta}^{C_{t}\mu}\|\leq\|S_{\rho,\delta}^{-1}\|\|S_{\rho,\delta}-S^{C_{t}\mu}_{\rho,\delta}\|\in
O\bigg{(}\frac{C_{t}\mu}{\delta(1+\rho\mu)}\bigg{)},$

and hence

$\|(S_{\rho,\delta}^{-1}-(S_{\rho,\delta}^{C_{t}\mu})^{-1})\mathbf{v}\|\leq\|(S_{\rho,\delta}^{C_{t}\mu})^{-1}\|\|S_{\rho,\delta}^{-1}\|\|(S_{\rho,\delta}-S^{C_{t}\mu}_{\rho,\delta})\mathbf{v}\|\in
O\bigg{(}\frac{C_{t}\mu}{\delta^{2}(1+\rho\mu)}\|\mathbf{v}\|\bigg{)}.$

∎

The next technical result is useful for the proof of Corollary 4.

######  Lemma 7.

$\|\bar{\boldsymbol{\xi}}_{p}^{j}\|$ is uniformly bounded, i.e. there exists a
constant $C_{9}>0$ such that for all $j\in\mathbb{N}$

$\|\bar{\boldsymbol{\xi}}_{p}^{j}\|\leq C_{9}$

###### Proof.

For the sake of clarity of notation, we do not include the index $j$ in the
proof.

To bound $\|\bar{\boldsymbol{\xi}}_{p}\|$, consider the following estimate

$\|\bar{\boldsymbol{\xi}}_{p}\|\leq\|{\boldsymbol{\xi}}_{p}\|+\|A\|\Big{(}\|(\Theta^{-1}+\rho
I)^{-1}X^{-1}{\boldsymbol{\xi}}_{\mu,\sigma}\|+\|(\Theta^{-1}+\rho
I)^{-1}\|\|{\boldsymbol{\xi}}_{d}\|\Big{)}.$

We already know the following estimates

$\|{\boldsymbol{\xi}}_{p}\|\leq\frac{\mu
n}{\gamma_{p}}\leq\frac{\mu^{0}n}{\gamma_{p}},\qquad\|{\boldsymbol{\xi}}_{d}\|\leq\frac{\mu
n}{\gamma_{d}}\leq\frac{\mu^{0}n}{\gamma_{d}},\qquad\|(\Theta^{-1}+\rho
I)^{-1}\|\leq\frac{1}{\rho}.$

To estimate $\|(\Theta^{-1}+\rho
I)^{-1}X^{-1}{\boldsymbol{\xi}}_{\mu,\sigma}\|$, we proceed as in (35):

$\displaystyle\|(\Theta^{-1}+\rho
I)^{-1}X^{-1}{\boldsymbol{\xi}}_{\mu,\sigma}\|=\|(\Theta^{-1}+\rho
I)^{-1}(S\mathbf{e}-\sigma\mu X^{-1}\mathbf{e})\|\leq$
$\displaystyle\leq\|(\Theta^{-1}+\rho I)^{-1/2}\|\|(\Theta^{-1}+\rho
I)^{-1/2}X^{-1/2}S^{1/2}\|\big{(}\|X^{1/2}S^{1/2}\mathbf{e}\|+\sigma\mu\|X^{-1/2}S^{-1/2}\mathbf{e}\|\big{)}.$

It is straightforward to prove that

$\|(\Theta^{-1}+\rho
I)^{-1/2}\|\leq\frac{1}{\rho^{1/2}},\qquad\|(\Theta^{-1}+\rho
I)^{-1/2}X^{-1/2}S^{1/2}\|\leq 1.$

The remaining terms can be bounded using the properties of the neighbourhood

$\|X^{1/2}S^{1/2}\mathbf{e}\|\leq\sqrt{\mu\bar{\gamma}n},\quad\sigma\mu\|X^{-1/2}S^{-1/2}\mathbf{e}\|\leq\sigma\sqrt{\frac{\mu
n}{\underline{\gamma}}}.$

Since $\mu\leq\mu^{0}$, we deduce that $\|\bar{\boldsymbol{\xi}}_{p}\|\leq
C_{9}$, for some positive constant $C_{9}$. ∎

Finally, we show that, for a small enough constant $C_{t}$, the inexactness
introduced by the sparsification strategy satisfies the Assumption (16).
Therefore, an algorithm that includes such a sparsification strategy retains
the polynomial complexity of the inexact IPM shown in Section 3.

######  Corollary 4.

If in Algorithm 2 we generate the search directions using
$(S^{C_{t}\mu}_{\rho,\delta})^{-1}$ with $C_{t}$ sufficiently small, i.e. if
we compute the search directions using (11), (12) and (13) where
$S^{C_{t}\mu}_{\rho,\delta}$ substitutes $S_{\rho,\delta}$, then Algorithm 2
is convergent.

###### Proof.

Using Theorem 7, we have

$(S^{C_{t}\mu}_{\rho,\delta})^{-1}\bar{\boldsymbol{\xi}}_{p}=S^{-1}_{\rho,\delta}\bar{\boldsymbol{\xi}}_{p}+\boldsymbol{\boldsymbol{\psi}}$

where $\|\boldsymbol{\boldsymbol{\psi}}\|\leq
C_{10}\frac{C_{t}\mu}{\delta^{2}(1+\rho\mu)}\|\bar{\boldsymbol{\xi}}_{p}\|$
for some constant $C_{10}>0$. Hence

$S_{\rho,\delta}\underbrace{(S^{C_{t}\mu}_{\rho,\delta})^{-1}\bar{\boldsymbol{\xi}}_{p}}_{=\Delta\mathbf{y}}=\bar{\boldsymbol{\xi}}_{p}+S_{\rho,\delta}\boldsymbol{\boldsymbol{\psi}}.$

Recall (16) and Lemma 5; the thesis follows observing that there exists a
constant $C_{11}>0$ s.t.

$\|S_{\rho,\delta}\boldsymbol{\boldsymbol{\psi}}\|\leq
C_{11}\frac{C_{t}\mu}{\delta^{2}(1+\rho\mu)}\leq
C_{\text{inexact}}\mathbf{x}^{T}\mathbf{s}$

where the last inequality holds if

$C_{t}<\frac{\delta^{2}(1+\rho\mu)nC_{\text{inexact}}}{C_{11}}.$

∎

##  6 Numerical Results

The proposed method is compared with Lemon (Library for Efficient Modelling
and Optimization on Networks) [26], an extremely efficient C++ implementation
of the network simplex method, that has been shown to significantly outperform
other popular implementations, like Cplex, see e.g. [5, 48]. The network
simplex method has been shown [41] to be very competitive against other
algorithms specifically developed for discrete OT, while remaining very robust
and adaptable to many types of problems. Let us highlight that the other
algorithms available in Lemon (cost scaling, capacity scaling, cycle
cancelling) produced worse results than the network simplex.

All the computational tests discussed in this section are performed using a
Dell PowerEdge R740 running Scientific Linux 7 with $4\times$ Intel Gold 6234
3.3G, 8C/16T, 10.4GT/s, 24.75M Cache, Turbo, HT (130W) DDR4-2933, with 500GB
of memory. The PS-IPM implementation closely follows the one from [6] and is
written in Matlab®. The software versions used for the numerical experiments
are as follows: Matlab R2022a, Lemon 1.3.1 and GCC 4.8.5 as the C++ compiler.

We stop Algorithm 1, when

$\|\mathbf{g}-A^{T}\mathbf{y}-\mathbf{s}\|_{\infty}\leq\;R\cdot
tol\;\wedge\;{\|\mathbf{b}-A\mathbf{x}\|_{1}}\leq\;R\cdot
tol\;\wedge\;C_{\mathbf{x},\mathbf{s}}\leq\;tol,$ (42)

where

$tol=10^{-10},\,\,\,\,R:=\max\\{\|A\|_{\infty},\|\mathbf{b}\|_{1},\|\mathbf{c}\|_{1}\\},$

and

$C_{\mathbf{x},\mathbf{s}}:=\max_{i}\\{\min\\{|(\mathbf{x}_{i}\mathbf{s}_{i})|,|\mathbf{x}_{i}|,|\mathbf{s}_{i}|\\}\\}.$

Concerning the choice of the parameters in Algorithm 1, we set
$\sigma_{r}=0.7$. Moreover, to prevent wasting time on finding excessively
accurate solutions in the early PPM sub-problems, we set $\tau_{1}=10^{-4}$,
i.e. we use as inexactness criterion for the PPM method

$\|\mathbf{r}_{k}(\mathbf{x}_{k+1},\mathbf{y}_{k+1}))\|<10^{4}\sigma_{r}^{k}\min\\{1,\|(\mathbf{x}_{k+1},\mathbf{y}_{k+1})-(\mathbf{x}_{k},\mathbf{y}_{k})\|.$

Indeed, in our computational experience, we have found that driving the IPM
solver to a high accuracy in the initial PPM iterations is unnecessary and
usually leads to a significant deterioration of the overall performance.

Concerning Algorithm 2, we set as regularization parameters $\rho=10^{-4}$ and
$\delta=10^{-6}$. Moreover, in order to find the search direction, we employ a
widely used predictor-corrector method [30]. This issue represents the main
point where practical implementation deviates from the theory in order to gain
computational efficiency.

Finally, concerning the test problems, in all the following experiments we
generate the load vector $\boldsymbol{\rho_{1}-\rho_{0}}$ in (1) randomly and
such that the sum of its entries is zero (to guarantee feasibility of the
optimization problem), with only $10\%$ of them being nonzeros. Moreover, we
fix the weight of each edge at $1$.

###  6.1 Analysis of the sparsification strategy

In this section, we compare three possible solution strategies inside the PS-
IPM: Cholesky factorization (using Matlab’s chol function) applied to the full
normal equations matrix (14); Cholesky factorization (always using Matlab’s
chol function) applied to the sparsified matrix (41); preconditioned conjugate
gradient (PCG) (using Matlab’s pcg function) applied to the sparsified matrix
(41) with incomplete Cholesky preconditioner (computed using Matlab’s ichol
function). More in particular, as sparsification parameter in (40) we use
$C_{t}=0.4$, $\texttt{`droptol'}=10^{-3}$ in ichol and
$\texttt{`tol'}=10^{-1}\mu$ in pcg.

We test the above mentioned solution strategies on various instances generated
with the CONTEST generator [43]; in particular, we considered the graphs pref,
kleinberg, smallw and erdrey, with a fixed number of $100,000$ nodes and
different densities (i.e. average number of edges per node). Therefore, for
these instances, $m=|V|=100,000$ and $n=|E|=m\cdot\text{density}$.

In the upper panels of Figures 1 and 2 we report the computational time of the
three approaches for various values of densities (chosen in relation to the
properties of the graph), whereas in the lower panels we report the total
number of IPM iterations. From the presented numerical results, it is clear
that the sparsification strategy, in conjunction with the iterative solution
of the linear systems, provides a clear advantage over the use of a direct
factorization. As can be expected, the iterative method and the sparsification
strategy become more advantageous when the size of the problem (number of
edges) increases. On the other hand, it is important to note that the use of
the sparsified Newton equations in conjunction with the full Cholesky
factorization presents only limited advantages in terms of computational time
when compared to the Cholesky factorization of the full Newton normal
equation. This is the case because the resulting inexact IPM requires,
generally, more iterations to converge (see lower panels of Figures 1 and 2).
Advantages of the proposed approach become clearer when the graphs are denser.

Figure 1: Comparison of sparsified and full normal equations approach, using
full Cholesky factorization or incomplete Cholesky as preconditioner for PCG,
in terms of IPM iterations and computational time, for problems pref and
kleinberg. Figure 2: Comparison of sparsified and full normal equations
approach, using full Cholesky factorization or incomplete Cholesky as
preconditioner for PCG, in terms of IPM iterations and computational time, for
problems erdrey and smallw.

###  6.2 Results on randomly generated graphs

In this section, we compare the PS-IPM algorithm, using the sparsified normal
equations matrix and the PCG, with the network simplex solver of Lemon. For
PS-IPM we use the same parameters as proposed in Section 6.1. The graphs used
in this section come from the generator developed in [44] and already used for
OT on graphs in [12]. This generator produces random connected graphs with a
number of nodes varying from $1,000$ to $10,000,000$ and degrees of each node
in the range $[1,10]$, with an average of $5$ edges per node. For each size,
$10$ graphs and load vectors are generated and tested. These parameters
closely resemble the ones used in [12].

Figure 3 shows the comparison of the computational time between PS-IPM and
Lemon: for each size of the problem (indicated by the total number of edges),
we report the summary statistics of the execution times using Matlab’s
boxplot.

Figure 3: Box Plots of the computational times of PS-IPM and Lemon, for
randomly generated graphs. The red and black intervals show the spread of the
measured computational times (red, on the left, is PS-IPM; black, on the
right, is Lemon), while the (blue) crosses indicate the outliers.

For small size problems, Lemon is the clear winner, by two orders of
magnitude; however, as the size increases, the performance difference between
the two methods reduces and for the largest instance considered, Lemon becomes
one order of magnitude slower than PS-IPM.

Figure 4 shows the average computational time against the number of edges
(from $5,000$ to $50M$) in a logarithmic scale, the corresponding regression
lines and their slopes. Both the proposed method and the network simplex (see
[35]) are known to have polynomial complexity in terms of number of iterations
(although the estimates are usually very pessimistic for IPMs); from the
computational results presented, we can estimate the practical time complexity
of both methods. Recall that, in a log-log plot, polynomials of the type
$x^{m}$ appear as straight lines with slope $m$. Using linear regression, we
can estimate that the time taken by Lemon grows with exponent approximately
$2.06$, while the time taken by PS-IPM grows with exponent approximately
$1.28$, providing a considerable advantage for large sizes.

Looking at the full set of results as reported in Figure 3, let us mention
finally the fact that the variance of the computational times over the 10 runs
for a given problem size is smaller when using PS-IPM, especially for large
sizes, indicating that the method is more robust and less dependent on the
specific problem being solved. This is a very desirable property.

Figure 4: Logarithmic plot of the computational time for randomly generated
graphs. The time taken by Lemon (red circles) grows as $(\text{number of
edges})^{2.06}$; the time taken by PS-IPM (blue triangles) grows as
$(\text{number of edges})^{1.28}$.

###  6.3 Results on SuiteSparse graphs

Results on randomly generated problems do not necessarily represent the
ability of an optimization method to tackle problems coming from real world
applications. Therefore, in this section, we show the results of applying PS-
IPM and Lemon to some sparse graphs from the SuiteSparse matrix collection
[10]. The characteristics of the graphs considered are shown in Table 1: the
number of nodes, edges and the average number of edges per node. All the
graphs are undirected and connected. Due to the fact that the considered
graphs are particularly sparse, in the numerical results presented in this
section, we solve the sparsified normal equations using the full Cholesky
factorization.

Table 1: Details of the graphs from the SuiteSparse matrix collection, ordered by increasing number of edges. Name | Nodes | Edges | Density
---|---|---|---
nc2010 | 288,987 | 1,416,620 | 4.9
NACA0015 | 1,039,183 | 6,229,636 | 6.0
great-britain-osm | 7,733,822 | 16,313,034 | 2.1
hugetric-00010 | 6,592,765 | 19,771,708 | 3.0
hugetric-00020 | 7,122,792 | 21,361,554 | 3.0
hugetrace-00010 | 12,057,441 | 36,164,358 | 3.0
hugetrace-00020 | 16,002,413 | 47,997,626 | 3.0
delaunay-n23 | 8,388,608 | 50,331,568 | 6.0

Figure 5 shows the computational times for the eight problems considered,
using PS-IPM and Lemon. Apart from the problem nc2010, which represents a
relatively small instance in out dataset, on all the other problems PS-IPM
consistently outperforms Lemon in terms of required computational time. In
particular, for the problems hugetric-00010, hugetric-00020, hugetrace-00010
and hugetrace-00020, which reach up to 16 million nodes and 48 million edges,
PS-IPM is one order of magnitude faster than Lemon.

Notice that graphs of these sizes (and larger) appear in many modern practical
applications, e.g. social networks, PageRank, analysis of rail/road networks,
energy models, to mention a few.

Looking at the regression lines and their slopes, we notice that the time
taken by Lemon grows with exponent approximately $2.07$ while the time taken
by PS-IPM grows with exponent approximately $1.40$. These values are very
close to the ones found previously for randomly generated graphs. The data of
Figure 5 however has a more erratic behaviour than the times shown in Figure
4, because the properties of each graph considered are different and because
we are not averaging over $10$ different instances of each problem.

Let us highlight also that the time taken by Lemon seems to be more problem
dependent, while PS-IPM looks more consistent and robust.

Figure 5: Logarithmic plot of the computational time for the SuiteSparse
problems. The time taken by Lemon (red circles) grows as $(\text{number of
edges})^{2.07}$; the time taken by PS-IPM (blue triangles) grows as
$(\text{number of edges})^{1.40}$. The instances are ordered as in Table 1.

##  7 Conclusion

An efficient computational framework for the solution of Optimal Transport
problems on graphs has been presented in this paper. Such framework relies on
Proximal-Stabilized Interior Point Method and clever sparsifications of the
normal Newton equations to compute the inexact search directions. The proposed
technique is sound and polynomial convergence guarantee has been established
for the inner inexact IPM. Extensive numerical experiments show that for large
scale problems, a simple prototype Matlab implementation is able to outperform
consistently a highly specialized and very efficient C++ implementation of the
network simplex method.

We highlight also that Interior Point Methods are more easily parallelizable
than simplex-like methods; for huge scale problems, for which high performance
computing resources need to be used, the use of IPMs with proper
parallelization may be the only viable strategy to solve these problems.

## References

  * [1] R. K. Ahuja, T. L. Magnanti, and J. B. Orlin, Network Flows, Prentice-Hall, Englewood Cliffs, NJ, 1993.
  * [2] A. Altman and J. Gondzio, Regularized symmetric indefinite systems in interior point methods for linear and quadratic optimization, Optim. Methods Softw., 11/12 (1999), pp. 275–302.
  * [3] P. Armand and J. Benoist, Uniform boundedness of the inverse of a Jacobian matrix arising in regularized interior-point methods, Math. Program., 137 (2013), pp. 587–592.
  * [4] S. Bellavia, Inexact interior-point method, J. Optim. Theory Appl., 96 (1998), pp. 109–121.
  * [5] J. Castro and S. Nasini, A specialized interior-point algorithm for huge minimum convex cost flows in bipartite networks, European J. Oper. Res., 290 (2021), pp. 857–869.
  * [6] S. Cipolla and J. Gondzio, Proximal stabilized interior point methods and low-frequency-update preconditioning techniques, J. Optim. Theory Appl., 197 (2023), p. 1061–1103.
  * [7] J. Cornelis and W. Vanroose, Convergence analysis of a regularized inexact interior-point method for linear programming problems, 2021.  10.48550/ARXIV.2105.01333[math.OC].
  * [8] M. Cuturi, Sinkhorn distances: Lightspeed computation of optimal transport, in Adv. in Neural Information Processing Systems, C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Weinberger, eds., vol. 26, Curran Associates, Inc., 2013.
  * [9] S. I. Daitch and D. A. Spielman, Faster approximate lossy generalized flow via interior point algorithms, in Proceedings of the Fortieth Annual ACM Symposium on Theory of Computing, STOC ’08, New York, NY, USA, 2008, Association for Computing Machinery, p. 451–460.
  * [10] T. A. Davis and Y. Hu, The University of Florida sparse matrix collection, ACM Trans. Math. Software, 38 (2011), pp. Art. 1, 25.
  * [11] R. S. Dembo, S. C. Eisenstat, and T. Steihaug, Inexact Newton methods, SIAM J. Numer. Anal., 19 (1982), pp. 400–408.
  * [12] M. Essid and J. Solomon, Quadratically regularized optimal transport on graphs, SIAM J. Sci. Comput., 40 (2018), pp. A1961–A1986.
  * [13] E. Facca and M. Benzi, Fast iterative solution of the optimal transport problem on graphs, SIAM J. Sci. Comput., 43 (2021), pp. A2295–A2319.
  * [14] A. Frangioni and C. Gentile, New preconditioners for KKT systems of network flow problems, SIAM J. Optim., 14 (2004), pp. 894–913.
  * [15] A. Frangioni and S. Serra Capizzano, Spectral analysis of (sequences of) graph matrices, SIAM J. Matrix Anal. Appl., 23 (2001), pp. 339–348.
  * [16] M. P. Friedlander and D. Orban, A primal-dual regularized interior-point method for convex quadratic programs, Math. Program. Comput., 4 (2012), pp. 71–107.
  * [17] J. Gondzio, Interior point methods 25 years later, European J. Oper. Res., 218 (2012), pp. 587–601.
  * [18] , Matrix-free interior point method, Comput. Optim. Appl., 51 (2012), pp. 457–480.
  * [19] C. Gottschlich and D. Schuhmacher, The shortlist method for fast computation of the Earth mover’s distance and finding optimal solutions to transportation problems, PLoS ONE, 9 (2014), p. e110214.
  * [20] I. Gulrajani, F. Ahmed, M. Arjovsky, V. Dumoulin, and A. Courville, Improved training of Wasserstein GANs, in Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, Red Hook, NY, USA, 2017, Curran Associates Inc., p. 5769–5779.
  * [21] S. Haker, L. Zhu, A. Tannenbaum, and S. Angenent, Optimal mass transport for registration and warping, Int. J. Comput. Vis, 60 (2004), pp. 225–240.
  * [22] L. Kantorovich, On the translocation of masses, Manag. Sci., 5 (1958), pp. 1–4.
  * [23] C. T. Kelly, Iterative Methods for Linear and Nonlinear Equations, SIAM, Philadelphia, USA, 1995.
  * [24] M. Kojima, N. Megiddo, and S. Mizuno, A primal-dual infeasible interior-point algorithm for linear programming, Math. Program., 61 (1993), pp. 263–280.
  * [25] J. Korzak, Convergence analysis of inexact infeasible-interior-point algorithms for solving linear programming problems, SIAM J. Optim., 11 (2000), pp. 133–148.
  * [26] P. Kovacs, Minimum-cost flow algorithms: an experimental evaluation, Optim. Methods Softw., 30 (2015), pp. 94–127.
  * [27] H. Ling and K. Okada, An efficient Earth mover’s distance algorithm for robust histogram comparison, IEEE Trans. Pattern Anal. Mach. Intell., 29 (2007), pp. 840–853.
  * [28] F. J. Luque, Asymptotic convergence analysis of the proximal point algorithm, SIAM J. Control Optim., 22 (1984), pp. 277–293.
  * [29] T. A. Manteuffel, An incomplete factorization technique for positive definite liner systems, Math. Comp., 34 (1980), pp. 473–497.
  * [30] S. Mehrotra, On the implementation of a primal-dual interior-point methods, SIAM J. Optim., 2 (1992), pp. 575–601.
  * [31] Q. Merigot, A multiscale approach to optimal transport, Computer Graphics Forum, 30 (2011), pp. 1583–1592.
  * [32] L. Métivier, R. Brossier, Q. Mérigot, E. Oudet, and J. Virieux, An optimal transport approach for seismic tomography: application to 3d full waveform inversion, Inverse Problems, 32 (2016), p. 115008.
  * [33] A. Natale and G. Todeschi, Computation of optimal transport with finite volumes, ESAIM: M2AN, 55 (2021), pp. 1847–1871.
  * [34] C.-C. Ni, Y.-Y. Lin, F. Luo, and J. Gao, Community detection on networks with Ricci flow, Scientific Reports, 9 (2019), p. 9984.
  * [35] J. B. Orlin, A polynomial time primal network simplex algorithm for minimum cost flows, Math. Program., 78 (1997), pp. 109–129.
  * [36] G. Peyre and M. Cuturi, Computational optimal transport: with applications to data science, Found. Trends and Trends in Machine Learning, 11 (2019), pp. 355–607.
  * [37] L. F. Portugal, M. G. C. Resende, G. Veiga, and J. J. Júdice, A truncated primal-infeasible dual-feasible network interior point method, Networks, 35 (2000), pp. 91–108.
  * [38] M. G. Resende and G. Veiga, An implementation of the dual affine scaling algorithm for minimum-cost flow on bipartite uncapacitated networks, SIAM J. Optim., 3 (1993), pp. 516–537.
  * [39] R. T. Rockafellar, Monotone operators associated with saddle-functions and minimax problems, in Nonlinear Functional Analysis (Proc. Sympos. Pure Math., Vol. XVIII, Part 1, Chicago, Ill., 1968), Amer. Math. Soc., Providence, R.I., 1970, pp. 241–250.
  * [40] M. Saunders and T. J. A., Solving regularized linear programs using barrier methods and KKT systems, Technical Report SOL 96-4, Systems Optimization Laboratory, Dept. of Operations Research, Stanford University, Stanford, CA 94305, USA,, (1996).
  * [41] J. Schrieber, D. Schuhmacher, and C. Gottschlich, DOTmark - a benchmark for discrete optimal transport, IEEE Access, 5 (2017), pp. 271–282.
  * [42] D. A. Spielman and S.-H. Teng, Nearly linear time algorithms for preconditioning and solving symmetric, diagonally dominant linear systems, SIAM J. Matrix Anal. Appl., 35 (2014), pp. 835–885.
  * [43] A. Taylor and D. J. Higham, CONTEST: A controllable test matrix toolbox for Matlab, ACM Trans. Math. Softw., 35 (2009).
  * [44] F. Viger and M. Latapy, Efficient and simple generation of random simple connected graphs with prescribed degree sequence, in International Computing and Combinatorics Conference, Springer, 2005, pp. 440–449.
  * [45] N. K. Vishnoi, $Lx=b$ Laplacian solvers and their algorithmic applications, Found. Trends Theor. Comput. Sci., 8 (2012), pp. front matter, 1–141 (2013).
  * [46] J. Wijesinghe and P. Chen, Matrix-free interior point methods for point set matching problems, 2022.  arXiv:2202.09763[math.OC].
  * [47] S. J. Wright, Primal-dual interior-point methods, Society for Industrial and Applied Mathematics, Philadelphia, PA, 1997.
  * [48] F. Zanetti and J. Gondzio, An interior-point-inspired algorithm for linear programs arising in discrete optimal transport, INFORMS J. Comput., (2023).