# General-purpose preconditioning for regularized interior point methods

ERGO Technical Report 21-004

Jacek Gondzio  Spyridon Pougkakiotis  John W. Pearson

###### Abstract

In this paper we present general-purpose preconditioners for regularized
augmented systems, and their corresponding normal equations, arising from
optimization problems. We discuss positive definite preconditioners, suitable
for CG and MINRES. We consider “sparsifications” which avoid situations in
which eigenvalues of the preconditioned matrix may become complex. Special
attention is given to systems arising from the application of regularized
interior point methods to linear or nonlinear convex programming problems.

##  1 Introduction

In this paper, we are concerned with applying Krylov-subspace methods for the
efficient solution of systems of the following form:

$\underbrace{\begin{bmatrix}-(Q+\rho I_{n})&A^{\top}\\\ A&\delta
I_{m}\end{bmatrix}}_{K}\begin{bmatrix}\Delta x\\\ \Delta
y\end{bmatrix}=\begin{bmatrix}\xi_{1}\\\ \xi_{2}\end{bmatrix},$ (1.1)

where $A\in\mathbb{R}^{m\times n}$ (with $m\leq n$), $Q\succeq
0\in\mathbb{R}^{n\times n}$, $I_{n}$ is the identity matrix of size $n$, and
$\delta,\ \rho>0$. Such systems arise in a plethora of applications [6], which
go far beyond optimization. However, in this paper we restrict the discussion
to the case of regularized systems arising in Interior Point Methods (IPMs)
for optimization [1, 3, 20, 34, 40, 43, 46]. Due to the potential large
dimensions of the systems, they are often solved by means of iterative
techniques, usually from the family of Krylov-subspace methods [22]. To
guarantee efficiency of such methods the possibly ill-conditioned system (1.1)
often needs to be appropriately preconditioned and, indeed, there exists a
rich literature which addresses the issue (see the discussions in [4, 6, 8, 9,
16, 17, 15], and the references therein).

Multiple preconditioning approaches have been developed in the literature,
used to accelerate the associated iterative methods. These can be divided into
positive definite (e.g. [6, 10, 21, 30, 32, 33, 45]) and indefinite ones (e.g.
[16, 25, 26, 30, 42]). The latter are often employed within long-recurrence
non-symmetric solvers (such as the Generalized Minimal RESidual method [42]),
while the former can be used within short-recurrence methods (such as the
MINimal RESidual method [35]). A comprehensive study of saddle point systems
and their associated “optimal” preconditioners can be found in [6]. Indefinite
preconditioners are significantly more difficult to analyze and a simple
spectral analysis is not sufficient to deduce their effectiveness (see [23]).
On the other hand, positive definite preconditioners are often easier to
analyze, and the eigenvalues of the preconditioned matrices allow one to
theoretically compare different preconditioning approaches.

In this paper, we are focused on systems arising from the application of
regularized interior point methods for the solution of an arbitrary convex
programming problem. In this case, $Q$ represents the Hessian of the primal
barrier problem’s objective function (or the Hessian of the Lagrangian in the
nonlinear programming case), $A$ represents the constraint matrix (or the
Jacobian of the constraints in the nonlinear programming case), while $\rho$
and $\delta$ are the primal and dual regularization parameters, respectively.
We note that the IPM may contribute a term to the $(1,1)$ or the $(2,2)$ block
of (1.1), depending on the form of the constraints and non-negativity
variables. Here we assume that the term is added in the $(1,1)$ block. For
example, the matrix $Q$ may be written as $Q\coloneqq H+\Theta^{-1}$, where
$H$ is the Hessian of the Lagrangian and $\Theta\coloneqq XZ^{-1}$ is a
diagonal IPM scaling matrix (assuming $x,\ z$ are the primal and dual non-
negative variables, respectively, while $X,\ Z$ denote the diagonal matrices
with diagonal entries taken from vectors $x,\ z$, respectively) which
originates from the use of the logarithmic barrier.

We present positive definite preconditioning approaches that can be used
within MINRES [35] or the Conjugate Gradient method [24], and we provide some
basic spectral analysis results for the associated preconditioned systems.
More specifically, we consider preconditioners which are derived by
“sparsifications” of system (1.1), that is, by dropping specific entries from
sparse matrices $Q$ and $A$, thus making them more sparse and hence easier to
factorize. Various such approaches have been proposed to date and include:
preconditioners which exploit an early guess of a basic–nonbasic partition of
variables to drop columns from $A$ [33], constraint preconditioners [9, 16,
17, 15], inexact constraint preconditioners [8] which drop specific entries in
the matrices $Q$ and $A$, and of course a plethora of preconditioners which
involve various levels of incomplete Cholesky factorizations of the matrix in
(1.1), see for example [10]. The literature on preconditioners is growing
rapidly and we refer the interested reader to [5, 6, 13, 37, 48] and the
references therein, for a detailed discussion.

We consider dropping off-diagonal entries of $Q$, but restrict the elimination
of entries in $A$ only to the removal of complete columns. Additionally, we
consider sparsifying parts of rows of the Schur complement corresponding to
system (1.1). Such a strategy guarantees that we avoid situations in which
eigenvalues of the preconditioned matrix may become complex (such as those
employed in [8]), which as a consequence would have required us to employ non-
symmetric Krylov methods.

In order to construct the preconditioners, by following [7], we take advantage
of the properties of the logarithmic barrier, that allow us to know in advance
which columns of the problem matrix are important and which are less
influential. In particular, assuming the aforementioned representation of $Q$
as $Q=H+\Theta^{-1}$, the logarithmic barrier indicates which variables of the
problem are likely to be inactive in the solution. More precisely, the
variables are naturally split into “basic”–$\mathcal{B}$ (not in the simplex
sense), “non-basic”–$\mathcal{N}$, and “undecided”–$\mathcal{U}$. Hence, as
IPMs progress towards optimality, we expect the following partition of the
diagonal barrier matrix $\Theta^{-1}$:

$\displaystyle\forall
j\in\mathcal{N}:\left(\Theta^{\left(\mathcal{N},\mathcal{N}\right)}\right)^{-1}=\bm{\Theta}\left(\mu^{-1}\right),\quad\forall
j\in\mathcal{B}:\left(\Theta^{\left(\mathcal{B},\mathcal{B}\right)}\right)^{-1}=\bm{\Theta}\left(\mu\right),\quad\forall
j\in\mathcal{U}:\left(\Theta^{\left(\mathcal{U},\mathcal{U}\right)}\right)^{-1}=\bm{\Theta}\left(1\right),$

where $\mu$ is the barrier parameter (and is such that $\mu\rightarrow 0$),
$\mathcal{N}$, $\mathcal{B}$ and $\mathcal{U}$ are mutually disjoint, and
$\mathcal{N}\cup\mathcal{B}\cup\mathcal{U}=\\{1,\ldots,n\\}$, while
$\bm{\Theta}(\cdot)$ denotes that two positive quantities are of the same
order of magnitude (see the notation section at the end of the introduction).
Given the large magnitude of the diagonal elements of $Q$ for any
$j\in\mathcal{N}$ (assuming $\mu$ is close to zero), we expect that the
corresponding columns of $A$ (i.e. $A^{(:,\mathcal{N})}$) will not contribute
important information, and thus can be set to zero when constructing a
preconditioner for (1.1). In [7], the Hessian was approximated by its
diagonal. In this paper, we extend this work by allowing the utilization of
non-diagonal Hessian information. More specifically, we showcase how to
analyze, construct, and apply the inverse of preconditioners in which we only
drop non-diagonal elements of $Q$ corresponding to diagonal elements in
$\mathcal{N}$. We should note at this point that such a splitting of $Q$
occurs in other second-order optimization methods as well, such as those based
on augmented Lagrangian strategies (e.g. see [27]).

Furthermore, we discuss some approaches for dealing with problems for which
the matrix $A$ may contain a subset of dense columns or rows. Any dense column
or row in $A$ can pose great difficulties when trying to factorize the
associated saddle-point matrix (or a preconditioner approximating it). Thus,
it is desirable to alleviate the dangers of such columns or rows, by
appropriate “sparsifications” of the preconditioner, allowing us to reduce the
memory requirements of applying its inverse.

All such “sparsifications” are captured in a general result presented in
Section 2 which provides the spectral analysis of an appropriate
preconditioned normal equations’ matrix. The main theorem sheds light on
consequences of sparsifying rows of the normal equations corresponding to
(1.1), or dropping columns of $A$, and demonstrates that the former might
produce a larger number of non-unit eigenvalues. In Section 3, these normal
equation approximations are utilized to construct positive definite block-
diagonal preconditioners for the saddle point system in (1.1), and the
spectral properties of the resulting preconditioned matrices are also
discussed. Additionally, an alternative saddle-point preconditioner based on a
$LDL^{\top}$ decomposition is presented.

All of the preconditioning approaches discussed are compared numerically on
saddle-point systems arising from the application of a regularized IPM for the
solution of real-life linear and convex quadratic programming problems. In
particular, we present some numerical results on certain test problems taken
from the Netlib (see [31]) and the Maros–Mészáros (see [29]) collections.
Then, we test the preconditioners on examples of Partial Differential Equation
(PDE) optimal control problems. All preconditioning approaches have been
implemented within an Interior Point-Proximal Method of Multipliers (IP-PMM)
framework, which is a polynomially convergent primal-dual regularized IPM,
based on the developments in [40, 41]. A robust implementation is provided.

It is worth stressing that the proposed preconditioners are general and do not
assume the knowledge of special structures which might be present in the
matrices $Q$ and $A$ (such as block-diagonal, block-angular, network, PDE-
induced, etc.). Therefore they may be applied within general-purpose IPM
solvers for linear and convex quadratic programming problems.

##### Notation.

Throughout this paper we use lowercase Roman and Greek letters to indicate
vectors and scalars. Capitalized Roman fonts are used to indicate matrices.
Superscripts are used to denote the components of a vector/matrix. Sets of
indices are denoted by caligraphic capital fonts. For example, given
$M\in\mathbb{R}^{m\times n}$, $v\in\mathbb{R}^{n}$,
$\mathcal{R}\subseteq\\{1,\ldots,m\\}$, and
$\mathcal{C}\subseteq\\{1,\ldots,n\\}$, we set
$v^{\mathcal{C}}:=(v^{i})_{i\in\mathcal{C}}$ and
$M^{(\mathcal{R},\mathcal{C})}:=(m^{(i,j)})_{i\in\mathcal{R},j\in\mathcal{C}}$,
where $v^{i}$ is the $i$-th entry of $v$ and $m^{(i,j)}$ the $(i,j)$-th entry
of $M$. Additionally, the full set of indices is denoted by a colon. In
particular, $M^{(:,\mathcal{C})}$ denotes all columns of $M$ with indices in
$\mathcal{C}$. Given a matrix $M$, we denote the diagonal matrix with the same
diagonal elements as $M$ by $\textnormal{Diag}(M)$. We use $\lambda_{\min}(B)$
($\lambda_{\max}(B)$, respectively) to denote the minimum (maximum) eigenvalue
of an arbitrary square matrix $B$ with real eigenvalues. Similarly,
$\sigma_{\min}(B)$ ($\sigma_{\max}(B)$, respectively) denotes the minimum
(maximum) singular value of an arbitrary rectangular matrix $B$. We use
$0_{m,n}$ to denote a matrix of size $m\times n$ with entries equal to $0$.
Furthermore, we use $I_{n}$ to indicate the identity matrix of size $n$. For
any finite set $\mathcal{A}$, we denote by $|\mathcal{A}|$ its cardinality.
Finally, given two positive functions $T,\
f\colon\mathbb{R}_{+}\mapsto\mathbb{R}_{+}$, we write $T(x)=\bm{\Theta}(f(x))$
if these functions are of the same order of magnitude, that is, there exist
constants $c_{1},\ c_{2}>0$ and some $x_{0}\geq 0$ such that $c_{1}f(x)\leq
T(x)\leq c_{2}f(x)$, for all $x\geq x_{0}$.

##### Structure of the article.

The rest of this paper is organized as follows. In Section 2 we present some
preconditioners suitable for the normal equations. Then, in Section 3, we
adapt these preconditioners to regularized saddle point systems. Subsequently,
in Section 4 we focus on saddle point systems arising from the application of
regularized IPMs to convex programming problems, and present some numerical
results. Finally, in Section 5, we deliver our conclusions.

##  2 Regularized normal equations

We begin by defining the regularized normal equations matrix (or Schur
complement) $M\coloneqq AGA^{\top}+\delta I_{m}\in\mathbb{R}^{m\times m}$,
corresponding to (1.1), where $G\coloneqq(Q+\rho I_{n})^{-1}\succ
0\in\mathbb{R}^{n\times n}$. In this section, we derive and analyze
preconditioning approaches for $M$. As we have already mentioned in the
introduction, we achieve a simplification of the preconditioner by setting to
zero certain columns of the matrix $A$ (and consequently sparsifying the
corresponding rows and columns of $Q$), as well as by sparsifying certain rows
of the matrix $M$.

More specifically, we define two integers, $k_{c}$ and $k_{r}$, such that
$0\leq k_{c}\leq n$ and $0\leq k_{r}\leq m$. The former counts the number of
columns of $A$ (and corresponding columns and rows of $Q$) that are set to
zero (that are sparsified, respectively), while the latter counts the number
of rows of the matrix $M$ that are sparsified. At this point, we assume that
we have been given these columns or rows, but we later specify how these can
be chosen (see Remark 1 and Section 4). In order to highlight these given
columns and rows, we assume that we are given two permutation matrices; a
column permutation $\mathscr{P}_{c}\in\mathbb{R}^{n\times n}$, and a row
permutation $\mathscr{P}_{r}\in\mathbb{R}^{m\times m}$. Applied to the
constraint matrix $A$ in (1.1), these permutations bring all the columns and
rows which will need to be treated specially to the leading positions of
columns and rows of $\mathscr{P}_{r}A\mathscr{P}_{c}$, respectively.

Given the previous permutation matrices, let us firstly define an
approximation to the matrix $Q$. In particular, we approximate $Q$ by the
following block-diagonal and positive semi-definite matrix:

$\widehat{Q}\coloneqq\mathscr{P}_{c}\begin{bmatrix}\widehat{Q}_{1}&0_{k_{c},(n-k_{c})}\\\
0_{n-k_{c},k_{c}}&\widehat{Q}_{2}\end{bmatrix}\mathscr{P}_{c}^{\top},\quad\text{assuming}\quad
Q\equiv\mathscr{P}_{c}\begin{bmatrix}Q_{1}&Q_{3}^{\top}\\\
Q_{3}&Q_{2}\end{bmatrix}\mathscr{P}_{c}^{\top},$ (2.1)

where $\widehat{Q}_{1}=Q_{1}$ or $\widehat{Q}_{1}=\textnormal{Diag}(Q_{1})$,
and $\widehat{Q}_{2}=Q_{2}$ or $\widehat{Q}_{2}=\textnormal{Diag}(Q_{2})$
(both cases are treated concurrently). The column permutation
$\mathscr{P}_{c}$ reorders symmetrically both rows and columns of the matrix
$Q$ in (1.1), and places the $k_{c}$ columns and rows which will be sparsified
at the leading $(1,1)$ block of the permuted version of the matrix $Q$. Given
this approximation of $Q$, let us define an approximated normal equations’
matrix that will be of interest when analyzing the spectral properties of the
preconditioned matrices derived in this paper:

$\widehat{M}\coloneqq A\widehat{G}A^{\top}+\delta
I_{m},\quad\widehat{G}\coloneqq\left(\widehat{Q}+\rho
I_{n}\right)^{-1}\equiv\mathscr{P}_{c}\begin{bmatrix}\left(\widehat{Q}_{1}+\rho
I_{k_{c}}\right)^{-1}&0_{k_{c},(n-k_{c})}\\\
0_{(n-k_{c}),k_{c}}&\left(\widehat{Q}_{2}+\rho
I_{n-k_{c}}\right)^{-1}\end{bmatrix}\mathscr{P}_{c}^{\top}.$ (2.2)

In what follows, we derive a preconditioner for the approximated normal
equations’ matrix $\widehat{M}$. We should note that system (1.1) is solved
using the normal equations only if $Q$ is diagonal (due to obvious numerical
considerations), in which case $Q\equiv\widehat{Q}$, and thus
$M\equiv\widehat{M}$. If this is not the case, we would like to derive a
preconditioner for the approximate normal equations’ matrix $\widehat{M}$.
Later on, and in particular in Section 3, this is utilized to derive and
analyze a preconditioner for the matrix $K$, defined in (1.1).

We proceed by introducing some notation, for convenience of exposition. Given
the definitions in (2.1), (2.2), we can write:

$B\coloneqq\mathscr{P}_{r}A\widehat{G}^{\frac{1}{2}}\mathscr{P}_{c}=\begin{bmatrix}B_{11}&B_{12}\\\
B_{21}&B_{22}\end{bmatrix},$

where $\mathscr{P}_{r}$ is a given row-permutation matrix,
$B_{11}\in\mathbb{R}^{k_{r}\times k_{c}}$,
$B_{12}\in\mathbb{R}^{k_{r}\times(n-k_{c})}$,
$B_{21}\in\mathbb{R}^{(m-k_{r})\times k_{c}}$, and
$B_{22}\in\mathbb{R}^{(m-k_{r})\times(n-k_{c})}$. Notice that the aim of the
row-permutation matrix $\mathscr{P}_{r}$, is to bring on top all rows of the
matrix $\widehat{M}$, that we are planning to sparsify in order to construct
the preconditioner. Let us further introduce the following notation:

$\mathscr{P}_{r}\widehat{M}\mathscr{P}_{r}^{\top}\equiv\begin{bmatrix}\widehat{M}_{11}&\widehat{M}_{21}^{\top}\\\
\widehat{M}_{21}&\widehat{M}_{22}\end{bmatrix},$

where $\widehat{M}_{11}$, $\widehat{M}_{21}$, and $\widehat{M}_{22}$ are
defined as:

$\displaystyle\widehat{M}_{11}\coloneqq
B_{11}B_{11}^{\top}+B_{12}B_{12}^{\top}+\delta I_{k_{r}}\
\in\mathbb{R}^{k_{r}\times k_{r}},$ $\displaystyle\widehat{M}_{21}\coloneqq
B_{21}B_{11}^{\top}+B_{22}B_{12}^{\top}\ \in\mathbb{R}^{(m-k_{r})\times
k_{r}},$ $\displaystyle\widehat{M}_{22}\coloneqq
B_{21}B_{21}^{\top}+B_{22}B_{22}^{\top}+\delta I_{m-k_{r}}\
\in\mathbb{R}^{(m-k_{r})\times(m-k_{r})}.$

In what follows, we present two preconditioning strategies for the matrix
$\widehat{M}$. Both approaches exploit the sparsification of the matrix
$\widehat{M}$. The first approach relies on a Cholesky decomposition of a
sparsified matrix, while the second approach is based on an $LDL^{\top}$
decomposition of a sparsified augmented system matrix, which is used to
implicitly derive a preconditioner for $\widehat{M}$.

###  2.1 A Cholesky-based preconditioner

Our first proposal is to consider preconditioning
$\mathscr{P}_{r}\widehat{M}\mathscr{P}_{r}^{\top}$ with the following matrix:

$P_{NE,(k_{c},k_{r})}\coloneqq\begin{bmatrix}\widehat{M}_{11}&0_{k_{r},(m-k_{r})}\\\
0_{(m-k_{r}),k_{r}}&\widetilde{M}_{22}\end{bmatrix},\qquad\widetilde{M}_{22}\coloneqq\widehat{M}_{22}-B_{21}B_{21}^{\top}.$
(2.3)

The notation $P_{NE,(k_{c},k_{r})}$ signifies that this is a preconditioner
for the normal equations, in which we drop $k_{c}$ columns from the matrix $A$
and sparsify $k_{r}$ rows of the matrix $\widehat{M}$. Notice that if
$k_{c}=0$ (that is, we only sparsify certain rows of $\widehat{M}$ to
construct the preconditioner), we can write
$\small{B\equiv\mathscr{P}_{r}A\widehat{G}^{\frac{1}{2}}=\begin{bmatrix}B_{12}\\\
B_{22}\end{bmatrix}}$, while $B_{11}$, $B_{21}$ are zero-dimensional, and
hence absent. In this case, we have

$P_{NE,(0,k_{r})}\coloneqq\begin{bmatrix}\widehat{M}_{11}&0_{k_{r},(m-k_{r})}\\\
0_{(m-k_{r}),k_{r}}&\widehat{M}_{22}\end{bmatrix}.$

On the other hand, if $k_{r}=0$ (that is, we only drop $k_{c}$ columns from
$A$ to construct the preconditioner), we have
$B\equiv\begin{bmatrix}B_{21}&B_{22}\end{bmatrix}$, and $B_{11},\ B_{12}$ are
absent. Then, we obtain

$P_{NE,(k_{c},0)}=\widetilde{M}_{22}.$

Notice that the latter is obtained since $\widehat{Q}$ is block-separable
(with respect to the permutation $\mathscr{P}_{c}$), and thus dropping the
respective $k_{c}$ columns of $A$, results in dropping $B_{21}B_{21}^{\top}$.
For simplicity of notation, for the rest of this subsection we set
$P_{NE,(k_{c},k_{r})}\equiv P_{NE}$.

In the following theorem, we analyze the spectrum of the preconditioned matrix
$P_{NE}^{-1}\mathscr{P}_{r}\widehat{M}\mathscr{P}_{r}^{\top}$, with respect to
the spectrum of the associated matrices.

######  Theorem 2.1.

The preconditioned matrix
$P_{NE}^{-1}\mathscr{P}_{r}\widehat{M}\mathscr{P}_{r}^{\top}$ has at least
$\max\\{m-(2k_{r}+k_{c}),0\\}$ eigenvalues at $1$. If $k_{c}>0$ and $k_{r}>0$,
all remaining eigenvalues lie in the following interval:

$I_{k_{c},k_{r}}\coloneqq\left[\frac{\delta}{\delta+\sigma_{\max}^{2}(B)},2+\frac{\lambda_{\max}(B_{21}B_{21}^{\top})}{\delta+\lambda_{\min}(B_{22}B_{22}^{\top})}\right].$

If $k_{c}>0$ and $k_{r}=0$, the previous interval reduces to

$I_{k_{c}}\coloneqq\left[1,1+\frac{\lambda_{\max}(B_{21}B_{21}^{\top})}{\delta+\lambda_{\min}(B_{22}B_{22}^{\top})}\right],$

while if $k_{r}>0$ and $k_{c}=0$, we obtain

$I_{k_{r}}\coloneqq\left[\frac{\delta}{\delta+\sigma_{\max}^{2}(B)},2\right].$

###### Proof.

Given an arbitrary eigenvalue $\lambda$ (which must be positive since
$P_{NE}\succ 0$ and $\widehat{M}\succ 0$) corresponding to a unit eigenvector
$v$, let us write the generalized eigenproblem as:

$\begin{bmatrix}\widehat{M}_{11}&\widehat{M}_{21}^{\top}\\\
\widehat{M}_{21}&\widehat{M}_{22}\end{bmatrix}\begin{bmatrix}v_{1}\\\
v_{2}\end{bmatrix}=\lambda\begin{bmatrix}\widehat{M}_{11}v_{1}\\\
\widetilde{M}_{22}v_{2}\end{bmatrix}.$ (2.4)

We separate the analysis into two cases.

Case 1: Let $v_{2}\in\textnormal{Null}(\widehat{M}_{21}^{\top})$. Firstly, we
notice that:

$\textnormal{dim}\left(\textnormal{Null}\left(\widehat{M}_{21}^{\top}\right)\right)=(m-k_{r})-\textnormal{rank}\left(\widehat{M}_{21}^{\top}\right)\geq\max\\{m-2k_{r},0\\}.$

Two sub-cases arise here. For the first sub-case, we notice that if $v_{1}\neq
0$, then from positive definiteness of $\widehat{M}_{11}$, combined with the
first block equation of (2.4), we obtain that $\lambda=1$. In turn, we claim
that this implies that $v_{2}\in\textnormal{Null}(B_{21}B_{21}^{\top})$ and
$v_{1}\in\textnormal{Null}(\widehat{M}_{21})$. To see this, assume that
$v_{2}\notin\textnormal{Null}(B_{21}B_{21}^{\top})$. Then from the second
block equation of (2.4) we obtain:

$\widehat{M}_{21}v_{1}+\widehat{M}_{22}v_{2}=\widetilde{M}_{22}v_{2}~{}~{}\Rightarrow~{}~{}\widehat{M}_{21}v_{1}=-B_{21}B_{21}^{\top}v_{2},$

where we used the definition of $\widetilde{M}_{22}$. If
$v_{2}\notin\textnormal{Null}(B_{21}B_{21}^{\top})$, this implies that
$v_{2}^{\top}B_{21}B_{21}^{\top}v_{2}>0$. The previous equation then yields
that

$v_{2}^{\top}\widehat{M}_{21}v_{1}=-v_{2}^{\top}B_{21}B_{21}^{\top}v_{2}~{}~{}\Rightarrow~{}~{}0=-v_{2}^{\top}B_{21}B_{21}^{\top}v_{2}<0,$

which follows from the base assumption (i.e.
$v_{2}\in\textnormal{Null}(\widehat{M}_{21}^{\top})$), and results in a
contradiction. Hence, $v_{2}\in\textnormal{Null}(B_{21}B_{21}^{\top})$. On the
other hand, if $v_{1}\notin\textnormal{Null}(\widehat{M}_{21})$ then the
second block equation yields directly a contradiction, since we have shown
that $v_{2}\in\textnormal{Null}(B_{21}B_{21}^{\top})$.

Next we consider the second sub-case, i.e. $v_{1}=0$. Combined with our base
assumption, the first block equation of (2.4) becomes redundant. From the
second block equation of the eigenproblem, and using $v_{1}=0$, we obtain:

$\begin{split}v_{2}^{\top}\widehat{M}_{22}v_{2}=&\ \lambda
v_{2}^{\top}\widetilde{M}_{22}v_{2}\\\
\Rightarrow~{}~{}v_{2}^{\top}\left(\widetilde{M}_{22}+B_{21}B_{21}^{\top}\right)v_{2}=&\
\lambda v_{2}^{\top}\widetilde{M}_{22}v_{2}.\end{split}$

Hence we have that:

$\lambda=1+\frac{v_{2}^{\top}(B_{21}B_{21}^{\top})v_{2}}{v_{2}^{\top}\widetilde{M}_{22}v_{2}}\leq
1+\frac{\lambda_{\max}(B_{21}B_{21}^{\top})}{\delta+\lambda_{\min}(B_{22}B_{22}^{\top})}.$

All eigenvalues in this case can be bounded by the previous inequality and
there will be at most $\textnormal{rank}(B_{21}B_{21}^{\top})$ non-unit
eigenvalues. On the other hand, if
$v_{2}\in\textnormal{Null}(B_{21}B_{21}^{\top})$, then trivially $\lambda=1$.
This concludes the first case. Notice that this case would occur necessarily
if $k_{r}=0$, and thus, we obtain the interval $I_{k_{c}}$.

Case 2: In this case, we assume that
$v_{2}\notin\textnormal{Null}(\widehat{M}_{21}^{\top})$. In what follows we
assume $\lambda\neq 1$ (noting that $\lambda=1$ would only occur if
$v_{1}\in\textnormal{Null}(\widehat{M}_{21})$ and
$v_{2}\in\textnormal{Null}(B_{21}B_{21}^{\top})$), and there are at most
$2k_{r}$ such eigenvalues. Given the previous assumption, and using the first
block equation in (2.4), we obtain:

$v_{1}=\frac{1}{\lambda-1}\widehat{M}_{11}^{-1}\widehat{M}_{21}^{\top}v_{2}.$

Substituting the previous into the second block equation of (2.4), yields the
following generalized eigenproblem:

$\begin{split}\left(\widehat{M}_{21}\widehat{M}_{11}^{-1}\widehat{M}_{21}^{\top}+(\lambda-1)B_{21}B_{21}^{\top}\right)v_{2}=(\lambda-1)^{2}\widetilde{M}_{22}v_{2},\end{split}$
(2.5)

where we used the definitions of $\widetilde{M}_{22}$ and $\widehat{M}_{22}$.
Premultiplying (2.5) by $v_{2}^{\top}$ and rearranging yields the following
quadratic algebraic equation that $\lambda$ must satisfy in this case:

$\lambda^{2}+\beta\lambda+\gamma=0,$ (2.6)

where

$\beta\coloneqq-2-\frac{v_{2}^{\top}B_{21}B_{21}^{\top}v_{2}}{v_{2}^{\top}\widetilde{M}_{22}v_{2}}$

and

$\gamma\coloneqq
1-\frac{v_{2}^{\top}\left(\widehat{M}_{21}\widehat{M}_{11}^{-1}\widehat{M}_{21}^{\top}-B_{21}B_{21}^{\top}\right)v_{2}}{v_{2}^{\top}\widetilde{M}_{22}v_{2}}.$

Let us notice that the smallest eigenvalue is at least as large as
$\frac{\delta}{\delta+\sigma_{\max}^{2}(B)}$. This follows from positive
definiteness of $P_{NE}$ and $M$, and the bound can be deduced by noticing
that
$\lambda_{\min}(P_{NE}^{-1}\widehat{M})\geq\frac{\lambda_{\min}(\widehat{M})}{\lambda_{\max}(P_{NE})}\geq\frac{\delta}{\delta+\sigma_{\max}^{2}(B)}$.
Still we need to find an upper bound for the largest eigenvalue. To that end,
notice that:

$\gamma=\frac{v_{2}^{\top}\left(\widehat{M}_{22}-\widehat{M}_{21}\widehat{M}_{11}^{-1}\widehat{M}_{21}^{\top}\right)v_{2}}{v_{2}^{\top}\widetilde{M}_{22}v_{2}},$

which follows from the definition of $\widetilde{M}_{22}$. Positive
definiteness of $\widehat{M}$ then implies that $\gamma>0$. From the last
relation we also have that:

$0<\gamma\leq
1+\frac{v_{2}^{\top}B_{21}B_{21}^{\top}v_{2}}{v_{2}^{\top}\widetilde{M}_{22}v_{2}}=\frac{v_{2}^{\top}\widehat{M}_{22}v_{2}}{v_{2}^{\top}\widetilde{M}_{22}v_{2}}\leq
1+\frac{\lambda_{\max}(B_{21}B_{21}^{\top})}{\delta+\lambda_{\min}(B_{22}B_{22}^{\top})}\eqqcolon\gamma_{u}.$

Furthermore,
$\beta_{l}\coloneqq-\left(2+\frac{\lambda_{\max}\left(B_{21}B_{21}^{\top}\right)}{\delta+\lambda_{\min}\left(B_{22}B_{22}^{\top}\right)}\right)\leq\beta\leq-2.$
From the previous inequality, one can also observe that $\gamma<-\beta-1$.

Returning to (2.6), we first consider the following solution:

$\lambda_{-}=\frac{1}{2}\left(-\beta-\sqrt{\beta^{2}-4\gamma}\right).$

It is easy to see that $\beta^{2}-4\gamma$ is always larger than $0$. Next, we
notice that the relation for $\lambda_{-}$ is increasing with respect to
$\gamma$. We omit finding a lower bound for $\lambda_{-}$ since this was
established earlier. For the upper bound, we use the fact that
$\gamma<-\beta-1$, to obtain:

$\lambda_{-}<\frac{1}{2}\left(-\beta-\sqrt{\beta^{2}+4(\beta+1)}\right)=\frac{1}{2}\left(|\beta|-|\beta+2|\right)=1,$

since $\beta\leq-2$ (also, in the beginning of this case, we have treated
$\lambda_{-}=1$ separately).

Finally, we consider the other solution of (2.6), which reads:

$\lambda_{+}=\frac{1}{2}\left(-\beta+\sqrt{\beta^{2}-4\gamma}\right).$

Firstly, we can easily notice that $\lambda_{+}>1$. Subsequently, upon
noticing that $\lambda_{+}$ is decreasing with respect to $\gamma$, we can
obtain the following obvious bound:

$\lambda_{+}\leq|\beta|\leq-\beta_{l}.$

Now, let us observe that dropping $\widehat{M}_{21}$ and
$\widehat{M}_{21}^{\top}$ yields at most
$k_{r}+\textnormal{rank}(\widehat{M}_{21}^{\top})\leq 2k_{r}$ eigenvalue
outliers. Similarly, dropping $B_{21}B_{21}^{\top}$ from the $(2,2)$ block of
$M$ yields at most $\textnormal{rank}\left(B_{21}\right)\leq k_{c}$ eigenvalue
outliers. Hence, there will be at least
$\max\left\\{m-(2k_{r}+k_{c}),0\right\\}$ eigenvalues of the preconditioned
matrix at $1$.

Finally, the case where $k_{c}=0$ and $k_{r}>0$ follows by a direct
generalization of [14, Theorem 4.1], and completes the proof. ∎∎

######  Remark 1.

Now that we have presented the spectral properties of the preconditioned
system, let us discuss the use of such a preconditioning strategy. In
practice, we solve system (1.1) using the normal equations, only when matrix
$Q$ is diagonal. As already mentioned, in this case $\widehat{M}=M$, and thus
$P_{NE}$ is a preconditioner for the normal equations’ matrix. In Section 3,
we discuss how $P_{NE}$ is utilized to construct preconditioners for the
saddle point matrix in (1.1), even if $Q$ is not diagonal (in which case
$\widehat{M}$ is an approximation of the normal equations’ matrix).

Let us now discuss how to choose which columns of $A$ (or rows of
$\widehat{M}$, respectively) to drop (sparsify, respectively).

  * $\bullet$

Firstly, it often happens in optimization, and especially when solving systems
arising from the application of an interior point method (as already mentioned
in the introduction), to have certain diagonal elements of $G$ that are very
small. In view of this property, and given the bound presented in Theorem 2.1,
we can observe that dropping all columns corresponding to small diagonal
elements in $G$ results in manageable and not too sizeable outliers. Such a
preconditioner was proposed in [7] for the case where
$\widehat{Q}=\textnormal{Diag}(Q)$, and arises as a special case of $P_{NE}$
in (2.3), by choosing $k_{r}=0$ and a suitable permutation matrix
$\mathscr{P}_{c}$, traversing first the $k_{c}$ indices corresponding to the
smallest diagonal elements of $G$.

  * $\bullet$

Secondly, it is common in many application areas to have a small number of
columns or rows of $A$ that are dense. Such columns (or rows) could pose
significant difficulties as they produce dense factors when one tries to
factorize the normal equations (e.g. using a Cholesky decomposition). This is
especially the case for dense columns. A single dense column of $A$ with $p$
non-zero entries induces a dense window of size $p\times p$ in the normal
equations (we refer the reader to the discussion in [2, Section 4]). The use
of a preconditioner like the one defined in (2.3) serves the purpose of
dropping (sparsifying, respectively) such columns of $A$ (rows of
$\widehat{M}$, respectively), thus making the Cholesky factors of $P_{NE}$
significantly more sparse. For example, we may find two permutation matrices
$\mathscr{P}_{r}$, $\mathscr{P}_{c}$ which sort the rows and columns,
respectively, of $A$ in descending order of their number of non-zeros, and
write $\widehat{A}=\mathscr{P}_{r}A\mathscr{P}_{c}$. Then, the resulting
normal equations read as
$\mathcal{P}_{r}^{\top}BB^{\top}\mathcal{P}_{r}+\delta I_{m}$. As long as the
number of dropped columns or rows is low (which is observed in several
applications), the number of outliers produced by this dropping strategy is
manageable. While some of these outliers will be dangerously close to zero
(given that the regularization parameter $\delta>0$ is small), they can be
dealt with efficiently. We should note, however, that if $\widehat{Q}$ is non-
diagonal, without a sparse representation of its inverse, this strategy would
be unattractive to employ, and thus in this work, we consider it only when
$\widehat{Q}$ is diagonal. Indeed, as we discuss in Section 2.2, in this case
we never explicitly form the normal equations. Instead, we appropriately
utilize a $LDL^{\top}$ decomposition, and the fill-in produced by few dense
rows or columns of $A$ can be prevented.

######  Remark 2.

As we discuss later, a case of interest would be to only drop certain $k_{c}$
columns of $A$, which results in introducing at most $k_{c}$ eigenvalue
outliers in the preconditioned matrix. Similarly, only sparsifying certain
($k_{r}$) rows of $\widehat{M}$, introduces at most $2k_{r}$ non-unit
eigenvalues. Notice that dropping columns is expected to be more useful in
general, resulting in fewer outliers and possibly in greater gains (either in
terms of processing time or memory requirements). Indeed, notice that, on the
one hand, dropping dense columns of $A$ can result in a significant reduction
of the fill-in within a factorization of the normal equations, while, on the
other hand, dropping any column of $A$ corresponding to a small diagonal
element of $G$ yields a not too sizeable outlier. However, in certain special
applications one has to resort to sparsifying “problematic” rows. Indeed, we
refer the reader to [14, Section 4], where such a row-sparsifying strategy was
key to the efficient solution of fMRI classification problems.

###  2.2 An $LDL^{\top}$-based preconditioner

Next, we present an alternative to the preconditioner in (2.3). This approach
offers an alternative for dealing with a small set of dense columns or rows of
the matrix $A$, while remaining efficient when the approximate Hessian
$\widehat{Q}$, given in (2.1), is non-diagonal. More specifically, let us
divide the columns of the matrix $A$ into two mutually exclusive sets
$\mathcal{B}$ and $\mathcal{N}$, based solely on the magnitude of the
respective diagonal elements of the matrix $G$, and not on the density of the
columns of $A$. Then, using the column-dropping strategy presented in the
previous section (assuming that the variables corresponding to $\mathcal{N}$
are such that $Q^{(j,j)}\geq Q^{(i,i)}$, for all $j\in\mathcal{N}$ and all
$i\in\mathcal{B}$), we propose approximating the matrix $\widehat{M}$ by the
following preconditioner:

${P}_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}=A^{(:,\mathcal{B})}\widehat{G}^{(\mathcal{B},\mathcal{B})}\left(A^{(:,\mathcal{B})}\right)^{\top}+\delta
I_{m},$ (2.7)

which was proposed in [7], for the special case where $\widehat{G}$ was
diagonal. Notice that the block-separable structure of $\widehat{Q}$, given in
(2.1), implies that
$(\widehat{G}^{(\mathcal{B},\mathcal{B})})^{-1}=\widehat{Q}^{(\mathcal{B},\mathcal{B})}+\rho
I_{|\mathcal{B}|}$. Given our previous discussion, we would like to avoid
applying the inverse of this preconditioner by means of a Cholesky
decomposition, as a single dense column of $A$ in $\mathcal{B}$ could result
in dense Cholesky factors, while the potential non-diagonal nature of
$\widehat{Q}^{\left(\mathcal{B},\mathcal{B}\right)}$ might prevent us from
even efficiently forming it. Instead, we form an appropriate saddle point
system to compute the action of the approximated normal equations. More
specifically, given an arbitrary vector $y\in\mathbb{R}^{m}$, instead of
computing ${P}_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}^{-1}y$ using a
Cholesky decomposition, we can compute

$\underbrace{\begin{bmatrix}-\left(Q^{(\mathcal{B},\mathcal{B})}+\rho
I_{|\mathcal{B}|}\right)&\left(A^{(:,\mathcal{B})}\right)^{\top}\\\
A^{(:,\mathcal{B})}&\delta
I_{m}\end{bmatrix}}_{\widetilde{P}_{NE}}\begin{bmatrix}w_{1}\\\
w_{2}\end{bmatrix}=\begin{bmatrix}0_{|\mathcal{B}|}\\\ y\end{bmatrix},$ (2.8)

by means of an $LDL^{\top}$ decomposition of the previous saddle point matrix.
Then, we notice that returning $w_{2}$ is equivalent to computing
${P}_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}^{-1}y$.

Following the discussion in [2, Section 4], we know that using an $LDL^{\top}$
structure to factorize the matrix in (2.8) can result in significant memory
savings compared to the Cholesky decomposition of
${P}_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}$. Notice that in view of the
regularized nature of the systems under consideration (indeed, we have assumed
that $G$ is positive definite), we can use the result in [46], stating that
matrices like the one in (2.8) are _quasi-definite_ ; any symmetric
permutation of such matrices admits an $LDL^{\top}$ decomposition.

While this approach might seem expensive, it can provide significant time and
memory savings, especially in cases where $A^{(:,\mathcal{B})}$ contains dense
columns. In the previous section we discussed a strategy for alleviating this
issue, noting however that such a strategy can only be used to deal with a
small number of dense columns. On the contrary, if we have a sizeable subset
of the columns of $A^{(:,\mathcal{B})}$ that are dense, we could delay their
pivot order within the $LDL^{\top}$, thus significantly reducing the overall
fill-in of the decomposition factors, without introducing any eigenvalue
outliers in the preconditioned system. Of course, finding the optimal
permutation for the $LDL^{\top}$ decomposition is an NP-hard problem, however,
there have been developed several permutation heuristics tailored to such
symmetric decompositions. Moreover, in the $LDL^{\top}$ factorization, the
pivots are computed dynamically to ensure both stability and sparsity. In view
of the previous, the preconditioner based on solving (2.8) is expected to be
more stable than its counterpart based on the Cholesky decomposition. Finally,
difficulties arising from dense rows or in general “problematic” rows can also
be alleviated using a heuristic proposed in [28].

On the other hand, by using this approach we avoid explicitly forming the
preconditioner $P_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}$. This is
especially important in cases where $\widehat{Q}$ is non-diagonal, and thus
forming $P_{NE}$ can be extremely expensive. Hence, the approach presented in
this subsection allows us to utilize non-diagonal information in a practical
way when constructing an approximation for the matrix $Q$.

##  3 Regularized saddle point matrices

Let us now consider the regularized saddle point system in (1.1). In what
follows, we discuss two families of preconditioning strategies, noting their
advantages and disadvantages. All presented preconditioners are positive
definite in order to be usable within the MINRES method, which is a short-
recurrence iterative solver, suitable for solving symmetric indefinite or
quasi-definite systems. This allows us to avoid non-symmetric long-recurrence
solvers like the GMRES method.

###  3.1 Block-diagonal preconditioners

The most common approach is to employ a block-diagonal preconditioner (see [6,
7, 32, 45]). To construct such a preconditioner, we need approximations for
the $(1,1)$ block $F\coloneqq Q+\rho I_{n}$ of the coefficient matrix in
(1.1), and for its associated Schur complement $M=A\left(Q+\rho
I_{n}\right)^{-1}A^{\top}+\delta I_{m}$.

In this section, we assume that $Q$ is approximated as shown in (2.1), and
thus can potentially contain non-diagonal blocks. Concerning the approximation
of the Schur complement matrix $M$, we can employ the preconditioner
$P_{NE,\left(k_{c},k_{r}\right)}$ given in (2.3). Then, we can define the
following preconditioner for the coefficient matrix $K$ in (1.1):

$P_{AS,(k_{c},k_{r})}\coloneqq\begin{bmatrix}\widehat{Q}+\rho I_{n}&0_{n,m}\\\
0_{m,n}&P_{NE,\left(k_{c},k_{r}\right)}\end{bmatrix}\equiv\begin{bmatrix}\widehat{F}&0_{n,m}\\\
0_{m,n}&P_{NE,\left(k_{c},k_{r}\right)}\end{bmatrix}.$ (3.1)

For the rest of this subsection, let $P_{AS,(k_{c},k_{r})}\equiv P_{AS}$ and
$P_{NE,\left(k_{c},k_{r}\right)}\equiv P_{NE}$.

In order to analyze the spectrum of the preconditioned matrix $P_{AS}^{-1}K$,
we introduce some notation for simplicity of exposition. We work with positive
definite similarity transformations of the associated matrices, defined as

$\widetilde{F}=\widehat{F}^{-1/2}F\widehat{F}^{-1/2},\quad\widetilde{M}_{NE}=P_{NE}^{-1/2}\widehat{M}P_{NE}^{-1/2},$
(3.2)

where $\widehat{M}\coloneqq A\widehat{F}^{-1}A^{\top}+\delta I_{m}$. Then, we
set

$\begin{array}[]{lcllcllcl}\alpha_{NE}&=&\lambda_{\min}\left(\widetilde{M}_{NE}\right),{}&~{}\beta_{NE}&=&\lambda_{\max}\left(\widetilde{M}_{NE}\right),{}&~{}\kappa_{NE}&=&\dfrac{\beta_{NE}}{\alpha_{NE}},\\\\[6.00006pt]
\alpha_{F}&=&\lambda_{\min}\left(\widetilde{F}\right),{}&~{}\beta_{F}&=&\lambda_{\max}\left(\widetilde{F}\right),{}&~{}\kappa_{F}&=&\dfrac{\beta_{F}}{\alpha_{F}}.\end{array}$

From the definition of $\widehat{Q}$ given in (2.1), we can observe that
$\alpha_{F}\leq 1\leq\beta_{F}$ as

$\frac{1}{n}\sum_{i=1}^{n}\lambda_{i}\left(\widehat{F}^{-1}F\right)=\frac{1}{n}\textnormal{Tr}\left(\widehat{F}^{-1}F\right)=1.$

On the other hand, notice that $\alpha_{NE}$ and $\beta_{NE}$ are bounded
directly from Theorem 2.1. Indeed, from (3.2), we observe that we need to
bound the spectrum of an approximate preconditioned Schur complement matrix,
since $M$ has been substituted by $\widehat{M}$. This is exactly what the
analysis in Section 2.1 does. Below we provide a theorem analyzing the
spectral properties of the matrix $P_{AS}^{-1}K$.

######  Theorem 3.1.

The eigenvalues of $P_{AS}^{-1}K$ lie in the union of the following intervals:

$I_{-}=\left[-\beta_{F}-\sqrt{\beta_{NE}},-\alpha_{F}\right];\quad
I_{+}=\left[\frac{1}{2}\left(-\beta_{F}+\sqrt{\beta_{F}^{2}+4\alpha_{NE}}\right),1+\sqrt{\beta_{NE}-1}\right].$

###### Proof.

The proof, which follows by trivially extending [7, Theorem 3], is summarized
in Appendix A for completeness. ∎

The authors of [7] make use of the approximation
$\widehat{Q}=\textnormal{Diag}(Q)$, and
${P}_{NE}={P}_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}$ (where the latter
is defined in (2.7), with $\widehat{G}=(\widehat{Q}+\rho I_{n})^{-1}$).
Approximating $Q$ by its diagonal allows one to form the normal equations’
preconditioner (i.e. ${P}_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}$), thus
enabling the efficient use of a Cholesky factorization. However, there might
exist problems for which a better approximation of the matrix $Q$ is required.
In this case, one could consider a sparsified version of $Q$ like the one
given in (2.1). While this could lead to reasonable approximations, the
problem of fill-in introduced by $(\widehat{Q}+\rho I_{n})^{-1}$ in the Schur
complement approximation, ${P}_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}$,
would (in general) remain.

In order to address the previous issue, we make use of the $LDL^{\top}$-based
preconditioner defined in Section 2.2. As in Section 2.2, we divide the
columns of $A$ using the sets $\mathcal{B}$ and $\mathcal{N}$, where
$\mathcal{N}$ contains all indices corresponding to the largest diagonal
elements of $Q$. Assume that $Q$ is approximated by $\widehat{Q}$ as given in
(2.1), where the permutation matrix $\mathscr{P}_{c}$ traverses first the
column indices in $\mathcal{N}$. For example, a reasonable approximation would
be the following

$\mathscr{P}_{c}^{\top}\widehat{Q}\mathscr{P}_{c}=\begin{bmatrix}\textnormal{Diag}\left(Q^{(\mathcal{N},\mathcal{N})}\right)&0_{\lvert\mathcal{N}\rvert,\lvert\mathcal{B}\rvert}\\\
0_{\lvert\mathcal{B}\rvert,\lvert\mathcal{N}\rvert}&Q^{(\mathcal{B},\mathcal{B})}\end{bmatrix}.$
(3.3)

The effect of this approximation in the context of regularized IPMs has been
analyzed in detail in [39]. We notice that
$(Q^{(\mathcal{B},\mathcal{B})}+\rho I_{|\mathcal{B}|})^{-1}$ does not
introduce significant fill-in in the $(2,2)$ block of the preconditioner in
(3.1), as we implicitly invert this block using the methodology presented in
Section 2.2.

######  Remark 3.

Notice that further approximations can be employed here. In particular, we
could define a banded approximation of $Q$ and then employ the approximation
proposed earlier. The implicit inversion of the Schur complement, outlined in
Section 2.2, gives us complete freedom on how to approximate $Q$, and hence we
no longer rely on diagonal approximations. We return to this point in the
numerical experiments.

###  3.2 Factorization-based preconditioners

Finally, given the regularized nature of the systems under consideration, we
can construct factorization-based preconditioners for MINRES. In particular,
we can compute $K=LDL^{\top}$ (with $K$ in (1.1)), where $D$ is a diagonal
matrix (since $K$ is quasi-definite [46]) with $n$ negative and $m$ positive
elements on its diagonal. Then, by defining $P_{K}\coloneqq
L|D|^{\frac{1}{2}}$, the preconditioned saddle point matrix reads:

$P_{K}^{-1}KP_{K}^{-\top}=|D|^{-1}D,$

and hence contains only two distinct eigenvalues: $-1$ and $1$ [21, 34]. As
before, let us assume that we have available a splitting of the columns of $A$
such that $A\mathscr{P}_{c}=[A^{\mathcal{B}}\ \ A^{\mathcal{N}}]$, where
$\mathcal{B}$ contains indices corresponding to the smallest diagonal elements
of $Q$. Then, we can precondition $K$, left and right, by
$\widehat{P}_{K}\coloneqq\widehat{L}|\widehat{D}|^{\frac{1}{2}}$, where
$\widehat{K}=\widehat{L}\widehat{D}\widehat{L}^{\top}$ and:

$\widehat{K}\coloneqq\begin{bmatrix}-\widehat{Q}&\widehat{A}^{\top}\\\
\widehat{A}&\delta I_{m},\end{bmatrix},$

with $\widehat{A}\coloneqq[A^{\mathcal{B}}\ \
0_{m,|\mathcal{N}|}]\mathscr{P}_{c}^{\top}$, and $\widehat{Q}$ defined as in
(3.3). Notice that by setting several columns of $A$ to zero, as well as by
sparsifying the respective rows and columns of $Q$, the cost of applying the
inverse of $\widehat{K}$ is significantly reduced when compared to that
required to apply the inverse of $K$.

Further limited-memory versions of this preconditioner can be employed, e.g.
by using the methodologies presented in [34, 44]. Other approximations of the
blocks of $\widehat{K}$, based on the structure of the problem at hand, could
also be possible, as already mentioned in the previous subsection.

We should note, however, that this approach is less stable than the approach
presented in Section 3.1. This is because we are required to use only diagonal
pivots during the $LDL^{\top}$ decomposition for this methodology to work
(indeed, notice that the presence of a non-diagonal matrix $D$ in the
factorization of $K$ would not allow the use of such a preconditioning
strategy). If the regularization parameters $\delta$ or $\rho$ have very small
values, the stability of the factorization could be compromised (since we
enforce the use of only diagonal $1\times 1$ pivots), and we would have to
heavily rely on stability introduced by means of uniform [43] or weighted
regularization [1]. On the other hand, the methodology presented in Section
3.1 would not be affected by the occasional use of $2\times 2$ pivots within
the $LDL^{\top}$ factorization for the implicit inversion of the approximate
Schur complement. Of course the latter is not the case if the “Analyze” phase
of the factorization (used to determine the pivot order) is performed
separately, however, the subset of columns in $\mathcal{B}$ may change
significantly from one iteration to the next, making this strategy less
attractive. Nevertheless, this factorization-based approach can be more
efficient than the approach presented in Section 3.1, when solving certain
non-separable convex programming problems. This is because the approach in
Section 3.1 requires the computation of an $LDL^{\top}$ decomposition of the
coefficient matrix $\widetilde{P}_{NE}$ in (2.8) (with potential $2\times 2$
pivots) as well as a Cholesky decomposition of $\widehat{Q}+\rho I_{n}$ (or
some iterative scheme which could be application-dependent, as in [38]).

##  4 Regularized IPMs: numerical results

Let us now focus on the case of the regularized saddle point systems (and
their respective normal equations) arising from the application of regularized
IPMs on convex programming problems. The MATLAB code, which is based on the
IP-PMM presented in [40, 41], can be found on
GitHub111https://github.com/spougkakiotis/IP-PMM_QP_Solver.

In all the presented experiments, a 6-digit accurate solution is requested.
The reader is referred to [7, Section 4], and [40, Section 5] for the
implementation details of the algorithm (such as termination criteria, the
employed predictor–corrector scheme for the solution of the Newton system, as
well as the tuning of the algorithmic regularization parameters). The
associated iterative methods (i.e. PCG or MINRES) are adaptively terminated if
the following accuracy is reached:
$\frac{\min\\{10^{-3},\max\\{10^{-1}\cdot\mu_{k},\texttt{tol}\\}\\}}{\max\\{1,\|\texttt{rhs}\|\\}}$,
where $\texttt{tol}=10^{-6}$, $\mu_{k}$ is the barrier parameter at iteration
$k$, and rhs is the right hand side of the system being solved. This adaptive
stopping criterion is based on the developments in [12]. When PCG is employed
we allow at most 100 iterations per linear system solved, while for MINRES up
to 200 iterations are allowed. If the maximum number of Krylov iterations is
reached, an inexact Newton direction is accepted if it is at least 3-digit
accurate. Any Cholesky decomposition is computed via the chol function of
MATLAB. When an $LDL^{\top}$ decomposition is employed, we utilize the ldl
function of MATLAB. In this case, the minimum pivot threshold is adaptively
set to $\text{pivot}_{\text{thr}}=0.1\cdot\min\\{\delta,\rho,10^{-4}\\}$. This
is done to ensure that no $2\times 2$ pivots are used during the
factorization, ensuring that the factorization remains efficient. However, in
the context of the preconditioner in Section 2.2, where $2\times 2$ pivots can
safely be used (unlike the preconditioner presented in Section 3.2, which
requires the use of $1\times 1$ pivots), this mechanism is turned off when
$\min\\{\delta,\rho\\}\leq 10^{-8}$, and we set
$\text{pivot}_{\text{thr}}=10^{-6}$ to ensure stability. All the presented
experiments were run on a PC with a 2.2GHz Intel Core i7 processor (hexa-
core), 16GB RAM, run under the Windows 10 operating system. The MATLAB version
used was 2019a.

###  4.1 Linear programming

Let us initially focus on Linear Programming (LP) problems of the following
form:

$\underset{x\in\mathbb{R}^{n}}{\text{min}}\ c^{\top}x,\ \ \text{s.t.}\ Ax=b,\
x^{\mathcal{I}}\geq 0,\ x^{\mathcal{F}}\ \textnormal{free},$ (LP)

where $A\in\mathbb{R}^{m\times n}$, $\mathcal{I}\cap\mathcal{F}=\emptyset$,
and $\mathcal{I}\cup\mathcal{F}=\\{1,\ldots,n\\}$. Applying regularized IPMs
to problems like (LP), one often solves a regularized normal equations system
at every iteration. Such systems have a coefficient matrix of the following
form:

$M=AGA^{\top}+\delta I_{m},\qquad G^{(i,i)}=\begin{cases}\frac{1}{\rho}&\
\textnormal{if}\ i\in\mathcal{F},\\\ \frac{1}{\rho+z^{i}/x^{i}}&\
\textnormal{if}\ i\in\mathcal{I},\end{cases}$

where $\delta,\ \rho>0$ and $z\in\mathbb{R}^{n}$ (where $z^{\mathcal{I}}\geq
0$, $z^{\mathcal{F}}=0$) are the dual slack variables. Notice that the IPM
penalty parameter $\mu$ is often tuned as
$\mu=\frac{(x^{\mathcal{I}})^{\top}z^{\mathcal{I}}}{n}$ and we expect that
$\mu\rightarrow 0$. As already mentioned in the introduction, the variables
are naturally split into “basic”–$\mathcal{B}$, “non-basic”–$\mathcal{N}$, and
“undecided”–$\mathcal{U}$. Hence, as IPMs progress towards optimality, we
expect the following partition of the quotient
$\frac{x^{\mathcal{I}}}{z^{\mathcal{I}}}$:

$\displaystyle\forall j\in\mathcal{N}:x^{j}\rightarrow 0,\quad
z^{j}\rightarrow\widehat{z}^{j}>0$ $\displaystyle\Rightarrow\quad$
$\displaystyle\frac{x^{j}}{z^{j}}=\frac{x^{j}z^{j}}{(z^{j})^{2}}=\bm{\Theta}(\mu),$
$\displaystyle\forall j\in\mathcal{B}:x^{j}\rightarrow\widehat{x}^{j}>0,\quad
z^{j}\rightarrow 0$ $\displaystyle\Rightarrow\quad$
$\displaystyle\frac{x^{j}}{z^{j}}=\frac{(x^{j})^{2}}{x^{j}z^{j}}=\bm{\Theta}(\mu^{-1}),$
$\displaystyle\forall j\in\mathcal{U}:x^{j}=\bm{\Theta}(1),\quad
z^{j}=\bm{\Theta}(1)$ $\displaystyle\Rightarrow\quad$
$\displaystyle\frac{x^{j}}{z^{j}}=\bm{\Theta}(1),$

where $\mathcal{N}$, $\mathcal{B}$, and $\mathcal{U}$ are mutually disjoint,
and $\mathcal{N}\cup\mathcal{B}\cup\mathcal{U}=\mathcal{I}$. For the rest of
this section, we assume that $\delta=\bm{\Theta}(\rho)=\bm{\Theta}(\mu)$. This
assumption is based on the developments in [40, 41], where a polynomially
convergent regularized IPM is derived for convex quadratic and linear positive
semi-definite programming problems, respectively. Following [7], we could
precondition the matrix $M$ using the following matrix:

$P_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}=A^{\left(\colon,\mathcal{R}\right)}G^{\left(\mathcal{R},\mathcal{R}\right)}\left(A^{\left(\colon,\mathcal{R}\right)}\right)^{\top}+\delta
I_{m},$ (4.1)

where $\mathcal{R}\coloneqq\mathcal{F}\cup\mathcal{B}\cup\mathcal{U}$. Then,
by [7, Theorem 1] (or by applying Theorem 2.1), we obtain:

$\lambda_{\max}\left(P_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}^{-1}M\right)\leq
1+\frac{\underset{j\in\mathcal{N}}{\max}\left(G^{(j,j)}\right)}{\delta}\sigma_{\max}^{2}(A),\qquad\lambda_{\min}\left(P_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}^{-1}M\right)\geq
1.$

The preconditioner in (4.1) is a special case of the preconditioner defined in
Section 2. Indeed, as already indicated by its notation, it can be derived by
setting $k_{r}=0$ and then by dropping all columns belonging to $\mathcal{N}$,
i.e. we set $k_{c}=|\mathcal{N}|$ and we drop the $k_{c}$ columns of $A$
corresponding to the smallest diagonal elements of $G$. Notice that in the
linear programming case, the matrix $\widehat{M}$ that is analyzed in Section
2 coincides with $M$, since $G$ is diagonal.

From our previous remarks, we notice that

$\underset{j\in\mathcal{N}}{\max}\left(G^{(j,j)}\right)=\bm{\Theta}(\mu)=\bm{\Theta}(\delta)$

implies that the spectrum of the preconditioned matrix remains bounded and is
asymptotically independent of $\mu$ (assuming that $\delta=\bm{\Theta}(\mu)$).
While this preconditioner performs very well in practice (see [7, Section 4]),
it can be expensive to compute in certain cases, as its inverse needs to be
applied by means of a Cholesky decomposition. To that end, we propose to
further approximate this matrix as indicated in Section 2. This idea is based
on the fact that the Preconditioned Conjugate Gradient (PCG) method is
expected to converge in a small number of iterations, if the preconditioned
system matrix can be written as:

$P^{-1}M=I_{m}+U+V,$

where $P$ is the preconditioner, $M$ is the normal equations matrix, $U$ is a
low-rank matrix, and $V$ is a matrix with small norm. In our case, dropping
the part of the normal equations corresponding to $\mathcal{N}$ contributes
the small-norm term (that is,
$V=A^{\left(\colon,\mathcal{N}\right)}G^{\left(\mathcal{N},\mathcal{N}\right)}(A^{\left(\colon,\mathcal{N}\right)})^{\top}$,
the norm of which is of the order of magnitude of $\mu$), and furthermore
dropping a few dense columns (or sparsifying certain rows) contributes the
low-rank term (indeed, as already shown in Theorem 2.1, dropping $k_{c}$ dense
columns of $A$ and sparsifying $k_{r}$ dense rows of $M$ yields at most
$2k_{r}+k_{c}$ outliers, and thus $\textnormal{rank}(U)\leq 2k_{c}+k_{r}$,
where, in this case, $k_{c}$ does not account for columns corresponding to the
index set $\mathcal{N}$).

To construct such a preconditioner, we first need to note that $\mathcal{R}$
will change at every IPM iteration. However, we can heuristically choose which
columns to drop (and/or rows to sparsify) based on the sparsity pattern of
$A$. To that end, at the beginning of the optimization procedure, we count the
number of non-zeros of each column and row of $A$, respectively. These can
then be used to sort the columns and rows of $A$ in descending order of their
number of non-zero entries. These sorted columns and rows can easily be
represented by means of two permutation matrices $\mathscr{P}_{c}$ and
$\mathscr{P}_{r}$. We note that this is a heuristic, and it is not guaranteed
to identify the most “problematic” columns or rows (which can be sources of
difficulty for IPMs). For a discussion on such heuristics, and alternatives,
the reader is referred to [2, Section 4], and the references therein.

####  4.1.1 Numerical results

Initially, we present some results to show the effect of dropping dense
columns of $A$ and then of sparsifying dense rows of $M$ using the strategy
outlined in Section 2. Then, we present a comparison between the
preconditioner in (2.3) (that is, $P_{NE,(k_{c},k_{r})}$), the preconditioner
given in (2.7) or (4.1) (denoted as
$P_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}$), noting that (2.7) and (4.1)
are equivalent, and the one in (2.8) (that is, $\widetilde{P}_{NE}$). Notice
that in the linear programming case, employing $\widetilde{P}_{NE}$ and
$P_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}$ should yield identical results
in exact arithmetic. The difference between these preconditioning strategies,
is that in the latter case (i.e. $\widetilde{P}_{NE}$) the action of the
former preconditioner (i.e. $P_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}$)
is computed implicitly by means of an $LDL^{\top}$ factorization of
$\widetilde{P}_{NE}$, as indicated in Section 2.2.

##### Dropping dense columns versus factorizing directly.

We run IP-PMM on all problems from the Netlib collection that have some dense
columns, where dense is defined in this case to be a column with at least 15%
non-zero elements. We note that these were the only problems within the Netlib
collection having any columns with such a density of non-zero elements. We
compare an IP-PMM using Cholesky factorization for the solution of the
associated Newton system (Chol.), with an IP-PMM that uses the preconditioner
$P_{NE,\left(k_{c},0\right)}$ presented in Section 2 alongside PCG. The latter
method is only allowed to drop dense columns (at most 30) to create the
preconditioner (and thus $k_{c}\leq 30$). The results are collected in Table
1, where $k_{c}$ denotes the number of dense columns that are dropped to
create the preconditioner, nnz denotes the number of non-zero elements present
in the Cholesky factor, Avg. Krylov It. denotes the average number of Krylov
iterations performed from the inexact approach, while Krylov Last denotes the
number of Krylov iterates performed in the last IP iteration of the inexact
approach.

Table 1: The effect of dropping dense ($>15\%$) columns of $A$ (Netlib
collection).

Name | $\mathbf{k_{c}}$ | nnz | time (s) | Avg. Krylov It. | Krylov Last
---|---|---|---|---|---
Chol. | $P_{NE,\left(k_{c},0\right)}$ | Chol. | $P_{NE,\left(k_{c},0\right)}$
BLEND | 5 | 1,006 | 736 | 0.05 | 0.05 | 14.50 | 14
FIT1P | $20$ | 197,676 | 26,706 | 1.16 | 0.32 | 18.67 | 35
FIT2P | 16 | 4,516,500 | 1,962,616 | 27.17 | 16.21 | 17.33 | 33
FORPLAN | $8$ | 3,810 | 2,918 | 0.11 | 0.11 | 2.81 | 3
ISRAEL | $30$ | 12,261 | 1,744 | 0.14 | 0.30 | 53.50 | 100
SEBA | $14$ | 55,937 | 2,238 | 0.35 | 0.18 | 13.83 | 15

From Table 1 we can immediately see that certain dense columns present in the
constraint matrix $A$ can have a significant impact on the sparsity pattern of
the Cholesky factors. This is a well-known fact (see for example the
discussion in [2, Section 4]). Notice that the Netlib collection contains only
small- to medium-scale instances. For such problems, memory is not an issue,
and hence direct methods tend to be faster than their iterative alternatives
(like PCG). Despite the small size of the presented problems, we can see
tremendous memory savings (and even a decrease in CPU time) for problems
FIT1P, FIT2P, and SEBA, by eliminating only a small number of dense columns.
On the other hand, for problems where we observe an increase in CPU time (e.g.
see ISRAEL), the memory savings can be significant, making this acceptable.

##### Sparsifying dense rows versus factorizing directly.

Next, we consider the case where the inexact version of IP-PMM is only allowed
to sparsify dense rows, where dense is defined in this case to be a row with
at least 25% non-zero elements.

Before moving to the numerical results, let us note some differences between
sparsifying rows of $M$ and dropping columns of $A$. Firstly, as we have shown
in Section 2, sparsifying $k$ rows can potentially introduce twice as many
outliers, while dropping $k$ columns introduces at most $k$ eigenvalue
outliers. Furthermore, the potential density induced in the Cholesky factors
by a single dense column is usually more significant than that introduced by a
single dense row. However, we cannot know in advance how effective the
dropping of a column will be. On the other hand, sparsifying dense rows of $M$
introduces a certain separability in the approximate normal equations matrix,
allowing us to estimate very well the memory savings.

In Table 2 we compare the direct IP-PMM, to its inexact version, the
preconditioner of which is only allowed to sparsify at most 30 dense rows of
$M$. Any problem with at least one dense row from the Netlib collection is
considered.

Table 2: The effect of dropping dense ($>25\%$) rows of $M$ (Netlib
collection).

Name | $\mathbf{k_{r}}$ | nnz | time (s) | Avg. Krylov It. | Krylov Last
---|---|---|---|---|---
Chol. | $P_{NE,\left(0,k_{r}\right)}$ | Chol. | $P_{NE,\left(0,k_{r}\right)}$
BEACONFD | $17$ | 2,903 | 1,475 | 0.07 | 0.09 | 13.59 | 23
BLEND | 1 | 1,006 | 959 | 0.05 | 0.04 | 2.88 | 3
D6CUBE | 6 | 55,179 | 52,757 | 0.27 | 0.44 | 12.73 | 18
FIT1D | $11$ | 14,726 | 4,973 | 0.22 | 0.38 | 23.50 | 46
FIT2D | $12$ | 139,843 | 49,513 | 1.51 | 2.47 | 21.30 | 44
ISRAEL | $3$ | $12,261$ | $11,758$ | 0.14 | 0.14 | 6.96 | 8
STANDATA | 1 | 3,416 | 3,168 | 0.07 | 0.07 | 2.94 | 3
STANDGUB | 1 | 3,418 | 3,170 | 0.08 | 0.08 | 2.98 | 3
STANDMPS | 1 | 5,529 | 5,185 | 0.09 | 0.09 | 2.91 | 3
WOOD1P | $27$ | 19,088 | 13,879 | 0.34 | 0.44 | 14.95 | 13

From Table 2 we can observe that the required memory to form the Cholesky
factors is consistently decreased but CPU time is often increased by the row-
dropping strategy. We should note that an increase in CPU time usually relates
to the size of the problems under consideration, and CPU time as well as a
memory advantages can be observed if the problem is sufficiently large with
sufficiently many dense rows. In particular, this row sparsifying strategy was
successfully used within IP-PMM in [14, Section 4] in order to tackle fMRI
sparse approximation problems (in which the constraint matrix contains
thousands of dense rows). Memory requirements were significantly lowered,
allowing this inexact version to outperform its exact counterpart, while being
competitive with standard state-of-the-art first-order methods used to solve
such problems.

##### The Cholesky versus the $LDL^{\top}$ approach.

Let us now provide some numerical evidence for the potential benefits and
drawbacks of the approach presented in Section 2.2 over that presented in
Section 2.1. To that end, we run three inexact versions of IP-PMM, on some of
the most challenging instances within the Netlib collection. The first
approach uses the preconditioner given in (2.3) (denoted as
$P_{NE,(k_{c},k_{r})}$, allowing at most 15 dense columns/rows to be
dropped/sparsified), the second uses the preconditioner given in (2.7)
(denoted as $P_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}$; notice that this
is the same as the former preconditioner, without employing the strategy of
dropping/sparsifying dense columns/rows), while the third version uses the
preconditioner in (2.8), denoted as $\widetilde{P}_{NE}$. In all three cases
the set $\mathcal{B}$, used to decide which columns are dropped irrespectively
of their density, is determined as indicated at the beginning of this section.

Table 3: Cholesky-based versus the $LDL^{\top}$-based preconditioner (Netlib
collection).

Name | Krylov Its. | max nnz | time (s)
---|---|---|---
$P_{NE,(k_{c},k_{r})}$ | $P_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}$ | $\widetilde{P}_{NE}$ | $P_{NE,(k_{c},k_{r})}$ | $P_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}$ | $\widetilde{P}_{NE}$ | $P_{NE,(k_{c},k_{r})}$ | $P_{NE,\left(\lvert\mathcal{N}\rvert,0\right)}$ | $\widetilde{P}_{NE}$
AGG | 2,966 | 2,966 | 2,268 | $1.60\cdot 10^{4}$ | $1.60\cdot 10^{4}$ | $1.40\cdot 10^{4}$ | 0.42 | 0.42 | 0.58
DFL001 | 4,778 | 4,778 | 4,969 | $1.55\cdot 10^{6}$ | $1.55\cdot 10^{6}$ | $1.54\cdot 10^{6}$ | 8.40 | 8.40 | 31.51
FIT1P | 1,166 | 636 | 491 | $1.20\cdot 10^{5}$ | $2.00\cdot 10^{5}$ | $1.36\cdot 10^{4}$ | 0.58 | 1.56 | 0.51
FIT2P | 2,555 | 1,880 | 1,867 | $4.19\cdot 10^{6}$ | $4.60\cdot 10^{6}$ | $9.4\cdot 10^{4}$ | 29.45 | 43.10 | 15.95
PILOT | 3,558 | 3,558 | 3,530 | $1.99\cdot 10^{5}$ | $1.99\cdot 10^{5}$ | $4.02\cdot 10^{5}$ | 3.79 | 3.79 | 9.16
QAP12 | 1,583 | 1,583 | 1,591 | $2.48\cdot 10^{6}$ | $2.48\cdot 10^{6}$ | $1.61\cdot 10^{6}$ | 2.09 | 2.09 | 5.24
QAP15 | 1,704 | 1,704 | 1,708 | $8.83\cdot 10^{6}$ | $8.83\cdot 10^{6}$ | $5.28\cdot 10^{6}$ | 20.56 | 20.56 | 25.42
SEBA | 2,678 | 2,396 | 2,313 | $2.18\cdot 10^{3}$ | $5.54\cdot 10^{4}$ | $7.5\cdot 10^{3}$ | 0.33 | 0.80 | 0.84

From Table 3 we can observe that the $LDL^{\top}$-based preconditioner can
provide substantial (memory and/or CPU time) benefits for certain problems
(e.g. see problems FIT1P, FIT2P, QAP12, QAP15, SEBA). Nevertheless, we should
note that this approach is usually slower, albeit more stable (as a pivot re-
ordering is computed at every iteration, and the pivots of the $LDL^{\top}$
factorization are chosen to ensure stability as well as efficiency). We
observe that instances AGG, DFL001, PILOT did not benefit from the use of this
strategy, neither in terms of efficiency nor memory requirements, despite a
comparable number of Krylov iterations. This comes in line with our
observations in Section 2.2, since neither of the aforementioned instances
contains any dense rows or columns. Notice that the stability and efficiency
of the preconditioner $\widetilde{P}_{NE}$, depends heavily on the choice for
the threshold for the ldl function of MATLAB. Larger values imply better
stability, however at the cost of efficiency, since more $2\times 2$ pivots
will be chosen during the $LDL^{\top}$ factorization. The stability of this
approach can be guaranteed by using a large-enough pivot threshold.
Additionally, there are instances without dense rows or columns (see QAP12,
QAP15), in which the $LDL^{\top}$-based preconditioner (i.e.
$\widetilde{P}_{NE}$) provides significant advantages in terms of memory
requirements. Finally, we note that for problems AGG, PILOT, QAP12, and QAP15,
the two Cholesky-based variants are exactly the same, as no dense columns or
rows were present.

There is a long-standing discussion on the comparison between the Cholesky and
the $LDL^{\top}$ decompositions. The former tend to be faster and usually
easier to implement, while the latter tend to be slower, more stable, and more
general. For more on this subject, the reader is referred to [2, Section 4]
and the references therein.

###  4.2 Convex quadratic programming

Next, we consider problems of the following form:

$\min_{x\in\mathbb{R}^{n}}c^{\top}x+\frac{1}{2}x^{\top}Hx,\quad\textnormal{s.t.}\
Ax=b,\ x^{\mathcal{I}}\geq 0,\ x^{\mathcal{F}}\textnormal{ free},$ (QP)

where $H\in\mathbb{R}^{n\times n}$ is the positive semi-definite Hessian
matrix. Let us notice that a similar partitioning of the variables as that
presented in Section 4.1 also holds in this case. Hence, the index set
$\mathcal{N}$ guides us on which columns of $A$ to drop. In the case where $H$
is either diagonal, or can be well-approximated by a diagonal, the discussion
of Section 4.1, about dropping dense columns of $A$ (or sparsifying dense rows
of the approximated Schur complement $\widehat{M}$ given in (2.2)), also
applies here.

In what follows we make use of three different preconditioners. We compare the
two block-diagonal preconditioners given in Section 3.1. The first is called
$P^{C}_{AS,\left(k_{c},k_{r}\right)}$ (where the superscript $C$ stands for
Cholesky, which is used to invert the $(2,2)$ block of this preconditioner),
and employs a diagonal approximation for $Q$, allowing one to drop dense
columns and/or sparsify dense rows as shown in Section 2.1, and the second is
called $P^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$ (where the
superscript $L$ stands for $LDL^{\top}$), and employs a block-diagonal
approximation of $Q$, using the implicit inversion of the Schur complement
proposed in Section 2.2. The block-diagonal preconditioners are also compared
against the factorization-based preconditioner presented in Section 3.2,
termed as $\widehat{P}_{K}$.

####  4.2.1 Numerical results

In the following experiments, we employ MINRES to solve the associated Newton
systems. Initially, we present the comparison of the three preconditioning
strategies over some problems from the Maros–Mészáros collection of convex
quadratic programming problems. Then, the two block-diagonal preconditioning
approaches are compared over some Partial Differential Equation (PDE)
optimization problems.

##### Maros–Mészáros collection.

In Table 4, we report on the runs of the three methods on a diverse set of
non-separable instances within the Maros–Mészáros test set.

Table 4: Comparison of QP preconditioners (Maros–Mészáros collection).

Name | Krylov Its. | max nnz | time (s)
---|---|---|---
$P^{C}_{AS,\left(k_{c},k_{r}\right)}$ | $P^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$ | $\widehat{P}_{K}$ | $P^{C}_{AS,\left(k_{c},k_{r}\right)}$ | $P^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$ | $\widehat{P}_{K}$ | $P^{C}_{AS,\left(k_{c},k_{r}\right)}$ | $P^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$ | $\widehat{P}_{K}$
CVXQP2_L | 3,241 | 2,057 | 3,078 | $5.06\cdot 10^{4}$ | $2.35\cdot 10^{6}$ | $1.47\cdot 10^{6}$ | 9.01 | 29.74 | 33.52
CVXQP2_M | 3,019 | 2,588 | 3,128 | $5.04\cdot 10^{3}$ | $1.15\cdot 10^{5}$ | $6.15\cdot 10^{4}$ | 1.02 | 1.54 | 1.44
DUAL3 | 911 | 543 | 992 | $1.12\cdot 10^{2}$ | $6.77\cdot 10^{3}$ | $6.77\cdot 10^{3}$ | 0.15 | 0.11 | 0.14
GOULDQP3 | 1,236 | 1,039 | 814 | $4.89\cdot 10^{3}$ | $1.05\cdot 10^{4}$ | $1.05\cdot 10^{4}$ | 0.31 | 0.37 | 0.27
MOSARQP2 | 752 | 730 | 803 | $1.71\cdot 10^{4}$ | $2.18\cdot 10^{4}$ | $2.25\cdot 10^{4}$ | 0.12 | 0.16 | 0.18
POWELL20 | 1,531 | 1,537 | 1,538 | $2.10\cdot 10^{4}$ | $6.20\cdot 10^{4}$ | $6.20\cdot 10^{4}$ | 1.64 | 2.41 | 5.24
Q25FV47 | 6,213 | 5,994 | 8,117 | $2.13\cdot 10^{4}$ | $6.68\cdot 10^{4}$ | $1.14\cdot 10^{5}$ | 2.12 | 3.75 | 4.61
QETAMACRO | 4,901 | 4,661 | 6,378 | $1.13\cdot 10^{4}$ | $2.62\cdot 10^{4}$ | $2.46\cdot 10^{4}$ | 0.71 | 1.40 | 1.60
QISRAEL | 4,516 | 3,367 | 5,920 | $2.17\cdot 10^{2}$ | $4.60\cdot 10^{3}$ | $4.84\cdot 10^{3}$ | 0.49 | 0.61 | 0.80
QSHIP12L | 5,664 | 5,070 | 5,803 | $1.09\cdot 10^{4}$ | $2.06\cdot 10^{5}$ | $2.10\cdot 10^{5}$ | 3.01 | 4.92 | 5.69
STCQP1 | 3,410 | 2,543 | 2,691 | $7.00\cdot 10^{5}$ | $1.25\cdot 10^{5}$ | $1.26\cdot 10^{5}$ | 7.48 | 5.49 | 5.21
STCQP2 | 3,074 | 2,704 | 1,918 | $5.44\cdot 10^{4}$ | $2.13\cdot 10^{5}$ | $2.13\cdot 10^{5}$ | 4.30 | 6.40 | 3.96
UBH1 | 681 | 679 | 700 | $8.40\cdot 10^{4}$ | $2.13\cdot 10^{5}$ | $2.13\cdot 10^{5}$ | 1.97 | 3.78 | 4.17

From Table 4, one can observe that most of the time
$P^{C}_{AS,\left(k_{c},k_{r}\right)}$ is rather inexpensive, and naturally
requires some additional Krylov iterations. On the other hand,
$P^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$ delivers faster
convergence of the Krylov solver, at the cost of additional memory (since we
utilize non-diagonal Hessian information). However, while the same is true for
most problems when employing $\widehat{P}_{K}$, the latter can be prone to
numerical inaccuracy (since we do not allow the use of $2\times 2$ pivots in
the $LDL^{\top}$ factorization). Whether the use of non-diagonal Hessian
information is beneficial should depend on the problem under consideration. In
the above experiments, this proved to be beneficial for only 4 out of the 13
instances tested (that is DUAL3, GOULDQP3, STCQP1, STCQP2). Nevertheless, we
can observe that all three approaches are competitive, while
$P^{C}_{AS,\left(k_{c},k_{r}\right)}$ and
$P^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$ are both very stable.

##### PDE-constrained optimization instances.

Next, we compare the preconditioning approaches on some Partial Differential
Equation (PDE) optimization problems. In particular, we consider the
$\text{L}^{1}/\text{L}^{2}$-regularized Poisson control problem, as well as
the $\text{L}^{1}/\text{L}^{2}$-regularized convection–diffusion control
problem with control bounds. We should emphasize at this point that while
beskope preconditioners have been created for PDE problems of this form, here
we treat the discretized problems as if we hardly knew anything about their
structure, to demonstrate the generality of the approaches presented in this
paper.

We consider problems of the following form:

$\begin{split}\min_{\mathrm{y},\mathrm{u}}\ &\
\mathrm{J}(\mathrm{y}(\bm{x}),\mathrm{u}(\bm{x}))\coloneqq\frac{1}{2}\|\rm{y}-\bar{\rm{y}}\|_{L^{2}(\Omega)}^{2}+\frac{\alpha_{1}}{2}\|\rm{u}\|_{L^{1}(\Omega)}+\frac{\alpha_{2}}{2}\|\rm{u}\|_{L^{2}(\Omega)}^{2},\\\
\text{s.t.}\ &\
\mathrm{D}\mathrm{y}(\bm{x})=\mathrm{u}(\bm{x})+\mathrm{g}(\bm{x}),\\\ &\
\mathrm{u_{a}}(\bm{x})\leq\mathrm{u}(\bm{x})\leq\mathrm{u_{b}}(\bm{x}),\end{split}$
(4.2)

where $(\rm{y},\rm{u})\in\text{H}^{1}(\Omega)\times\text{L}^{2}(\Omega)$,
$\mathrm{D}$ denotes some linear differential operator associated with the
differential equation, $\bm{x}$ is a $2$-dimensional spatial variable, and
$\alpha_{1},\ \alpha_{2}\geq 0$ denote the regularization parameters of the
control variable. We note that other variants for $\mathrm{J}({\rm y},{\rm
u})$ are possible, including measuring the state misfit and/or the control
variable in other norms, as well as alternative weightings within the cost
functionals. In particular, the methods tested here also work well for
$\text{L}^{2}$-norm problems (e.g. see [36]). We consider problems of the form
of (4.2) to create an extra level of difficulty for our solvers.

The problem is considered on a given compact spatial domain $\Omega$, where
$\Omega\subset\mathbb{R}^{2}$ has boundary $\partial\Omega$, and is equipped
with Dirichlet boundary conditions. The algebraic inequality constraints are
assumed to hold a.e. on $\Omega$. We further note that ${\rm u_{a}}$ and ${\rm
u_{b}}$ may take the form of constants, or functions in spatial variables,
however we restrict our attention to the case where these represent constants.

Problems in the form of (4.2) are often solved numerically, by means of a
discretization method. In the following experiments we employ the Q1 finite
element discretization implemented in
IFISS222https://personalpages.manchester.ac.uk/staff/david.silvester/ifiss/default.htm
(see [18, 19]). Applying the latter yields a sequence of non-smooth convex
programming problems, which can be transformed to the smooth form of (QP), by
introducing some auxiliary variables to deal with the $\ell_{1}$ terms
appearing in the objective (see [38, Section 2]). In order to restrict the
memory requirements of the approach, we consider an additional approximation
of $H$ in the preconditioner
$P^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$. In the cases under
consideration, the resulting Hessian matrix takes the following form:

$H=\begin{bmatrix}J_{M}&0_{d,d}&0_{d,d}\\\
0_{d,d}&\alpha_{2}J_{M}&-\alpha_{2}J_{M}\\\
0_{d,d}&-\alpha_{2}J_{M}&\alpha_{2}J_{M}\end{bmatrix},$

where $J_{M}$ is the mass matrix of size $d$. When non-diagonal Hessian
information is utilized within the preconditioner, we approximate each block
of $H$ by its diagonal (i.e. $\widetilde{J}_{M}=\textnormal{Diag}(J_{M})$; an
approximation which is known to be optimal [47]). The resulting matrix is then
further approximated as discussed in Section 3. From now on, the $LDL^{\top}$
preconditioner, which is based on an approximation of
$P^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$, is referred to as
$\widehat{P}^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$, in order to
stress the additional level of approximation employed within the Hessian
matrix. For these examples, the preconditioning strategy based on
$\widehat{P}_{K}$ (given in Section 3.2) behaved significantly worse, and
hence was not included in the numerical results. The preconditioner
$\widehat{P}^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$ can be useful in
that it allows us to employ block-diagonal preconditioners of which the Schur
complement takes into account non-diagonal information of the Hessian matrix
$H$. In certain cases, this can result in a faster convergence of IP-PMM, as
compared to ${P}^{C}_{AS,\left(k_{c},k_{r}\right)}$ (see Table 6).

The first problem that we consider is the two-dimensional
$\text{L}^{1}/\text{L}^{2}$-regularized Poisson optimal control problem, with
bound constraints on the control and free state, posed on the domain
$\Omega=(0,1)^{2}$. Following [38, Section 5.1], we consider the constant
control bounds $\rm{u_{a}}=-2$, $\rm{u_{b}}=1.5$, and the desired state
$\bar{\rm{y}}=\sin(\pi x_{1})\sin(\pi x_{2})$. In Table 5, we fix
$\alpha_{2}=10^{-2}$ (which we find to be the most numerically interesting
case), and we present the runs of the method using the different
preconditioning approaches, with increasing problem size, and varying
$\text{L}^{1}$ regularization parameter (that is $\alpha_{1}$). To reflect the
change in the grid size, we report the number of variables of the optimization
problem after transforming it to the IP-PMM format. We also report the overall
number of Krylov iterations required for IP-PMM to converge (and the number of
IP-PMM iterations in brackets), the maximum number of non-zeros stored in
order to apply the inverses of the associated preconditioners, as well as the
required CPU time.

Table 5: Comparison of QP preconditioners (Poisson Control: problem size and
varying regularization).

$\bm{n}$ | $\bm{\alpha_{1}}$ | Krylov (IP) Its. | max nnz | time (s)
---|---|---|---|---
${P}^{C}_{AS,\left(k_{c},k_{r}\right)}$ | $\widehat{P}^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$ | ${P}^{C}_{AS,\left(k_{c},k_{r}\right)}$ | $\widehat{P}^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$ | ${P}^{C}_{AS,\left(k_{c},k_{r}\right)}$ | $\widehat{P}^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$
$2.11\cdot 10^{4}$ | $10^{-2}$ | 1,353 (13) | 868 (13) | $3.88\cdot 10^{5}$ | $4.65\cdot 10^{5}$ | 5.81 | 5.91
$10^{-4}$ | 1,586 (14) | 1,015 (14) | $3.88\cdot 10^{5}$ | $4.65\cdot 10^{5}$ | 6.53 | 6.58
$10^{-6}$ | 1,586 (14) | 1,013 (14) | $3.88\cdot 10^{5}$ | $4.65\cdot 10^{5}$ | 6.75 | 6.61
$8.32\cdot 10^{4}$ | $10^{-2}$ | 1,327 (14) | 887 (14) | $2.12\cdot 10^{6}$ | $2.21\cdot 10^{6}$ | 22.67 | 39.60
$10^{-4}$ | 1,759 (15) | 1,054 (15) | $2.12\cdot 10^{6}$ | $2.21\cdot 10^{6}$ | 27.58 | 44.14
$10^{-6}$ | 1,575 (14) | 934 (14) | $2.12\cdot 10^{6}$ | $2.21\cdot 10^{6}$ | 24.83 | 40.23
$3.30\cdot 10^{5}$ | $10^{-2}$ | 246 (8) | 203 (8) | $1.06\cdot 10^{7}$ | $1.09\cdot 10^{7}$ | 27.52 | 75.89
$10^{-4}$ | 246 (8) | 204 (8) | $1.06\cdot 10^{7}$ | $1.09\cdot 10^{7}$ | 27.20 | 76.73
$10^{-6}$ | 246 (8) | 204 (8) | $1.06\cdot 10^{7}$ | $1.09\cdot 10^{7}$ | 27.38 | 77.15
$1.32\cdot 10^{6}$ | $10^{-2}$ | 193 (7) | 158 (7) | $5.51\cdot 10^{7}$ | $5.38\cdot 10^{7}$ | 99.14 | 308.79
$10^{-4}$ | 193 (7) | 158 (7) | $5.51\cdot 10^{7}$ | $5.38\cdot 10^{7}$ | 101.78 | 318.35
$10^{-6}$ | 193 (7) | 158 (7) | $5.51\cdot 10^{7}$ | $5.38\cdot 10^{7}$ | 99.71 | 318.01

We can draw several observations from the results in Table 5. Firstly, one can
observe that in this case, a diagonal approximation of $H$ is sufficiently
good to deliver very fast convergence of MINRES. The block-diagonal
preconditioner using non-diagonal Hessian information (i.e.
$\widehat{P}^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$) required
consistently fewer MINRES iterations (and not necessarily more memory; see the
three largest experiments), however, this did not result in a reduction in CPU
time. There are several reasons for this. Firstly, the Hessian of the problem
becomes “no less” diagonally dominant as the problem size is increased. As a
result, the diagonal approximation of it remains robust with respect to the
problem size for the problem under consideration. On the other hand, the
algorithm uses the built-in MATLAB function ldl to factorize the
preconditioner $\widehat{P}^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$.
While this implementation is very stable, it employs a dynamic permutation at
each IP-PMM iteration, which slows down the algorithm. In this case, a
specialized method using preconditioner
$\widehat{P}^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$ should employ a
separate symbolic factorization step, that could be used in subsequent IP-PMM
iterations, thus significantly reducing the CPU time. This is not done here,
however, as we treat these PDE-optimization problems as black-box (notice that
the implementation allows the user to feed an approximation of the Hessian,
but does not allow the user to use a different $LDL^{\top}$ decomposition). In
all the previous runs, the reported Krylov iterations include both the
predictor and the corrector steps of IP-PMM. Thus, the systems solved in each
case are twice the number of IP iterations. We should note that for the
problem under consideration employing a predictor–corrector scheme is not
necessary, however, we wanted to keep the implementation as general and robust
as possible, without tailoring it to specific applications. For this problem,
we can also observe that IP-PMM was robust with respect to the problem size
(i.e. IP-PMM convergence was not significantly affected by the size of the
problem). This is often observed when employing an IPM for the solution of PDE
optimization problems (e.g. see [36]), however, in theory one should expect
dependence of IPM on the problem size.

Next we consider the optimal control of the convection–diffusion equation,
i.e. $-\epsilon\rm{\Delta y}+\rm{w}\cdot\nabla y=u$, on the domain
$\Omega=(0,1)^{2}$, where $\rm{w}$ is the wind vector given by
$\rm{w}=[2x_{2}(1-x_{1})^{2},-2x_{1}(1-x_{2}^{2})]^{\top}$, with control
bounds $\rm{u_{a}}=-2$, $\rm{u_{b}}=1.5$ and free state (e.g. see [38, Section
5.2]). Once again, the problem is discretized using Q1 finite elements,
employing the Streamline Upwind Petrov-Galerkin (SUPG) upwinding scheme
implemented in [11]. We define the desired state as
$\rm{\bar{y}}=exp(-64((x_{1}-0.5)^{2}+(x_{2}-0.5)^{2}))$ with zero boundary
conditions. The diffusion coefficient $\epsilon$ is set as $\epsilon=0.02$.
The $\text{L}^{2}$ regularization parameter $\alpha_{2}$ is set as
$\alpha_{2}=10^{-2}$. We run IP-PMM with the two different preconditioning
approaches on the aforementioned problem, with different $\text{L}^{1}$
regularization values (i.e. $\alpha_{1}$) and with increasing problem size.
The results are collected in Table 6.

Table 6: Comparison of QP preconditioners (Convection–Diffusion Control:
problem size and varying regularization).

$\bm{n}$ | $\bm{\alpha_{1}}$ | Krylov (IP) Its. | max nnz | time (s)
---|---|---|---|---
${P}^{C}_{AS,\left(k_{c},k_{r}\right)}$ | $\widehat{P}^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$ | ${P}^{C}_{AS,\left(k_{c},k_{r}\right)}$ | $\widehat{P}^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$ | ${P}^{C}_{AS,\left(k_{c},k_{r}\right)}$ | $\widehat{P}^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$
$2.11\cdot 10^{4}$ | $10^{-2}$ | 3,947 (21) | 1,903 (19) | $3.88\cdot 10^{5}$ | $4.65\cdot 10^{5}$ | 15.94 | 11.67
$10^{-4}$ | 7,546 (25) | 2,721 (22) | $3.88\cdot 10^{5}$ | $4.65\cdot 10^{5}$ | 31.24 | 15.31
$10^{-6}$ | 7,489 (25) | 2,962 (23) | $3.88\cdot 10^{5}$ | $4.65\cdot 10^{5}$ | 31.29 | 16.42
$8.32\cdot 10^{4}$ | $10^{-2}$ | 3,464 (19) | 1,937 (19) | $2.12\cdot 10^{6}$ | $2.21\cdot 10^{6}$ | 49.05 | 65.98
$10^{-4}$ | 7,198 (25) | 2,976 (23) | $2.12\cdot 10^{6}$ | $2.21\cdot 10^{6}$ | 99.49 | 93.19
$10^{-6}$ | 7,150 (25) | 3,178 (24) | $2.12\cdot 10^{6}$ | $2.21\cdot 10^{6}$ | 98.82 | 98.90
$3.30\cdot 10^{5}$ | $10^{-2}$ | 4,037 (21) | 1,667 (18) | $1.07\cdot 10^{7}$ | $1.09\cdot 10^{7}$ | 297.96 | 285.77
$10^{-4}$ | 4,971 (22) | 2,542 (22) | $1.07\cdot 10^{7}$ | $1.09\cdot 10^{7}$ | 357.43 | 334.89
$10^{-6}$ | 5,418 (23) | 2,530 (22) | $1.07\cdot 10^{7}$ | $1.09\cdot 10^{7}$ | 372.40 | 354.23
$1.32\cdot 10^{6}$ | $10^{-2}$ | 384 (9) | 284 (9) | $5.51\cdot 10^{7}$ | $5.38\cdot 10^{7}$ | 158.17 | 436.49
$10^{-4}$ | 385 (9) | 286 (9) | $5.51\cdot 10^{7}$ | $5.38\cdot 10^{7}$ | 161.12 | 438.44
$10^{-6}$ | 385 (9) | 286 (9) | $5.51\cdot 10^{7}$ | $5.38\cdot 10^{7}$ | 165.57 | 443.67

In Table 6 we can observe that the IP-PMM convergence is improved when the
problem size is increased, which relates to the good conditioning of the
Hessian. On the other hand, IP-PMM convergence is affected by the
$\text{L}^{1}$ regularization parameter $\alpha_{1}$. Unlike in the Poisson
control problem, we can see clear advantages of using
$\widehat{P}^{L}_{AS,\left(\lvert\mathcal{N}\rvert,0\right)}$ instead of
${P}^{C}_{AS,\left(k_{c},k_{r}\right)}$ in this case. We can observe that in
this problem using non-diagonal Hessian information within the preconditioner
is significantly more important, and the reduced number of Krylov iterations
often translates into a reduction of the CPU time. As before, we should
mention that the reported number of Krylov iterations includes the solution of
both the predictor and the corrector steps for each IP-PMM iteration.

Overall, we observe that each of the presented approaches can be very
successful on a wide range of problems, including those of very large scale.
Although we have treated every problem as if we knew nothing about its
structure for these numerical tests, our a priori knowledge of the
preconditioners and of the problem’s structure could in principle aid us in
selecting a preconditioner without compromising their “general purpose”
nature.

##  5 Conclusions

In this paper we have presented several general-purpose preconditioning
methodologies suitable for primal-dual regularized interior point methods,
applied to convex optimization problems. All presented preconditioners are
positive definite and hence can be used within symmetric solvers such as PCG
or MINRES. After analyzing and discussing the different preconditioning
approaches, we have presented extensive numerical results, showcasing their
use and potential benefits for different types of practical applications of
convex optimization. A robust and general IP-PMM implementation, using the
proposed preconditioners, has been provided for the solution of convex
quadratic programming problems, and one can readily observe its ability of
reliably and efficiently solving general large-scale problems, with minimal
input from the user.

As a future research direction, we would like to include certain matrix-free
preconditioning methodologies that could be used as alternatives for huge-
scale instances that cannot be solved by means of factorization-based
preconditioners, due to memory requirements.

## Acknowledgments

JG and SP were supported by the Google project “Fast (1 + $x$)-order Methods
for Linear Programming”. JWP gratefully acknowledges support from the
Engineering and Physical Sciences Research Council (UK) grant EP/S027785/1. We
wish to remark that this study does not have any conflict of interest to
disclose.

##  Appendix A Proof of Theorem 3.1

For simplicity of notation, let $P_{AS,(k_{c},k_{r})}\equiv P_{AS}$. In order
to provide an outline of the proof of Theorem 3.1, which follows trivially by
extending the result in [7, Theorem 3], we need to introduce the notion of the
Rayleigh quotient for symmetric matrices. The numerical range of a symmetric
matrix $U$, denoted as $q(U)$, is defined as

$q(U)\coloneqq\left\\{z\in\mathbb{R},\textnormal{ s.t.}\
z=\frac{x^{\top}Ux}{x^{\top}x},\textnormal{ for some }x\in\mathbb{R}^{n},\
x\neq 0\right\\}.$

Given the notation of Section 3.1, an element of the Rayleigh quotient of
these matrices is represented as:

$\begin{array}[]{lcllcl}\gamma_{NE}&\in&q\left(\widetilde{M}_{NE}\right)=[\alpha_{NE},\beta_{NE}],{}{}&~{}~{}\gamma_{F}&\in&q\left(\widetilde{F}\right)=[\alpha_{F},\beta_{F}].\end{array}$

Similarly, an arbitrary element of $q(P_{NE})$ is denoted by

$\begin{array}[]{lcllclcl}\gamma_{p}&\in&[\lambda_{\min}(P_{NE}),\lambda_{\max}(P_{NE})]&\subseteq&\left[\delta,\dfrac{\sigma_{\max}^{2}(A)}{\rho}+\delta\right].\end{array}$

The eigenvalues of $P_{AS}^{-1}K$ are the same as those of

$P_{AS}^{-1/2}KP_{AS}^{-1/2}=\begin{bmatrix}\widehat{F}^{-1/2}&0_{n,m}\\\
0_{m,n}&P_{NE}^{-1/2}\end{bmatrix}\begin{bmatrix}-F&A^{\top}\\\ A&\delta
I_{m}\end{bmatrix}\begin{bmatrix}\widehat{F}^{-1/2}&0_{n,m}\\\
0_{m,n}&P_{NE}^{-1/2}\end{bmatrix}=\begin{bmatrix}-\widetilde{F}&R^{\top}\\\
R&\delta P_{NE}^{-1}\end{bmatrix},$

where $\widetilde{F}$ is defined in (3.2) and $R\coloneqq
P_{NE}^{-1/2}A\widehat{F}^{-1/2}.$ Any eigenvalue $\lambda$ of the
preconditioned matrix $P_{AS}^{-1/2}KP_{AS}^{-1/2}$ must therefore satisfy

$\displaystyle-\widetilde{F}w_{1}\ +$ $\displaystyle R^{\top}w_{2}$
$\displaystyle=\lambda w_{1},$ (A.1) $\displaystyle Rw_{1}\ +$
$\displaystyle\delta P_{NE}^{-1}w_{2}$ $\displaystyle=\lambda w_{2}.$ (A.2)

First, note that

$\displaystyle
RR^{\top}=P_{NE}^{-1/2}A\widehat{F}^{-1}A^{\top}P_{NE}^{-1/2}=P_{NE}^{-1/2}\left(\widehat{M}-\delta
I_{m}\right)P_{NE}^{-1/2}=\widetilde{M}_{NE}-\delta P_{NE}^{-1}.$

If $\lambda\not\in[-\beta_{F},-\alpha_{F}]$ then $\widetilde{F}+\lambda I_{n}$
is symmetric positive (or negative) definite; moreover $R^{\top}w_{2}\neq
0_{n}$. Then from (A.1) we obtain

$w_{1}=\left(\widetilde{F}+\lambda I_{n}\right)^{-1}R^{\top}w_{2},$

which, after substituting in (A.2), yields

$R\left(\widetilde{F}+\lambda I_{n}\right)^{-1}R^{\top}w_{2}+\delta
P_{NE}^{-1}w_{2}=\lambda w_{2}.$

Premultiplying by $w_{2}^{\top}$ and dividing by $\|w_{2}\|^{2}$, we obtain
the following equation, where we have set $z=R^{\top}w_{2}$:

$\lambda=\frac{z^{\top}\left(\widetilde{F}+\lambda
I_{n}\right)^{-1}z}{z^{\top}z}\frac{w_{2}^{\top}RR^{\top}w_{2}}{w_{2}^{\top}w_{2}}+\delta\frac{w_{2}^{\top}P_{NE}^{-1}w_{2}}{w_{2}^{\top}w_{2}}=\frac{1}{\gamma_{F}+\lambda}\left(\gamma_{NE}-\frac{\delta}{\gamma_{p}}\right)+\frac{\delta}{\gamma_{p}}.$

Hence, $\lambda$ must satisfy the following second-order algebraic equation:

$\lambda^{2}+\left(\gamma_{F}-\omega\right)\lambda-\left(\omega(\gamma_{F}-1)+\gamma_{NE}\right)=0,$

where we have set $\omega=\dfrac{\delta}{\gamma_{p}}$ satisfying $\omega\leq
1$. Notice that $\gamma_{NE}-\omega\geq 0$ by construction.

All bounds, except for the lower bound of $I_{+}$, follow directly by
following the developments in [7, Theorem 3]. Thus, we only derive the lower
bound for the positive eigenvalues of the preconditioned matrix, which can be
obtained by computing a lower bound for the positive eigenvalue solution of
the previous algebraic equation. In particular, we have

$\begin{split}\lambda_{+}=&\frac{1}{2}\left[\omega-\gamma_{F}+\sqrt{(\gamma_{F}-\omega)^{2}+4(\omega\gamma_{F}-\omega+\gamma_{NE})}\right]\\\
=&\frac{1}{2}\left[\omega-\gamma_{F}+\sqrt{(\gamma_{F}+\omega)^{2}+4(\gamma_{NE}-\omega)}\right].\end{split}$

We notice that $\lambda_{+}$ is a decreasing function with respect to the
variable $\gamma_{F}$ and increasing with respect to $\gamma_{NE}$. Hence, we
have that:

$\begin{split}\lambda_{+}\geq\
&\frac{1}{2}\left[\omega-\beta_{F}+\sqrt{(\beta_{F}+\omega)^{2}+4(\alpha_{NE}-\omega)}\right]\\\
\geq\
&\frac{1}{2}\left[-\beta_{F}+\sqrt{\beta_{F}^{2}+4\alpha_{NE}}\right],\end{split}$

where the last inequality follows because the penultimate expression is
increasing with respect to $\omega$.

## References

  * [1] A. Altman and J. Gondzio, Regularized Symmetric Indefinite Systems in Interior Point Methods for Linear and Quadratic Optimization, Optimization Methods and Software, 11 (1999), pp. 275–302, https://doi.org/10.1080/10556789908805754.
  * [2] E. D. Andersen, J. Gondzio, C. Mészáros, and X. Xu, Implementation of Interior Point Methods for Large Scale Linear Programming, in Interior Point Methods in Mathematical Programming, T. Terlaky, ed., Kluwer Academic Publishers, 1996, pp. 189–252, https://doi.org/10.1007/978-1-4613-3449-1_6.
  * [3] P. Armand and R. Omheni, A Mixed Logarithmic Barrier-Augmented Lagrangian Method for Nonlinear Optimization, Journal of Optimization Theory and Applications, 173 (2017), pp. 523–547, https://doi.org/10.1007/s10957-017-1071-x.
  * [4] S. Bellavia, V. De Simone, D. di Serafino, and B. Morini, On the Update of Constraint Preconditioners for Regularized KKT Systems, Computational Optimization and Applications, 65 (2016), pp. 339–360, https://doi.org/10.1007/s10589-016-9830-4.
  * [5] M. Benzi, Preconditioning Techniques for Large Linear Systems: A Survey, Journal of Computational Physics, 182 (2022), pp. 418–477, https://doi.org/10.1006/jcph.2002.7176.
  * [6] M. Benzi, G. H. Golub, and J. Liesen, Numerical Solution of Saddle Point Systems, Acta Numerica, 14 (2005), pp. 1–137, https://doi.org/10.1017/S0962492904000212.
  * [7] L. Bergamaschi, J. Gondzio, A. Martínez, J. W. Pearson, and S. Pougkakiotis, A New Preconditioning Approach for an Interior Point-Proximal Method of Multipliers for Linear and Convex Quadratic Programming, Numerical Linear Algebra with Applications, 28 (2020), p. e2361, https://doi.org/10.1002/nla.2361.
  * [8] L. Bergamaschi, J. Gondzio, M. Venturin, and G. Zilli, Inexact Constraint Preconditioners for Linear Systems Arising in Interior Point Methods, Computational Optimization and Applications, 36 (2007), pp. 137–147, https://doi.org/10.1007/s10589-006-9001-0.  see also Erratum, Computational Optimization and Applications 49, 401–406 (2011).
  * [9] L. Bergamaschi, J. Gondzio, and G. Zilli, Preconditioning Indefinite Systems in Interior Point Methods for Optimization, Computational Optimization and Applications, 28 (2004), pp. 149–171, https://doi.org/10.1023/B:COAP.0000026882.34332.1b.
  * [10] S. Bocanegra, F. Campos, and A. Oliveira, Using a Hybrid Preconditioner for Solving Large-Scale Linear Systems arising from Interior Point Methods, Computational Optimization and Applications, 36 (2007), pp. 149–164, https://doi.org/10.1007/s10589-006-9009-5.
  * [11] A. N. Brooks and T. J. R. Hughes, Streamline Upwind/Petrov-Galerkin Formulations for Convection Dominated Flows with Particular Emphasis on the Incompressible Navier-Stokes Equations, Computer Methods in Applied Mechanics and Engineering, 32 (1982), pp. 199–259, https://doi.org/10.1016/0045-7825(82)90071-8.
  * [12] S. Cafieri, M. D’Appuzo, V. De Simone, and D. di Serafino, Stopping Criteria for Inner Iterations in Inexact Potential Reduction Methods: A Computational Study, Computational Optimization and Applications, 36 (2007), pp. 165–193, https://doi.org/10.1007/s10589-006-9007-7.
  * [13] M. D’Apuzzo, 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, 45 (2010), pp. 283–310, https://doi.org/10.1007/s10589-008-9226-1.
  * [14] V. De Simone, D. di Serafino, J. Gondzio, S. Pougkakiotis, and M. Viola, Sparse Approximations with Interior Point Methods, arXiv:2102.13608, (2021).  Accepted for publication in SIAM Review.
  * [15] D. di Serafino and D. Orban, Constraint-Preconditioned Krylov Solvers for Regularized Saddle-Point Systems, SIAM Journal on Scientific Computing, 43 (2021), pp. A1001–A1026, https://doi.org/10.1137/19M1291753.
  * [16] H. S. Dollar, Constraint-Style Preconditioners for Regularized Saddle Point Problems, SIAM Journal on Matrix Analysis and Applications, 29 (2007), pp. 672–684, https://doi.org/10.1137/050626168.
  * [17] H. S. Dollar, N. I. M. Gould, W. H. A. Schilders, and A. J. Wathen, Using Constraint Preconditioners with Regularized Saddle-Point Problems, Computational Optimization and Applications, 36 (2007), pp. 249–270, https://doi.org/10.1007/s10589-006-9004-x.
  * [18] H. C. Elman, A. Ramage, and D. J. Silvester, Algorithm 866: IFISS, a Matlab Toolbox for Modelling Incompressible Flow, ACM Transactions on Mathematical Software, 33 (2007), p. 14, https://doi.org/10.1145/1236463.1236469.
  * [19] H. C. Elman, A. Ramage, and D. J. Silvester, IFISS: A Computational Laboratory for Investigating Incompressible Flow Problems, SIAM Review, 52 (2014), pp. 261–273, https://doi.org/10.1137/120891393.
  * [20] M. P. Friedlander and D. Orban, A Primal-Dual Regularized Interior-Point Method for Convex Quadratic Programs, Mathematical Programming Computation, 4 (2012), pp. 71–107, https://doi.org/10.1007/s12532-012-0035-2.
  * [21] P. E. Gill, W. Murray, D. B. Ponceleón, and M. A. Saunders, Preconditioners for Indefinite Systems Arising in Optimization, SIAM Journal on Matrix Analysis and Applications, 13 (1992), pp. 292–311, https://doi.org/10.1137/0613022.
  * [22] A. Greenbaum, Iterative Methods for Solving Linear Systems, Frontiers in Applied Mathematics. SIAM, Philadelphia, USA, 1997, https://doi.org/10.1137/1.9781611970937.
  * [23] A. Greenbaum, V. Pták, and Z. Strakoš, Any Nonincreasing Convergence Curve is Possible for GMRES, SIAM Journal on Matrix Analysis and Applications, 17 (1996), pp. 465–469, https://doi.org/10.1137/S0895479894275030.
  * [24] M. R. Hestenes and E. Stiefel, Method of Conjugate Gradients for Solving Linear Systems, Journal of Research of the National Bureau of Standards, 49 (1952), pp. 409–436, https://doi.org/10.6028/jres.049.044.
  * [25] I. C. F. Ipsen, A Note on Preconditioning Non-Symmetric Matrices, SIAM Journal on Scientific Computing, 23 (2001), pp. 1050–1051, https://doi.org/10.1137/S1064827500377435.
  * [26] C. Keller, N. I. M. Gould, and A. J. Wathen, Constraint Preconditioning for Symmetric Indefinite Linear Systems, SIAM Journal on Matrix Analysis and Applications, 21 (2000), pp. 1300–1317, https://doi.org/10.1137/S0895479899351805.
  * [27] X. Li, D. Sun, and K. C. Toh, An Asymptotically Superilinearly Convergent Semismooth Newton Augmented Lagrangian Method for Linear Programming, SIAM Journal on Optimization, 30 (2020), pp. 2410–2440, https://doi.org/10.1137/19M1251795.
  * [28] I. Maros and C. Mészáros, The Role of the Augmented System in Interior Point Methods, European Journal of Operational Research, 107 (1998), pp. 720–736, https://doi.org/10.1016/S0377-2217(97)00074-X.
  * [29] I. Maros and C. Mészáros, A Repository of Convex Quadratic Programming Problems, Optimization Methods and Software, 11 (1999), pp. 671–681, https://doi.org/10.1080/10556789908805768.
  * [30] M. F. Murphy, G. H. Golub, and A. J. Wathen, A Note on Preconditioning for Indefinite Linear Systems, SIAM Journal on Scientific Computing, 21 (2000), pp. 1969–1972, https://doi.org/10.1137/S1064827599355153.
  * [31] Netlib.  http://netlib.org/lp, 2011.
  * [32] Y. Notay, A New Analysis of Block Preconditioners for Saddle Point Systems, SIAM Journal on Matrix Analysis and Applications, 35 (2014), pp. 143–173, https://doi.org/10.1137/130911962.
  * [33] A. R. L. Oliveira and D. C. Sorensen, A New Class of Preconditioners for Large-Scale Linear Systems from Interior Point Methods for Linear Programming, Linear Algebra and its Applications, 394 (2005), pp. 1–24, https://doi.org/10.1016/j.laa.2004.08.019.
  * [34] D. Orban, Limited-Memory LD$L^{T}$ Factorization of Symmetric Quasi-Definite Matrices with Application to Constrained Optimization, Numerical Algorithms, 70 (2015), pp. 9–41, https://doi.org/10.1007/s11075-014-9933-x.
  * [35] C. C. Paige and M. A. Saunders, Solution of Sparse Indefinite Systems of Linear Equations, SIAM Journal on Numerical Analysis, 12 (1975), pp. 617–629, https://doi.org/10.1137/0712047.
  * [36] J. W. Pearson and J. Gondzio, Fast Interior Point Solution of Quadratic Programming Problems Arising from PDE-Constrained Optimization, Numerische Mathematik, 137 (2017), pp. 959–999, https://doi.org/10.1007/s00211-017-0892-8.
  * [37] J. W. Pearson and J. Pestana, Preconditioners for Krylov Subspace Methods: An Overview, GAMM-Mitteilungen, 43 (2020), p. e202000015, https://doi.org/10.1002/gamm.202000015.
  * [38] J. W. Pearson, M. Porcelli, and M. Stoll, Interior-Point Methods and Preconditioning for PDE-Constrained Optimization Problems Involving Sparsity Terms, Numerical Linear Algebra with Applications, 27 (2019), p. e2276, https://doi.org/10.1002/nla.2276.
  * [39] S. Pougkakiotis and J. Gondzio, Dynamic Non-Diagonal Regularization in Interior Point Methods for Linear and Convex Quadratic Programming, Journal of Optimization Theory and Applications, 181 (2019), pp. 905–945, https://doi.org/10.1007/s10957-019-01491-1.
  * [40] S. Pougkakiotis and J. Gondzio, An Interior Point-Proximal Method of Multipliers for Convex Quadratic Programming, Computational Optimization and Applications, 78 (2021), pp. 307–351, https://doi.org/10.1007/s10589-020-00240-9.
  * [41] S. Pougkakiotis and J. Gondzio, An Interior Point-Proximal Method of Multipliers for Linear Positive Semi-Definite Programming, Journal of Optimization Theory and Applications, 192 (2022), pp. 97–129, https://doi.org/10.1007/s10957-021-01954-4.
  * [42] Y. Saad and M. H. Schultz, GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems, SIAM Journal on Scientific Computing, 7 (1986), pp. 856–869, https://doi.org/10.1137/09070583.
  * [43] M. Saunders and J. A. Tomlin, Solving Regularized Linear Programs Using Barrier Methods and KKT Systems, Tech. Report SOL 96-4, Systems Optimization Laboratory, Department of Operations Research, Stanford University, Stanford, CA 94305, USA, December 1996.
  * [44] J. Scott and M. Tůma, A Schur Complement Approach to Preconditioning Sparse Linear Least-Squares Problems with some Dense Rows, Numerical Algorithms, 79 (2018), pp. 1147–1168, https://doi.org/10.1007/s11075-018-0478-2.
  * [45] D. J. Silvester and A. J. Wathen, Fast Iterative Solution of Stabilized Stokes Systems Part II: Using General Block Preconditioners, SIAM Journal on Numerical Analysis, 31 (1994), pp. 1352–1367, https://doi.org/10.1137/0731070.
  * [46] R. J. Vanderbei, Symmetric Quasidefinite Matrices, SIAM Journal on Optimization, 5 (1995), pp. 100–113, https://doi.org/10.1137/0805005.
  * [47] A. J. Wathen, Realistic Eigenvalue Bounds for the Galerkin Mass Matrix, IMA Journal of Numerical Analysis, 7 (1987), pp. 449–457, https://doi.org/10.1093/imanum/7.4.449.
  * [48] A. J. Wathen, Preconditioning, Acta Numerica, 24 (2015), pp. 329–376, https://doi.org/10.1017/S0962492915000021.