# A First-Order Numerical Algorithm
without Matrix Operations

Muhammad Adil, Ramtin Madani, Sasan Tavakkol, and Ali Davoudi Muhammad Adil is
with Palo Alto Research Center (PARC), Ramtin Madani, and Ali Davoudi are with
the University of Texas at Arlington. Sasan Tavakkol is with Google Research.
Ramtin Madani<EMAIL_ADDRESS>is the corresponding author. This work
is funded, in part, by the Office of Naval Research under award
N00014-18-1-2186, and approved for public release under DCN# 43-9295-22.

###### Abstract

This paper offers a matrix-free first-order numerical method to solve large-
scale conic optimization problems. Solving systems of linear equations pose
the most computationally challenging part in both first-order and second-order
numerical algorithms. Existing direct and indirect methods are either
computationally expensive or compromise on solution accuracy. Alternatively,
we propose an easy-to-compute decomposition method to solve sparse linear
systems that arise in conic optimization problems. Its iterations are
tractable, highly parallelizable, with closed-form solutions. This algorithm
can be easily implemented on distributed platforms, such as graphics
processing units, with orders-of-magnitude time improvement. The performance
of the proposed solver is demonstrated on large-scale conic optimization
problems and is compared with the state-of-the-art first-order solvers.

###### Index Terms:

Cone programming, numerical algorithms, optimization, parallel computing.

##  I Introduction

Conic optimization is used in various areas such as operation research,
machine learning, signal processing, and optimal control. Common solutions to
conic optimization problems are based on the interior point method (IPM) [1,
2, 3, 4, 5, 6], that are suitable for small- to medium-sized problems. At each
iteration, IPM solves a linear system of equations, mainly using Gauss-Jordan
[7], Gaussian elimination [8], LU decomposition [9], Cholesky decomposition
[10], QR decomposition, or Monte Carlo Methods [11]. These direct approaches
become prohibitively expensive by scale, and, rendering IPM-based methods
impractical for large-scale problems.

Matrix-free interior-point methods, including indirect or iterative methods,
are among the most popular methods for solving large-scale problem [12, 13].
However, direct methods become prohibitive [12, 14] at larger scale and Krylov
subspace iterative methods, such as preconditioned conjugate gradient (PCG),
become attractive alternatives [14]. In conjugate gradient method [15], an
iterative technique solves the Newton step instead of factorizing the Hessian
matrix directly. [16] proposed a Lagrangian dual predictor-corrector
algorithm, and applied the conjugate gradient method to solving linear
systems. In [17], an iterative solver is applied to the modified barrier
method for large-scale semidefinite programming optimization problems. The
conjugate gradient method, with a simple preconditioner, reduces the
computation time of solving semidefinite programming problems [18], and an
inexact semismooth Newton conjugate gradient approach in [19] improves the
solution accuracy.

The performance of iterative methods depends on the spectral properties of the
linear system and the condition number of the matrix involved [20]. Existing
preconditioners aim at reducing the condition number to the extend possible
[21]. A matrix-free algorithm for equality-constraint nonlinear programming in
[22] remedies ill-conditioned problems with rank-deficient Jacobian matrices.
A matrix-free solver PDCO [23] uses least squares minimal residual (LSMR) to
solve linear systems. [24] proposes a matrix-free IPM for quadratic programs,
where the Karush–Kuhn–Tucker (KKT) system is regularized to bound the
condition number and, then, a preconditioner is designed for the regularized
system. This approach substantially decreases the computational cost in each
iteration, albeit by trading off more iterations and lower accuracy.

Collectively, IPMs are inherently computationally expensive and hard to scale
for larger conic optimization problems. Alternatively, first-order methods
scale gracefully with moderate accuracy [25, 26, 27, 28]. The computational
complexity of each iteration is significantly less than IPM. They are suitable
when high accuracy is not needed [29, 30]. The convergence of first-order
methods in a limited number of iterations is an active area of research [31,
32, 33, 34, 35, 36, 37, 38, 39, 40].

First-order methods are required to solve a one time linear system [28].
Thanks to factorizing caching, in practice, the linear system is solved only
in the first iteration, factors are stored, and, then, reused in subsequent
iteration [28, 27, 30]. The most common factorizing caching approaches are LDL
and QR decompositions [41]. Direct methods or factorizing caching become
impractical at a very large scale, and point matrix-free or indirect
approaches become viable [27]. The first-order methods [41, 25, 27, 28] apply
conjugate gradient method to the resulting large-scale problems. First-order
methods struggle with accuracy, and matrix-free indirect conjugate gradient
approaches further hamper their ability. The computational complexity of
direct methods and low accuracy solution of indirect methods to tackle large
conic problems are the motivations behind the proposed matrix-free approach
for solving very large problems with a modest accuracy.

###  I-A Contributions

We develop a matrix-free first-order numerical algorithm to solve very large-
scale sparse conic optimization problems. The basic idea is to decompose the
constraint matrix into sparse factors, such that iterative steps do not
involve any matrix operations. These factors are easy-to-compute and require
minimal memory storage. The matrix inversion lemma makes the operations matrix
division free. We reformulate the standard conic optimization problem by
introducing auxiliary variables and, then, apply the proposed matrix-free
algorithm in conjunction with the well-known two-block alternating direction
method of multipliers (ADMM) [29, 42, 43]. Therefore, the computational burden
of solving the linear system in each iteration is taken out of the iterative
loop. The proposed algorithm admits parallel implementation, and is amenable
to graphics processing units (GPUs). We demonstrate the performance gain and
computational speedup of the proposed algorithm by conducting a range of
experiments and compare the results with other first-order numerical solvers.

###  I-B Paper Structure

The rest of this paper is organized as follows. We introduce the cone
programming and a brief description of operator splitting methods for such
problems in Section II. We analyze the existing direct and indirect methods
for solving the linear system, and provide a motivation for our numerical
algorithm in Section III. Section IV presents the proposed algorithm.
Numerical experiments and comparison with competing solvers are presented in
V. Section VI concludes the paper.

###  I-C Notations

Symbols $\mathbb{R}$ and $\mathbb{N}$ denote the set of real and natural
numbers, respectively. Matrices and vectors are represented by bold uppercase
and bold lowercase letters, respectively. Notation $\lVert\cdot{\rVert}_{2}$
refers to $\ell_{2}$ norm of either matrix or vector depending on the context,
and $\lvert\cdot\rvert$ represents the absolute value. The symbol
$\\!(\cdot)^{\\!\top}\\!$ represent the transpose operators. The notations
$\bm{I}_{n}$ refer to the $n\times n$ identity matrix. The symbol
$\mathcal{K}$ describes different types of cones. The superscript
$\\!(\cdot)^{\\!\mathrm{opt}}\\!$ refers to the optimal solution of
optimization problem. $o$ denotes the number of nonzero entries in matrix. The
symbol $\mathcal{L}$ represent the augmented Lagrangian function. The
notations “$\mathcal{P}_{1}$”, “$\mathcal{P}_{2}$”, represent the primal,
whereas “$\mathcal{D}$” denote the dual block of two-block ADMM. The symbols
$\varepsilon_{\rm prim},\varepsilon_{\rm dual},\varepsilon_{\mathrm{abs}}$,
and $\varepsilon_{\mathrm{rel}}$ are used for primal, dual, absolute, and
relative tolerance, respectively. $\varepsilon_{\rm gap}$ is the tolerance for
difference between primal and dual objective values.

##  II Preliminaries

We consider the class of convex optimization problems with linear objective
function, linear constraints, and conic constraints, in the form of:

$\displaystyle\underset{\begin{subarray}{c}\\!\\!\\!\\!\\!\\!\\!\\!\bm{x}\in\,\mathbb{R}^{n}\end{subarray}}{\text{minimize~{}~{}~{}}}$
$\displaystyle\bm{c}^{\top}\bm{x}$ (1a) subject to
$\displaystyle\bm{A}\bm{x}=\bm{b}$ (1b) $\displaystyle\bm{x}\in\mathcal{K}$
(1c)

where $\bm{x}\in\mathbb{R}^{n}$ is the primal decision variable and
$\bm{c}\in\mathbb{R}^{n}$, $\bm{A}\in\mathbb{R}^{m\times n}$, and
$\bm{b}\in\mathbb{R}^{m}$ are given. Additionally,
$\mathcal{K}\triangleq\mathcal{K}_{n_{1}}\times\mathcal{K}_{n_{2}}\times\cdots\times\mathcal{K}_{n_{k}}\subseteq\mathbb{R}^{n}$
is a non-empty, closed, convex cone, where each
$\mathcal{K}_{n_{i}}\subseteq\mathbb{R}^{n_{i}}$ is a Lorentz cone of size
$n_{i}$, i.e.,

$\displaystyle\mathcal{K}_{n_{i}}\triangleq\big{\\{}\bm{w}\in\mathbb{R}^{n_{i}}\,|\,w_{1}\geq\big{\|}[w_{2},\ldots,w_{n_{i}}]\big{\|}_{2}\big{\\}},$

and $n_{1}+n_{2}+\ldots+n_{k}=n$.

In order to solve 1, various first-order operator splitting methods have been
proposed in the past decade. First-order methods are particularly interesting
for the cases where iterative steps can be solved efficiently through explicit
formula, and a large number of iterations can be executed in a short amount of
time. Among most popular methods are the Douglas-Rachford Splitting (DRS)
technique and ADMM.

####  II-1 Douglas-Rachford Splitting

DRS was originally proposed in [44] to find numerical solutions of
differential equations for heat conduction problems, and it has been widely
used to solve separable convex optimization problems. Rather than operating on
the whole problem directly, DRS works on a splitting scheme to address each
component of the problem separately. In order to implement the DR splitting
method, one casts the problem (1) in the form of

minimize $\displaystyle f(\bm{x})+g(\bm{x})$ (2)

where $f,g:\mathbb{R}^{n}\to\mathbb{R}\cup\\{\infty\\}$ are the indicator
functions:

$\displaystyle
f(\bm{x})\\!\triangleq\\!\left\\{\begin{matrix}0&\\!\\!\mathrm{if}\
\bm{x}\in\mathcal{K}\\\ \infty&\\!\\!\mathrm{\
otherwise}\end{matrix}\right.\quad\mathrm{and}\quad
g(\bm{x})\\!\triangleq\\!\left\\{\begin{matrix}\bm{c}^{\top}\bm{x}&\\!\\!\mathrm{if}\
\bm{A}\bm{x}=\bm{b}\\\ \infty&\\!\\!\mathrm{otherwise}\end{matrix}\right.$

leading to the following iterative steps:

$\displaystyle\bm{x}\leftarrow\mathrm{prox}_{f}(\bm{z})$ (3a)
$\displaystyle\bm{z}\leftarrow\bm{z}+\mathrm{prox}_{\mu^{-1}g}(2\bm{x}-\bm{z})-\bm{x}$
(3b)

where $\mu$ is a fixed tuning parameter, and for every
$h:\mathbb{R}^{n}\to\mathbb{R}\cup\\{\infty\\}$ and
$\bm{w}_{0}\in\mathbb{R}^{n}$, the operator $\mathrm{prox}_{h}(\bm{w}_{0})$
returns the unique solution to the following optimization problem:

$\displaystyle\underset{\begin{subarray}{c}\\!\\!\\!\\!\\!\\!\\!\\!\bm{w}\in\mathbb{R}^{n}\end{subarray}}{\text{minimize~{}~{}~{}}}$
$\displaystyle h(\bm{w})+\frac{1}{2}\|\bm{w}-\bm{w}_{0}\|_{2}^{2}.$ (4)

Each iteration of DRS requires the evaluation of the proximal operators
$\mathrm{prox}_{f}(\cdot)$ and $\mathrm{prox}_{\mu^{-1}g}(\cdot)$. While the
evaluation of $\mathrm{prox}_{f}(\cdot)$ is parallelizable and enjoys a
closed-form solution, the evaluation of
$\mathrm{prox}_{\mu^{-1}g}(\bm{w}_{0})$ is the main bottleneck which requires
solving the following system of linear equations:

$\displaystyle\begin{bmatrix}\bm{I}_{n}&\bm{A}^{\top}\\\
\bm{A}&\bm{0}\end{bmatrix}\begin{bmatrix}\bm{w}\\\
\bm{\nu}\end{bmatrix}=\begin{bmatrix}\bm{w}_{0}-\mu^{-1}\bm{c}\\\
\bm{b}\end{bmatrix},$ (5)

where $\bm{w}=\mathrm{prox}_{\mu^{-1}g}(\bm{w}_{0})$.

####  II-2 Alternating Direction Method of Multipliers

ADMM is one of the most commonly used first-order methods for solving large-
scale optimization problems. ADMM can be analyzed as a special case of DRS, as
the former applied to the primal problem is equivalent to the latter applied
to the dual problem. A standard way of solving problem (1) via ADMM is through
the following reformulation

$\displaystyle\underset{\begin{subarray}{c}\\!\\!\\!\\!\\!\\!\\!\\!\bm{x}_{1},\bm{x}_{2}\in\mathbb{R}^{n}\end{subarray}}{\text{minimize~{}~{}~{}}}$
$\displaystyle f(\bm{x}_{1})+g(\bm{x}_{2})$ (6a) subject to
$\displaystyle\bm{x}_{1}=\bm{x}_{2}$ (6b)

leading to the steps

$\displaystyle\bm{x}_{1}\leftarrow\mathrm{prox}_{f}(\bm{x}_{2}-\mu^{-1}\bm{\lambda_{1}})$
(7a)
$\displaystyle\bm{x}_{2}\leftarrow\mathrm{prox}_{\mu^{-1}g}(\bm{x}_{1}+\mu^{-1}\bm{\lambda_{1}})$
(7b)
$\displaystyle\bm{\lambda}_{12}\leftarrow\bm{\lambda}_{12}+\mu(\bm{x}_{1}-\bm{x}_{2}).$
(7c)

where $\bm{x}_{1},\bm{x}_{2}$ are primal variables, $\bm{\lambda}_{12}$ is the
dual variable associated with (6b), and $\mu$ is a fixed tuning parameter.
Each iteration requires solving the system of linear equations (5).

##  III Solving the Linear System

While first-order methods are considered highly efficient for solving large-
scale problems with modest accuracy, evaluating the projection operator
$\mathrm{prox}_{\mu^{-1}g}$ at each iteration can become computationally
prohibitive. Depending on the structure and size of the constraint matrix
$\bm{A}$, there are different ways to solve the linear system (5). The common
approaches for solving (5) in several iterations are classified as direct and
indirect methods [27, 28].

####  III-1 Direct Methods

The linear system (5) can be solved by first factoring the matrix in (5) and,
then, performing forward or backward substitutions. The common practice is to
compute the factors at the beginning and, then, reuse them in subsequent
iterations. This process is known as factorization caching, which is widely
used by first-order solvers [27, 28, 25, 41]. When $\bm{A}$ is dense, the
common practice is to solve the linear system (5) by performing either
Cholesky or QR decompositions and, then, cache the factors for subsequent
iterations. The complexity in both approaches is $O(mn^{2})$ in first
iteration and $O(mn)$ in subsequent iterations. When $\bm{A}$ is sparse and
not too large, the direct sparse $\bm{L}\bm{D}\bm{L}^{\top}$ decomposition of
the matrix in (5) can be used with $\bm{L}$ being a lower triangular matrix
and $\bm{D}$ being a diagonal matrix [45]. In this case, the computational
complexity depends upon the number of nonzero entries of the factor $\bm{L}$.

####  III-2 Indirect Methods

Instead of solving (5) by factorization, indirect methods apply an iterative
procedure such as conjugate gradient (CG) [15, 46, 47] or LSQR [48, 49]. In
contrast to direct methods, indirect methods offer lower complexity and more
freedom for adaptive parameter selection since there is no factorization. In
each step of conjugate gradient, a vector is multiplied by $\bm{A}$ and
$\bm{A}^{\top}$, thus, the computational complexity of each step is
significantly lower.

##  IV UV Decomposition

In contrast to existing direct and indirect methods, herein, we propose an
efficient matrix-free decomposition to handle large-scale sparse optimization
problems of the form (1).

######  Definition 1.

Assume that $\bm{A}=\bm{U}\bm{V}^{\top}$, where $\bm{U}\in\mathbb{R}^{m\times
o}$ and $\bm{V}\in\mathbb{R}^{n\times o}$. The pair of matrices
$(\bm{U},\bm{V})$ are regarded as a UV decomposition of the matrix $\bm{A}$,
if $\bm{U}\bm{U}^{\top}$ and $\bm{V}\bm{V}^{\top}$ are both diagonal.

Such decomposition can be readily constructed as follows. Let
$\\{\bm{A}_{i_{k},j_{k}}\\}^{o}_{k=1}$ represent the non-zero elements of
$\bm{A}$ in an arbitrary order. The pair

$\displaystyle\bm{U}\triangleq\sum_{k=1}^{o}{\bm{A}_{i_{k},j_{k}}\bm{e}_{i_{k}}\dot{\bm{e}}_{k}^{\top}}$
(8a)
$\displaystyle\bm{V}\triangleq\sum_{k=1}^{o}{\dot{\bm{e}}_{k}^{\top}\ddot{\bm{e}}_{j_{k}}}$
(8b)

is a possible candidate for a UV decomposition of $\bm{A}$, where
$\\{\bm{e}_{k}\\}^{m}_{k=1}$, $\\{\dot{\bm{e}}_{k}\\}^{o}_{k=1}$, and
$\\{\ddot{\bm{e}}_{k}\\}^{n}_{k=1}$ represent the standard basis for
$\mathbb{R}^{m}$, $\mathbb{R}^{o}$, and $\mathbb{R}^{n}$, respectively. The
factor $\bm{U}$ in (8a) contains the nonzero elements of the matrix $\bm{A}$,
whereas $\bm{V}$ in (8b) encodes the locations. Next, we further elucidate
this decomposition by providing a simple illustrative example.

Example 1: Consider the following sparse matrix with $8$ non-zero elements:

$\displaystyle\bm{A}$ $\displaystyle:=\begin{bmatrix}1&0&4&6&8\\\ 0&0&5&0&0\\\
2&3&0&7&0\\\ \end{bmatrix}.$ (9a) The pair $\displaystyle\bm{U}$
$\displaystyle:=\begin{bmatrix}1&0&0&4&0&6&0&8\\\ 0&0&0&0&5&0&0&0\\\
0&2&3&0&0&0&7&0\\\ \end{bmatrix},$ (9b) $\displaystyle\bm{V}$
$\displaystyle:=\begin{bmatrix}1&1&0&0&0&0&0&0\\\ 0&0&1&0&0&0&0&0\\\
0&0&0&1&1&0&0&0\\\ 0&0&0&0&0&1&1&0\\\ 0&0&0&0&0&0&0&1\\\ \end{bmatrix}$ (9c)

satisfies $\bm{A}=\bm{U}\bm{V}^{\top}$. Additionally, each column of $\bm{U}$
and each column of $\bm{V}$ has exactly one non-zero element. Hence,
$(\bm{U},\bm{V})$ is a UV decomposition of the matrix $\bm{A}$.

Next, using the notion of UV decomposition, we seek to develop a matrix-free
numerical algorithm for solving problem 1. To this end, we reformulate the
problem as follows:

$\displaystyle\underset{\begin{subarray}{c}\\!\\!\\!\\!\\!\\!\\!\\!\bm{x}\in\,\mathbb{R}^{n}\\\
\\!\\!\\!\\!\\!\\!\\!\bm{y}\in\,\mathbb{R}^{o}\\\
\\!\\!\\!\\!\\!\\!\\!\\!\bm{z}\in\,\mathbb{R}^{n}\end{subarray}}{\text{minimize~{}~{}~{}}}$
$\displaystyle\bm{c}^{\top}\bm{x}$ (10a) subject to
$\displaystyle\bm{\lambda}:\;\bm{U}\bm{y}=\bm{b}$ (10b)
$\displaystyle\bm{\gamma}:\;\bm{y}=\bm{V}^{\top}\bm{x}$ (10c)
$\displaystyle\bm{\delta}:\;\bm{z}=\bm{x}$ (10d)
$\displaystyle\bm{z}\in\mathcal{K}\hskip 56.9055pt$ (10e)

where $\bm{y}\in\mathbb{R}^{o}$ and $\bm{z}\in\mathbb{R}^{n}$ are auxiliary
variables, and $(\bm{U},\bm{V})$ is a UV decomposition of the matrix $\bm{A}$.
It is straightforward to verify that the two formulations are equivalent. To
solve problem (10) via ADMM, we form the augmented Lagrangian function:

$\displaystyle\\!\\!\\!\mathcal{L}(\bm{x},\bm{y},\bm{z},\bm{\lambda},\bm{\gamma},\bm{\delta})\triangleq\bm{c}^{\top}\bm{x}$
$\displaystyle+\bm{\lambda}^{\top}(\bm{U}\bm{y}-b)+\dfrac{\mu}{2}\big{\|}\bm{U}\bm{y}-\bm{b}\big{\|}_{2}^{2}$
$\displaystyle+\bm{\gamma}^{\top}(\bm{y}\\!-\\!\bm{V}^{\\!\top}\bm{x})+\dfrac{\mu}{2}\big{\|}\bm{y}\\!-\\!\bm{V}^{\\!\top}\bm{x}\big{\|}_{2}^{2}$
$\displaystyle+\bm{\delta}^{\top}(\bm{z}-\bm{x})+\dfrac{\mu}{2}\big{\|}\bm{z}-\bm{x}\big{\|}_{2}^{2}$
(11)

where $\mu>0$ is a fixed parameter, and
$\bm{\lambda}\in\mathbb{R}^{m},\bm{\gamma}\in\mathbb{R}^{o}$, and
$\bm{\delta}\in\mathbb{R}^{n}$ are the Lagrange multipliers associated with
constraints 10b,10c and 10d, respectively.

We perform two block ADMM and regroup the primal and dual variables as follows

$\displaystyle\text{(Block 1)}\hskip 28.45274pt\mathcal{P}_{1}$
$\displaystyle=\\{\bm{x}\\}$ $\displaystyle\text{(Block 2)}\hskip
28.45274pt\mathcal{P}_{2}$ $\displaystyle=\\{\bm{y},\bm{z}\\}$
$\displaystyle\text{(Dual)}\hskip 39.83368pt\mathcal{D}$
$\displaystyle=\\{\bm{\lambda},\bm{\gamma},\bm{\delta}\\}$

Given the above partitioning of variables, each iteration of two block ADMM
for problem 10 involves the following steps:

####  IV-1 Block 1

This step consist of minimizing the Lagrangian function (11) with respect to
the variable $\bm{x}$, and freezing the other variables at their previous
values, i.e.,

$\displaystyle\bm{x}^{k+1}=\underset{\begin{subarray}{c}\\!\\!\\!\\!\\!\\!\\!\\!\bm{x}\in\mathbb{R}^{n}\end{subarray}}{\text{\rm{arg\
min~{}~{}~{}}}}\mathcal{L}(\bm{x},\bm{y}^{k},\bm{z}^{k},\bm{\lambda}^{k},\bm{\gamma}^{k},\bm{\delta}^{k})$

which enjoys the following closed-form solution

$\displaystyle\\!\\!\\!\\!\bm{x}^{k+1}\\!=\\!(\bm{I}\\!+\\!\bm{V}\bm{V}^{\top})^{-1}\Big{[}\bm{V}\big{(}\bm{y}^{k}\\!+\\!\frac{\bm{\gamma}^{k}}{\mu}\big{)}+\bm{z}^{k}+\frac{\bm{\delta}^{k}}{\mu}-\frac{\bm{c}}{\mu}\Big{]}.\\!\\!\\!$
(12)

Observe that $(\bm{I}+\bm{V}\bm{V}^{\top})$ is a diagonal positive definite
matrix and, once the factor $\bm{V}$ is known, computing the solution of 12
can be made matrix-free by storing the diagonal elements of
$(\bm{I}+\bm{V}\bm{V}^{\top})^{-1}$.

####  IV-2 Block 2

The next block consists of minimizing the Lagrangian function (11) with
respect to the variables $\bm{y}$ and $\bm{z}$, while $\bm{x}$ is fixed:

$\displaystyle(\bm{y}^{k+1},\bm{z}^{k+1})=\underset{\begin{subarray}{c}\\!\\!\\!\\!\\!\\!\\!\\!\bm{y}\in\mathbb{R}^{o}\\\
\\!\\!\\!\\!\\!\\!\\!\bm{z}\in\mathbb{R}^{n}\end{subarray}}{\text{\rm{arg\
min~{}~{}~{}}}}\mathcal{L}(\bm{x}^{k+1},\bm{y},\bm{z},\bm{\lambda}^{k},\bm{\gamma}^{k},\bm{\delta}^{k})$

which involves two parallel steps.

Minimization with respect to $\bm{y}$: leads to the following closed-form
solution:

$\displaystyle\\!\\!\\!\bm{y}^{k+1}\\!=\\!(\bm{I}\\!+\\!\bm{U}^{\top}\bm{U})^{-1}\Big{[}\bm{U}^{\top}\big{(}\bm{b}-\frac{\bm{\lambda}^{k}}{\mu}\big{)}\\!+\\!\bm{V}^{\top}\bm{x}^{k+1}-\frac{\bm{\gamma}^{k}}{\mu}\Big{]}.$
(13)

Now, according to the matrix inversion lemma, we have

$\displaystyle(\bm{I}+\bm{U}^{\top}\bm{U})^{-1}=(\bm{I}-\bm{U}^{\top}(\bm{I}+\bm{U}\bm{U}^{\top})^{-1}\bm{U})$
(14)

and since $(\bm{I}+\bm{U}\bm{U}^{\top})$ is a diagonal positive-definite
matrix, this step can be made matrix free as well, by storing the diagonal
elements of $(\bm{I}+\bm{U}\bm{U}^{\top})^{-1}$.

Minimization with respect to $\bm{z}$: involves projection onto the associated
Lorentz cones which is parallelizable and enjoys a closed-form solution as
well:

$\displaystyle\bm{z}^{k+1}=\mathrm{proj}_{\mathcal{K}}\Big{\\{}\bm{x}^{k+1}-\frac{\bm{\delta}^{k}}{\mu}\Big{\\}}.$
(15)

Observe that for each $i\in\\{1,\ldots,k\\}$,

$\displaystyle\mathrm{proj}_{\mathcal{K}_{i}}(\bm{w})\triangleq\left\\{\begin{matrix}\mathbf{0}^{n_{i}}&\\!\\!\alpha\leq-
w_{1}\\\ \bm{w}&\\!\\!\alpha\leq+w_{1}\\\
\frac{\bm{w}}{2}\\!+\\![\frac{\alpha}{2},\frac{w_{1}w_{2}}{2\alpha},\ldots,\frac{w_{1}w_{n_{i}}}{2\alpha}]&\\!\\!\mathrm{\
otherwise}\end{matrix}\right.$ (16)

where $\alpha=\|[w_{2},\ldots,w_{n_{i}}]\|_{2}$.

####  IV-3 Dual variables update

This steps involves the update of dual variable as follows:

$\displaystyle\bm{\lambda}^{k+1}$
$\displaystyle=\bm{\lambda}^{k}+\mu(\bm{U}\bm{y}^{k+1}-\bm{b})$ (17a)
$\displaystyle\bm{\gamma}^{k+1}$
$\displaystyle=\bm{\gamma}^{k}+\mu(\bm{y}^{k+1}-\bm{V}^{\top}\bm{x}^{k+1})$
(17b) $\displaystyle\bm{\delta}^{k+1}$
$\displaystyle=\bm{\delta}^{k}+\mu(\bm{z}^{k+1}-\bm{x}^{k+1})$ (17c)

Algorithm 1

1: $(\bm{A},\bm{b},\bm{c},\mathcal{K})$, fixed $\mu>0$, and initial points
$\bm{x},\bm{z},\bm{\delta}\in\mathbb{R}^{n},\bm{y},\bm{\gamma}\in\mathbb{R}^{o}$,
and $\bm{\lambda}\in\mathbb{R}^{m}$

2:Construct a UV decomposition of the matrix $\bm{A}$

3:$\bm{F}_{\bm{U}}:=\big{(}\bm{I}+\bm{U}\bm{U}^{\top}\big{)}^{-1}$

4:$\bm{F}_{\bm{V}}:=\big{(}\bm{I}+\bm{V}\bm{V}^{\top}\big{)}^{-1}$

5:repeat

6:
$\bm{x}\leftarrow\bm{F}_{\bm{V}}\Big{(}\bm{V}\big{(}\bm{y}+\frac{\bm{\gamma}}{\mu}\big{)}+\bm{z}+\frac{\bm{\delta}}{\mu}-\frac{\bm{c}}{\mu}\Big{)}$

7:
$\bm{y}\leftarrow(\bm{I}-\bm{U}^{\top}\bm{F}_{\bm{U}}\bm{U})\Big{(}\bm{U}^{\top}\big{(}\bm{b}-\frac{\bm{\lambda}}{\mu}\big{)}+\bm{V}^{\top}\bm{x}-\frac{\bm{\gamma}}{\mu}\Big{)}$

8:
$\bm{z}\leftarrow\mathrm{proj}_{\mathcal{K}}\Big{\\{}\bm{x}-\frac{\bm{\delta}}{\mu}\Big{\\}}$

9: $\bm{\lambda}\leftarrow\bm{\lambda}+\mu(\bm{U}\bm{y}-\bm{b})$

10: $\bm{\gamma}\leftarrow\bm{\gamma}+\mu(\bm{y}-\bm{V}^{\top}\bm{x})$

11: $\bm{\delta}\leftarrow\bm{\delta}+\mu(\bm{z}-\bm{x})$

12:until stopping criteria is met.

13: $\bm{x}^{\mathrm{opt}},\bm{\lambda}^{\mathrm{opt}}$

The above steps are summarized in Algorithm 1.

  * •

Steps 1: The equality constraint matrix $\bm{A}$ is decomposed into
$\bm{A}=\bm{U}\bm{V}^{\top}$ such that $\bm{U}\bm{U}^{\top}$ and
$\bm{V}\bm{V}^{\top}$ are both diagonal.

  * •

Steps 2 and 3: To make iterative steps matrix-free, we compute the
multiplication factors in these steps. These factors are easy-to-compute since
$(\bm{I}+\bm{U}\bm{U}^{\top})$ and $(\bm{I}+\bm{V}\bm{V}^{\top})$ are positive
definite diagonal matrices.

  * •

Steps 5 to 10: These steps are the ADMM primal and dual variables updates, as
well as the projection onto specific cones performed in Step 7.

##  V Numerical Experiments

We highlight the computational strength and scalability of the proposed
algorithm by testing it on a variety of randomly generated large-scale and
sparse linear programming (LP) and second-order cone programming (SOCP)
problems. We compare our methods with the first-order solvers POGS (Proximal
Graph Solver) [27], OSQP (Operator Splitting Solver for Quadratic Programs
[28], and SCS (Splitting Conic Solver) [25]. For each experiment, we stop
immediately after a better solution than the competing solver is obtained in
terms of residual norms, constraint violations, or objective gap. The
performance gain of the proposed algorithm is consistent across all problem
instances and different hardware architectures.

The proposed algorithm and competing solvers are implemented in MATLAB R2020a,
and all the experiments are carried out on a Linux-based DGX station with 20
2.2 GHz cores, Intel Xeon E5-2698 v4 CPU, with NVIDIA Tesla V100-DGXS-32GB
(128 GB total) GPU processor and 256 GB of RAM. The parallel nature of the
proposed algorithm helps us to take advantage of multi-core CPU. Note that our
implementation in MATLAB utilizes only a single GPU and does not benefit from
multiple GPUs. Moreover, no experiment is bounded by RAM or GPU memory of the
DGX station. We have used the MATLAB interface of POGS, OSQP v0.6.0 and SCS
v2.1.2.

Problem instances: All data is generated in such a way that the linear and
second order cone programming problems are feasible and bounded. The number of
nonzero entries of $\bm{A}$ are in the range of $10^{4}$ to $10^{7}$.

  * •

We generate $\bm{A}\in\mathbb{R}^{m\times n}$ to be a sparse random matrix
with $0.1\%$, $0.5\%$, and $1.0\%$ nonzero elements, which are drawn i.i.d.
(independently and identically distributed) from $\mathcal{N}(0,1)$.

  * •

$\bm{b}:=\bm{A}\times\mathrm{proj}_{\mathcal{K}}(\dot{\bm{x}})$, where the
elements of $\dot{\bm{x}}\in\mathbb{R}^{n}$ have i.i.d standard normal
distribution.

  * •

The elements of $\bm{c}\in\mathbb{R}^{n}$ have i.i.d standard normal
distribution.

  * •

For linear programming instances, we have $\mathcal{K}=\mathbb{R}_{+}^{n}$.
For SOCP instance, we have $\mathcal{K}=(\mathcal{K}_{4})^{\frac{n}{4}}$,
where $\mathcal{K}_{4}$ is the standard Lorentz cone of size $4$.

###  V-A Linear Programming

We consider randomly-generated linear programming problems, and compare the
performance of algorithm 1 with POGS and OSQP on a variety of sparse problems.
Experiments are continued until the run time of the competing solver reaches a
maximum time of 1200 seconds. The maximum time is chosen in such a way that
the experiments provide sufficient information to compare the computational
time for all solvers.

(a)  (b)

(c)

Figure 1: The performance of Algorithm 1 for linear programming in comparison
with OSQP for (a) $0.1\%$ (b) $0.5\%$ and (c) $1.0\%$ nonzero entries in
matrix $\bm{A}$.

(a)  (b)

(c)

Figure 2: The performance of Algorithm 1 for linear programming in comparison
with POGS for (a) $0.1\%$ (b) $0.5\%$ and (c) $1.0\%$ nonzero entries in
matrix $\bm{A}$.

####  V-A1 Comparison with OSQP

OSQP is a first-order general purpose open-source solver based on the
alternating direction method of multipliers. In Fig. 1, we compare our method
with OSQP by running several sparse linear programming instances. For each
instance, we use OSQP solver in its default settings and, then, provide these
parameters as input to our algorithm to achieve the same tolerance values.
Termination criteria: We stop when both primal and dual residuals are smaller
than some predefined tolerance limits $\varepsilon_{\rm prim}>0$ and
$\varepsilon_{\rm dual}>0$, i.e.,

$\displaystyle\|\bm{A}\bm{x}-\bm{b}\|_{\infty}<\varepsilon_{\rm prim}$ (18a)
$\displaystyle\|\bm{A}^{\top}\bm{\lambda}+\bm{c}\|_{\infty}<\varepsilon_{\rm
dual}$ (18b) $\displaystyle\varepsilon_{\rm prim}=\varepsilon_{\rm
abs}+\varepsilon_{\rm
rel}\max\\{\|\bm{A}\bm{x}\|_{\infty},\|\bm{b}\|_{\infty}\\}$ (18c)
$\displaystyle\varepsilon_{\rm dual}=\varepsilon_{\rm abs}+\varepsilon_{\rm
rel}\max\\{\|\bm{A}^{\top}\bm{\lambda}\|_{\infty},\|\bm{c}\|_{\infty}\\}$
(18d)

where $\varepsilon_{\rm abs}=10^{-4}$ and $\varepsilon_{\rm rel}=10^{-3}$ are
default absolute and relative tolerance values, respectively. We also satisfy
the following condition to terminate our algorithm

$\displaystyle|\bm{c}^{\top}\bm{x}+\bm{b}^{\top}\bm{\lambda}|<|\bm{c}^{\top}\bm{x}^{\rm
OSQP}+\bm{b}^{\top}\bm{y}^{\rm OSQP}|$

Figure 1a demonstrates that the CPU implementation of the proposed algorithm
is, at least, ten times faster than OSQP, whereas the GPU implementation shows
two orders-of-magnitude improvement. Similar results are depicted in Figures
1b and 2c, where we have used $0.5\%$ and $1.0\%$ nonzero values for the
constraint matrix $\bm{A}$. The CPU and GPU implementations of the proposed
algorithm can easily achieve ten times and hundred times improvements,
respectively, while satisfying similar, or even stricter, stopping criteria as
compared with OSQP.

####  V-A2 Comparison with POGS

We compare the performance of the proposed algorithm with POGS on both CPU and
GPU system architectures in figure 2. POGS is an open-source implementation of
graph projection splitting method that targets multi-core and GPU-based
systems for solving convex optimization problems. Herein, POGS would return
inaccurate solution in its default parameter settings. Hence, we made a slight
change in its default parameters by setting
$\varepsilon_{\mathrm{abs}}=10^{-5}$ and $\varepsilon_{\mathrm{rel}}=10^{-4}$
to obtain reasonable accuracy. The same tolerance parameters, along with
primal and dual solutions returned by POGS, are used as input parameters for
the proposed algorithm to meet the same stopping criteria.
Termination criteria: The stopping criteria of Algorithm 1 is when we exceed
both primal and dual feasibility of the solution produced by the competing
solver, i.e.,, when the following two criteria are met:

$\displaystyle\|\bm{A}\bm{x}-\bm{b}\|_{2}<\|\bm{A}\bm{x}^{\rm
POGS}-\bm{b}\|_{2}$ (19a)
$\displaystyle|\bm{c}^{\top}\bm{x}+\bm{b}^{\top}\bm{\lambda}|<|\bm{c}^{\top}\bm{x}^{\rm
POGS}+\bm{b}^{\top}\bm{\lambda}^{\rm POGS}|$ (19b)

where $\bm{x}^{\text{POGS}}$ and $\bm{\lambda}^{\text{POGS}}$ are primal and
dual solutions produced by the competing solver POGS under slightly-modified
default settings.

Figure 2a, 2b, and 2c, show notable time improvements offered by the proposed
algorithm in both cases. Our CPU and GPU implementation of Algorithm 1 shows
an order of magnitude time improvements compared to POGS, and this improvement
is further pronounced for larger problems.

###  V-B Second-order Cone Programming

We compare the performance of Algorithm 1 in with splitting conic solver (SCS)
on a variety of sparse conic problems.

####  V-B1 Comparison with SCS

Primarily written in C, SCS a first-order numerical optimization solver for
large-scale cone programs. This solver returns both primal and dual solutions
along with infeasibility certificate when applies.
Termination criteria: We use the following default stopping criteria of the
competing solver

$\displaystyle\|\bm{A}\bm{x}-\bm{b}\|_{2}\leq\varepsilon_{\rm
prim}(1+\|\bm{b}\|_{2}),$ (20a)
$\displaystyle\|\bm{A}^{\top}\bm{\lambda}+\bm{c}\|_{2}\leq\varepsilon_{\rm
dual}(1+\|\bm{c}\|_{2}),$ (20b)
$\displaystyle\bm{c}^{\top}\bm{x}+\bm{b}^{\top}\bm{\lambda}\leq\varepsilon_{\rm
gap}(1+|\bm{c}^{\top}\bm{x}|+|\bm{b}^{\top}\bm{\lambda}|)$ (20c)

where $\varepsilon_{\rm prim}=\varepsilon_{\rm dual}=\varepsilon_{\rm
gap}=10^{-3}$.

We have solved multiple SOCP problems for different sparse density and satisfy
the termination criteria of SCS solver and compare the computational time in
Figure 3.

(a)  (b)

(c)

Figure 3: The performance of Algorithm 1 for second order cone programming in
comparison with SCS direct for (a) $0.1\%$ (b) $0.5\%$ and (c) $1.0\%$ nonzero
entries in matrix $\bm{A}$.

##  VI Conclusion

We propose a computationally-efficient matrix-free first-order method for
solving large-scale sparse conic optimization problems. The computational
burden is managed by decomposing the equality constraint into sparse factors
that are easy-to-compute. A highly parallelizable and computationally-cheap
numerical algorithm is developed based on the alternating direction method of
multipliers. Each iteration of the proposed algorithm has a closed-form
solution with simple arithmetic operations. The parallel nature of the
proposed algorithm allows graphics processing unit (GPU) implementation and
speeds up the computational gains by an order-of-magnitude. The proposed
algorithm is applied to several linear and conic optimization problems. The
numerical experiments show that the proposed matrix-free algorithm
significantly achieves approximately an order-of-magnitude time improvement in
comparison with POGS, OSQP, or SCS solvers.

## References

  * [1] S. Mehrotra, “On the implementation of a primal-dual interior point method,” _SIAM Journal on Optimization_ , vol. 2, no. 4, pp. 575–601, 1992. [Online]. Available: https://doi.org/10.1137/0802028
  * [2] A. Domahidi, E. Chu, and S. Boyd, “ECOS: An SOCP solver for embedded systems,” in _2013 European Control Conference (ECC)_ , 2013, pp. 3071–3076.
  * [3] R. H. Tütüncü, K. C. Toh, and M. J. Todd, “Solving semidefinite-quadratic-linear programs using SDPT3,” _MATHEMATICAL PROGRAMMING_ , vol. 95, pp. 189–217, 2003.
  * [4] J. Nocedal and S. J. Wright, _Numerical Optimization_ , 2nd ed. New York, NY, USA: Springer, 2006.
  * [5] Y. Nesterov, _Lectures on Convex Optimization_ , 2nd ed. Springer Publishing Company, Incorporated, 2018.
  * [6] S. J. Wright, _Primal-Dual Interior-Point Methods_. USA: Society for Industrial and Applied Mathematics, 1997.
  * [7] S. C. Althoen and R. Mclaughlin, “Gauss-jordan reduction: A brief history,” _The American Mathematical Monthly_ , vol. 94, no. 2, pp. 130–142, 1987. [Online]. Available: https://doi.org/10.1080/00029890.1987.12000605
  * [8] P. S. Stanimirović and M. D. Petković, “Gauss–jordan elimination method for computing outer inverses,” _Applied Mathematics and Computation_ , vol. 219, no. 9, pp. 4667–4679, 2013. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0096300312010971
  * [9] R. Mittal and A. Al-Kurdi, “Lu-decomposition and numerical structure for solving large sparse nonsymmetric linear systems,” _Computers & Mathematics with Applications_, vol. 43, no. 1, pp. 131–155, 2002. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0898122101002796
  * [10] I. P. Stanimirović and M. B. Tasić, “Computation of generalized inverses by using the ldl decomposition,” _Applied Mathematics Letters_ , vol. 25, no. 3, pp. 526–531, 2012. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0893965911004630
  * [11] B. Fathi Vajargah, “A way to obtain monte carlo matrix inversion with minimal error,” _Applied Mathematics and Computation_ , vol. 191, no. 1, pp. 225–233, 2007. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0096300307002494
  * [12] S. Diamond and S. Boyd, “Matrix-free convex optimization modeling,” in _Optimization and its applications in control and data sciences_. Springer, 2016, pp. 221–264.
  * [13] J. Al-Jeiroudi, G. Gondzio, “Convergence analysis of the inexact infeasible interior-point method for linear optimization,” _Journal of Optimization Theory and Applications_ , vol. 141, no. 2, pp. 231–247, 2009. [Online]. Available: https://doi.org/10.1007/s10957-008-9500-5
  * [14] K.-C. Toh, “Solving large scale semidefinite programs via an iterative solver on the augmented systems,” _SIAM Journal on Optimization_ , vol. 14, no. 3, pp. 670–698, 2004. [Online]. Available: https://doi.org/10.1137/S1052623402419819
  * [15] M. R. Hestenes and E. Stiefel, “Methods of conjugate gradients for solving linear systems,” _Journal of research of the National Bureau of Standards_ , vol. 49, pp. 409–436, 1952.
  * [16] M. Fukuda, M. Kojima, and M. Shida, “Lagrangian dual interior-point methods for semidefinite programs,” _SIAM Journal on Optimization_ , vol. 12, no. 4, pp. 1007–1031, 2002. [Online]. Available: https://doi.org/10.1137/S1052623401387349
  * [17] M. KoăundefinedVara and M. Stingl, “On the solution of large-scale sdp problems by the modified barrier method using iterative solvers,” _Math. Program._ , vol. 109, no. 2–3, p. 413–444, Mar. 2007.
  * [18] C. Choi and Y. Ye, “Solving sparse semidefinite programs using the dual scaling algorithm with an iterative solver,” _Working paper, Department of Management Sciences, University of Iowa_ , 2000.
  * [19] X.-Y. Zhao, D. Sun, and K.-C. Toh, “A newton-cg augmented lagrangian method for semidefinite programming,” _SIAM Journal on Optimization_ , vol. 20, no. 4, pp. 1737–1765, 2010. [Online]. Available: https://doi.org/10.1137/080718206
  * [20] C. T. Kelley, _Iterative Methods for Linear and Nonlinear Equations_ , ser. Frontiers in Applied Mathematics. SIAM, 1995, no. 16. [Online]. Available: http://www.siam.org/books/textbooks/fr16_book.pdf
  * [21] M. DApuzzo, V. De Simone, and D. di Serafino, “On mutual impact of numerical linear algebra and large-scale optimization with focus on interior point methods,” _Computational Optimization and Applications_ , vol. 45, no. 2, pp. 283–310, 2010. [Online]. Available: https://doi.org/10.1007/s10589-008-9226-1
  * [22] F. E. Curtis, J. Nocedal, and A. Wächter, “A matrix-free algorithm for equality constrained optimization problems with rank-deficient jacobians,” _SIAM J. on Optimization_ , vol. 20, no. 3, p. 1224–1249, Sep. 2009.
  * [23] M. Saunders, B. Kim, C. Maes, A. Santiago, and M. Zahr, “PDCO: primal-dual interior method for convex Objectives,” _http://www.stanford.edu/group/SOL/ software/pdco.html_ , 2018.
  * [24] J. Gondzio, “Matrix-free interior point method,” _Comput. Optim. Appl._ , vol. 51, no. 2, p. 457–480, Mar. 2012. [Online]. Available: https://doi.org/10.1007/s10589-010-9361-3
  * [25] B. O’Donoghue, E. Chu, N. Parikh, and S. Boyd, “Conic optimization via operator splitting and homogeneous self-dual embedding,” _Journal of Optimization Theory and Applications_ , vol. 169, no. 3, pp. 1042–1068, Jun 2016\. [Online]. Available: https://doi.org/10.1007/s10957-016-0892-3
  * [26] Y. Zheng, G. Fantuzzi, A. Papachristodoulou, P. Goulart, and A. Wynn, “Chordal decomposition in operator-splitting methods for sparse semidefinite programs,” _Mathematical Programming_ , vol. 180, pp. 489–532, Mar. 2020.
  * [27] C. Fougner and S. Boyd, _Parameter selection and preconditioning for a graph form solver_. Cham: Springer International Publishing, 2018, pp. 41–61. [Online]. Available: https://doi.org/10.1007/978-3-319-67068-3_4
  * [28] B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd, “OSQP: An operator splitting solver for quadratic programs,” _ArXiv e-prints_ , Jan. 2018. [Online]. Available: https://arxiv.org/abs/1711.08013
  * [29] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” _Found. Trends Mach. Learn._ , vol. 3, no. 1, pp. 1–122, Jan. 2011. [Online]. Available: http://dx.doi.org/10.1561/2200000016
  * [30] P. Sopasakis, K. Menounou, and P. Patrinos, “SuperSCS: fast and accurate large-scale conic optimization,” in _18th European Control Conference (ECC)_ , 2019, pp. 1500–1505.
  * [31] M. Adil, S. Tavakkol, and R. Madani, “Rapid convergence of first-order numerical algorithms via adaptive conditioning,” 2021.
  * [32] R. Nishihara, L. Lessard, B. Recht, A. Packard, and M. Jordan, “A general analysis of the convergence of ADMM,” in _Proceedings of the 32nd International Conference on Machine Learning_ , vol. 37. ICML’15, 2015, pp. 343–352.
  * [33] P. Giselsson and S. Boyd, “Linear convergence and metric selection for Douglas-Rachford splitting and ADMM,” _IRE Transactions on Automatic Control_ , vol. 62, no. 2, pp. 532–544, 2 2017.
  * [34] J. Eckstein and W. Yao, “Understanding the convergence of the alternating direction method of multipliers: Theoretical and computational perspectives,” 2015.
  * [35] G. Banjac and P. Goulart, “Tight global linear convergence rate bounds for operator splitting methods,” _IEEE Transactions on Automatic Control_ , vol. 63, pp. 4126–4139, 2018.
  * [36] G. Banjac and P. J. Goulart, “Global linear convergence in operator splitting methods,” in _IEEE 55th Conference on Decision and Control (CDC)_ , 2016, pp. 233–238.
  * [37] A. Fu, J. Zhang, and S. P. Boyd, “Anderson accelerated Douglas-Rachford splitting,” _arXiv: Optimization and Control_ , 2019.
  * [38] M. Hong and Z. Luo, “On the linear convergence of the alternating direction method of multipliers,” _Mathematical Programming_ , vol. 162, no. 1-2, pp. 165–199, 2017.
  * [39] W. Ouyang, Y. Peng, Y. Yao, J. Zhang, and B. Deng, “Anderson acceleration for nonconvex ADMM based on Douglas-Rachford splitting,” _Computer Graphics Forum_ , vol. 39, no. 5, pp. 221–239, 2020. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1111/cgf.14081
  * [40] T. Pock and A. Chambolle, “Diagonal preconditioning for first order primal-dual algorithms in convex optimization,” in _International Conference on Computer Vision_ , Nov 2011, pp. 1762–1769.
  * [41] E. Chu, B. O’Donoghue, N. Parikh, and S. P. Boyd, “A primal-dual operator splitting method for conic optimization,” 2013.
  * [42] R. Madani, A. Kalbat, and J. Lavaei, “ADMM for sparse semidefinite programming with applications to optimal power flow problem,” in _2015 54th IEEE Conference on Decision and Control (CDC)_ , 2015, pp. 5932–5939.
  * [43] ——, “A low-complexity parallelizable numerical algorithm for sparse semidefinite programming,” _IEEE Transactions on Control of Network Systems_ , vol. 5, no. 4, pp. 1898–1909, 2018.
  * [44] J. Douglas and H. H. Rachford, “On the numerical solution of heat conduction problems in two and three space variables,” _Transactions of the American Mathematical Society_ , vol. 82, no. 2, pp. 421–439, 1956. [Online]. Available: http://www.jstor.org/stable/1993056
  * [45] R. J. Vanderbei, “Symmetric quasidefinite matrices,” _SIAM Journal on Optimization_ , vol. 5, no. 1, pp. 100–113, 1995. [Online]. Available: https://doi.org/10.1137/0805005
  * [46] G. H. Golub and C. F. Van Loan, _Matrix Computations (3rd Ed.)_. USA: Johns Hopkins University Press, 1996.
  * [47] J. Nocedal and S. Wright, _Numerical Optimization: Springer Series in Operations Research and Financial Engineering_. Springer, 2006.
  * [48] C. C. Paige and M. A. Saunders, “Lsqr: An algorithm for sparse linear equations and sparse least squares,” _ACM Trans. Math. Software_ , pp. 43–71, 1982.
  * [49] H. Huang, J. M. Dennis, L. Wang, and P. Chen, “A scalable parallel lsqr algorithm for solving large-scale linear system for tomographic problems: A case study in seismic tomography,” _Procedia Computer Science_ , vol. 18, pp. 581–590, 2013, 2013 International Conference on Computational Science. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1877050913003657