
\section{Introduction}

In semidefinite programming the ever increasing number of applications
\cite{AnjosLasserre2012,WolkowiczSaigalVandenberghe2000,BoydVandenberghe2007,HenrionKordaLasserre2020}
leads to a corresponding increase in demand for reliable and efficient
solvers. In general, interior point methods are the method of choice.
Yet, if the order of the matrices gets large and the affine matrix
functions involved do not allow to use decomposition or factorization
approaches such as proposed in
\cite{KimKojimaMevissenYamashita2011,BensonYeZhang2000,BurerMonteiro2003},
general interior point methods are no longer applicable. The limiting
factors are typically memory requirements and computation times
connected with forming and factorizing a Schur complemented system
matrix of the interior point KKT system.  In this case, the spectral
bundle method of \cite{HelmbergRendl2000} offers a viable alternative.

The spectral bundle method reformulates the semidefiniteness condition
via a penalty term on the extremal eigenvalues of a corresponding
affine matrix function and assumes these eigenvalues to be efficiently
computable by iterative methods. In each step it selects a subspace
close to the current active eigenspace. Then the next candidate point
is determined as the proximal point with respect to the extremal
eigenvalues of the affine matrix function projected onto this
subspace. The proximal point is the optimal solution to a quadratic
semidefinite subproblem whose matrix variable is of the order of the
dimension of the approximating subspace. If the subspace is kept
small, this allows to find approximately optimal solutions in
reasonable time. In order to reach solutions of higher precision it
seems unavoidable to go beyond the full active eigenspace
\cite{HelmbergOvertonRendl2014,DingGrimmer2020}. In the current
implementation within the callable library ConicBundle
\cite{ConicBundle2021} the quadratic subproblem is solved by an
interior point approach. Again for each of its KKT systems the
limiting work consists in collecting and factorizing a Schur
complement matrix whose order is typically the square of the dimension
of the active eigenspace. The main question addressed here is whether
it is possible to stretch these limits by developing a suitably
preconditioned iterative solver that allows to circumvent the
collection and factorization of this Schur complement.  The focus is
thus not on the spectral bundle method itself but on solving KKT
systems of related quadratic semidefinite programs by iterative
methods. Even though the methodology will be developed and discussed
for low rank properties that arise in the ConicBundle setting, some of
the considerations and ideas should be transferable to general conic
quadratic optimization problems whose quadratic term consists of
a positive diagonal plus a low rank Gram matrix or maybe even to
general positive definite systems of this form.

Here is an outline of the paper and its main
contributions. Section \ref{S:ConicBundleKKT} provides
the necessary background on the bundle philosophy underlying
ConicBundle and derives the KKT system of the bundle subproblem. The
core of the work is presented in Section \ref{S:lowrankprecond} on low
rank preconditioning for a Gram-matrix plus positive diagonal. For
slightly greater generality, denote the cone of positive
(semi)definite matrices of order $m$ by $\mathbb{S}^m_{++}\ (\mathbb{S}^m_+)$ and
let the system matrix be given in the form
\begin{displaymath}
  H=D+VV^\top\qquad\text{ with }D\in\mathbb{S}^m_{++}, V\in\mathbb{R}^{m\times n}, 
\end{displaymath}
where it is tacitly assumed that $D^{-1}$ times vector and $V$ times
vector are efficiently computable. Typically $n\le m$ but whenever $n$
is sizable one would like to approximate $V$ by a matrix
$\hat V\in\mathbb{R}^{m\times k}$ with significantly smaller $k<n$ to obtain a
preconditioner $\hat H=D+ \hat V\hat V^\top$ with low rank inverse
$\hat H^{-1}=D^{-1}-D^{-1}\hat V(I_k+\hat V^\top D^{-1}\hat
V)^{-1}\hat V^\top D^{-1}$. Comparing this to the inverse of $H$, the
goal is to capture the large eigenvalues of $V^\top D^{-1}V$, more
precisely the directions belonging to large singular values of
$D^{-\frac12}V$. By the singular value decomposition (SVD) this can be
achieved by the projection onto a subspace, say $D^{-\frac12}VP$ for a
suitably chosen orthogonal $P\in\mathbb{R}^{n\times k}$. Because the full SVD
is computationally too expensive, two other approaches will be
developed and analyzed here. In the first, Section
\ref{S:randprecond}, the orthogonal $P$ is generated by a Gaussian
matrix $\Omega\in\mathbb{R}^{n\times k}$. In the second, Section
\ref{S:detprecond}, some knowledge about the interior point method
leading to $V$ will be exploited in order to form $P$
deterministically.

The projection onto a random subspace may be motivated geometrically
by interpreting the Gram matrix $VV^\top$ as the inner products of the
row vectors of $V$. The result of Johnson-Lindenstrauss, cf.~\cite{Achlioptas2001,DasguptaGupta2002}, allows to approximate this
with low distortion by a projection onto a low dimensional
subspace. In matrix approximations this idea seems to have first
appeared in \cite{PapadimitriouRaghavanTamakiVempala2000}.  In
connection with preconditioning a recent probabilistic approach is
described in \cite{HighamMary2019} in the context of controlling the
error of a LU preconditioner. \cite{HalkoMartinssonTropp2011} gives an
excellent introduction to probabilistic algorithms for constructing
approximate matrix decompositions and provides useful
bounds. Building directly on their techniques we provide
deterministic and probabilistic bounds on the condition number of the
random subspace preconditioned system in theorems
\ref{th:deterministicbound} and \ref{th:probabilisticbound}. In
comparison to the moment analysis of the Ritz values of the
preconditioned matrix presented in Theorem~\ref{T:moments}, the bounds seem
to fall below expectation and are maybe still improvable.  Random
projections do not require any problem specific structural insights,
but it remains open how to choose the subspace dimension in order to
obtain an efficient preconditioner.

In contrast, identifying the correct subspace seems to work well for the
deterministic preconditioning routine. It exploits structural
properties of the KKT system's origin in interior point methods.
Within interior point methods iterative approaches have been
investigated in quite a number of works, in conjunction with
semidefinite optimization see \emph{e.g.}\  \cite{Toh2008,KocvaraStingl2007}.
These methods were
mostly designed for exploiting sparsity rather than low rank
structure.  During the last months of this work an approach closely
related to ours appeared in \cite{HabibiKavandKocvaraStingl2021}. It significantly
extends ideas of \cite{ZhangLavaei2017} for a deterministic
preconditioning variant. It assumes the rank of the optimal solution to
be known in advance and provides a detailed analysis for this case. Their
ideas and arguments heavily influenced the condition number analysis
of our approach presented in theorems \ref{T:deterministiccond} and
\ref{T:speccondbound}. In contrast to \cite{HabibiKavandKocvaraStingl2021}, our
algorithmic approach does not require any a priori knowledge on the
rank of the optimal solution. Rather, Theorem~\ref{T:speccondbound}
and Lemma~\ref{L:TTTPSC} motivate an
estimate on the singular value induced by certain directions
associated with active interior point variables, that seems to offer a good indicator for the relevance of the corresponding subspace.

In Section \ref{S:experiments} the performance of the preconditioning
approaches is illustrated relative to the direct solver on sequences
of KKT systems that arise in solving three large
scale instances within ConicBundle.  The deterministic approach turns out to
be surprisingly effective in identifying a suitable subspace.  It
provides good control on the condition number and reduces the number
of matrix vector multiplications significantly.  The selected instances are
also intended to demonstrate the differences in the potential of the
methods depending on the problem characteristics.  Roughly, the direct
solver is mainly attractive if the model is tiny, if significant parts of the Schur
complement can be precomputed for all KKT systems of the same
subproblem or if precision requirements get exceedingly high with the entire
bundle model being strongly active. In general, however, the iterative
approach with deteriministic preconditioner can be expected to lead to
significant savings in computation time in large scale applications.
In Section \ref{S:concl} the paper ends with some concluding remarks.

\textbf{Notation.} For matrices or vectors $A,B\in\mathbb{R}^{m\times n}$ the
(trace) inner product is denoted by $\ip{A}{B}=\mathop{\rm tr} B^\top A=\sum_{ij}A_{ij}B_{ij}$. $A\circ B=(A_{ij}B_{ij})$ denotes the elementwise or Hadamard product. 
Consider symmetric matrices $A,B\in\mathbb{S}^n$ of order $n$. For representing
these as vectors, 
the operator $\mathop{\rm svec}
A=(A_{11},\sqrt{2}A_{21},\dots,\sqrt{2}A_{n1},A_{22
},\sqrt{2}A_{32},\dots,A_{nn})^\top$ stacks the columns of the lower
triangle with offdiagonal elements multiplied by $\sqrt{2}$ so that $\ip{A}{B}=\mathop{\rm svec}(A)^\top\mathop{\rm svec}(B)$. For matrices $F,G\in \mathbb{R}^{k\times n}$ the symmetric Kronecker product $\otimes\unkern_{\unkern s}$ is defined by $(F\otimes\unkern_{\unkern s} G)\mathop{\rm svec}(A)=\frac12\mathop{\rm svec}(FAG^\top+GAF^\top)$.  The Loewner partial
order $A\succeq B$ ($A\succ B$) refers to $A-B\in\mathbb{S}^n_{+}$ ($A-B\in\mathbb{S}^n_{++})$ being positive semidefinite (positive definite). The eigenvalues of $A$ are denoted by $\lambda_{\max}(A)=\lambda_1(A)\ge \dots\ge \lambda_n(A)=\lambda_{\min}(A)$. The norm $\|\cdot\|$ refers to the
Euclidean norm for vectors and to the spectral norm for matrices. 
$I_n$ ($I$) denotes the identity matrix of order $n$ (or of appropriate size), the canonical unit vectors $e_i$ refer to
the $i$-th column of $I$. Unless stated explicitly otherwise, $\mathbbm{1}$ denotes the vector of all ones of appropriate size.
$\mathbb{E}$ refers to the expected value of a random variable, $\mathop{\mathrm{Var}}$ to its variance and
 ${\mathcal{N}}(\mu,\sigma^2)$ to the normal or Gaussian distribution with mean $\mu$ and standard deviation $\sigma$. 





\section{The KKT system of the ConicBundle Subproblem}\label{S:ConicBundleKKT}

The general setting of bundle methods deals with minimizing a (mostly closed) convex 
function $f\colon\mathbb{R}^m\to\overline{\mathbb{R}}:=\mathbb{R}\cup\{\infty\}$ over a closed convex
ground set $C\subseteq\mathop{\mathrm{dom}} f$ of simple structure like $\mathbb{R}^m$, a
box or a polyhedron,
\begin{displaymath}
\text{minimize}~f(y)~~\text{ subject to }~y\inC.
\end{displaymath}
Typically, $f$ is given by a first order oracle, \mbox{i.\,e.,}\  a routine that returns
for a given $\bar y\inC$ the function value $f(\bar y)$ and an
arbitrary subgradient $g\in\partial f(\bar y)$ from the
subdifferential of $f$ in $\bar y$. Value $f(\bar y)$ and subgradient
$g$ give rise to a
supporting hyperplane to the epigraph of $f$ in $\bar y$. The
algorithm collects these affine minorants in the \emph{bundle} to form a
cutting model of $f$. It will be convenient to arrange the value at zero
and the gradient in a pair $\omega=(\gamma=f(\bar y)-\ip{g}{\bar
  y},g)$ and to denote, for $y\in\mathbb{R}^m$, the minorant's value in $y$ by
$\omega(y):=\gamma+\ip{g}{y}$.

Let $\mathcal{W}_f=\{\omega=(\gamma,g)\in\mathbb{R}^{1+m}\colon \gamma+\ip{g}{y}\le f(y),
y\in \mathbb{R}^m\}$ denote the set of all affine minorants of $f$. For closed
$f$ we have $f(y)=\sup_{\omega\in\mathcal{W}_f}\omega(y)$. Any
compact subset $W\subseteq\mathcal{W}_f$ gives rise to a minorizing
cutting model
of $f$,
\begin{displaymath}
  W(y):=\max_{\omega\in W} \omega(y)\le f(y),\quad y\in\mathbb{R}^m. 
\end{displaymath}
At the beginning of iteration $k=0,1,\dots$ the bundle method's state is described by a
current \emph{stability center} $\hat y_k\in \mathbb{R}^m$, a compact cutting
model $W_k\subseteq \mathcal{W}_f$, and a proximity term, here the square 
of a norm $\|\cdot\|_{H_k}^2:=\ip{\cdot}{H_k\cdot}$ with positive definite $H_k$. The method
determines the next candidate $y_{k+1}\in\mathbb{R}^m$ as minimizer of the
\emph{augmented model} or \emph{bundle subproblem}
\begin{equation}\label{eq:candidate}
  y_{k+1}=\mathop {\rm argmin}_{y\inC} W_k(y)+\tfrac12\|y-\hat y_k\|_{H_{k}}^2.
\end{equation}
Solving this bundle subproblem may be viewed as determining a
saddle point $(y_{k+1},\bar\omega_{k+1}=(\bar\gamma_{k+1},\bar g_{k+1}))\in C\times\hbox{\rm conv} W_k$,
which exists for any closed convex $C$ by \cite{Rockafellar70},
Theorems 37.3 and 37.6, due to the strong convexity in $y$ and the compactness of
$W_k$, 
\begin{eqnarray*}
\bar\omega_{k+1}(y_{k+1})+\tfrac12\|y_{k+1}-\hat
  y_k\|_{H_{k}}^2&=&
  \inf_{y\inC}\sup_{\omega=(\gamma,g)\in
    W_k}\gamma+\ip{g}{y}+\tfrac12\|y-\hat
  y_k\|_{H_{k}}^2\\
&=&\sup_{\omega\in\hbox{\rm conv} W_k}\inf_{y\inC}\omega(y)+\tfrac12\|y-\hat
  y_k\|_{H_{k}}^2.
\end{eqnarray*}
Strong convexity in $y$ ensures uniqueness of $y_{k+1}$. First
order optimality with respect to $y$ implies
\begin{displaymath}
 \bar g_{k+1}+H_k(y_{k+1}-\hat y_k)\in -N_C(y_{k+1}),
\end{displaymath}
where $N_C(y)$ denotes the normal cone to $C$ at $y\in C$.
In the unconstrained case of $C=\mathbb{R}^m$ the \emph{aggregate}
$\bar\omega_{k+1}$ is also unique. Whether unique or not, the \emph{aggregate}
will refer to the solution $\bar\omega_{k+1}\in\hbox{\rm conv} W_k$  produced by the algorithmic approach for
solving \eqref{eq:candidate}. The progress predicted by the model is 
$f(\hat y_k)-\bar\omega_{k+1}(y_{k+1})=f(\hat
y_k)-W_{k}(y_{k+1})$. Actual
progress will be compared to a threshold value which arises from
damping the progress predicted by the model by some  $\kappa
\in(0,1)$ 
\begin{displaymath}
  \vartheta_{k+1}=\kappa[f(\hat y_k)-\bar\omega_{k+1}(y_{k+1})].
\end{displaymath}
Next, $f$ is evaluated at $y_{k+1}$ by calling the oracle which returns $f(y_{k+1})$ and a
new minorant $\omega_{k+1}$ with
$\omega_{k+1}(y_{k+1})=f(y_{k+1})$. If progress in objective value
is sufficiently large in comparison to the
progress predicted by the model, \mbox{i.\,e.,}\  
\begin{displaymath}
f(\hat
y_k)-f(y_{k+1})\ge \vartheta_{k+1},
\end{displaymath}
the method executes a \emph{descent step} which moves the center of
stability to the new point, $\hat y_{k+1}=y_{k+1}$. Otherwise, in a
\emph{null step}, the center remains unchanged,
$\hat y_{k+1}=\hat y_{k}$, but the new minorant $\omega_{k+1}$ is used
to improve the model. In fact, the requirement
$\{\bar\omega_{k+1},\omega_{k+1}\}\subseteq W_{k+1}$ ensures
convergence of the function values $f(\hat y_k)$ to $f^*=\inf_{y\in C} f(y)$ under mild technical conditions on $H_{k+1}$. For these it suffices, \emph{e.g.}\  to
fix $0<\underline{\lambda}\le\bar\lambda$ and to choose
$\underline{\lambda}I\preceq H_{k+1} \preceq \bar\lambda I$ following
a descent step and $H_k\preceq H_{k+1} \preceq \bar\lambda I$
following a null step, see \cite{BonnansGilbertLemarechalSagastizabal2006}. 

The decisive elements for an efficient implementation are the following:
\begin{itemize}
\item the choice of the cutting model $W_k$,
\item the choice of the proximal term, in our case of $H_k$,
\item the solution method for the bundle subproblem
  \eqref{eq:candidate} with corresponding structural requirements on supported ground
  sets $C$. 
\end{itemize}
and their interplay. While most bundle implementations employ
polyhedral cutting models combined with a suitable active set QP
approach, the ConicBundle callable library \cite{ConicBundle2021} is
primarily designed for (nonnegative combinations of) conic cutting
models built from symmetric cones. In particular the cone of positive
semidefinite matrices and the second order cone lead to
nonpolyhedral models that change significantly in each step. In
solving \eqref{eq:candidate} for these models, interior point methods
are currently the best option available, but with these methods the
cost of assembling the coefficients and solving the subproblem
dominates the work per iteration for most applications. This paper
explores possibilities to replace the classical Schur complement
approach for computing the Newton step by an iterative approach
in order to improve applicability to large scale problems.


In order to describe the main structure of the primal dual KKT-system
for \eqref{eq:candidate}, let us briefly sketch the
conic cutting models employed. These build on combinations of the
cone of nonnegative vectors, the second order cone and the cone
of positive semidefinite matrices, each with a specific trace
vector for measuring the ``size'' of its elements, 
\begin{align*}
  x\in\mathbb{R}^n_+&:=\{x\in\mathbb{R}^n\colon x_i\ge 0, i=1,\dots,n\} &&\Leftrightarrow\colon x\ge 0,&&\text{trace vector }\mathbbm{1}_n=\begin{sbmatrix}1\\[-1ex]\vdots\\1\end{sbmatrix},\\  
x\in {\ensuremath{\calQ^n}}&:=\left\{\vectwo{x_1}{\bar x}\in\mathbb{R}^n\colon x_1\ge\|\bar x\|\right\}&&
\Leftrightarrow\colon x\unrhd 0,&&\text{trace vector }e_1=\begin{sbmatrix}1\\0\\[-1ex]\vdots\\0\end{sbmatrix},\\
X\in\mathbb{S}^n_+&:=\left\{X\in\mathbb{S}^n\colon \lambda_{\min}(X)\ge 0\right\}&&\Leftrightarrow\colon X\succeq 0, &&\text{trace matrix }I_n.
\end{align*}
Cartesian products of these are described
by a triple $t=(n,q,s)$ with $n\in \mathbb{N}_0$, $q\in\mathbb{N}^{n_q}$, $s\in\mathbb{N}^{n_s}$ for some  $n_q\in\mathbb{N}_0$, $n_s\in\mathbb{N}_0$ specifying the cone
\begin{equation}\label{eq:St}
  \mathcal{S}^t_+:=\mathbb{R}_+^n\times \bigtimes_{i=1}^{n_q} \soc{q_i}\times \bigtimes_{i=1}^{n_s} \mathbb{S}_+^{s_{i}}.
\end{equation}
The cone $\mathcal{S}^t_+$ will be regarded as a full dimensional cone in $\mathbb{R}^t:=\mathbb{R}^{n(t)}$ with $n(t)=n+\mathbbm{1}^\top q+\sum_{i=1}^{n_s}{s_i+1\choose 2}$. Whenever convenient,
an $x\in \mathcal{S}^t_+$ or $\mathbb{R}^t$ will be split into
\begin{equation}\label{eq:xt}
x=(\xi^\top,(x^{(1)})^\top,\dots,(x^{(n_q)})^\top,\mathop{\rm svec}(X^{(1)})^\top,\dots,\mathop{\rm svec}(X^{(n_s)})^\top)^\top
\end{equation}
in the natural way. Indices or elements may be omitted if the corresponding counters $n$, $n_q$, $n_s$ happen to be 0 or 1. 
The \emph{trace} $\mathop{\rm tr}(\cdot)$ of an element $x\in S_t$ is defined to be 
\begin{displaymath}
  \mathop{\rm tr} x := \ip{\mathbbm{1}_t}{x}:= \ip{\mathbbm{1}}{\xi}+\sum_{i=1}^{n_q}\ip{e_1}{x^{(i)}}+\sum_{i=1}^{n_s} \ip{I}{X^{(i)}}.
\end{displaymath}
In ConicBundle each cutting model may be considered to be specified by a tuple  $M=(t,\alpha,K,{\mathcal{B}},\underline\omega)$ via
\begin{displaymath}
  W_M(y)=\maxset{x\in \mathcal{S}^t_+\\ \tau-\ip{\mathbbm{1}_t}{x}\in K}[\underline\omega(y) +\ip{{\mathcal{B}}(y)}{x}]
\end{displaymath}
where
\begin{align*}
  &t=(n,q,s)&&\text{ specifies the cone as above,}\\
  & \tau>0 &&\text{ gives the trace value or trace upper bound,}\\
  & K\in\{\{0\},\mathbb{R}_+\}&&\text{ specifies constant or bounded trace,}\\
  &\begin{array}{@{}r@{\;}r@{\,}c@{\,}l}\mathcal{B}\colon&\mathbb{R}^m&\to&\mathbb{R}^{n(t)}\\
     &y&\mapsto&\mathcal{B}(y)=B_0+By\end{array} &&\text{ represents the bundle as affine function,}\\
  & \underline{\omega}=(\underline{\gamma},\underline{g}) &&\text{ provides a constant offset subgradient}.
\end{align*}
For example, the standard polyhedral model for $h$ subgradients $\omega_i=(\gamma_i,g_i)$, $i=1,\dots,h$ is obtained for $M=(t=(h,0,0),\tau=1,K=\{0\},
(B_0=\begin{sbmatrix}
  \gamma_1\\[-1ex]\vdots\\\gamma_h
\end{sbmatrix},B=
\begin{sbmatrix}
  g_1^\top\\[-1ex]\vdots\\g_h^\top
\end{sbmatrix}),\underline{\omega}=0)$. Indeed, maximizing $\ip{{\mathcal{B}}(y)}{x}$ over $\xi\in\mathbb{R}^h_+$ with $\mathbbm{1}^\top\xi=1$ finds the best convex combination of the subgradients at $y$. In polyhedral models ${\mathcal{B}}$ may well be sparse, but in combination with positive semidefinite models large parts of ${\mathcal{B}}$ will be dense. Therefore we will not assume any specific structure in the bundle ${\mathcal{B}}$. 


Depending on the variable metric heuristic in use, see \cite{HelmbergOvertonRendl2014,HelmbergPichler2017} for typical choices in ConicBundle, the proximal term may either be a positive multiple of $I_m$ or it may be of the form $H=D_H+V_HV_H^\top\in\mathbb{S}^n_+$, where $D_H$ is a diagonal matrix with strictly positive diagonal entries and $V_H\in\mathbb{R}^{m\times h_H}$ specifies a rank $h_H$ contribution. If $V_H$ is present, it typically consists of dense orthogonal columns.

The final ingredient is the basic optimization set $C$, which may have a polyhedral description of the form
\begin{displaymath}
  C=\{y\colon \underline{a} \le A y \le \overline{a}, \underline{y}\le y\le\overline{y}\},
\end{displaymath}
where $A\in\mathbb{R}^{h_A\times m}$, $\underline{a}\in
(\mathbb{R}\cup\{-\infty\})^{h_A},\overline{a}\in(\mathbb{R}\cup\{\infty\})^{h_A}$,
$\underline{y}\in(\mathbb{R}\cup\{-\infty\})^m,\overline{y}\in(\mathbb{R}\cup\{\infty\})^m$
are given data.
If $A$ is employed in applications, we expect the number of rows $h_A$ of $A$ to be small in comparison to $m$. The set $C$ is tested for feasibility in advance, but no preprocessing is applied.

In order to reduce the indexing load in this presentation, we consider problem \eqref{eq:candidate} in $\hat y$ only for a single iteration and one model 
$(t,\alpha,K,\mathcal{B},\underline\omega=(\underline\gamma,\underline
g))$. Putting $b=-H\hat y+\underline{g}$ and $\delta=\frac12\ip{H\hat
  y}{\hat y}+\underline\gamma$ the bundle problem may be written in
the form 
\begin{equation}\label{eq:saddle}
  \minset{Ay-w=0\\\underline{a}\le w\le \overline{a}\\\underline{y}\le y \le \overline{y}}
  \maxset{x\in \mathcal{S}^t_+\\ \alpha-\ip{\mathbbm{1}_t}{x}-\sigma=0\\\sigma \in K}
  \left[\frac12\ip{Hy}{y}+\ip{b}{y}+\ip{B_0}{x}+\ip{By}{x}+\delta\right].
\end{equation}
The existence of saddle points is guaranteed by compactness of the $x$-set and
by strong convexity for $y$ due to $H\succ 0$, see \emph{e.g.}\  \cite{HiriartUrrutyLemarechal93a}. 
For the purpose of presentation, the primal dual KKT system
for solving the saddle point problem \eqref{eq:saddle} is built for $\underline{a}<\overline{a}$, $\underline{y}<\overline{y}$ and
$K=\mathbb{R}_+$. The extension to the equality cases follows quite naturally
and will be commented on at appropriate places.  Throughout we will
assume that $A,A^\top$ and $B,B^\top$ are given by matrix-vector multiplication
oracles. $A$ is assumed to have few rows, $B$ may actually have a
large number of columns, $H$ is a positive diagonal plus low rank, but
no further structural information is assumed to be available.

The original spectral bundle approach of \cite{HelmbergRendl2000} was designed for the
unconstrained case $C=\mathbb{R}^m$ which allows direct elimination of $y$ by
convex optimality. Setting up the maximization problem for $x$ then
requires forming the typically dense Schur complement
$BH^{-1}B^\top$. For increasing bundle sizes this is in fact the
limiting operation within each bundle iteration. The aim of
developing an iterative approach for \eqref{eq:saddle} is therefore
not only to allow for general $C$ but also to circumvent the explicit
computation of this Schur complement.

For setting up a primal-dual interior point approach for solving \eqref{eq:saddle}, the dual variables to constraints on the minimizing side
will be denoted by $s\in\mathbb{R}^{h_A}$, $s_{\underline{a}},s_{\overline{a}}\in\mathbb{R}^{h_A}_+$, $s_{\underline{y}},s_{\overline{y}}\in\mathbb{R}^m_+$, the dual variables to the constraints on the maximizing side will be 
$z\in\mathcal{S}^t_+$, $z_\alpha\in\mathbb{R}$ and $\zeta\in K^*=\mathbb{R}_+$.



With barrier parameter $\mu>0$ for the conic constraints the usual primal-dual KKT-System may be arranged in the form
\begin{equation}\label{eq:origKKT}
  \begin{array}{c}
  \begin{array}{ccccccccccc}
    Hy & + &A^\top s & + &B^\top x & &    &+&s_{\overline{y}} -s_{\underline{y}} &= & -b \\
    Ay &  &        &  &          & &                 &-&w &=&0 \\
    By &   &       &  &          & -&\mathbbm{1}_t\zeta     & +&z & =&-B_0 \\
       & &     &-& \ip{\mathbbm{1}_t}{x} & &                 &-&\sigma&=&-\alpha\\
       &- &s     & & & &                 &+& s_{\overline{a}}-s_{\underline{a}} &=& 0\\
  \end{array}\\\\
  \begin{array}{ccc@{\qquad\qquad}ccc}
    (w-\underline{a})\circ s_{\underline{a}} &=& \mu\mathbbm{1} & (\overline{a}-w)\circ s_{\overline{a}} &=& \mu\mathbbm{1} \\
    (y-\underline{y})\circ s_{\underline{y}} &=& \mu\mathbbm{1} & (\overline{y}-y)\circ s_{\overline{y}} &=& \mu\mathbbm{1} \\
    x\circ_tz &=& \mu \mathbbm{1}_t & \sigma\zeta &=&\mu.
  \end{array}
    \end{array}
\end{equation}
In this, ``$\circ$'' denotes the componentwise Hadamard product and ``$\circ_t$'' a canonical
generalization to the cone $\mathcal{S}^t_+$, employing the arrow operator for second order
cone parts and (typically symmetrized) matrix products for semidefinite parts.

In solving this by Newton's method, the linearization of the first perturbed complementarity
line yields
\begin{displaymath}
  \begin{array}{ccl}
    \Delta s_{\underline a}&=&\mu(w-\underline{a})^{-1}-s_{\underline{a}}-\Delta w\circ s_{\underline{a}}\circ(w-\underline{a})^{-1}\\
    \Delta s_{\overline a}&=&\mu(\overline{a}-w)^{-1}-s_{\overline{a}}+\Delta w\circ s_{\overline{a}}\circ(\overline{a}-w)^{-1}\\ \hline
    \Delta s_{\overline a}-\Delta s_{\underline a}& =&\Delta w\circ\left((s_{\underline{a}}\circ(w-\underline{a})^{-1})+(s_{\overline{a}}\circ(\overline{a}-w)^{-1})\right)\\
                           &&+s_{\underline{a}}-s_{\overline{a}}+\mu\left((\overline{a}-w)^{-1}-(w-\underline{a})^{-1}\right)\\
                           &=:& \Delta w\circ d_w+s_{\underline{a}}-s_{\overline{a}}-\mu c_w.
  \end{array}
\end{displaymath}
Likewise, the second perturbed complementarity line gives rise to
\begin{displaymath}
  \begin{array}{ccl}
    \Delta s_{\overline b}-\Delta s_{\underline b}& =&\Delta y\circ\left((s_{\underline{y}}\circ(y-\underline{y})^{-1})+(s_{\overline{y}}\circ(\overline{y}-y)^{-1})\right)\\
                                                  &&+s_{\underline{y}}-s_{\overline{y}}+\mu\left((\overline{y}-y)^{-1}-(y-\underline{y})^{-1}\right)\\
                                                  &=:& \Delta y\circ d_y+s_{\underline{y}}-s_{\overline{y}}-\mu c_y.
  \end{array}
\end{displaymath}
For dealing with the conic complementarity ``$\circ_t$'' we employ the
symmetrization operators ${\mathcal{E}}_t$ and ${\mathcal{F}}_t$ of \cite{ToddTohTuetuencue98} in diagonal block form corresponding to ${\mathcal{S}}^t_+$, which give rise to a symmetric positive definite 
${\mathfrak{X}}_t={\mathcal{E}}_t^{-1}{\mathcal{F}}_t\succ 0$ with diagonal block structure according to ${\mathcal{S}}^t_{++}$. With this the last perturbed complementarity line results in
\begin{displaymath}
  \begin{array}{ccl}
    \Delta z&=&-{\mathfrak{X}}_t^{-1}\Delta x + \mu x^{-1}-z,\\
    \Delta \sigma & = & -\zeta^{-1}\sigma\Delta \zeta + \mu \zeta^{-1}-\sigma.
  \end{array}
\end{displaymath}
Employing the linearization of the defining equation for $s$, 
\begin{displaymath}
  \Delta s_{\overline a}-\Delta s_{\underline a}=\Delta s + s +s_{\underline{a}}- s_{\overline{a}},
\end{displaymath}
the variable $\Delta w$ may now be eliminated via
\begin{displaymath}
  \Delta w=d_w^{-1}\circ \Delta s+d_w^{-1}\circ s+\mu d_w^{-1}\circ c_w.
\end{displaymath}
Put $D_y:=\mathop{\rm Diag}(d_y)>0$ and $D_w:=\mathop{\rm Diag}(d_w)>0$, then the Newton step is obtained by solving the system
\begin{equation}\label{eq:fullKKT}
  \left[
  \begin{array}{cccc}
    H+D_y&A^\top&B^\top & 0 \\
    A & -D_w^{-1} &  0  & 0 \\
    B & 0 & -{\mathfrak{X}}_t^{-1} & -\mathbbm{1}_t \\
    0 & 0 & -\mathbbm{1}_t^\top & \zeta^{-1}\sigma
  \end{array}\right]
  \left(
  \begin{array}{c}
    \Delta y\\ \Delta s \\ \Delta x \\ \Delta \zeta
  \end{array}\right)
  =
  \left(
  \begin{array}{c}
    r_y\\ r_s \\ r_x \\ r_\zeta
  \end{array}\right),
\end{equation}
with right hand side 
\begin{displaymath}
  \left(
  \begin{array}{c}
    r_y\\ r_s \\ r_x \\ r_\zeta
  \end{array}\right)
  =
  \left(
  \begin{array}{ccl}
    -(Hy+B^\top x+A^\top s+b)&+&\mu c_y\\
    -(Ay-w)+d_w^{-1}\circ s&+&\mu d_w^{-1}\circ c_w \\
    -(By+B_0-\mathbbm{1}_t\zeta)&-&\mu x^{-1}\\
    -(\alpha-\ip{\mathbbm{1}_t}{x})&+&\mu \zeta^{-1}
  \end{array}\right).
\end{displaymath}
The right hand side may be modified as usual to obtain predictor and corrector right hand sides, but this will not be elaborated on here. 
Note,  for $K=\{0\}$ the same system works with $\sigma=0$ and without the centering
terms associated with $\zeta$ in $r_\zeta$. Likewise, whenever line $i$ of $A$ corresponds to an equation, \mbox{i.\,e.,}\  $\underline{a}_i=\overline{a}_i$, the respective entry of $d_w^{-1}$ has to be replaced by zero.

Eq.~\eqref{eq:fullKKT} is a symmetric indefinite system, that could be solved
by appropriate iterative methods directly. So far, however,
we were not able to conceive suitable general preconditioning approaches
for exploiting the given structural properties in the full system.
Surprisingly, a viable path seems to be offered by the traditional
Schur complement approach after all. The resulting system allows to perform matrix vector
multiplications at minimal additional cost and is frequently positive
 definite. 

To see this, first take the Schur complement with respect to the ${\mathfrak{X}}_t^{-1}$ block,
\begin{displaymath}\label{eq:XiZSchur}
  \left[
  \begin{array}{@{}c@{\quad}c@{\quad}c@{\ }}
    H+D_y+B^\top{\mathfrak{X}}_tB &A^\top& -B^\top{\mathfrak{X}}_t\mathbbm{1}_{t} \\
    A &  -D_w^{-1}  & 0 \\
    -(B^\top{\mathfrak{X}}_t\mathbbm{1}_{t})^\top  & 0 & \zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t \\
  \end{array}\right].
\end{displaymath}
Assuming $d_w>0$ (no equality constraint rows in $A$), eliminate the second and third block with further Schur complements and
split $B^\top {\mathfrak{X}}_t B=B^\top{\mathfrak{X}}_t^{\frac12}{\mathfrak{X}}_t^{\frac12}B$,
\begin{equation}\label{eq:SchurH}
  \tilde H= \left[
    H+D_y+B^\top{\mathfrak{X}}_t^{\frac12}(I-\tfrac{{\mathfrak{X}}_t^{\frac12}\mathbbm{1}_{t}({\mathfrak{X}}_t^{\frac12}\mathbbm{1}_{t})^\top}{\zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t}){\mathfrak{X}}_t^{\frac12}B +A^\top D_wA  
  \right].
\end{equation}
Also in the equality case of $\sigma=0$ the matrix $I-\tfrac{{\mathfrak{X}}^{\frac12}\mathbbm{1}_{t}({\mathfrak{X}}^{\frac12}\mathbbm{1}_{t})^\top}{\zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t}\succeq 0$ is positive semidefinite, so the resulting system is positive definite. Equality constraints in $A$ induce zero diagonal elements in $D_w$ (or in $-D_w^{-1}$). In this case the corresponding rows will not be eliminated and give rise to an indefinite system
of the form $\begin{sbmatrix}
  \tilde H & \tilde A^\top \\ \tilde A & 0
\end{sbmatrix}$ with a large positive definite block $\tilde H$ and hopefully few
further rows in $\tilde A$. For such systems it is well studied how to
employ a preconditioner for $\tilde H$ to solve the full indefinite
system with \emph{e.g.}\ \ MINRES, see \cite{ElmanSilvesterWathen2006}.

The cost of multiplying the full KKT-matrix of \eqref{eq:fullKKT} by a
vector is roughly the same as that of multiplying $\tilde H$ by a vector. Indeed, the same
multiplications arise for $H+D_y,A,A^\top,B,B^\top$. So it remains to compare
the cost of a multiplication by $\begin{sbmatrix}
 -{\mathfrak{X}}^{-1} &-\mathbbm{1}_t\\-\mathbbm{1}_t^\top &\zeta^{-1}\sigma 
\end{sbmatrix}$
to a multiplication by ${\mathfrak{X}}_t^{\frac12}(I-\tfrac{{\mathfrak{X}}_t^{\frac12}\mathbbm{1}_{t}({\mathfrak{X}}_t^{\frac12}\mathbbm{1}_{t})^\top}{\zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t}){\mathfrak{X}}_t^{\frac12}={\mathfrak{X}}_t-\tfrac{{\mathfrak{X}}_t\mathbbm{1}_{t}({\mathfrak{X}}_t\mathbbm{1}_{t})^\top}{\zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t}$. Recall that ${\mathfrak{X}}_t$ is a block diagonal matrix with a separate block for each cone $\mathbb{R}_+$, $\soc{q_i}$, $\mathbb{S}^{s_i}_+$ specified by $t$ and the cost of
multiplying by ${\mathfrak{X}}_t$ or ${\mathfrak{X}}_t^{-1}$ is identical. Thus the only difference
are the multiplications by $\mathbbm{1}_t$ in the first case and by the precomputed  vector ${\mathfrak{X}}_t\mathbbm{1}_{t}$ in the second. The vector ${\mathfrak{X}}_t\mathbbm{1}_{t}$ may be formed at almost
negligible cost along with setting up ${\mathfrak{X}}_t$. So there is no noteworthy difference in the cost of matrix vector multiplications between the two systems and no
structural advantages are lost when working with $\tilde H$ instead of the full system. We will therefore concentrate on developing a preconditioner for $\tilde H$.

For this note that $\tilde H$ of \eqref{eq:SchurH} arises from adding a
Gram matrix to a positive diagonal,
\begin{displaymath}
  \tilde H=D+VV^\top,
\end{displaymath}
where (recall $H=D_H+V_HV_H^\top$)
\begin{equation}\label{eq:tildeV}
  D=D_H+D_y \quad\text{and}\quad V=\left[V_H, A^\top D_w^{\frac12}, B^\top{\mathfrak{X}}_t^{\frac12}(I-\tfrac{{\mathfrak{X}}_t^{\frac12}\mathbbm{1}_{t}({\mathfrak{X}}_t^{\frac12}\mathbbm{1}_{t})^\top}{\zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t})^{\frac12}\right].
\end{equation}
Note that the multiplication of $V$ with a vector requires only a little bit more than half the number of operations of multiplying $\tilde H$ (or the full KKT matrix) with a vector. This suggests to explore possibilities of finding low rank approximations of $V$ for preconditioning. 

\section{Low rank preconditioning a Gram-matrix plus positive diagonal}\label{S:lowrankprecond}

Consider a matrix
\begin{displaymath}
  H=D+VV^\top
\end{displaymath}
with a positive definite matrix $D\in\mathbb{S}^m_{++}$ and
$V\in\mathbb{R}^{m\times n}$. In our application $D$ is diagonal, but the
results apply for general $D\succ 0$. This is applicable in practice
as long as $D^{-1}$ can be applied efficiently to vectors.
Matrix $V$ is assumed to be given by a matrix-vector multiplication
oracle, \mbox{i.\,e.,}\  $V$ and $V^\top$ may be multiplied by vectors but the
matrix does not have to be available explicitly.  

For motivating the following preconditioning approaches, first
consider computing the singular value decomposition of
$$D^{-\frac12}V=Q_H\left[\begin{array}{c}\Sigma\\
    0\end{array}\right]P_H^\top$$ with orthogonal
$Q_H\in\mathbb{R}^{m\times m}$, diagonal
$\Sigma=\mathop{\rm Diag}(\sigma_1,\dots,\sigma_n)$ ordered nonincreasingly by
$\sigma_1\ge\cdots\ge\sigma_n\ge 0$ and orthogonal
$P_H\in\mathbb{R}^{n\times n}$ (for convenience, it is assumed that $n\le
m$). Then
\begin{displaymath}
  H=D^{\frac12}
  Q_H\left[\begin{array}{cc}I_{ n}+\Sigma^2& 0 \\ 0 & I_{m- n}\end{array}\right]
   Q_H^\top D^{\frac12},
\end{displaymath}
\begin{displaymath}
  H^{-1}=D^{-\frac12}
Q_H\left[\begin{array}{cc}(I_{ n}+\Sigma^2)^{-1}& 0 \\ 0 &
                                                                      I_{m- n}\end{array}\right]
Q_H^\top D^{-\frac12}.
\end{displaymath}
When $\Sigma$ is replaced by the $k$ largest singular values, this gives rise to a
good ``low rank'' preconditioner, see Theorem~\ref{T:kappa} below.  Computing the full matrix
$D^{-\frac12}V$ and its singular value decomposition will in general
be too costly or even impossible.  Instead the general idea is to work
with $D^{-\frac12}V\Omega$ for some random or deterministic choice of
$\Omega\in\mathbb{R}^{n\times k}$. 

Multiplying by a random $\Omega$ may be thought of as giving rise to a
subspace approximation in the style of Johnson-Lindenstrauss,
cf.~\cite{Achlioptas2001,DasguptaGupta2002}, and this formed the
starting point of this investigation. The actual randomized approach
and analysis, however, mainly builds on
\cite{HalkoMartinssonTropp2011} and the bounding techniques presented
there. For the deterministic preconditioning variant the
recent work \cite{HabibiKavandKocvaraStingl2021} provided strong
guidance for analysing the condition number.

Here, $\Omega$ will mostly consist of orthonormal columns. Yet it is
instructive to consider more general cases, as well.
An arbitrary $\Omega\in\mathbb{R}^{n\times k}$ gives rise to the preconditioner
\begin{equation}\label{eq:HOmega}
  H_\Omega=D+V\Omega\Omega^\top V=D^{\frac12}
  Q_H\left[\begin{array}{cc}I_{ n}+\Sigma P_H^\top\Omega\Omega^\top P_H\Sigma& 0 \\ 0 & I_{m- n}\end{array}\right]
   Q_H^\top D^{\frac12}.
\end{equation}
Putting $G_{\Omega}=D^{\frac12}Q_H
\begin{sbmatrix}
  I_{ n}+\Sigma P_H^\top\Omega\Omega^\top P_H\Sigma & 0  \\ 0 & I_{m-n}
\end{sbmatrix}^{\frac12}$ we have $H_{\Omega}=G_{\Omega}G_{\Omega}^\top$. The preconditioner is the better the closer $G_{\Omega}^{-1}HG_{\Omega}^{-T}$ is to the identity. In the analysis of convergence rates, see \emph{e.g.}\  \cite{ElmanSilvesterWathen2006}, this enters via the condition number 
\begin{displaymath}
  \kappa_{\Omega}:=\frac{\lambda_{\max}(G_{\Omega}^{-1}HG_{\Omega}^{-T})}{\lambda_{\min}(G_{\Omega}^{-1}HG_{\Omega}^{-T})}=\frac{\lambda_{\max}(H^{\frac12}H_{\Omega}^{-1}H^{\frac12})}{\lambda_{\min}(H^{\frac12}H_{\Omega}^{-1}H^{\frac12})}=\frac{\lambda_{\max}(H_{\Omega}^{-\frac12}HH_{\Omega}^{-\frac12})}{\lambda_{\min}(H_{\Omega}^{-\frac12}HH_{\Omega}^{-\frac12})}.
\end{displaymath}
In this, the equations follow from $BB^\top$ and $B^\top B$ having the same eigenvalues for $B\in\mathbb{R}^{n\times n}$. 
\begin{theorem}\label{T:kappa} Let $H=D+VV^\top\in \mathbb{S}^m_{++}$ with positive
  definite $D\in\mathbb{S}^m_{++}$ and $V\in \mathbb{R}^{m\times n}$ with
  $n<m$ and singular value decomposition
  $D^{-\frac12}V=Q_H\Sigma P_H^\top$, $Q_H^\top Q_H=I_m$, $P_H^\top
  P_H=I_n$, $\Sigma=\mathop{\rm Diag}(\sigma_1\ge\dots\ge\sigma_n)\in\mathbb{S}^n_+$. For  $\Omega\in\mathbb{R}^{n\times k}$
  the preconditioner $H_\Omega$ of \eqref{eq:HOmega} results in condition number  
  \begin{eqnarray*}
    \kappa_\Omega&=&\frac{\max\{1,\lambda_{\max}((I_{ n}+\Sigma^2)^{\frac12}(I_{ n}+\Sigma P_H^\top \Omega\Omega^\top P_H\Sigma)^{-1}(I_{ n}+\Sigma^2)^{\frac12})\}}{\min\{1,\lambda_{\min}((I_{ n}+\Sigma^2)^{\frac12}(I_{ n}+\Sigma P_H^\top\Omega\Omega^\top P_H\Sigma)^{-1}(I_{ n}+\Sigma^2)^{\frac12})\}}\\
    &=&\frac{\max\{1,\lambda_{\max}((I_{ n}+\Sigma^2)^{-\frac12}(I_{ n}+\Sigma P_H^\top \Omega\Omega^\top P_H\Sigma)(I_{ n}+\Sigma^2)^{-\frac12})\}}{\min\{1,\lambda_{\min}((I_{ n}+\Sigma^2)^{-\frac12}(I_{ n}+\Sigma P_H^\top\Omega\Omega^\top P_H\Sigma)(I_{ n}+\Sigma^2)^{-\frac12})\}}
  \end{eqnarray*}
  In particular, for $0\le k<n$ and $\Omega=(P_H)_{\bullet,[1,\dots,k]}$ the condition number's value is $1+\sigma_{k+1}^2$.
\end{theorem}
\begin{proof} For $G_\Omega$ as above direct computation yields
\begin{eqnarray*}
  G_{\Omega}^{-1}HG_{\Omega}^{-T}&=&\begin{sbmatrix}
  I_{ n}+\Sigma\Omega\Omega^\top\Sigma & 0  \\ 0 & I_{m-n}
\end{sbmatrix}^{-\frac12}\begin{sbmatrix}
  I_{ n}+\Sigma^2 & 0  \\ 0 & I_{m-n}
\end{sbmatrix}\begin{sbmatrix}
  I_{ n}+\Sigma\Omega\Omega^\top\Sigma & 0  \\ 0 & I_{m-n}
\end{sbmatrix}^{-\frac12}\\
 &=&\begin{sbmatrix}
  (I_{ n}+\Sigma\Omega\Omega^\top\Sigma)^{-\frac12}(I_{ n}+\Sigma^2)(I_{ n}+\Sigma\Omega\Omega^\top\Sigma)^{-\frac12} & 0  \\ 0 & I_{m-n}
\end{sbmatrix}.                      
\end{eqnarray*}
The eigenvalues of $(I_{ n}+\Sigma\Omega\Omega^\top\Sigma)^{-\frac12}(I_{ n}+\Sigma^2)(I_{ n}+\Sigma\Omega\Omega^\top\Sigma)^{-\frac12}$ coincide with
those of $(I_{ n}+\Sigma^2)^{\frac12}(I_{ n}+\Sigma\Omega\Omega^\top\Sigma)^{-1}(I_{ n}+\Sigma^2)^{\frac12}$, because for $B=(I_{ n}+\Sigma\Omega\Omega^\top\Sigma)^{-\frac12}(I_{ n}+\Sigma^2)^{\frac12}$ the two matrices are $BB^\top$ and $B^\top B$. This gives rise to the first line. The second follows because for positive definite $A$ there holds $\lambda_{\max}(A)=1/\lambda_{\min}(A^{-1})$ and $\lambda_{\min}(A)=1/\lambda_{\max}(A^{-1})$. 
\end{proof}
Consider now a subspace spanned by $k$ orthonormal columns collected in
some matrix $P_ \ell\in \mathbb{R}^{n\times k}$ which hopefully generates most of
the large directions of  $D^{-\frac12}V$. In this orthonormal case a simpler bound on the condition number may be obtained by following the argument of Th.~5.1 in \cite{HabibiKavandKocvaraStingl2021}.
\begin{theorem}\label{T:deterministiccond}
  Let $H=D+VV^\top$ with positive definite $D\in\mathbb{S}^m_{++}$ and
  general $V\in\mathbb{R}^{m\times n}$, let $P=[P_\ell,P_s]\in\mathbb{R}^{n\times n}$,
  $PP^\top=I_n$. Preconditioner $H_{P_\ell}=D+VP_\ell P_\ell^\top
  V^\top$ has condition number $\kappa_{P_\ell}\le
  1+\lambda_{\max}(H_{P_\ell}^{-\frac12}VP_sP_s^\top V^\top
  H_{P_\ell}^{-\frac12})$. Equality holds if and only if $\mathop{\rm rank}(VP_s)<m$. 
\end{theorem}
\begin{proof}
  Because $H=D+[VP_\ell,VP_s][VP_\ell,VP_s]^\top=H_{P_\ell}+VP_sP_s^\top V^\top$ we have
  \begin{displaymath}
    H_{P_\ell}^{-\frac12}H H_{P_\ell}^{-\frac12}= I_m+H_{P_\ell}^{-\frac12}VP_sP_s^\top V^\top H_{P_\ell}^{-\frac12}.
  \end{displaymath}
  The second summand is positive semidefinite with minimum eigenvalue
  0 if and only if $\mathop{\rm rank}(VP_s)<m$. Thus,
  \begin{eqnarray*}
    \lambda_{\min}(H_{P_\ell}^{-\frac12}HH_{P_\ell}^{-\frac12})&\ge &1,\\
    \lambda_{\max}(H_{P_\ell}^{-\frac12}HH_{P_\ell}^{-\frac12})&= & 1+\lambda_{\max}(H_{P_\ell}^{-\frac12}VP_sP_s^\top V^\top H_{P_\ell}^{-\frac12}). 
  \end{eqnarray*}
  By $\kappa(H_{P_\ell}^{-\frac12}HH_{P_\ell}^{-\frac12})=\frac{\lambda_{\max}(H_{P_\ell}^{-\frac12}HH_{P_\ell}^{-\frac12})}{\lambda_{\min}(H_{P_\ell}^{-\frac12}HH_{P_\ell})}$ the result is proved.
\end{proof}
Building on these two theorems we first analyze randomized approaches
that do not make any assumptions on structural properties of
$D^{-\frac12}V$ but only require a multiplication oracle. Afterwards
we present a deterministic approach that exploits some knowledge of
the bundle subproblem and the interior point algorithm. The
corresponding routines supply a $\hat V=V\Omega$. The actual
preconditioning routine, Alg.~\ref{alg:precond} below, does not use
$H_\Omega$ directly, but a truncated preconditioner $H_{\Omega\hat P}$
that drops all singular values of $D^{-\frac12}V\Omega$ that are
less than one. The inverse is then formed via a Woodbury-formula, see
\cite[\S0.7.4]{HornJohnson85}.  Note, depending on the expected number
of calls to the routine and the structure preserved in $\hat V$, it may
or may not pay off to also precompute $\hat V\hat P$. For diagonal $D$
and dense $\hat V$ the cost of applying this preconditioner is
$O(m+mk+k\hat k)$.

\begin{algorithm}[Preconditioning by truncated
  $H_\Omega=D+V\Omega(V\Omega)^\top$]\label{alg:precond}\quad\\
\textbf{Input:} $v\in\mathbb{R}^m$, $D\in\mathbb{S}^n_{++}$, precomputed $\hat V=V\Omega\in\mathbb{R}^{n\times k}$
and, for $\hat V^\top D^{-1}\hat V=P\mathop{\rm Diag}(\hat
\lambda_1\ge\dots\ge\hat \lambda_k)P^\top$, 
$\hat k=\max\{0,i\colon\hat\lambda_i\ge 1\}$,
$\hat \Lambda=\mathop{\rm Diag}(\hat\lambda_1,\dots,\hat\lambda_{\hat k})$, $\hat
P=P_{\bullet,[1,\dots,\hat k]}$.\\
  \textbf{Output:} $(H_{\Omega\hat P})^{-1}v$.
  \begin{compactenum}
  \item $v\gets D^{-1}v$.
  \item If $\hat k>0$ set $v \gets v - D^{-1}\hat
    V\hat P(I+\hat\Lambda)^{-1}\hat P^\top\hat V^\top v$.
  \item return $v$.
  \end{compactenum}
\end{algorithm}




\subsection{Preconditioning by Random Subspaces}\label{S:randprecond}

For the random subspace approach fix some
$k\in\mathbb{N}$ with $2\le k<n$, let $\Omega$ be a $n\times k$ random matrix whose elements are independently
identically distributed by the normal distribution $\mathcal{N}(0,\frac1{k})$. For this
$\Omega$ consider first the low rank approximation $D^{-\frac12} V\Omega= 
Q_H\mat{c}{ \Sigma\\ 0}
P_H^\top\Omega$. Because the normal distribution is invariant under orthogonal
transformations, we may assume $ P_H=I$ and analyze the setting  $
Q_H\mat{c}{ \Sigma\\ 0}\Omega$ giving rise to the low rank approximation by the random
matrix
\begin{displaymath}
  H_\Omega=D^{\frac12}
  Q_H\left[\begin{array}{cc}I_{ n}+\Sigma\Omega\Omega^\top\Sigma& 0 \\ 0 & I_{m- n}\end{array}\right]
   Q_H^\top D^{\frac12}.
\end{displaymath}
In view of Theorem~\ref{T:kappa} such a preconditioner is good if $(I+\Sigma^2)^{-\frac12}(I+\Sigma
\Omega\Omega^\top\Sigma)(I+\Sigma^2)^{-\frac12}$ is close to the
identity. Based on the Johnson-Lindenstrauss interpretation,
it seems likely that large portions of the spectrum will be close to
one. This can be justified to some extent by studying the moments of
the Ritz values.
\begin{theorem}\label{T:moments}
  Let $\Omega\in\mathbb{R}^{ n\times k}$ have its elements i.i.d.\ according to
  the normal distribution $\mathcal{N}(0,\frac1{k})$, then for any $x\in\mathbb{R}^{ n}$ the
quadratic form
\begin{displaymath}
  q(x)=x^\top(I+\Sigma^2)^{-\frac12}(I+\Sigma
  \Omega\Omega^\top\Sigma)(I+\Sigma^2)^{-\frac12}x
\end{displaymath}
has expected value $\mathbb{E}(q(x))=\|x\|^2$ and variance $\mathop{\mathrm{Var}}(q(x))=\frac{2}{k}\big(\sum_{i=1}^{n}\frac{\sigma_i^2}{1+\sigma_i^2}x_i^2\big)^2$.
\end{theorem}
\begin{proof} Let $\Omega=(\omega_{ij})$ with i.i.d.\ elements
  $\omega_{ij}$ from $\mathcal{N}(0,\frac{1}{k})$. Recall that
  $\mathbb{E}(\omega_{ij})=0$, $\mathbb{E}(\omega_{ij}^2)=\frac1k$,
  $\mathbb{E}(\omega_{ij}^3)=0$, $\mathbb{E}(\omega_{ij}^4)=3/k^2$ and that for
  independent random variables $X,Y$ there holds $\mathbb{E}(XY)=\mathbb{E}(X)\mathbb{E}(Y)$.

  The expected value of the quadratic form evaluates to
\begin{eqnarray*}
  \lefteqn{\mathbb{E}\bigg(x^\top(I+\Sigma^2)^{-\frac12}(I+\Sigma
  \Omega\Omega^\top\Sigma)(I+\Sigma^2)^{-\frac12}x\bigg)=}\\
  &=&\sum_{i=1}^{
      n}\frac{1}{1+\sigma_i^2}x_i^2+\mathbb{E}\bigg(\sum_{h=1}^k\big(\sum_{i=1}^{
      n}\frac{\sigma_i}{\sqrt{1+\sigma_i^2}} \omega_{hi} x_{i}\big)^2\bigg)\\
  &=&\sum_{i=1}^{
      n}\frac{1}{1+\sigma_i^2} x_i^2+\sum_{h=1}^k\sum_{i=1}^{
      n}\frac{\sigma_i^2}{1+\sigma_i^2}x_{i}^2\mathbb{E} \omega_{hi}^2 \\
  &=& \sum_{i=1}^{ n} x_i^2\bigg(\frac{1}{1+\sigma_i^2} +
      \frac{\sigma_i^2}{1+\sigma_i^2}\sum_{h=1}^k\mathbb{E} \omega_{hi}^2\bigg)\\
  &=& \sum_{i=1}^{ n} x_i^2\bigg(\frac{1}{1+\sigma_i^2} +
      \frac{\sigma_i^2}{1+\sigma_i^2}\sum_{h=1}^k\frac1k\bigg) =
      \sum_{i=1}^{ n} x_i^2= \|x\|^2.
\end{eqnarray*}

For determining the variance, the second moment may be computed as follows.
\begin{eqnarray*}
\lefteqn{\mathbb{E}\bigg(\bigg[\sum_{h=1}^k\big(\sum_{i=1}^{
      n}\frac{\sigma_i}{\sqrt{1+\sigma_i^2}} \omega_{hi} x_{i}\big)^2\bigg]^2\bigg)=}\\
  &=&\mathbb{E}\bigg(\bigg[\sum_{h=1}^k\sum_{i=1}^{
      n}\sum_{i'=1}^{
      n}\frac{\sigma_i}{\sqrt{1+\sigma_i^2}}\frac{\sigma_{i'}}{\sqrt{1+\sigma_{i'}^2}}\omega_{hi} x_{i}\omega_{hi'} x_{i'}\bigg]^2\bigg)\\
  &=&\mathbb{E}\bigg(\sum_{h,h'=1}^k~
      \sum_{i,i',j,j'=1}^{ n} \\
  &&\phantom{\mathbb{E}\bigg(}\frac{\sigma_i}{\sqrt{1+\sigma_i^2}}\frac{\sigma_{i'}}{\sqrt{1+\sigma_{i'}^2}}\frac{\sigma_j}{\sqrt{1+\sigma_j^2}}\frac{\sigma_{j'}}{\sqrt{1+\sigma_{j'}^2}}x_{i} x_{i'}x_{j} x_{j'}\omega_{hi}\omega_{hi'}\omega_{h'j}\omega_{h'j'}\bigg).
\end{eqnarray*}
In the cases of $h\neq h'$ only terms with $i=i'$ and $j=j'$ are not zero. These
evaluate to $\mathbb{E} \omega_{hi}^2 \mathbb{E} \omega_{h'j}^2=\frac1{k^2}$ giving
\begin{eqnarray*}
  \frac{k(k-1)}{k^2}\sum_{i=1}^{
      n}\sum_{j=1}^{
      n}\frac{\sigma_i^2}{1+\sigma_i^2}x_{i}^2
     \frac{\sigma_j^2}{1+\sigma_j^2}x_{j}^2.
\end{eqnarray*}
For each $h=h'$ there remain $(i=i'=j=j')$ with value $\mathbb{E} \omega_{hi}^4=\frac3{k^2}$,
\begin{eqnarray*}
  \frac{3k}{k^2}\sum_{i=1}^{ n}\bigg(\frac{\sigma_i^2}{1+\sigma_i^2}\bigg)^2x_{i}^4,
\end{eqnarray*}
and the three pairings $(i=i',j=j')$,
$(i=j,i'=j')$ and $(i=j',i'=j)$ each with value $\frac1{k^2}$,
\begin{eqnarray*}
  \frac{3k}{k^2}\sum_{i\neq j}\frac{\sigma_i^2}{1+\sigma_i^2}x_{i}^2
     \frac{\sigma_j^2}{1+\sigma_j^2}x_{j}^2.
\end{eqnarray*}
Summing up these three expressions yields
\begin{displaymath}
  \mathbb{E}\bigg(\bigg[\sum_{h=1}^k\big(\sum_{i=1}^{
      n}\frac{\sigma_i}{\sqrt{1+\sigma_i^2}} \omega_{hi} x_{i}\big)^2\bigg]^2\bigg)=\frac{k^2+2k}{k^2}\bigg(\sum_{i=1}^{ n}\frac{\sigma_i^2}{1+\sigma_i^2}x_i^2\bigg)^2.
  \end{displaymath}
The result now follows from the usual $\mathop{\mathrm{Var}} X=\mathbb{E}(X^2)-(\mathbb{E} X)^2$ for any random variable $X$. 
\end{proof}
This suggests that even for relatively small $k$ the behavior of the
preconditioned system may be expected to be reasonably close to the
identity for a large portion of the directions.  The result, however,
does not seem to open a path towards good bounds on the condition
number.


A first possibility is offered by Theorem~\ref{T:deterministiccond}. 
Recall that for an arbitrary matrix $A\in\mathbb{R}^{m\times n}$ the projector
\begin{displaymath}
  \mathbf{P}_A=A(A^\top A)^\dagger A^\top
\end{displaymath}
projects any vector of $\mathbb{R}^m$ onto the range space of $A$ and
$\mathbf{P}$ depends only on this range space. Computationally
it may be determined by a $QR$-decomposition of $A=Q_AR_A$  with orthogonal $Q_A\in\mathbb{R}^{m\times n'}$ for $n'=\mathop{\rm rank}(A)$ via $\mathbf{P}_A=Q_AQ_A^\top$. The formula allows to verify $\mathbf{P}_A=\mathbf{P}_A^\top$, $\mathbf{P}_A\mathbf{P}_A=\mathbf{P}_A$ and $\mathbf{P}_AA=A$ by direct computation. Furthermore, for any $B\in\mathbb{R}^{n\times h}$ there holds $\mathbf{P}_{AB}\preceq \mathbf{P}_A\preceq I_m$ because of the containment relations between the ranges. 

In the following $\Omega\Omega^\top$ in $H_{\Omega}$ will be replaced by
the projector $\mathbf{P}_\Omega$. The random low rank approximation to be considered reads
\begin{displaymath}
  H_{P_\Omega}=D^{\frac12}
  Q_H\left[\begin{array}{cc}I_{ n}+\Sigma \mathbf{P}_\Omega\Sigma& 0 \\ 0 & I_{m- n}\end{array}\right]
   Q_H^\top D^{\frac12}.
\end{displaymath}
\begin{corollary}
   Let $H=D+VV^\top\in \mathbb{R}^m$ with positive
   definite $D\in\mathbb{S}^m_{++}$ and $V\in \mathbb{R}^{m\times n}$.
 
 
 
 
  Given $\Omega\in\mathbb{R}^{n\times k}$, let
  $\mathbf{P}_\Omega=\Omega(\Omega^\top \Omega)^\dagger\Omega^\top$. For the
  preconditioner $H_{\mathbf{P}_\Omega}$ the condition number satisfies
  $\kappa_{\mathbf{P}_\Omega}\le 1+\|D^{-\frac12}V(I-\mathbf{P}_{\Omega})\|^2$,
  where $\|\cdot\|$ denotes the spectral norm. 
\end{corollary}
\begin{proof}
  Let $\mathbf{P}_\Omega=P_\ell P_\ell^\top$ with $P_\ell\in\mathbb{R}^{n\times
    k'}$ for some $k'\le k$ and $P_\ell^\top P_\ell=I_{k'}$. Add
  orthonormal columns $P_s$ so that $P=[P_\ell,P_s]$ satisfies
  $PP^\top=I_n$. Note, $I-\mathbf{P}_\Omega=P_sP_s^\top$ is the projector onto
  the orthogonal complement. Use this choice in Theorem~\ref{T:deterministiccond},
  then $H_{\mathbf{P}_\Omega}=H_{P_\ell}\succeq D$. Observe that by $D\in \mathbb{S}^m_{++}$
  for any $\lambda\ge 0$, $A\in \mathbb{S}^m$  the relation $D^{-\frac12}AD^{-\frac12}\preceq
  \lambda I_m$ is equivalent to $A\preceq \lambda D$ and implies $A\preceq
  \lambda (D+VP_\ell P_\ell^\top V^\top)=\lambda H_{P\ell}$ which is
  equivalent to
  $H_{P_\ell}^{-\frac12}AH_{P_\ell}^{-\frac12}\preceq
  \lambda I_m$. Therefore  
  \begin{displaymath}
    \lambda_{\max}(H_{P_\ell}^{-\frac12}VP_sP_s^\top V^\top
    H_{P_\ell}^{-\frac12})\le\lambda_{\max}(D^{-\frac12}VP_sP_s^\top V^\top
    D^{-\frac12})=\|D^{-\frac12}V(I-\mathbf{P}_\Omega)\|^2.
  \end{displaymath}
\end{proof}
While this bound is rather straight forward to derive, it does not
seem strong enough to observe a reduced influence of the largest
singular values of $D^{-\frac12}V$.  Indeed, in its derivation only the
diagonal of $H_{\mathbf{P}_\Omega}$ was considered and the
influence of $V\Omega$ is lost.

In order to obtain stronger bounds, the rather involved techniques layed out in
\cite{HalkoMartinssonTropp2011} seem to be required.
The next steps and results follow their arguments closely. This time
the $H_{P_\Omega}$-part is kept inverted in the analysis of the
condition number.

Because $I+\Sigma P_\Omega\Sigma\preceq I+\Sigma^2$ there holds
$$ (I+\Sigma^2)^{\frac12}(I+\Sigma
P_{\Omega}\Sigma)^{-1}(I+\Sigma^2)^{\frac12}\succeq I.$$ By Theorem~\ref{T:kappa} the condition number is bounded by
\begin{displaymath}
  \kappa_{P_\Omega}\le \lambda_{\max}( (I+\Sigma^2)^{\frac12}(I+\Sigma P_{\Omega}\Sigma)^{-1}(I+\Sigma^2)^{\frac12})
\end{displaymath}
and will attain it, whenever $n<m$.  In terms of $\Omega$, the best
possible outcome is an event resulting in $P_\Omega=\mat{c}{I_{k}\\0}$ (see, \emph{e.g.}\ 
\cite[7.4.52]{HornJohnson85}). It corresponds to the truncated SVD and
gives $\kappa_{
  \begin{sbmatrix}
I_{k}\\0
\end{sbmatrix}
} = 1+\sigma_{k+1}^2$. Aiming for something more
realistic, one hopes for a good coverage of the first $k$ singular
values when oversampling with $p$ additional columns.

The first step in the analysis is to obtain a deterministic bound for
a fixed $\Omega\in\mathbb{R}^{ n\times(k+p)}$ as outlined in
\cite[\S9.2]{HalkoMartinssonTropp2011}.
\begin{theorem}\label{th:deterministicbound}
  Given $\sigma_1\ge\dots\ge\sigma_{ n} \ge 0$ and a matrix $\Omega\in\mathbb{R}^{ n\times (k+p)}$ with $k\le  n$ so that the
  first $k$ rows of $\Omega$ are linearly independent, split $\Sigma=\mat{cc}{\Sigma_1& 0\\ 0&\Sigma_2}$ into blocks $\Sigma_1=\mathop{\rm Diag}(\sigma_1,\dots,\sigma_k)$ and $\Sigma_2=\mathop{\rm Diag}(\sigma_{k+1},\dots,\sigma_{ n})$ and
$\Omega=\mat{c}{\Omega_1\\\Omega_2}$ into the first $k$ rows
$\Omega_1\in\mathbb{R}^{k\times(k+p)}$ and the last $ n-k$ rows
$\Omega_2\in\mathbb{R}^{( n-k)\times k+p}$. Then
\begin{displaymath}
  \lambda_{\max}( (I+\Sigma^2)^{\frac12}(I+\Sigma P_{\Omega}\Sigma)^{-1}(I+\Sigma^2)^{\frac12})\le 2+ \sigma_{k+1}^2+\|(I+\Sigma_2^2)^{\frac12}\Omega_2\Omega_1^\dagger\|^2.
\end{displaymath}
\end{theorem}
\begin{proof}
By assumption $\Omega_1$ has full row rank and the range space of the matrix
\begin{displaymath}
  Z=\Omega\cdot\Omega_1^\dagger= \mat{c}{I_k\\ F=\Omega_2\Omega_1^\dagger} \in\mathbb{R}^{ n\times k}
\end{displaymath}
is contained in the range space of $\Omega$. Hence $P_Z\preceq P_\Omega$ and
\begin{equation}\label{eq:lmaxle}
\lambda_{\max}( (I+\Sigma^2)^{\frac12}(I+\Sigma P_{\Omega}\Sigma)^{-1}(I+\Sigma^2)^{\frac12})\le 
\lambda_{\max}( (I+\Sigma^2)^{\frac12}(I+\Sigma P_{Z}\Sigma)^{-1}(I+\Sigma^2)^{\frac12}).
\end{equation}
The projector $P_Z$ computes to
\begin{displaymath}
  P_Z=\mat{c}{I\\ F}\mat{c}{I+F^\top F}^{-1}\mat{c}{I\\mathbb{F}}^\top.
\end{displaymath}
Use this in the Woodbury-formula for inverses of rank adjustments \cite[\S0.7.4]{HornJohnson85} for $(I+\Sigma P_{Z}\Sigma)^{-1}$ to obtain
\begin{eqnarray}
(I+\Sigma P_{Z}\Sigma)^{-1}&=&I-
                                           \mat{c}{\Sigma_1\\\Sigma_2
  F}\mat{c}{I+F^\top
  F+\Sigma_1^2+F^\top\Sigma_2^2F}^{-1}\mat{c}{\Sigma_1\\ \Sigma_2F
  }^\top\nonumber\\
                                       &=& I-\mat{c}{\Sigma_1\\\Sigma_1F}\mat{c}{I+\Sigma_1^2+F^\top(I+\Sigma_2^2)F}^{-1}\mat{c}{\Sigma_1\\\Sigma_2F}^\top\nonumber\\
  &\preceq& I- \mat{c}{\Sigma_1\\\Sigma_2F}\mat{c}{I+\Sigma_1^2+\|(I+\Sigma_2^2)^{\frac12}F\|^2I}^{-1}\mat{c}{\Sigma_1\\\Sigma_2F}^\top\hspace*{-1.5ex}.\label{eq:invpreceq}
\end{eqnarray}
The last line follows, because
$\|(I+\Sigma_2^2)^{\frac12}F\|^2=\lambda_{\max}(F^\top(I+\Sigma_2^2)F)=:\bar\lambda$
and therefore $F^\top(I+\Sigma_2^2)F\preceq \bar\lambda I$ giving
\begin{displaymath}
  I+\Sigma_1^2+F^\top(I+\Sigma_2^2)F\preceq
  I+\Sigma_1^2+\bar\lambda I \quad \Leftrightarrow\quad
  (I+\Sigma_1^2+F^\top(I+\Sigma_2^2)F)^{-1}\succeq
  (I+\Sigma_1^2+\bar\lambda I)^{-1},
\end{displaymath}
so the relation is implied by semidefinite scaling. Put
$\Lambda=(I+\Sigma_1^2+\bar\lambda I)$ and note that the second
diagonal block of \eqref{eq:invpreceq} asserts
$\Sigma_2F\Lambda^{-1}F^\top\Sigma_2\preceq I$, then
\begin{eqnarray*}
  \lefteqn{S:=(I+\Sigma^2)^{\frac12}(I+\Sigma
  P_{\Omega}\Sigma)^{-1}(I+\Sigma^2)^{\frac12}\overset{\eqref{eq:lmaxle}\eqref{eq:invpreceq}}{\preceq}}\\
&\preceq& (I+\Sigma^2)^{\frac12}\mat{cc}{I-\Sigma_1\Lambda^{-1}\Sigma_1 & -\Sigma_1\Lambda^{-1}F^\top\Sigma_2\\-\Sigma_2F\Lambda^{-1}\Sigma_1 &I-\Sigma_2F\Lambda^{-1}F^\top\Sigma_2} (I+\Sigma^2)^{\frac12}\\
  &\preceq& (I+\Sigma^2)^{\frac12}\mat{cc}{(\Lambda-\Sigma_1^2)\Lambda^{-1} & -\Sigma_1\Lambda^{-1}F^\top\Sigma_2\\-\Sigma_2F\Lambda^{-1}\Sigma_1 &I} (I+\Sigma^2)^{\frac12}\quad :=\bar S.
\end{eqnarray*}
Employing \cite[Prop.\ 8.3]{HalkoMartinssonTropp2011} now results in
  \begin{displaymath}
  \lambda_{\max}(S)\le\lambda_{\max}(\bar S)\le\|(I+\Sigma_1^2)(\Lambda-\Sigma_1^2)\Lambda^{-1}\|+\|I+\Sigma_2^2\|.
\end{displaymath}
The last term evaluates to $\lambda_{\max}(I+\Sigma_2^2)=1+\sigma_{k+1}^2$. For the second last term substituting in the definitions of $\Lambda$ and $\bar\lambda$ yields 
\begin{eqnarray*}
  \lefteqn{\|(I+\Sigma_1^2)(\Lambda-\Sigma_1^2)\Lambda^{-1}\|=}\\
  &=&(1+\|(I+\Sigma_2^2)^{\frac12}\Omega_2 \Omega_1^\dagger\|^2)\cdot \max_{i\in\{1,\dots,k\}} \frac{1+\sigma_i^2}{1+\sigma_i^2+\|(I+\Sigma_2^2)^{\frac12}\Omega_2 \Omega_1^\dagger\|^2}\\
  &=& \max_{i\in\{1,\dots,k\}} \frac{1+\sigma_i^2+\|(I+\Sigma_2^2)^{\frac12}\Omega_2 \Omega_1^\dagger\|^2+\sigma_i^2\|(I+\Sigma_2^2)^{\frac12}\Omega_2 \Omega_1^\dagger\|^2}{1+\sigma_i^2+\|(I+\Sigma_2^2)^{\frac12}\Omega_2 \Omega_1^\dagger\|^2}\\
  &\le& 1+\min\{\sigma_1^2,\|(I+\Sigma_2^2)^{\frac12}\Omega_2 \Omega_1^\dagger\|^2\}.
\end{eqnarray*}
\end{proof}
The current bound falls somewhat short of expectation because of
the identity in $\|(I+\Sigma_2^2)^{\frac12}\Omega_2\Omega_1^\dagger\|^2$.
By Theorem~\ref{T:kappa} and $I_n\preceq I_n+\Sigma P_\Omega \Sigma\preceq I_n+\Sigma^2$,
the use of projectors will never result in condition numbers larger than $1+\sigma_1^2$, so the
influence of the dimension seems to be too dominant in this.
Maybe a better bound is achievable by a more
sophisticated argument. 

The deterministic bound allows to also make use of the
probabilistic bounds on $\|(I+\Sigma_2^2)^{\frac12}\Omega_2 \Omega_1^\dagger\|$ given in \cite{HalkoMartinssonTropp2011}. These shed some light on the advantage of employing oversampling by $p$ additional random vectors in $\Omega$. In our application, oversampling
corresponds to computing the singular values of
$D^{-\frac12}V\mathbf{P}_\Omega$ for $k+p$ columns in order to get better
control on the $k$ largest singular values of $D^{-\frac12}V$ by the preconditioner $H_{\mathbf{P}_\Omega}$.
\begin{theorem}\label{th:probabilisticbound}
  In the setting of Theorem~\ref{th:deterministicbound} let
  $\Omega$ be drawn as a standard Gaussian matrix. Then
  \begin{eqnarray*}
    \lefteqn{\mathbb{E}\,
    \lambda_{\max}\big((I+\Sigma^2)^{\frac12}(I+\Sigma
    P_{\Omega}\Sigma)^{-1}(I+\Sigma^2)^{\frac12}\big)\le}\\
    &\le&2+
          \sigma_{k+1}^2+\bigg(\sqrt{\frac{k}{p-1}}(1+\sigma_{k+1}^2)
          +\frac{e\sqrt{k+p}}{p}\big(\sum_{i=k+1}^{
          n}(1+\sigma_i^2)\big)^{\frac12}\bigg)^2.
  \end{eqnarray*}
  Furthermore, if $p\ge 4$ then for all $u,t\ge 1$ the probability for
  \begin{eqnarray*}
    \lefteqn{\lambda_{\max}\big((I+\Sigma^2)^{\frac12}(I+\Sigma
    P_{\Omega}\Sigma)^{-1}(I+\Sigma^2)^{\frac12}\big)>}\\
    &>&2+
          \sigma_{k+1}^2+\bigg(t\big[\sqrt{\tfrac{3k}{p+1}}+u\tfrac{e\sqrt{k+p}}{p+1}\big](1+\sigma_{k+1}^2)
          +t\tfrac{e\sqrt{k+p}}{p+1}\big(\sum_{i=k+1}^{
          n}(1+\sigma_i^2)\big)^{\frac12}\bigg)^2    
      \end{eqnarray*}
      is at most $2t^{-p}+e^{-u^2/2}$. \\
      The same bounds hold for the condition number
      $\kappa_{P_\Omega}$.
\end{theorem}
\begin{proof}
A central and complex step in \cite[proof of Th.\ 10.2]{HalkoMartinssonTropp2011} is to establish the relation
  \begin{displaymath}
 \mathbb{E} \|(I+\Sigma_2^2)^{\frac12}\Omega_2 \Omega_1^\dagger\|\le\sqrt{\frac{k}{p-1}}\|(1+\Sigma_2^2)^{\frac12}\|+\frac{e\sqrt{k+p}}{p}\|(I+\Sigma_2^2)^{\frac12}\|_F
\end{displaymath}
which directly yields the bound on the expected value via
Theorem~\ref{th:deterministicbound}. 

Likewise, in \cite[proof of Th.\ 10.8]{HalkoMartinssonTropp2011} the
authors derive for $p\ge 4$ and $u,t\ge 1$ 
\begin{displaymath}
\begin{array}{r@{\;}l}  
  \mathop{\rm Prob}\bigg\{\|(I+\Sigma_2^2)^{\frac12}\Omega_2
    \Omega_1^\dagger\|>&\|(I+\Sigma_2^2)^{\frac12}\|\sqrt{\frac{3k}{p+1}}\cdot
    t+\|(I+\Sigma_2^2)^{\frac12}\|_F\frac{e\sqrt{k+p}}{p+1}\cdot
                         t\\
  &+\|(I+\Sigma_2^2)^{\frac12}\|\frac{e\sqrt{k+p}}{p+1}\cdot
    ut\bigg\}\le 2t^{-p}+e^{-u^2/2}.
    \end{array}
\end{displaymath}
\end{proof}
Again, the presence of the identity in the deterministic bound of
Theorem~\ref{th:deterministicbound} has a major
impact also in these probabilistic bounds. Indeed, one would hope
that a better deterministic bound helps to prove stronger decay. 


Without some a priori knowledge of the singular values of
$D^{-\frac12}V\in\mathbb{R}^{m\times n}$ it is difficult to determine a
suitable number of columns for $\Omega$, \mbox{i.\,e.,}\  a suitable dimension of
the random subspace. For huge $m$ the Johnson-Lindenstrauss result as
presented in \cite{DasguptaGupta2002} suggests that for $k$ at most
$4(\varepsilon^2/2-\varepsilon^3/3)^{-1}\ln m$ a suitably chosen random $\Omega\in\mathbb{R}^{n\times k}$ results
in a distortion of $1\pm\varepsilon$ with sufficiently high
probability. This roughly translates to that
each matrix element of $D^{-\frac12}VV^\top D^{-\frac12}$ and
$D^{-\frac12}V\Omega\Omega^\top V^\top D^{-\frac12}$ differs by at
most this factor. When considering the sizes of $m$ aimed for here ---
the dimension of the design space will be a few thousands to a few
hundred thousands --- this number is still too large for efficient
computations even for a moderate $\varepsilon=0.1$. Indeed, the burden
of forming the preconditioner and of applying it would exceed the gain
by far. \cite{HalkoMartinssonTropp2011} propose an algorithmic variant
for identifying a significant drop in singular values, but this requires
successive matrix vector multiplications and these are quite costly in
practice. In preliminary experiments with a number of tentative
randomized variants, those relating the number of columns to the number
of matrix-vector multiplications of the previous solve seemed
reasonable. It will turn out, however, that even the cost of this is
too high and the gain too small in comparison to the deterministic
approach of the next section. The latter appears to capture the
important directions quite well and offers better possibilities to
exploit structural properties of the data. Due to the rather clear
superiority of the deterministic approach, the numerical experiments of
Section \ref{S:experiments} will only present results for one
particular randomized variant that performed best among the tentative
versions. It attempts to identify the
most relevant subspace by storing and extending the important
directions generated in the previous round. For completeness and
reproducability, its details are given in Alg.~\ref{alg:randprecond},
but in view of the rather disencouraging results we refrain from further
discussions.

\begin{algorithm}[Randomized subspace selection forming $\hat V=VP_\Omega$]\label{alg:randprecond}\quad\\
  \textbf{Input:} $V\in\mathbb{R}^{m\times n}$, $D\in\mathbb{S}^n_{++}$,
  previous relevant subspace $P_{old}\in\mathbb{R}^{n\times
    \underline{k}}$ (initially $\underline{k}=0$),
  previous number of multiplications $n_{mult}$, previous $\hat k$ of Alg.~\ref{alg:precond}\\
  \textbf{Output:} $\hat V$ (and stores $P_{old}$)
  \begin{compactenum}
  \item If ($\underline{k}=0$) then
    \begin{compactenum}[(a)]
    \item set $k=\min\{n,3+2\hat k,\lceil\sqrt{n_{mult}\frac{n_{mult}+n}4}-\frac{n_{mult}}2\rceil\}$,
    \item generate a standard Gaussian $\Omega\in\mathbb{R}^{n\times k}$ and set $P_\Omega\gets\Omega$,
    \end{compactenum}
    else
    \begin{compactenum}[(a)]
    \item set $k_+=\max\{3,\lfloor\frac{\sqrt{n_{mult}}}2\rfloor-\underline{k}\}$,
    \item generate a standard Gaussian $\Omega\in\mathbb{R}^{n\times k_+}$ and set $P_\Omega\gets[P_{old},\Omega]$.
    \end{compactenum}
  \item Orthonormalize $P_\Omega$, reset $k$ to its number of columns, set $\hat V=V P_\Omega$.
  \item Compute eigenvalue decomposition $\hat V^\top D^{-1}\hat V=P\mathop{\rm Diag}(\hat
    \lambda_1\ge\dots\ge\hat \lambda_k)P^\top$.
  \item Compute threshold $\bar\lambda=\max\{10,e^{\frac1{10}\ln
      \hat\lambda_1-\frac9{10}\ln\hat\lambda_k}\}$ (enforce $\hat\lambda_k>0$).
  \item Set $\underline{k}\gets\min\big\{k,\max\{3,i>3\colon
    \hat\lambda_i>\bar\lambda\}\big\}$ and set $P_{old}\gets
    P_{\Omega}P_{\bullet,[1,\dots,\underline{k}]}$.
    \item Return $\hat V$. 
  \end{compactenum}
\end{algorithm}


\subsection{A Deterministic Subspace Selection Approach}\label{S:detprecond}

In the conic bundle method, $H=D+VV^\top$ of
Theorem~\ref{T:deterministiccond} is of the form described in
\eqref{eq:tildeV}. An inspection of the column blocks of this $V$
suggests to concretize the bound of Theorem~\ref{T:deterministiccond} for
interior point related applications as follows. 
\begin{theorem}\label{T:speccondbound}
  Given $D\in\mathbb{S}^{m}_{++}$ and $B\in\mathbb{R}^{n\times m}$, let $X\in\mathbb{S}^n_+$ have eigenvalue decomposition $X=[P_\ell,P_s]
  \begin{sbmatrix}
    \Lambda_\ell & 0\\
    0 & \Lambda_s
  \end{sbmatrix}[P_\ell,P_s]^\top$ with $[P_\ell,P_s]^\top[P_\ell,P_s]=I_n$ and diagonal $\Lambda_\ell\in\mathbb{S}^k_+$, $\Lambda_s\in\mathbb{S}^{n-k}_+$. Put $V=B^\top X^{\frac12}$.
  For $H=D+VV^\top$ and preconditioner $H_{P_\ell}=D+VP_\ell P_\ell^\top
  V^\top$ the condition number is bounded by 
  \begin{displaymath}
  \kappa_{P_\ell}\le 1+\sum_{i=1}^{n-k}(\Lambda_s)_{ii}\|B^\top (P_s)_{\bullet,i}\|^2_{D^{-1}}\le 1+(n-k)\bar\rho\bar\beta^2.
\end{displaymath}
where $\bar\rho=\max_{i=1,\dots,n-k}(\Lambda_s)_{ii}$ and $\bar\beta=\max_{i=1,\dots,n}\|B_{i,\bullet}\|_{D^{-1}}$.

\end{theorem}
\begin{proof}
  We show $\lambda_{\max}(H_{P_\ell}^{-\frac12}VP_sP_s^\top V^\top H_{P_\ell}^{-\frac12})\le \sum_{i=1}^{n-k}(\Lambda_s)_{ii}\|B^\top (P_s)_{\bullet,i}\|^2_{D^{-1}}$, then the statement follows by Theorem~\ref{T:deterministiccond}. Note that
$VP_sP_s^\top V^\top=B^\top P_s\Lambda_sP_s^\top B$.
  Furthermore, if $\lambda\ge 0$ satisfies $B^\top B\preceq \lambda D$, it also satisfies
  $B^\top B\preceq \lambda (D+VP_\ell P_\ell^\top V^\top)$, therefore
  \begin{eqnarray*}
    \lambda_{\max}(H_{P_\ell}^{-\frac12}VP_sP_s^\top V^\top H_{P_\ell}^{-\frac12})
    &\le& \lambda_{\max}(D^{-\frac12}B^\top P_s\Lambda_sP_s^\top BD^{-\frac12})\\
    &\le& \mathop{\rm tr}(D^{-\frac12}B^\top P_s\Lambda_sP_s^\top  BD^{-\frac12})\\
    &=&\mathop{\rm tr} (\Lambda_sP_s^\top B D^{-1}B^\top P_s)\\
    &=& \sum_{i=1}^{n-k}(\Lambda_s)_{ii}\|B^\top (P_s)_{\bullet,i}\|^2_{D^{-1}}\le 1+(n-k)\bar\rho\bar\beta^2.
  \end{eqnarray*}
 
\end{proof}
Note, the proof weakens $H_{P_\ell}$ to $D$, so the bound cannot be expected to be strong. Yet, it provides a good rule of thumb on which columns of $D^{-\frac12}BX^{\frac12}$ should be included, namely those with large value
$\Lambda_{ii}\|B^\top (P_s)_{\bullet,i}\|^2_{D^{-1}}$. 

In interior point methods the spectral decomposition and the size of the eigenvalues of $X$ of Theorem~\ref{T:speccondbound} strongly depend on the current iteration, in particular on the value of the barrier parameter $\mu$. Therefore it is worth to set up a new preconditioner for each new KKT system. In order to do so in a computationally efficient way, the following dynamic selection heuristic for $P_\ell$ with respect to $V$ of \eqref{eq:tildeV} tries to either pick columns of $V$ directly by including
unit vectors in $P_\ell$ or to at least take
linear combinations of few columns of $V$ in order to reduce the cost of matrix-vector multiplications and to 
preserve potential structural proprieties. So instead of forming $P_{\ell}$,
the heuristic builds $\hat V=VP_\ell$ directly by appending (linear combinations
of selected) columns of $V$ to $\hat V$. Also, it will often only
employ approximate eigenvalues $\lambda_i$ and eigenvectors $p_i$ of
the $X$ described in Theorem~\ref{T:speccondbound}. Generally, it will include those in $\hat V$ for which an estimate of $\lambda_i\|B^\top p_i\|_{D^{-1}}^2$ exceeds
a given bound $\underline\rho$. In order to reduce the number of
matrix vector multiplications, $\|B^\top p_i\|^2_{D^{-1}}$ will only
be computed for those $p_i$ with
$\big(\sum_{j=1}^n(p_i)_j^2\|(B^\top)_{\bullet,j}\|_{D^{-1}}\big)^2\ge
\underline\rho$ where the column norms of $B^\top$ are precomputed for each KKT systems.  The
implementation uses $\underline{\rho}=10$. Next the selections are
explained by going through $V$ of \eqref{eq:tildeV} step by
step for each of its three column groups $V_H$, $A^\top D_w^{\frac12}$ and
$B^\top{\mathfrak{X}}_t^{\frac12}(I-\tfrac{{\mathfrak{X}}_t^{\frac12}\mathbbm{1}_{t}({\mathfrak{X}}_t^{\frac12}\mathbbm{1}_{t})^\top}{\zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t})^{\frac12}$. Concerning the third group, it will become clear in the discussion of the semidefinite part that in practice it is advantageous to replace the square root ${\mathfrak{X}}_t^{\frac12}$ in the factorization of ${\mathfrak{X}}_t$ by a more general, possibly nonsymmetric factorization ${\mathfrak{X}}_t={\mathfrak{F}}_t{\mathfrak{F}}_t^\top$. The
matrix ${\mathfrak{F}}_t$ will have the same block structure and leads to a similar rank one correction by the transformed trace vector ${\mathfrak{F}}_t^\top\mathbbm{1}_t$,
\begin{displaymath}
  {\mathfrak{X}}_t-\tfrac{{\mathfrak{X}}_t\mathbbm{1}_{t}({\mathfrak{X}}_t\mathbbm{1}_{t})^\top}{\zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t}= {\mathfrak{F}}_t(I-\tfrac{{\mathfrak{F}}_t^{\top}\mathbbm{1}_{t}({\mathfrak{F}}_t^\top\mathbbm{1}_{t})^\top}{\zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{F}}_t{\mathfrak{F}}_t^\top\mathbbm{1}_t})^{\frac12}(I-\tfrac{{\mathfrak{F}}_t^{\top}\mathbbm{1}_{t}({\mathfrak{F}}_t^\top\mathbbm{1}_{t})^\top}{\zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{F}}_t{\mathfrak{F}}_t^\top\mathbbm{1}_t})^{\frac12}{\mathfrak{F}}_t^\top.
\end{displaymath}


\begin{algorithm}[Deterministic column selection heuristic forming $\hat V$]\label{alg:deterministicprecond}\quad\\
  \textbf{Input:} $D_H$, $V_H$, $D_y$, $A$, $D_w$, $B$, ${\mathfrak{X}}_t$, $\zeta$, $\sigma$ specifying $D$ and $V\in\mathbb{R}^{m\times n}$ of \eqref{eq:tildeV}\\
  \textbf{Output:} $\hat V\in\mathbb{R}^{m\times n'}$ for some $n'\le n$ with $\hat V= VP_{\ell}$, $P_\ell^\top P_\ell=I_{n'}$.
  \begin{compactenum}
  \item Initialize $\hat V\gets 0\in\mathbb{R}^{m \times 0}$, $\underline{\rho}:=\frac{10}m\mathop{\rm tr} D_H$
    \item Find ${\mathcal{J}}_{V_h}=\{j\colon \|(V_H)_{\bullet,j}\|_{D^{-1}}^2\ge\underline\rho\}$ and set $\hat V\gets [\hat V,(V_H)_{\bullet,{\mathcal{J}}_{V_H}}]$.
    \item Find ${\mathcal{J}}_A=\{j\colon (D_w)_{jj}\|(A^\top)_{\bullet,j}\|_{D^{-1}}^2\ge\underline\rho\}$ and set $\hat V\gets [\hat V,(A^\top D_w)_{\bullet,{\mathcal{J}}_A}]$.
    \item Compute ${\mathfrak{F}}_t^{\top}\mathbbm{1}_t$,
      $\eta=\zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{F}}_t{\mathfrak{F}}_t^\top\mathbbm{1}_t$, $B^\top{\mathfrak{F}}_t^{\top}\mathbbm{1}_t$ and for each
      conic diagonal block of ${\mathfrak{X}}_t$ call
      append\_``cone''\_columns$(\hat V)$ with corresponding parameters.
    \item Return $\hat V$. 
  \end{compactenum}
\end{algorithm}

The first group of columns $V_H\in\mathbb{R}^{m\times h_H}$ matches, in the notation of Theorem~\ref{T:speccondbound}, (a subblock of) $B^\top=V_H$ and (a diagonal block) $X=I_{h_H}$. The heuristic appends those columns $j$ to $\hat V$ that satisfy $\|(V_H)_{\bullet,j}\|^2_{D^{-1}}\ge \underline\rho$.

For the second group of columns $A^\top D_w^{\frac12}$, Theorem~\ref{T:speccondbound} applies
to $B^\top=A^\top$ and $X=D_w=\mathop{\rm Diag}(d_1,\dots,d_{h_A})$. Thus, column $j$ is appended to $\hat V$ if $d_j\|(A^\top)_{\bullet,j}\|_{D^{-1}}^2\ge \underline\rho$. 

With the comment above regarding ${\mathfrak{F}}_t$, the third column group is
formed by a term $B^\top{\mathfrak{F}}_t(I-\tfrac{{\mathfrak{F}}_t^{\top}\mathbbm{1}_{t}({\mathfrak{F}}_t^{\top}\mathbbm{1}_{t})^\top}{\zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t})^{\frac12}$ for each cutting model (we assume just one here). With respect to Theorem~\ref{T:speccondbound}, $B$ is just right
and $X$ is the  positive \mbox{(semi-)}definite matrix ${\mathfrak{X}}_t-\tfrac{{\mathfrak{X}}_t\mathbbm{1}_{t}({\mathfrak{X}}_t\mathbbm{1}_{t})^\top}{\zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t}={\mathfrak{F}}_t(I-\tfrac{{\mathfrak{F}}_t\mathbbm{1}_{t}({\mathfrak{F}}_t\mathbbm{1}_{t})^\top}{\zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{F}}_t{\mathfrak{F}}_t^\top\mathbbm{1}_t}){\mathfrak{F}}^\top$. Recall  that ${\mathfrak{X}}_t$ is a block diagonal matrix with the structure of the diagonal blocks
governed by the linearization of the perturbed complementarity conditions of the
various cones.
The overarching rank one modification by
${\mathfrak{F}}_t\mathbbm{1}_{t}$ couples the blocks within the same cutting model
and reappears in some form in each block together with $\eta=\zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t$. Observe that with $\|{\mathfrak{F}}_t^{\top}\mathbbm{1}_{t}\|^2=\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t$
\begin{displaymath}
  (I-\tfrac{{\mathfrak{F}}_t^{\top}\mathbbm{1}_{t}({\mathfrak{F}}_t^{\top}\mathbbm{1}_{t})^\top}{\eta})^{\frac12}=I-\big(1-\tfrac{\sqrt{\zeta^{-1}\sigma}}{\sqrt{\eta}}\big)\tfrac{{\mathfrak{F}}_t^{\top}\mathbbm{1}_{t}}{\|{\mathfrak{F}}_t^{\top}\mathbbm{1}_{t}\|}\tfrac{({\mathfrak{F}}_t^{\top}\mathbbm{1}_{t})^\top}{\|{\mathfrak{F}}_t^{\top}\mathbbm{1}_{t}\|}.
\end{displaymath}
For each column $p$ of $P_\ell$ computing $B^\top{\mathfrak{F}}_t(I-\tfrac{{\mathfrak{F}}_t^{\top}\mathbbm{1}_{t}({\mathfrak{F}}_t^{\top}\mathbbm{1}_{t})^\top}{\eta})^{\frac12}p$ splits into
\begin{displaymath}
 B^\top{\mathfrak{F}}_tp- \ip{{\mathfrak{F}}_t^{\top}\mathbbm{1}_{t}}{p}\tfrac{1}{\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_{t}} \big(1-\tfrac{\sqrt{\zeta^{-1}\sigma}}{\sqrt{\eta}}\big)B^\top{\mathfrak{X}}_t\mathbbm{1}_{t}.
\end{displaymath}
Thus, by keeping the support of $p$ restricted to single blocks, the proper
column computations can be kept restricted to the respective
block. This also holds for the coefficient $\ip{{\mathfrak{F}}_t^{\top}\mathbbm{1}_{t}}{p}$. The overarching vector $B^\top{\mathfrak{X}}_t\mathbbm{1}_{t}$
needs to be evaluated only once and can be added to the columns afterwards. The latter step  only requires the respective coefficients but not the vectors of $P_\ell$.
This allows to speed up the process of forming $\hat V$
considerably. Therefore, when forming the conceptional $P_\ell$ in the
heuristic, the influence of ${\mathfrak{F}}_t^{\top}\mathbbm{1}_{t}$ on eigenvalues
and eigenvectors of the blocks will mostly be considered as restricted
to each single block. Next the actual selection procedure is described
for ${\mathfrak{X}}_t$ blocks corresponding to cones $\mathbb{R}^h_+$
(Alg.~\ref{alg:nncselect}) and $\mathbb{S}^h_+$ (Alg.~\ref{alg:pscselect})
with Nesterov-Todd-scaling \cite{NesterovTodd98,ToddTohTuetuencue98}.

\begin{algorithm}[append\_$\mathbb{R}^h_+$\_columns$(\hat V)$]\label{alg:nncselect}\quad\\
  \textbf{Input:} column indices $J\in\mathbb{N}^h$ and $x\circ z^ {-1}$
  of this block in ${\mathfrak{X}}_t$, $B^\top_{\bullet,J}$,  $\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t$, $B^\top{\mathfrak{X}}_t\mathbbm{1}_t$,
  $\eta=\zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t$, $D$, threshold $\underline{\rho}$. \\
  \textbf{Output:} updated $\hat V$.
  \begin{compactenum}
    \item For each $i=1,\dots,h$ with $(\frac{x_i}{z_i}-\frac1\eta
      \frac{x_i^2}{z_i^2})\|(B^\top)_{\bullet,J(i)}\|_{D^{-1}}^2\ge\underline{\rho}$ 
      set
      $$\alpha\gets\sqrt{\tfrac{x_i}{z_i}}\tfrac{1}{\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t} \big(1-\tfrac{\sqrt{\zeta^{-1}\sigma}}{\sqrt{\eta}}\big),$$
      $$\hat b_i=\sqrt{\tfrac{x_i}{z_i}}(B^\top)_{\bullet,J(i)}- \alpha B^\top{\mathfrak{X}}_t\mathbbm{1}_t,$$
      and if $\|\hat b_i\|_{D^{-1}}^2>\underline\rho$ set 
      $\hat V\gets [\hat
      V,\hat b_i].$
    \end{compactenum}
\end{algorithm}
For Alg.~\ref{alg:nncselect} consider, within the cone specified by
$t$, a block with indices $J\in\mathbb{N}^h$ representing a nonnegative cone
$\mathbb{R}^h_+$ with primal dual pair $(x,z)$. The corresponding ``diagonal
block'' in ${\mathfrak{X}}_t$ is of the form $\mathop{\rm Diag}(x\circ z^{-1})$ and for ${\mathfrak{F}}_t$ it is $\mathop{\rm Diag}(x\circ z^{-1})^{\frac12}$. The
relevant part of the trace vector ${\mathfrak{X}}_t\mathbbm{1}_t$ reads
$\mathop{\rm Diag}(x\circ z^{-1})\mathbbm{1}=x\circ z^{-1}$. Considering the influence of
the trace vector as restricted to this block alone gives
$\mathop{\rm Diag}(x\circ z^{-1})-\frac{1}{\eta}(x\circ z^{-1})(x\circ
z^{-1})^\top$ with the correct overarching
$\eta=\zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t$. The eigenvectors to
large eigenvalues of this matrix have their most important coordinates
associated with the largest diagonal entries. The heuristic appends
the columns
\mbox{$B^\top{\mathfrak{X}}_t^{\frac12}(I-\tfrac{{\mathfrak{X}}_t^{\frac12}\mathbbm{1}_{t}({\mathfrak{X}}_t^{\frac12}\mathbbm{1}_{t})^\top}{\eta})^{\frac12}e_{J(i)}$}
to $\hat V$ for those $e_{J(i)}$ with
$(\frac{x_i}{z_i}-\frac1\eta
\frac{x_i^2}{z_i^2})\|(B^\top)_{\bullet,J(i)}\|_{D^{-1}}^2\ge\underline{\rho}$.



Note that in interior point methods $x_iz_i\approx\mu$ for barrier
parameter \mbox{$\mu\searrow 0$} and $x_i\to x_i^{opt}$, $z_i\to z_i^{opt}$.
Due to $\eta\ge\frac{x_i}{z_i}$ with $\eta$ mostly much larger,
the estimated value roughly behaves like
$\frac{x_i^2}\mu\|(B^\top)_{\bullet,J(i)}\|_{D^{-1}}^2$ and, indeed, by
experience it seems that columns are almost exclusively included only
for active $x_i^{opt}>0$ and only as $\mu$ gets small enough. When
computing high precision solutions with small $\mu$, the rank of
the preconditioner can thus be expected to match the number of
active subgradients in the cutting model. Theorem~\ref{T:speccondbound}
suggests that in iterative methods these columns have to be
included in some form in order to obtain reliable convergence
behavior.

Alg.~\ref{alg:pscselect} below deals with a positive semidefinite cone
$\mathbb{S}^h_+$ with Nesterov-Todd-scaling. For the current purposes it
suffices to know that the diagonal block of ${\mathfrak{X}}_t$ indexed by
appropriate $J\in\mathbb{N}^{h+1\choose 2}$ is of the form
$W\otimes\unkern_{\unkern s} W$ for a positive definite $W\in\mathbb{S}^h_{++}$; see
\cite{ToddTohTuetuencue98} for its efficient computation and for an
appendix of convenient rules for computing with symmetric Kronecker
products. The next result derives the eigenvectors and eigenvalues
when considering the rank one correction restricted to this block.
\begin{lemma}\label{L:TTTPSC}
  Let $W=P_W\Lambda_W P_W^\top$ with $\Lambda_W=\mathop{\rm Diag}(\lambda^W_1\ge \dots\ge\lambda^W_h>0)$ and $P_W^\top P_W=I_h$, $P_W=[w_1,\dots,w_h]$. Furthermore let $U=\Lambda_W^2-\frac1\eta(\Lambda_W^2\mathbbm{1})(\Lambda_W^2\mathbbm{1})^\top$ have eigenvalue decomposition
  $U=P_U\Lambda_UP_U^\top$ with $P_U^\top P_U=I_h$.  
  The eigenvalues of $W\otimes\unkern_{\unkern s} W-\frac1\eta\big((W\otimes\unkern_{\unkern s} W)\mathop{\rm svec}(I_n)\big)\big((W\otimes\unkern_{\unkern s} W)\mathop{\rm svec}(I_h)\big)^ \top$ are $\lambda^U_i=(\Lambda_U)_{ii}$ with eigenvectors $\sum_{j=1}^h (P_U)_{ji}\mathop{\rm svec}(w_jw_j^\top)$ for $i=1,\dots,h$ and $\lambda_i^W\lambda^W_j$ with eigenvectors $\frac1{\sqrt{2}}\mathop{\rm svec}(w_iw_j^\top+w_jw_j^\top)$ for $1\le i< j\le h$. 
\end{lemma}
\begin{proof}
  By \cite{AlizadehHaeberlyOverton98} the eigenvalues of $(W\otimes\unkern_{\unkern s} W)$ are $\lambda^W_i\lambda^W_j$
  for $1\le i\le j\le h$ with orthonormal eigenvectors
  \begin{equation}\label{eq:wij}
    \begin{array}{r@{\;\,}l@{\quad\text{for }}l}
      w_{ii}:=&\mathop{\rm svec}(w_iw_i^\top)&1\le i=j\le h,\\
      w_{ij}:=&\frac1{\sqrt2}\mathop{\rm svec}(w_iw_j^\top+w_jw_i^\top)& 1\le i<j\le h.
    \end{array}
  \end{equation}
  To see this \emph{e.g.}\  for $i<j$ observe $w_{ij}^\top w_{ij}=\frac12[2\ip{w_iw_j^\top}{w_iw_j^\top}+2\ip{w_iw_j^\top}{w_jw_i^\top}]$ and $(W\otimes\unkern_{\unkern s} W)w_{ij}=\frac1{\sqrt2}\mathop{\rm svec}(W w_iw_j^\top W+ W w_jw_i^\top W)= \lambda^W_i\lambda^W_j w_{ij}$.

  From $(W\otimes\unkern_{\unkern s} W)\mathop{\rm svec}(I_h)=\mathop{\rm svec}(W^2)=\sum_{i=1}^h(\lambda^W_i)^2w_{ii}$ one obtains
  \begin{displaymath}
    \begin{array}{r@{\;}l@{\quad\text{for }}l}
      \mathop{\rm svec}(W^2)^\top w_{ii}=&(\lambda_i^W)^2&i=1,\dots,h,\\
      \mathop{\rm svec}(W^2)^\top w_{ij}=&0&1\le i<j\le h.
    \end{array}
\end{displaymath}
The eigenvector-sorting $P_{\otimes\unkern_{\unkern s}}=[w_{11},w_{22},\dots,w_{hh},w_{12},w_{13},\dots,w_{h-1,h}]$ gives
\begin{displaymath}
  W\otimes\unkern_{\unkern s} W-\frac1\eta\mathop{\rm svec}(W^2)\mathop{\rm svec}(W^2)^\top=P_{\otimes\unkern_{\unkern s}}
  \begin{sbmatrix}
    U & 0 &\cdots & 0 \\
    0 & \lambda^W_1\lambda^W_2 & \ddots &\vdots\\
    \vdots &\ddots &\ddots& 0\\
    0 &     \cdots & 0 & \lambda^W_{h-1}\lambda^W_h
  \end{sbmatrix}P_{\otimes\unkern_{\unkern s}}^\top
\end{displaymath}
The result now follows by direct computation. 
\end{proof}
For semidefinite blocks, numerical experience indicates that it is
indeed worth to determine the eigenvalue decomposition of $U$ as in
Lemma \ref{L:TTTPSC}.  Finding the eigenvalues and eigenvectors
roughly requires the same amount of work as forming $W$ and is of no
concern.  With $J\in \mathbb{N}^{h+1\choose 2}$ denoting the column indices of
this block within $B^\top$, columns to corresponding eigenvectors are
computed by
$(\Lambda_U^{\frac12})_{ii}\cdot (B^\top)_{\bullet,J}\sum_{j=1}^h
(P_U)_{ji}w_{jj}$ or
$\sqrt{\lambda^W_i\lambda^W_j}(B^\top)_{\bullet,J}w_{ij}$. This
involves linear combinations of ${h+1\choose 2}$ columns and is
computationally expensive if the order $h$ of $W$ gets large. Indeed, when testing
all columns by their correct norms $\|(B^\top)_{\bullet,J}w_{ij}\|_{D^{-1}}^2$, too much time is spent in forming the preconditioner. Therefore  the
heuristic Alg.~\ref{alg:pscselect} first selects candidate eigenvectors to use for
$P_\ell$ via the rough estimate $\sum_{{\ensuremath{\hat\imath}}=1}^{h\choose 2}(w_{ij})_{\ensuremath{\hat\imath}}^2\|(B^\top)_{\bullet,J({\ensuremath{\hat\imath}})}\|_{D^{-1}}^2=\|w_{ij}\|_{\mathop{\rm Diag}(BD^{-1}B^\top)_J}^2$. 
For the selected eigenvectors it then
computes the precise values after the following transformation that is only seemingly involved.

In order to also account for the possibly overarching contribution of ${\mathfrak{F}}_t\mathbbm{1}_{t}$ it is advantageous to find a representation equivalent to $B^\top X^{\frac12} P_\ell$ with orthonormal
columns in $P_\ell$ as in Theorem~\ref{T:speccondbound} for a suitable
factorization of $X$ other than its square root. For this, let $V_W=P_W\Lambda_W^{\frac12}$, then $W\otimes\unkern_{\unkern s} W=(V_W\otimes\unkern_{\unkern s} V_W)(V_W^\top\otimes\unkern_{\unkern s} V_W^\top)$. Because $(V_W^\top\otimes\unkern_{\unkern s} V_W^\top)\mathop{\rm svec} I=\mathop{\rm svec}(V_W^\top V_W)=\mathop{\rm svec} \Lambda_W$ and $(V_W\otimes\unkern_{\unkern s} V_W)=(P_W\otimes\unkern_{\unkern s} P_W)(\Lambda_W^{\frac12}\otimes\unkern_{\unkern s}\Lambda_W^{\frac12})$, the notation of Lemma~\ref{L:TTTPSC} and its proof allows to rephrase the semidefinite block of ${\mathfrak{X}}_t-\tfrac{{\mathfrak{X}}_t\mathbbm{1}_{t}({\mathfrak{X}}_t\mathbbm{1}_{t})^\top}{\eta}$ as
\begin{eqnarray*}
  \lefteqn{W\otimes\unkern_{\unkern s} W-\frac{(W\otimes\unkern_{\unkern s} W)\mathop{\rm svec}(I)\mathop{\rm svec}(I)^\top(W\otimes\unkern_{\unkern s} W)}{\eta}=}\\
  &=&(V_W\otimes\unkern_{\unkern s} V_W)(I-\frac{\mathop{\rm svec}(\Lambda_W)\mathop{\rm svec}(\Lambda_W)^\top}{\eta})(V_W\otimes\unkern_{\unkern s} V_W)^\top\\
  &=&(P_W\otimes\unkern_{\unkern s} P_W)(\Lambda_W^{\frac12}\otimes\unkern_{\unkern s}\Lambda_W^{\frac12})(I-\frac{\mathop{\rm svec}(\Lambda_W)\mathop{\rm svec}(\Lambda_W)^\top}{\eta})(\Lambda_W^{\frac12}\otimes\unkern_{\unkern s}\Lambda_W^{\frac12})(P_W\otimes\unkern_{\unkern s} P_W)^\top\\
  &=&P_{\otimes\unkern_{\unkern s}}\begin{sbmatrix}
    U & 0 &\cdots & 0 \\
    0 & \lambda^W_1\lambda^W_2 & \ddots &\vdots\\
    \vdots &\ddots &\ddots& 0\\
    0 &     \cdots & 0 & \lambda^W_{h-1}\lambda^W_h
  \end{sbmatrix}P_{\otimes\unkern_{\unkern s}}^\top\\
  &=&(P_W\otimes\unkern_{\unkern s} P_W)FF^\top(P_W\otimes\unkern_{\unkern s} P_W)^\top,      
\end{eqnarray*}
where
\begin{displaymath}
F=(\Lambda_W^{\frac12}\otimes\unkern_{\unkern s}\Lambda_W^{\frac12})(I-\frac{\mathop{\rm svec}(\Lambda_W)\mathop{\rm svec}(\Lambda_W)^\top}{\eta})^{\frac12}.
\end{displaymath}
This suggests to put $V=B^\top P_{\otimes\unkern_{\unkern s}} F$ and to derive the columns
corresponding to $P_\ell$ via the singular value decomposition of $F=Q_F\Sigma_FP_F^\top$. Lemma~\ref{L:TTTPSC} provides the squared singular values $\Sigma_F^2$ and the eigenvectors give the left-singular vectors in $Q_F$. For $e_{ij}:=\frac1{\sqrt2}\mathop{\rm svec}(e_ie_j^\top+e_je_i^\top)$ there holds $\mathop{\rm svec}(\lambda_W)^\top e_{ij}=0$, so the right-singular vectors corresponding to $\sqrt{\lambda_i\lambda_j}$ read $(P_F)_{\bullet,ij}=e_{ij}$. The remaining right-singular vectors of $P_F$ may be computed via $P_F=F^\top Q_F\Sigma_F^{-1}$. In this it is sufficient and convenient to consider only the $U$ block, \mbox{i.\,e.,}\  the support restricted to the $ii$-coordinates. Denote the columns of $P_U=[u_1,\dots,u_h]$ in Lemma~\ref{L:TTTPSC}
by $u_j$ for $j=1,\dots,h$, then the corresponding right-singular vectors $u^F_j\in\mathbb{R}^h$  read for $\Lambda_W=\mathop{\rm Diag}(\lambda^W)$ and $\Lambda_U=\mathop{\rm Diag}(\lambda^U_1,\dots,\lambda^U_h)$
\begin{eqnarray}
  u^F_j&=&(I-\tfrac{\lambda^W(\lambda^W)^\top}{\eta})^{\frac12}\Lambda_W\cdot u_j\cdot \tfrac1{\sqrt{\lambda^U_j}}\label{eq:ujF}\\
       &=& \tfrac1{\sqrt{\lambda^U_j}}\bigg(\Lambda_W u_j-\frac{\sqrt{\eta}-\sqrt{\eta-\|\lambda^W\|^2}}{\sqrt{\eta}\|\lambda^W\|^2}\ip{\lambda^W}{\Lambda_Wu_j}\lambda^W\bigg).\nonumber
\end{eqnarray}
By expanding the $U$ block to the correct positions, the right-singular vector to singular value $\sqrt{\lambda^U_j}$ is $(P_F)_{\bullet,{jj}}=\mathop{\rm svec}(\mathop{\rm Diag}(u_j^F))$  for $j=1,\dots,h$.

With these preparations the selected semidefinite columns are appended to $\hat V$ as follows. First note that the semidefinite block with
coordinates $J$ of the factor ${\mathfrak{F}}_t$ is $(V_W \otimes\unkern_{\unkern s} V_W)$, which is nonsymmetric in general. The transformed trace vector ${\mathfrak{F}}_t^\top\mathbbm{1}_t$ reads $({\mathfrak{F}}_t^\top\mathbbm{1}_t)_J=(V_W^\top\otimes\unkern_{\unkern s} V_W^\top)\mathop{\rm svec} I=\mathop{\rm svec}(\Lambda_W)$.
If column $p_{ij}^F$ of $P^F$ with $1\le i\le j\le h$ is selected for $P_\ell$ by the heuristic, the column to be appended to $\hat V$ reads
\begin{displaymath}
 (B^\top)_{\bullet,J} (V_W \otimes\unkern_{\unkern s} V_W)p_{ij}^F- \ip{\mathop{\rm svec}{\Lambda_W}}{p_{ij}^F}\tfrac{1}{\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_{t}} \big(1-\tfrac{\sqrt{\zeta^{-1}\sigma}}{\sqrt{\eta}}\big)B^\top{\mathfrak{X}}_t\mathbbm{1}_{t}.
\end{displaymath}
If the selected indices satisfy $i<j$, the vector $p_{ij}^F$ is just
$e_{ij}=\frac1{\sqrt2}\mathop{\rm svec}(e_ie_j^\top+e_je_i^\top)$. By
$(V_W \otimes\unkern_{\unkern s} V_W)e_{ij}=\frac{\sqrt{\lambda_i^W\lambda_j^W}}{\sqrt2}\mathop{\rm svec}(w_iw_j^\top+w_jw_i^\top)=\sqrt{\lambda_i^W\lambda_j^W}w_{ij}$ and $\ip{\mathop{\rm svec}{\Lambda_W}}{e_{ij}}=0$ the column computation simplifies to
\begin{displaymath}
  \sqrt{\lambda_i^W\lambda_j^W}(B^\top)_{\bullet,J}w_{ij}=\sqrt{2\lambda_i^W\lambda_j^W}\big[w_i^\top{\mathop{\rm svec}}^{-1}([B^\top]_{k,J})w_j\big]_{k=1,\dots,n}.
\end{displaymath}
Typically, several mixed eigenvectors $w_{ij}$ have the same index $i$
corresponding to a large value $\lambda_i^W$, so it
quickly pays off to precompute $w_i^\top\mathop{\rm svec}^{-1}([B^\top]_{k,J})$
and to use these $h$-vectors for each $w_j$. For ease of presentation
this implementational detail is not described in
Alg.~\ref{alg:pscselect}. Also, this is
not helpful for the non-mixed vectors
$p_{jj}^F=\mathop{\rm svec}(\mathop{\rm Diag}(u_j^F))$, because
\begin{displaymath}
  (V_W\otimes\unkern_{\unkern s} V_W)p_{jj}^F=(V_W\otimes\unkern_{\unkern s} V_W)\sum_{i=1}^h (u_j^F)_i\mathop{\rm svec}{e_ie_i^\top}=\sum_{i=1}^h (u_j^F)_i\lambda^W_iw_{ii} 
\end{displaymath}
consists of a linear combination over all $w_{ii}$.
Fortunately only few of the non-mixed vectors are among those selected
for preconditioning. The transformed trace vector coefficient for
$p_{jj}^F$ evaluates to $\ip{\Lambda_W}{\mathop{\rm Diag}(u_j^F)}=\ip{\lambda^W}{u_j^F}$.
With this, the algorithm for appending semidefinite columns reads as follows.
\begin{algorithm}[append\_$\mathbb{S}^h_+$\_columns$(\hat V)$]\label{alg:pscselect}\quad\\
  \textbf{Input:}  column indices $J\in\mathbb{N}^{h+1\choose 2}$
  and Nesterov-Todd scaling matrix $W\succ 0$ of this block in
  ${\mathfrak{X}}_t$, $B^\top_{\bullet,J}$,  $\mathbbm{1}_t{\mathfrak{X}}_t\mathbbm{1}_t$, $B^\top{\mathfrak{X}}_t\mathbbm{1}_t$,
  $\eta=\zeta^{-1}\sigma+\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t$, $D$, threshold $\underline{\rho}$ \\
  \textbf{Output:} updated $\hat V$.
  \begin{compactenum}
  \item Compute norms $\|(B^\top)_{\bullet,J(i)}\|_{D^{-1}}$, set $\hat\rho=\underline{\rho}/\max_{i=1,\dots,{h+1\choose 2}}\|(B^\top)_{\bullet,J(i)}\|_{D^{-1}}^2$,\\
    compute eigenvalue decomposition $W=P_W\Lambda_W
    P_W^\top$, $\Lambda_W=\mathop{\rm Diag}(\lambda^W)$ with $\lambda^W_{1}\ge \dots\ge\lambda^W_{h}$,
    let $w_{ij}$ be defined by \eqref{eq:wij}.
  \item If $(\lambda^W_{1})^2<\hat\rho$ do nothing and return $\hat V$.  
  \item Compute
    $U=\Lambda_W^2-\frac1\eta(\lambda^W)(\lambda^W)^\top$,
    eigenvalue decomposition $U=P_U\Lambda_U P_U^\top$,\\
    with $P=[u_1,\dots,u_h]$.
  \item For each ${\ensuremath{\hat\imath}}=1,\dots,h$ with $(\Lambda_U)_{{\ensuremath{\hat\imath}}\hati}\ge\hat\rho$ do:\\
    Compute $\hat w_{{\ensuremath{\hat\imath}}\hati}=\sum_{i=1}^h (u_{\ensuremath{\hat\imath}})_iw_{ii}$ ($\in\mathbb{R}^{h+1\choose 2}$).\\
    If $(\Lambda_U)_{{\ensuremath{\hat\imath}}\hati}\sum_{j=1}^{h+1\choose 2}(
    \hat w_{{\ensuremath{\hat\imath}}\hati})_j^2\|(B^\top)_{\bullet,J(j)}\|_{D^{-1}}^2\ge\underline{\rho}$ then:
    \begin{compactenum}[(a)]
    \item Compute $u_{{\ensuremath{\hat\imath}}}^F$ according to \eqref{eq:ujF} and set 
      \begin{displaymath}
        \alpha\gets\ip{\lambda^W}{u_{\ensuremath{\hat\imath}}^F} \tfrac{1}{\mathbbm{1}_t^\top{\mathfrak{X}}_t\mathbbm{1}_t} \big(1-\tfrac{\sqrt{\zeta^{-1}\sigma}}{\sqrt{\eta}}\big),
      \end{displaymath}
      \begin{displaymath}
        \hat b_{{\ensuremath{\hat\imath}}\hati}=(B^\top)_{\bullet,J}\sum_{i=0}^h (u_{\ensuremath{\hat\imath}}^F)_i\lambda^W_iw_{ii}-\alpha B^\top {\mathfrak{X}}_t\mathbbm{1}_t.
      \end{displaymath}
    \item If $\|\hat b_{{\ensuremath{\hat\imath}}\hati}\|_{D^{-1}}^2\ge\underline\rho$ set $\hat V\gets [\hat
      V,\hat b_{{\ensuremath{\hat\imath}}\hati}]$.
    \end{compactenum}
  \item For each $1\le{\ensuremath{\hat\imath}}<{\ensuremath{\hat\jmath}}\le h$ with $\lambda^W_{{\ensuremath{\hat\imath}}}\lambda^W_{{\ensuremath{\hat\jmath}}}>\hat\rho$ do:\\
    If $\sqrt{\lambda^W_{{\ensuremath{\hat\imath}}}\lambda^W_{{\ensuremath{\hat\jmath}}}}\sum_{j=1}^{h+1\choose 2}(
    \hat w_{{\ensuremath{\hat\imath}}{\ensuremath{\hat\jmath}}})_j^2\|(B^\top)_{\bullet,J(j)}\|_{D^{-1}}^2\ge\underline{\rho}$ set \\    
    $$\hat b_{{\ensuremath{\hat\imath}}{\ensuremath{\hat\jmath}}}=\sqrt{\lambda^W_{{\ensuremath{\hat\imath}}}\lambda^W_{{\ensuremath{\hat\jmath}}}}(B^\top)_{\bullet,J}w_{{\ensuremath{\hat\imath}}{\ensuremath{\hat\jmath}}}$$
    and if $\|\hat b_{{\ensuremath{\hat\imath}}{\ensuremath{\hat\jmath}}}\|_{D^{-1}}^2\ge\underline{\rho}$ set $\hat V\gets [\hat
      V,\hat b_{{\ensuremath{\hat\imath}}{\ensuremath{\hat\jmath}}}]$.
      \item Return $\hat V$.
    \end{compactenum}
\end{algorithm}
As for the
linear case it can be argued that for small barrier parameter $\mu$
the number of selected columns corresponds at least to the order of the
active submatrix in the cutting model. Thus if $\hat h\le h$ eigenvalues of $X\in\mathbb{S}^h_+$
converge to positive values in the optimum, the heuristic will end up
with selecting at least ${\hat h+1\choose 2}$ columns once $\mu$ gets small.


For second order cones $\soc{h}$ the structural properties of the arrow operator and
the Nesterov-Todd-direction allow to restrict considerations to just
two directions per cone for preconditioning, but as the computational
experiments do not involve second order cones this will not be discussed here.















 


















\section{Numerical Experiments}\label{S:experiments}

The purpose of the numerical experiments is to explore and compare the
behavior and performance of the pure and preconditioned iterative variants to
the original direct solver on KKT instances that arise in the course
of solving large scale instances by the conic bundle method.

It has to be emphasized that the experiments are by no means designed
and intended to investigate the efficiency of the conic bundle method
with internal iterative solver. Indeed, many aspects of the
ConicBundle code \cite{ConicBundle2021} such as the cutting model
selection routines, the path following predictor-corrector approach
and the internal termination criteria have been tuned to work
reasonably well with the direct solver. As the theory suggests and the
results support, the performance of iterative methods depends more on
the size of the active set than on the size of the model. Thus
somewhat larger models might be better in connection with iterative
solvers. Also, the predictor-corrector approach is particularly
efficient if setting up the KKT system is expensive. For iterative
methods with deterministic preconditioning this hinges on the cost of
forming the preconditioner which gets expensive once the barrier
parameter gets small. Furthermore iterative methods might actually
profit from staying in a rather narrow neighborhood of the central
path. Therefore many implementational decisions need to be reevaluated
for iterative solvers. This is out of scope for this paper. Hence, the
experiments only aim to highlight the relative performance of the
solvers on sequences of KKT systems that currently arise in
ConicBundle.

The experiments will report on the performance for three different
instances: the first, denoted by MC, is a classical semidefinite
relaxation of Max-Cut on a graph with 20000 nodes as described in
\cite{GoemansWilliamson95,HelmbergRendl2000}, the second, BIS, is a
semidefinite Minimum-Bisection relaxation improved by dynamic
separation of odd cycle cutting planes on the support of the Boeing
instance KKT\_traj33 giving a graph on 20006 nodes explained in
\cite{Helmberg2004}, and the third, MMBIS, refers to a
min-max-bisection problem problem shifting the edge weights so as to
minimize a restricted maximum cut on a graph of 12600 nodes. All three
have a single semidefinite cutting model which consists of a
semidefinite cone with up to one nonnegative variable, so the model
cone ${\mathcal{S}}^t_+$ of \eqref{eq:St} typically has $t=(1,[],[h])$ for
some $h\in\mathbb{N}$.  In the Max-Cut instance the design variables are
unconstrained, in the Bisection instance the design variables
corresponding to the cutting planes are sign constrained ($D_y$ is
needed) and in the min-max-bisection problem some design variables
have bounds and there are linear equality and inequality constraints
($D_y$, $D_w$ and $A$ appear). Throughout, the proximal term is a
multiple of the identity for a dynamic weight, \mbox{i.\,e.,}\  $H_k=u_kI$ with
$u_k>0$ controlled as in \cite{HelmbergKiwiel2002}.

In each case ConicBundle is run with default settings for the
internal constrained QP solver with direct
KKT solver for the bundle subproblems. Whenever a new KKT-system arises, it is
solved consecutively but independently on the same machine by
\begin{compactitem}
\item (DS) the original direct solver,
\item (IT) MINRES without preconditioning (the implementation follows \cite{ElmanSilvesterWathen2006}),
\item (RP) MINRES with randomized preconditioning
  (Alg.~\ref{alg:precond} with  Alg.~\ref{alg:randprecond}),
\item (DP) MINRES with deterministic preconditioning (Alg.~\ref{alg:precond} with
Alg.~\ref{alg:deterministicprecond}).
\end{compactitem}
Only the results of the direct solver are then used to continue the
algorithm. Note, for nonsmooth optimization problems tiny deviations
in the solution of the subproblem may lead to huge differences in the
subsequent path of the algorithm. Therefore running the bundle method with
different solvers would quickly lead to incomparable KKT systems. 

The details of the direct solver DS are of little relevance at this
point. Suffice it to say that its main work consists in Schur
complementing the $H$ and $\zeta^{-1}\sigma$ blocks of the KKT system \eqref{eq:fullKKT} into the joined $\mathop{\rm Diag}(D_w^{-1},{\mathfrak{X}}_t^{-1})$
block and factorizing this. In the Max-Cut setting (no $D_y$), the
$H$ block is constant throughout each bundle subproblem. In this case the
Schur complement is precomputed once for each bundle subproblem ---
thus for several KKT systems --- and this make this approach extremely
efficient as long as the order $h$ of the semidefinite model is
small. Precomputation is no longer possible if $D_y$ is needed which
is the case in the two other instances.  Finally, if $A$ is also
present, the system to be factorized in every iteration gets
significantly larger. These differences motivated the choice of the
instances and explain part of the strong differences in the performance
of the solvers.

For Max-Cut and Bisection the iterative solver could exploit the positive
definiteness of the system by employing conjugate gradients instead of
MINRES. The min-max-bisection problem comprises equality constraints
in $A$, so the system is no longer positive definite and conjugate
gradients are not applicable. Employing MINRES for all three
facilitates the comparison, in particular as MINRES seemed to perform
numerically better on the other instances as well. MINRES computes the residual norm with respect to the inverse of the preconditioner and
the implementation uses this norm for termination. To safeguard against
effects due to the changes in this norm, the relative precision requirement
is multiplied, in the notation of Alg.~\ref{alg:precond}, by the
factor $\big(\sqrt[m]{\prod_{i=1}^{\hat k}(1+\hat\lambda_i)^{-1}}\cdot\min_i (D^{-1})_i\big)^{\frac12}$.

The results on the three instances will be presented in eight plots per
instance. The first four compare all four solvers, the last four plots
are devoted to information that is only relevant for iterative
solvers, so DS will not appear in these. 
\begin{compactenum}
\item Plot ``time per subproblem (seconds)'' gives for each of the
  four methods a box plot on the seconds (in logarithmic scale)
  required to solve the subproblems.  For each subproblem this is the
  sum of the time required for initializing/forming and solving all
  KKT systems of this subproblem.  This is needed, because in the case
  of Max-Cut instance MC, the direct solver DS forms the Schur complement of the
  $H$-block only once per subproblem and this is also accounted for
  here.
\item Plot ``subproblem time (seconds) per iteration'' displays
  the same cumulative time per subproblem in seconds (in logarithmic
  scale) for each successive iteration so that the development
  in solution time is aligned to the progress of the bundle method. 
\item Plot ``time per subproblem vs.\ bundle size'' serves to highlight
  the dependence of the solution time on the size of the cutting model
  (number of rows of $B$). For this the subproblems are grouped in the bundle size ranges $(0,50]$, $(50,500]$, $(500,1500]$, $(1500,\infty]$. Instead of these numbers, the actual maxima of the groups are shown in the bottom line of the plot.   
\item Plot ``time per subproblem vs.\ last $\mu$'' illustrates the
  dependence of the solution time on the last barrier parameter $\mu$
  for which the subproblem has to be solved. Roughly this corresponds
  to the precision required for the subproblem. Results are presented
  for a subdivision of the subproblems into four groups of equal
  cardinality (up to integer division) sorted according to the $\mu$
  value of their respective last KKT system. The minimum $\mu$
  of each group is given in the bottom line of the plot.
\item Plot ``time per KKT system (seconds)'' compares exclusively the
  iterative methods on the KKT systems belonging to the four different
  ranges of the barrier parameter $\mu$ as collected over all
  subproblems. The first three box plots give
  the box plot statistics on the seconds (in logarithmic scale)
  spent in solving KKT systems for barrier parameter values
  $\mu\ge 100$, the next three for $100>\mu\ge 1$, etc. Note, DS
  would require the same time for all KKT systems of the same
  subproblem, because its solution time does not depend on $\mu$
  or the associated required relative precision which was set to
  $\min\{0.01\cdot \mu,10^{-6}\}$. 
\item Plot ``matrix vector multiplications per KKT system'' shows box plots on
  the number of matrix-vector multiplications (in logarithmic scale)
  needed by MINRES, again subdivided into the same ranges of barrier
  parameter values.
\item Plot ``KKT system condition number estimate'' presents the box plot
  statistics of an estimate of the condition number (in logarithmic
  scale) for the same ranges of the barrier parameter. The estimate
  is obtained by a limited number of Lanczos iterations on the respective
  (non-)preconditioned system of the $H$ block; a possibly remaining
  equality part of $A$ is ignored in this. Computation times for the
  condition number are not included in the time measurement listed above.  
\item Plot ``preconditioning columns per KKT system'' gives the box
  plot statistics of the number of columns $\hat k$ in
  Alg.~\ref{alg:precond} for RP and DP for the usual ranges of the
  barrier parameter.
\end{compactenum}
In all box plots, the width of the boxes indicates the relative size of
the number of instances in the group, the horizontal lines of the boxes give the
values of the upper quartile, the median and the lower quartile. The
upper whisker shows the largest value below upper
quartile$ + 1.5\cdot$IQR, where IQR$ = $(upper quartile $-$ lower
quartile) is the interquartile range. The lower whisker displays the
smallest value above lower quartile$ - 1.5\cdot$IQR. The stars show
maximum and minimum value.

Computation times refer to a virtualized compute server of 40 Intel Xeon Processor (Cascadelake) cores with 600 GB RAM under Ubuntu 18.04. This virtual machine
is hosted on hardware consisting of two processors Intel(R) Xeon(R)
Gold 6240R CPU with 2.40GHz with 24 cores and 768 GB RAM. The code,
however, is purely sequential and does not exploit any parallel
computation possibilities. 

\subsection{Max-Cut (Instance MC, Figure \ref{fig:MC})}

\begin{figure}[htbp]
\centerline{
\begin{tabular}{@{}cc@{}}
\includegraphics[width=.48\textwidth]{MC_subp.pdf}
&
\includegraphics[width=.48\textwidth]{MC_iter.pdf}
\\[1ex]
\includegraphics[width=.48\textwidth]{MC_bund.pdf}
&
\includegraphics[width=.48\textwidth]{MC_prec.pdf}
\\[1ex]
\includegraphics[width=.48\textwidth]{MC_time.pdf}
&
\includegraphics[width=.48\textwidth]{MC_mult.pdf}
\\[1ex]
\includegraphics[width=.48\textwidth]{MC_cond.pdf}
&
\includegraphics[width=.48\textwidth]{MC_rank.pdf}
\end{tabular}
}\vspace*{-.5ex}
\caption{\label{fig:MC} Instance MC, 338 subproblems, 2832 KKT systems. In order to highlight the dependence on the barrier parameter $\mu$ with corresponding precision requirements the results for KKT systems are grouped into $\mu$-value ranges $(\infty,100]$, $(100,1]$, $(1,0.01]$, $(0.01,0)$.}
\end{figure}

The graph was randomly generated (\cite{rudy}, call \mbox{\texttt{rudy
    -rnd\_graph 20000 1 1}} for 20000 nodes, edge density one percent,
seed value 1). The semidefinite relaxation gives rise to an
unconstrained problem with $20000$ variables. Each variable influences
one of the diagonal elements of the Laplace matrix of the graph with
cost one and the task is to minimize the maximum eigenvalue of the
Laplacian times the number of nodes, see \cite{HelmbergRendl2000} for the
general problem description.

For graphs of this type but smaller size like 5000 or 10000 nodes the
direct solver DS still seemed to perform better, so rather large sizes
are needed to see some advantage of iterative methods. Other than that
the relative behavior of the solvers was similar also for the smaller
sizes. The jaggies within subproblem time in the second plot are due
to the reduction of the model to its active part after each descent
step while the model typically increases in size during null steps.
During the very first iterations the bundle is tiny and DS is the best
choice.  Once the bundle size increases sufficiently, the iterative
methods dominate. Over time, as precision requirements get higher and
the choice of the bundle subspace converges, the advantage of
iterative methods decreases. In the final phase of high precision the
direct solver may well be more attractive again.

The plots also show that for this instance (and presumably for most
instances of this random type) the performance of IT (MINRES without
preconditioning) is almost as good as DP (deterministic
preconditioning) while RP (randomized preconditioning) is not
competitive. Note that the condition number does not grow excessively
for IT in this instance. Deterministic preconditioning succeeds
in keeping the condition number almost exactly at the intended
value 10. For smaller values of $\mu$, so for higher precision
requirements, DP requires distinctly fewer matrix-vector
multiplications, but it then also selects a large number of columns. In
comparison to no preconditioning DP helps to improve stability but
does not lead to significantly better computation times except maybe
for the very last phase of the algorithm with high precision
requirements.


\subsection{Minimum Bisection (Instance BIS, Figure \ref{fig:BIS})}

\begin{figure}[htbp]
\centerline{
\begin{tabular}{@{}cc@{}}
\includegraphics[width=.48\textwidth]{BIS_subp.pdf}
&
\includegraphics[width=.48\textwidth]{BIS_iter.pdf}
\\[1ex]
\includegraphics[width=.48\textwidth]{BIS_bund.pdf}
&
\includegraphics[width=.48\textwidth]{BIS_prec.pdf}
\\[1ex]
\includegraphics[width=.48\textwidth]{BIS_time.pdf}
&
\includegraphics[width=.48\textwidth]{BIS_mult.pdf}
\\[1ex]
\includegraphics[width=.48\textwidth]{BIS_cond.pdf}
&
\includegraphics[width=.48\textwidth]{BIS_rank.pdf}
\end{tabular}
}\vspace*{-.5ex}
\caption{\label{fig:BIS} Instance BIS, 180 subproblems, 5380 KKT systems. In order to highlight the dependence on the barrier parameter $\mu$ with corresponding precision requirements the results for KKT systems are grouped into $\mu$-value ranges $(\infty,100]$, $(100,1]$, $(1,0.01]$, $(0.01,0)$.}
\end{figure}

The semidefinite relaxation of minimum bisection is similar in nature
to max-cut, but in addition to the single diagonal elements there is a
variable with coefficient matrix of all ones. Furthermore, variables
with sparse coefficient matrices corresponding to odd cycles in the
underlying graph are added dynamically in rounds, see
\cite{Helmberg2004} for the general framework and also for the origin
of the instance KKT\_traj33 with 20006 nodes and roughly 260000 edges.

Again, after the very first iterations the iterative methods turn out
to perform distinctly better in the initial phase of the
algorithm. Iterative methods get less attractive as precision
requirements increase.  The model size is often rather small (a bit
larger than the active set of about 150 columns) which is favorable for
DS. Indeed, additional output information of the log file indicates
that the performance of DS drops off whenever the cutting model is
significantly larger than that.

While for this instance RP is better than IT, the advantage of DP over
the other iterative variants is quite apparent and its superiority
also increases with precision requirements and smaller $\mu$. In fact,
for DP the condition number and the number of matrix-vector multiplications
decrease again for smaller $\mu$. Possible causes might be that the
active set is easier to identify correctly.  Due to the reduction in
matrix-vector multiplications, computation time does not increase for
DP in spite of a growing number of columns in the preconditioner.


\subsection{A Min-Max-Bisection Problem (Instance MMBIS, Figure \ref{fig:MMBIS})}

This problem arose in the context of an unpublished
attempt\footnote{together with B. Filipecki (TU Chemnitz), S. Heyder
  (TU Ilmenau), Th.\ Hotz (TU Ilmenau) within BMBF-project grant
  05M18OCA.} to optimize vaccination rates for five population groups
$N_1\dot\cup\dots\dot\cup N_5=N$ in a virtual town of $n=|N|$
inhabitants. Briefly, within the town $k$ anonymous people are assumed to be
infectious. There is vaccine for at most $\underline n$ people. The
aim is to reduce the spreading rate of the disease by vaccinating each
person with the respective group's probability. The task of
determining these vaccination rates motivated the following ad hoc model
which would be hard to justify rigorously. In a graph
$G=(N,E)$ each edge $ij=\{i,j\}\in E$ with $i\in N_{\ensuremath{\hat\imath}}$,
$j\in N_{\ensuremath{\hat\jmath}}$ has a weight $\hat w_{{\ensuremath{\hat\imath}}{\ensuremath{\hat\jmath}}}$ representing the infectiousness
of the typical contact for these two persons of the respective groups.
It will be convenient
to define the weighted Laplacians $L_{{\ensuremath{\hat\imath}}{\ensuremath{\hat\jmath}}}=\hat w_{{\ensuremath{\hat\imath}}{\ensuremath{\hat\jmath}}}\sum_{ij\in E, i\in N_{\ensuremath{\hat\imath}},j\in N_{\ensuremath{\hat\jmath}}}(e_i-e_j)(e_i-e_j)^\top$. In
this simplified approach, vaccination rates $v_{\ensuremath{\hat\imath}},v_{\ensuremath{\hat\jmath}}$ of the node
groups reduce a nominal infectiousness $\hat w_{{\ensuremath{\hat\imath}}{\ensuremath{\hat\jmath}}}$ between
these groups by the factor
$y_{{\ensuremath{\hat\imath}}{\ensuremath{\hat\jmath}}}\ge\max\{0,1-v_{\ensuremath{\hat\imath}}-v_{\ensuremath{\hat\jmath}}\}$.
The spreading probability to be minimized is considered proportional to the restricted
max-cut value
\begin{eqnarray*}
  &\displaystyle\max_{{\mathcal{I}}\subset N,|{\mathcal{I}}|=k}\sum y_{{\ensuremath{\hat\imath}}{\ensuremath{\hat\jmath}}}\hat w_{{\ensuremath{\hat\imath}}{\ensuremath{\hat\jmath}}}\cdot|\{ij\in E\colon i\in{\mathcal{I}}\cap N_{\ensuremath{\hat\imath}},j\in N_{\ensuremath{\hat\jmath}}\setminus{\mathcal{I}}\}|=&\\
  &=\displaystyle\max_{x\in\{-1,1\}^n\atop (\mathbbm{1}^\top x)^2=(n-2k)^2}\frac14 x^\top\bigg(\sum_{{\ensuremath{\hat\imath}}\le{\ensuremath{\hat\jmath}}}y_{{\ensuremath{\hat\imath}}{\ensuremath{\hat\jmath}}}L_{{\ensuremath{\hat\imath}}{\ensuremath{\hat\jmath}}}\bigg)x,& 
\end{eqnarray*}
For determining the vaccination rates the combinatorial problem is replaced by the usual (dual) semidefinite relaxation
\begin{displaymath}
  \begin{array}{rl}
    \text{minimize} &\frac{n}4\lambda_{\max}(\sum_{{\ensuremath{\hat\imath}}\le{\ensuremath{\hat\jmath}}}y_{{\ensuremath{\hat\imath}}{\ensuremath{\hat\jmath}}}L_{{\ensuremath{\hat\imath}}{\ensuremath{\hat\jmath}}}-\mathop{\rm Diag}(d)-u\mathbbm{1}\one^\top)+\mathbbm{1}^\top d+(n-2k)^2u\\
    \text{subject to} & y_{{\ensuremath{\hat\imath}}{\ensuremath{\hat\jmath}}}\ge 1-v_{\ensuremath{\hat\imath}}-v_{\ensuremath{\hat\jmath}},\qquad{\ensuremath{\hat\imath}}\le{\ensuremath{\hat\jmath}},\\
                    &\sum_{{\ensuremath{\hat\imath}}} |N_{\ensuremath{\hat\imath}}|v_{\ensuremath{\hat\imath}}=\underline n,\\
    & d\in\mathbb{R}^n,u\in\mathbb{R},y\ge 0, v\ge 0.
  \end{array}
\end{displaymath}

\begin{figure}[htbp]
\centerline{
\begin{tabular}{@{}cc@{}}
\includegraphics[width=.48\textwidth]{MMBIS_subp.pdf}
&
\includegraphics[width=.48\textwidth]{MMBIS_iter.pdf}
\\[1ex]
\includegraphics[width=.48\textwidth]{MMBIS_bund.pdf}
&
\includegraphics[width=.48\textwidth]{MMBIS_prec.pdf}
\\[1ex]
\includegraphics[width=.48\textwidth]{MMBIS_time.pdf}
&
\includegraphics[width=.48\textwidth]{MMBIS_mult.pdf}
\\[1ex]
\includegraphics[width=.48\textwidth]{MMBIS_cond.pdf}
&
\includegraphics[width=.48\textwidth]{MMBIS_rank.pdf}
\end{tabular}
}\vspace*{-.5ex}
\caption{\label{fig:MMBIS} Instance MMBIS, 28 subproblems, 642 KKT systems. In order to highlight the dependence on the barrier parameter $\mu$ with corresponding precision requirements the results for KKT systems are grouped into $\mu$-value ranges $(\infty,100]$, $(100,1]$, $(1,0.01]$, $(0.01,0)$.}
\end{figure}

In this case the resulting KKT system also has an equality and several inequality constraints in the block $A$. Preconditioning results are presented for the KKT systems of an instance with $n=12600$ inhabitants splitting into groups of sizes $5770, 6000, 600, 30, 200$, with $k=126$ infectious persons and $\underline n=1260$ available vaccinations.

In the actual computations the bundle size grows surprisingly
fast. This not only entails enormous memory requirements but also
excessive computation times for DS; indeed, computations of DS may
exceed those of DP by a factor of 70. In consequence comparative
results can only be reported for a very limited number of subproblem
evaluations. In particular, the precision requirements remain rather
moderate throughout these iterations. Still, the same initial behavior
can be observed as for the previous two instances. For very small
bundle sizes DS is best. Once the bundle size grows, the iterative
methods take over. Among the iterative solvers RP is better than IT,
but DP is the method of choice. It succeeds in tightly controlling the
condition number by selecting rather few columns. With this DP
requires the fewest matrix vector multiplications which seems to pay
of quickly on this instance.

\section{Conclusions}\label{S:concl}

In search for efficient low rank preconditioning techniques for the
iterative solution of the internal KKT system of the quadratic bundle
subproblem two subspace selection heuristics --- a randomized and a
deterministic variant --- were proposed.  For the randomized approach
the results are ambivalent in theory and in practice; obtaining a good
subspace this way seems to be difficult and the cost of exploratory
matrix-vector multiplications quickly dominates. In contrast, the
deterministic subspace selection approach allows to control the
condition number (and with it the number of matrix vector
multiplications) at a desired level without the need to tune any
parameters in theory as well as on the test instances. On these
instances, for low precision requirements (large barrier parameter)
the selected subspace is neglectably small. For high precision
requirements (small barrier parameter) the subspace grows to the
active model subspace. If the bundle size is close to this active
dimension, the work in forming the preconditioner may be comparable to
forming the Schur complement for the direct solver. Still, for large
scale instances the deterministically preconditioned iterative
approach seems to be preferable.

Conceivably it is possible to profit in ConicBundle from the
advantages of the deterministic iterative and the direct solver by
switching dynamically between both. The current experiments relied on
a predictor-corrector approach that was tuned for the direct
solver. In view of the properties of the iterative approach it may
well be worth to devise a different path following strategy for the
iterative approach, in particular for the initial phase of the
interior point method when the barrier parameter is still
comparatively large and the work invested in forming the
preconditioner is still neglectable.  Similar ideas should be
applicable to interior point solvers for solving convex quadratic
problems with low rank structure.

\textbf{Acknowledgments} I have profited a lot from discussions with
many colleagues, in part years back. In particular I have to thank
K.-C.~Toh as well as my colleagues O.~Ernst, R.~Herzog, A.~Pichler and
M.~Stoll in Chemnitz. Much of the preparatory restructuring of
ConicBundle was done during my sabbatical at the University of
Klagenfurt, thank you to F.~Rendl and A.~Wiegele for making this
possible.  The support of research grants 05M18OCA of the German
Federal Ministry of Education and Research and the DFG CRC 1410 is
gratefully acknowledged.


