%% \CharacterTable
%%  {Upper-case    \A\B\C\D\E\F\G\H\I\J\K\L\M\N\O\P\Q\R\S\T\U\V\W\X\Y\Z
%%   Lower-case    \a\b\c\d\e\f\g\h\i\j\k\l\m\n\o\p\q\r\s\t\u\v\w\x\y\z
%%   Digits        \0\1\2\3\4\5\6\7\8\9
%%   Exclamation   \!     Double quote  \"     Hash (number) \#
%%   Dollar        \$     Percent       \%     Ampersand     \&
%%   Acute accent  \'     Left paren    \(     Right paren   \)
%%   Asterisk      \*     Plus          \+     Comma         \,
%%   Minus         \-     Point         \.     Solidus       \/
%%   Colon         \:     Semicolon     \;     Less than     \<
%%   Equals        \=     Greater than  \>     Question mark \?
%%   Commercial at \@     Left bracket  \[     Backslash     \\
%%   Right bracket \]     Circumflex    \^     Underscore    \_
%%   Grave accent  \`     Left brace    \{     Vertical bar  \|
%%   Right brace   \}     Tilde         \~}

\documentclass[final,oneside,runningheads,10pt]{svjour2}
\usepackage{mathptmx,amsmath,amssymb,xspace,citesort}
\usepackage[ruled,section]{algorithm}
\usepackage{algorithmic}
\usepackage[pdftex]{graphicx}
%\DeclareGraphicsExtensions{.eps}

% Mathematical Programming macros

\smartqed  % flush right qed marks, e.g. at end of proof
\journalname{Mathematical Programming}

% Macros and definitions
\def\labelitemi{--} \def\ie{\textit{i.e.\xspace}}
\def\eg{\textit{e.g.\xspace}}
\def\etc{\textit{etc.}}
\def\st{\ensuremath{\mathrm{\;\;s.t.\;\;}}}

\def\BHplus  {\ensuremath{\mathrm{BH}^{+}}}
\def\Btwo    {\ensuremath{\mathrm{B}_2}}
\def\BeO     {\ensuremath{\mathrm{BeO}}}
\def\Ctwo    {\ensuremath{\mathrm{C}_2}}
\def\Ctwoplus{\ensuremath{\mathrm{C}_{2}^{+}}}
\def\Litwo   {\ensuremath{\mathrm{Li}_2}}
\def\LiF     {\ensuremath{\mathrm{LiF}}}
\def\NaH     {\ensuremath{\mathrm{NaH}}}

\def\half{\ensuremath{\frac{1}{2}}}
\def\svec{\ensuremath{\mathrm{svec}}}
\def\LAMBDA{{\sc Lambda}\xspace}
\def\reals{\ensuremath{\mathbb{R}}}
\def\E{\ensuremath{\mathbb{E}}}

\newcommand{\junk}[1]{}
\newcommand{\secref}[1]{\S\ref{#1}}
\newcommand{\tabref}[1]{Table~\ref{#1}}
\newcommand{\figref}[1]{Figure~\ref{#1}}
\newcommand{\algoref}[1]{Algorithm~\ref{#1}}

\newcommand{\supth}[1]{\ensuremath{#1^{\mathrm{th}}}}
\newcommand{\norm}[2]{\ensuremath{\left\|#1\right\|_{#2}}}
\newcommand{\maps}[2]{\ensuremath{: #1 \rightarrow #2}}
\newcommand{\mapselt}[2]{\ensuremath{: #1 \mapsto #2}}
\newcommand{\ip}[2]{\ensuremath{\left\langle #1,\, #2\right\rangle}}
\newcommand{\subdiff}[2]{\ensuremath{\partial #1(#2)}}
\newcommand{\epssubdiff}[3]{\ensuremath{\partial_{#3} #1(#2)}}
\newcommand{\tr}[1]{\ensuremath{\mathrm{tr}(#1)}}
\newcommand{\adj}[1]{\ensuremath{{#1}^{T}}}
\newcommand{\setof}[1]{\ensuremath{\left\{ #1 \right\} }}
\newcommand{\sym}[1]{\ensuremath{\mathbb{S}^{#1}}}
\newcommand{\Span}[1]{\ensuremath{\textrm{Span}\setof{#1}}}
\newcommand{\abs}[1]{\ensuremath{\left\vert #1\right\vert}}

\begin{document}
\title{Solving Large-Scale Semidefinite Programs in
Parallel\thanks{This work was supported in part by NSF grants
DMS-0215373 and DMS-0238008.}}

\author{Madhu~V.~Nayakkankuppam}

\authorrunning{M.~V.~Nayakkankuppam}

\titlerunning{Solving Large-Scale SDP's in Parallel}

\institute{
Computational Finance Group, 
Bloomberg LP, 
731 Lexington Avenue, 
New York, NY 10022, 
USA.  
Email: Madhu.Nayakkankuppam@gmail.com.
}

\dedication{Dedicated to the memory of Jos Sturm.}

% \date{Received March 06 2005; revised March 2006;
% accepted July 20, 2006.}

\maketitle

\begin{abstract}
We describe an approach to the parallel and distributed
solution of large-scale, block structured semidefinite
programs using the spectral bundle method.  Various
elements of this approach (such as data distribution, an
implicitly restarted Lanczos method tailored to handle
block diagonal structure, a mixed polyhedral-semidefinite
subdifferential model, and other aspects related to
parallelism) are combined in an implementation called
\LAMBDA, which delivers faster solution times than
previously possible, and acceptable parallel scalability
on sufficiently large problems.
\end{abstract}

\keywords{Semidefinite programming, eigenvalue
optimization, subgradient bundle methods, Lanczos method,
parallel computing.}

\subclass{
90C22, % Semidefinite Programming
90C06, % Large-scale problems
65F15, % Eigenvalues, eigenvectors
65Y05. % Parallel computation
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}
\label{sec:intro}

\paragraph{\textbf{Background.}} There has been a recent
renewed interest in first order subgradient bundle methods for
semidefinite programming (SDP) and eigenvalue optimization owing to
the fact that storage requirements and solution time of interior-point
methods scale poorly with increasing problem size.  In contrast, even
for large-scale problems, subgradient bundle methods can produce a
solution of low accuracy relatively quickly using only a modest amount
of memory.  To see this difference clearly, denote by $\sym{n}$ the
space of real, symmetric $n\times n$ matrices equipped with the trace
inner product, and consider the standard primal form of an SDP
\begin{equation}
\label{eq:primal-sdp}
\min\limits_{X\in\sym{n}}\;\;\ip{C}{X}\quad\st\;\;AX = b;
\;\; X \succeq 0,
\end{equation}
where $b \in \reals^m$, $A:\sym{n}\rightarrow \reals^m : X
\mapsto \adj{\left[
\ip{A_1}{X},\ldots,\ip{A_m}{X}\right]}$ is a linear
operator, the matrices $C,\,X,\,A_i\;\;(i=1,\ldots,m)\in
\sym{n}$, and $X \succeq 0$ means that $X$ is constrained
to be positive semidefinite.  The subproblem of finding a
search direction in every iteration of a typical
primal-dual interior-point method requires the solution of
an $m\times m$ linear system involving the so-called
\emph{Schur complement} matrix.  It requires $O(m^2n^2 +
mn^3)$ operations to form, and $O(m^3)$ operations to
factorize, the Schur complement matrix, while storing this
matrix and the primal variable $X$ requires $O(m^2)$ and
$O(n^2)$ memory respectively.  In practice, substantial
savings in both storage requirements and computational
time may be realized by exploiting
sparsity~\cite{FKMN00,NFFKM03}, resorting to dual-only
methods~\cite{BYZ00}, employing a low-rank factorization
technique~\cite{BM03a,BM03b}, or using carefully
preconditioned iterative
methods~\cite{TK02,Bur03,Toh04,KS05} for this linear
system.  Nevertheless, as a rule of thumb, serial
implementations of primal-dual interior-point methods are
limited to problems approximately as large as $n \approx
2,000$ and $m\approx 5,000$, with such an instance
requiring a few hundred megabytes of memory and a few
hours of solution time on a typical workstation; see
\cite{Mit03} for a recent benchmarking of SDP codes.

On the other hand, Helmberg and Rendl~\cite{HR99} observed
that whenever the primal feasible matrices $X$ are
constrained to have constant trace $a > 0$, the dual
of~\eqref{eq:primal-sdp} is equivalent to the nonsmooth,
convex, unconstrained, eigenvalue optimization problem
\begin{equation}
\label{eq:evo}
\min\limits_{y\in\reals^m}\;\; a\,\lambda_{\max}(\adj{A}y - C) -
\ip{b}{y},
\end{equation}
where 
\[
\adj{A}:\reals^m \rightarrow \sym{n}:y\mapsto
\sum\limits_{i=1}^m y_i A_i
\]
is the adjoint of $A$, and $\lambda_{\max}(\cdot)$ is the
maximal eigenvalue function.  Subgradient bundle methods,
pioneered by C.~Lemar\'{e}chal and further developed by
K.~C.~Kiwiel and others, approximate this nonsmooth
objective with a small number of cutting planes, each
derived from a subgradient of the objective function.
Hence, the storage requirements are $O(m)$.  The subproblem
of determining a search direction typically requires the
solution of a quadratic program whose size depends on the
number of cutting planes retained.  The \emph{spectral
bundle method} of~\cite{HR99}, obtained by specializing
the \emph{proximal bundle method} of Kiwiel~\cite{Kiw90},
has been remarkably successful in solving SDP relaxations
--- as large as $n\approx 6,000$ and $m\approx 20,000$ ---
of combinatorial optimization problems, within tens of
minutes on a typical workstation.

\paragraph{\textbf{Summary.}}  In this paper, we describe \LAMBDA, a
\emph{general-purpose}, \emph{portable}, \emph{parallel}
implementation of the spectral bundle method, \ie a code
that handles most types of SDP's, particularly those with
block diagonal structure (general-purpose), that is based
on the Message Passing Interface (MPI) for execution on a
variety of parallel platforms (portable), and that
operates on distributed data (parallel).  Thus \LAMBDA is
suited for problems which cannot be accommodated within
the memory of a single computer, or which cannot be solved
within a reasonable time on a single computer, or both.
Some preliminary studies to this end were conducted with
an earlier prototype~\cite{NT03} based partly on Matlab.
Numerical results from the much improved present implementation
show that the techniques proposed here significantly
extend current limits on problem sizes, with acceptable
scalability of solution times on sufficiently large
problems.

Beginning with the proximal bundle method
(\secref{sec:algorithm}), we describe how the problem is
formulated (\secref{sec:problem-formulation}) and how
problem data are distributed
(\secref{sec:data-distribution}).  We then explain how
each step of the algorithm is implemented
(\secref{sec:subgradient-computation}
--~\secref{sec:update}), focusing on block structure and
parallelism.  Numerical results illustrating various
aspects of the code are supplied
in~\secref{sec:numerical-results} prior to concluding
perspectives in~\secref{sec:perspectives}.

\paragraph{\textbf{Notation.}} As usual, lowercase letters denote
vectors and uppercase letters denote matrices.  The $\supth{j}$
coordinate vector is denoted $e_j$, and $I_n$ is the $n\times n$
identity matrix.  Parentheses are used to refer to diagonal blocks of
a block diagonal matrix (\eg $C(k)$ for the $\supth{k}$ diagonal block
of $C$) or of a block vector.  A superscript may denote an iteration
index or an exponent, with usage to be inferred from context.
Subscripts index an element in a set (\eg $A_i \; (i=1.\ldots,m))$, an
entry in a vector, a column vector, or an entry in a matrix named with
the same uppercase letter (\eg $P = [p_1,\ldots,p_r]$, or $c_{ij}$ for
the $(i,j)$ entry of $C$, or $c_{ij}(k)$ for the $(i,j)$ entry of
$C(k)$); again, usage should be evident from context.  Additionally,
some Matlab-like array indexing notation (\eg $1:n$ for indices $1$
through $n$), as well as some Matlab built-in function names are
employed in pseudocode for algorithms.  The $\svec(\cdot)$ operator
isometrically maps a symmetric matrix to a vector.  Throughout, we use
the $2$-norm on vectors.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{The Algorithm}
\label{sec:algorithm}
In this section, we describe the basic ideas behind the
proximal bundle method of Kiwiel~\cite{Kiw90}, and its
adaptation to the formulation~\eqref{eq:evo} by the
spectral bundle method.

\subsection{The Proximal Bundle Method}
\label{subsec:prox-bundle-method}
Consider the minimization of a finite-valued, convex
function $f\maps{\reals^m}{\reals}$.  Although $f$ may be
nondifferentiable, it is well known that its Moreau-Yosida
regularization
\begin{equation}
\label{eq:moreau-yosida}
f_{\rho}(x) = \min\limits_{y \in \reals^m} \left\{ f(y) + \frac{\rho}{2}
\norm{y-x}{}^2 \right\} \quad (\rho > 0)
\end{equation}
is a Lipschitz continuously differentiable convex function
with the same set of minimizers as $f$~\cite{HL93}.  The
unique minimizer $\tilde{x}$ in~\eqref{eq:moreau-yosida},
called the proximal point of $x$, satisfies
\[
\nabla f_{\rho}(x) = \rho\left(x - \tilde{x}\right),
\]
and hence represents a displacement of $x$ along the
negative gradient of $f_{\rho}$, with $\rho^{-1}$ being a
step length.  The proximal point algorithm for minimizing
$f$ is then to minimize $f_{\rho}$ by updating the current
point $x$ to its proximal point $\tilde{x}$.

The proximal bundle algorithm approximates $\tilde{x}$ by
replacing the $f(y)$ term in the right hand side
of~\eqref{eq:moreau-yosida} with a \emph{cutting plane
model} based on the subdifferential~\cite{Roc70,HL93} of
$f$.  Given a current iterate $y^k$ and some subset of
past iterations $J^k \subseteq \setof{1,\ldots,k}$, the
so-called ``bundle'', namely function values $f(y^j)
\;\;(j\in J^k)$ and a set $G^k = \setof{g^j \in
\subdiff{f}{y^j}\;:\;j \in J^k}$ of subgradients,
furnishes a polyhedral cutting plane model
\begin{equation}
\label{eq:cutting-plane-model}
\hat{f}^k(y) = \max\limits_{j \in J^k} f(y^j) + \ip{g^j}{y-y^j},
\end{equation}
which minorizes $f$.  Using $\hat{f}^k(y)$ from~\eqref{eq:cutting-plane-model} 
in lieu of $f(x)$ and setting $x = y^k$ in~\eqref{eq:moreau-yosida} yields
\begin{equation}
\label{eq:subproblem}
\min\limits_{y \in \reals^m}\;\; \left\{ \max\limits_{j \in
J^k} f(y^j) + \ip{g^j}{y-y^j} \right\} + \frac{\rho}{2}
\norm{y - y^k}{}^2
\end{equation}
as the subproblem, which is readily formulated as a
quadratic program.  By solving this program, one obtains
the optimal subgradient $g^j$ in~\eqref{eq:subproblem},
call it $\bar{g}$, as well as the unique minimizer
$\bar{y}$ which approximates the proximal point of $y^k$.
%%By interchanging the (unconstrained) minimization and the maximization
%%in~\eqref{eq:subproblem}~\cite{Roc70}, this solution
%%may be obtained as
%%\begin{equation}
%%\label{eq:subproblem}
%%\bar{y} = y^k - \frac{1}{\rho}\bar{g}, \text{ with }
%%\bar{g} \in \arg\min\limits_{g \in \cch{G^k}\;\; \frac{1}{2\rho} \ip{\bar{g}}{\bar{g}}.
%%\end{equation}
If $\bar{y}$ meets the 'sufficient descent' criterion
\begin{equation}
\label{eq:serious-step}
f(\bar{y}) \leq f(y^k) + m_L \left(\hat{f}^k(\bar{y})
- f(y^k)\right)
\end{equation}
for a given constant $m_L \in (0,1/2)$, it is accepted as
the new iterate by setting $y^{k+1} = \bar{y}$; this is
termed a \emph{serious step}.  Otherwise, $y^{k+1} = y^k$,
and we have a \emph{null step}.  In either case, the
cutting plane model $G^k$ is improved by adding $\bar{g}$
and a new subgradient $g \in \subdiff{f}{\bar{y}}$ to
yield $G^{k+1}$, thus making a serious step more likely
when~\eqref{eq:subproblem} is solved in the next
iteration.  The stopping criterion
\begin{equation}
\label{eq:stop-test}
f(y^k) - \hat{f}^k(\bar{y}) \leq \delta\left(\mid
  f(y^k)\mid + 1\right)
\end{equation}
tests if the upper bound on the maximum attainable descent
in this step (left hand side of~\eqref{eq:stop-test})
falls below a small fraction $\delta$ of the magnitude of
the objective value (right hand side
of~\eqref{eq:stop-test}), in which case $y^k$ is returned
as the computed solution.

Two more steps need to be addressed to complete the
description of the algorithm.  First, to keep storage
bounded, less significant subgradients in $G^k$ may be
lumped together into an \emph{aggregate subgradient}.
Second, for good practical performance, the penalty
parameter $\rho$ has to be judiciously varied at every
iteration, as was already recognized in the early
computational works of Lemar\'{e}chal.  Both issues have
been addressed by Kiwiel~\cite{Kiw90,Kiw83}, but we defer
their discussion to~\secref{sec:update}, and outline the
essence of the method in~\algoref{algo:PBM}.

\begin{algorithm}
\caption{(Proximal Bundle Method)}
\label{algo:PBM}
\begin{algorithmic}[1]
\STATE \textbf{(Input)} Given a starting point $y^0$, an
improvement parameter $m_L \in (0, 1/2)$, a termination
parameter $\delta > 0$, a penalty parameter $\rho^0 \geq
\rho_{\min} > 0$, where $\rho_{\min}$ is a lower bound for
the $\rho^k$ sequence.

\STATE \textbf{(Output)} An approximate minimizer $\tilde{y}$
to~\eqref{eq:evo}. 

\STATE \textbf{(Initialization)} Set the iteration
counter $k = 0$.  Compute a subgradient $g^0 \in
\subdiff{f}{y_0}$, and set $G^0 = \setof{g^0}$.

\FOR{$k = 0,1,\ldots$ until termination}

  \STATE \textbf{(Subproblem)}
  Solve~\eqref{eq:subproblem} to obtain $\bar{g}$ and
  $\bar{y}$.

  \IF{\eqref{eq:serious-step} holds} 

    \STATE Set $y^{k+1} = \bar{y}$. \COMMENT{serious step}

  \ELSE
    \STATE Set $y^{k+1} = y^k$. \COMMENT{null step}

  \ENDIF

  \IF{\eqref{eq:stop-test} is satisfied} 

     \STATE Stop.  Return $\tilde{y} = y^k$ as the
     computed solution.  

  \ENDIF

  \STATE \textbf{(Update)} Add $\bar{g}$ and a new
  subgradient $g \in \subdiff{f}{\bar{y}}$ to $G^k$, and
  perform aggregation (if necessary) to obtain $G^{k+1}$.
  Update $\rho^k$ to its new value $\rho^{k+1}$. (Details
  are deferred to~\secref{sec:update}.) 

\ENDFOR \COMMENT{---- $k$ loop ----}

\end{algorithmic}
\end{algorithm}

The convergence of the algorithm, which requires only that
the updated bundle contain $\bar{g}$ and one new
subgradient from $\subdiff{f}{\bar{y}}$, may be summarized
as follows.
\begin{theorem}[Kiwiel~\cite{Kiw83,Kiw90}]
Assume that~\algoref{algo:PBM} is executed with $\delta =
0$ in the stopping test~\eqref{eq:stop-test}.  If the
algorithm performs a finite number of serious steps, \ie
performs only null steps from some iteration $k$ onwards,
then $y^k \in \arg\min\,f$.  If the algorithm generates an
infinite sequence of serious steps, then $\setof{y^k}$ is
a minimizing sequence for $f$, and is convergent to a
minimizer of $f$, if one exists.
\end{theorem}

\subsection{The Spectral Bundle Method}
\label{sec:spectral-bundle-method}
From an algorithmic standpoint, it is convenient and more
general to treat the exact subgradients $g^j\;\;(j\in
J^k)$ as $\epsilon$-subgradients of $f$ at $y^k$ for an
appropriate $\epsilon \geq 0$, and to then view the
cutting plane model~\eqref{eq:cutting-plane-model} as
being derived from a particular inner approximation of the
$\epsilon$-subdifferential of $f$ at $y^k$; this opens up
the possibility of choosing other inner approximations
designed to better exploit the structure of $f$ and its
$\epsilon$-subdifferential.  Each choice for such an inner
approximation yields a different model for the
$\epsilon$-subdifferential of $f$ with its own 'bundling
process', yet the resulting algorithms can all be analyzed
within the same unified framework.

The \emph{spectral bundle method} of Helmberg and
Rendl~\cite{HR99} specializes the proximal bundle method
to the objective function in~\eqref{eq:evo} using an
$\epsilon$-subdifferential model tailored to exploit the
structure of the maximal eigenvalue function.  There is no
loss of generality in assuming $a = 1$, so let us consider
the objective function
\begin{equation}
\label{eq:simplified-objective}
f(y) = \lambda_{\max}(\adj{A}y - C) -\ip{b}{y}.
\end{equation}
Rayleigh's variational formulation of the maximal
eigenvalue
\[
\lambda_{\max}(Z) = \max\limits_{W\succeq 0;\; \tr{W} = 1}
\ip{W}{Z}
\]
as the support function of the compact, convex set
$\setof{W\succeq 0\;: \;\tr{W} = 1}$ leads to the well
known~\cite{Fle85,Ove88,HY95,HL93}
$\epsilon$-subdifferential formula
\begin{equation}
\label{eq:epssubdiff-lambda-max}
\epssubdiff{\lambda_{\max}}{Z}{\epsilon} =
\setof{W \succeq 0\;:\; \tr{W} = 1; \;
\ip{W}{Z} \geq \lambda_{\max}(Z) - \epsilon}
\end{equation}
for $\epsilon \geq 0$.  Exact and $\epsilon$-subgradients
$W \in \epssubdiff{\lambda_{\max}}{Z}{\epsilon}$ are
obtained by taking $W = p\adj{p}$, with $p$ being a
normalized eigenvector to the maximal eigenvalue (exact
subgradient) or a nearly maximal eigenvalue
($\epsilon$-subgradient) of $Z$.  A standard chain rule
from nonsmooth calculus~\cite{HL93} then provides
\begin{equation}
\label{eq:epssubdiff-f}
\epssubdiff{f}{y}{\epsilon} = \setof{AW - b \;:\;W \in
  \epssubdiff{\lambda_{\max}}{\adj{A}y - C}{\epsilon}}.
\end{equation}
Thus eigenvectors corresponding to the maximum eigenvalue
of $Z = \adj{A}y - C$ give rise to exact subgradients
$(\epsilon = 0)$ of $f$, while eigenvectors corresponding
to nearly maximum eigenvalues yield
$\epsilon$-subgradients of $f$.  By a harmless abuse of
terminology, we refer to both exact and
$\epsilon$-subgradients as simply 'subgradients', dropping
the $\epsilon$ prefix in the sequel.  In similar vein, we
loosely use the term 'maximal eigenpairs' to denote the
largest few eigenvalues of a matrix and their
corresponding eigenvectors, the latter being called
'maximal eigenvectors'.

At a given iteration $k$, suppose $P^k =
[p_1,\ldots,p_r]$ is an orthonormal matrix whose $r$
columns are maximal eigenvectors of $\adj{A}y - C$ at the
current iterate $y = y^k$, and possibly at some past
iterates $y = y^j \; (j \in J^k)$.  Whereas the
traditional proximal bundle method employs a polyhedral
cutting plane model $G^k$ (see
above~\eqref{eq:cutting-plane-model}) with a typical
subgradient of the form $A(p_i \adj{p_i}) - b$, a
novel feature of the spectral bundle method is the
\emph{semidefinite cutting plane model}
\[
G^k = \setof{AW - b\; :\; W = \alpha \overline{W} + P^k V
\adj{(P^k)}, \; \alpha \geq 0, \;V \succeq 0,\; \alpha +
\tr{V} = 1},
\]
where $\overline{W}$ represents an aggregate subgradient.  

In terms of computation, the evaluation of the objective function and
its subgradients amounts to computing maximal eigenpairs of the matrix
$Z = \adj{A}y - C$.  The Lanczos method (further details
in~\secref{sec:subgradient-computation}) efficiently computes the
largest eigenpair when intrinsic structure in problem data allows fast
multiplications by $Z$.  In fact, incurring only a slightly larger
computational expense, the method can extract several leading
eigenpairs, provided these eigenvalues are not too closely clustered.
Although such a clustering is inevitable as the iterates $y^k$
approach a minimizer of $f$, the Lanczos method provides an effective
way to add several $\epsilon$-subgradients in Step~14 (Update) of the
algorithm at least in the early stages, when the maximal eigenvalue of
$Z$ is typically well separated from the rest of the spectrum.  The
net result is an enhanced subdifferential model and convergence in
significantly fewer iterations, with a slightly increased cost per
iteration.

Furthermore, the dominant expense of computing maximal
eigenvectors of $Z$ is better justified by using them
within a rich semidefinite cutting plane model, rather
than within a crude (by comparison) polyhedral model.
Consequently, the subproblem --- traditionally a quadratic
program --- now becomes a conic program involving linear,
convex quadratic, and semidefinite constraints.  However,
this conic program is much smaller than the original
semidefinite program, and can be solved quickly with an
interior-point method.  The additional computational
expense incurred in solving this conic program is amply
compensated by convergence in many fewer iterations.

Beginning with the problem formulation in the next
section, we describe how each step of the algorithm may be
parallelized and adapted to block structured semidefinite
programs.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Problem Formulation}
\label{sec:problem-formulation}
\LAMBDA treats SDP's with the usual block diagonal
structure, \ie the data matrices $C,\, A_i
\;\;(i=1,\ldots,m)$ are block diagonal with block sizes
$n_1,\ldots,n_s$, and $\sym{n}$ in~\eqref{eq:primal-sdp}
now denotes by $\sym{n_1}\times \ldots \times \sym{n_s}$.
For SDP's whose primal feasible matrices may not have
constant trace, \LAMBDA allows the user to specify an
upper bound on the trace of some optimal primal solution.
This bound has to be inferred by the user based on problem
data, or must be guessed large enough so that at least one
primal solution satisfies the bound.  In this case, the primal
problem~\eqref{eq:primal-sdp} may be rewritten as
\[
\min\limits_{X\in\sym{n}}\;\;\ip{C}{X}\quad\st\;\;AX =
b;\;\;\tr{X} \leq a; \;\; X \succeq 0,
\]
and its dual as
\begin{equation}
\label{eq:modified-evo}
\min\limits_{y\in\reals^m}\;\;
a\,\max\left\{\lambda_{\max}(\adj{A}y - C),\;0\right\} - \ip{b}{y}.
\end{equation}
This is transformed without loss of sparsity to the form
in~\eqref{eq:evo} by adding an extra $1\times 1$ block to the
$X, C, \,A_i\;(i=1,\ldots,m)$ matrices:
\begin{align*}
X' = \left[
\begin{array}{cc}
X & 0 \\
0 & x
\end{array}
\right], & \quad
C' = \left[
\begin{array}{cc}
C & 0 \\
0 & 0
\end{array}
\right], & \quad
A_i' = \left[
\begin{array}{cc}
A_i & 0 \\
0   & 0 \\
\end{array}
\right]\;\;(i=1,\ldots,m),
\end{align*}
resulting in the final form in which the SDP is stored
internally.  We assume that $a = 1$ and use the simplified
form~\eqref{eq:simplified-objective} in the remainder of
the paper.

\section{Data Distribution}
\label{sec:data-distribution}

\LAMBDA distributes the matrices $C,\, A_i
\;\;(i=1,\ldots,m)$ and the entries of the vector $b_i$ to
multiple processors in the following simple way.  If $p$
processors are available, the index set
$\setof{1,\ldots,m}$ is partitioned into $p$ mutually
disjoint subsets:
\begin{equation}
\label{eq:partition}
\setof{1,\ldots,m} = M_0 \cup \ldots \cup M_{p-1},
\end{equation}
and processor $k$ is allocated the matrices $A_i \;\; (i
\in M_k)$, and that portion of the vector $b$ with entries
$b_i\;\;(i\in M_k)$.  Processor $0$, designated the root
(master) processor, also holds the matrix $C$ in addition
to its share of the $A_i\;(i\in M_0)$ matrices.  This
admittedly simple scheme offers some key advantages.
First, distributing entire matrices (rather than portions
thereof) facilitates a simple implementation.  The scheme
may be consistently applied to all problems, whether block
structured or not, whereas partitioning within a matrix
involves some detailed book-keeping, especially when such
partitions straddle block boundaries.  Second, scalability
relies on the requirement that $m$ (the number of primal
constraints) is a large multiple of the number of
available processors, while $n_i$ (the block sizes) may be
modest in comparison; this is usually the case in many
applications.  Third, based on sparsity or presence of
other structure in the data, it is easy to determine
\textit{a priori} a suitable partitioning
in~\eqref{eq:partition} that ensures a good load balance
among the processors. (Some limitations in this regard are
mentioned in~\secref{sec:perspectives}.)  Fourth, by
parallelizing the application of the $A(\cdot)$ and the
$\adj{A}(\cdot)$ operators, the computationally most
intensive parts of the algorithm (subgradient computation,
forming the data defining the subproblem) may be
effectively parallelized, leaving the serially executed
portion (solution of the subproblem, bundle update \etc)
at an acceptable level.  The root processor bears the
additional overhead of executing the serial portion.
Finally, as a fringe benefit, this scheme conserves some
memory by allowing distributed storage of the polyhedral
part (details in~\secref{sec:subproblem}) of the
subgradient bundle.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Subgradient Computation}
\label{sec:subgradient-computation}

Since maximal eigenvectors of $Z = \adj{A}y - C$ furnish
$\epsilon$-subgradients of $f$ at $y$, the Lanczos method may be
used to compute maximal eigenpairs of $Z$ --- a matrix
which inherits the combined sparsity of the $C,\,A_i\;
(i=1,\ldots,m)$ matrices.  \LAMBDA implements an
\emph{implicitly restarted block structured Lanczos
method} with an \emph{active block strategy}.  We provide
only a brief overview of the implicitly restarted Lanczos
method, referring to the excellent exposition
in~\cite{Par98,ARPACK} for full details.  Our focus here
is rather on efficiently handling block structure.

\subsection{The Implicitly Restarted Lanczos Method}
\label{subsec:lanczos-method}
The Lanczos method is essentially a Rayleigh-Ritz
procedure that approximates invariant subspaces of $Z =
\adj{A}y - C$ from a sequence $K^j(Z,v^1)\;(j=0,1,\ldots)$
of expanding Krylov subspaces
\[
K^j(Z,v^1) = \Span{v^1,Zv^1,Z^2 v^1,\ldots,Z^{j-1} v^1},
\]
where $v^1$ is a vector of unit norm with a nontrivial
component in the desired eigenspace.  The method uses $Z$
only via matrix-vector products, and is particularly
effective if the structure in $Z$ allows these products to be
quickly computed.  Each step of the method involves
augmenting a \emph{Lanczos factorization}, which consumes
progressively increasing storage and computing time.  The speed
of convergence improves if $v^1$ has a large
component in the maximal eigenspace and the maximal
eigenvalue is well separated from the rest of the
spectrum, but in practice, some form of restarting must be
employed.  If $p$ maximal eigenpairs are desired, a simple
scheme would be to restart the Lanczos process every $p$
steps, using the last Lanczos vector $v^p$ as the starting
vector $v^1$ for the next round of iterations.  At this
stage, spectral transformations (such as Chebyshev
acceleration) may be employed to enhance the component of
$v^1$ in the space spanned by the desired eigenvectors,
and hence to speed up convergence.

In our present setting, the matrix $Z = \adj{A}y - C$ is
distributed.  Inter-processor communication adds overhead
to matrix-vector products, making spectral transformations
via explicit use of matrix polynomials quite expensive.
Hence we have adopted \emph{implicit
restarting}~\cite{Sor92}, a technique by which a Lanczos
factorization of length $p+d$ is compressed to one of
length $p$ in a way that retains the information most
relevant to the desired eigenvectors of $Z$.  The ARPACK
Users' Guide~\cite{ARPACK} clearly explains how this
procedure is tantamount to a spectral transformation with
a filter polynomial of high degree.

In exact arithmetic, the method generates an expanding
orthogonal basis --- the Lanczos basis --- for the
nested Krylov subspaces $K^j(Z,v^1)\;(j=0,1,\ldots)$.
With rounding errors, periodic reorthogonalization is
necessary~\cite{Par98}.  \LAMBDA uses \emph{complete
reorthogonalization} of the Lanczos basis prior to each
restart.

\subsection{Handling Block Diagonal Structure}
\label{subsec:block-structure}
When the matrices $C,\,A_i \;\;(i=1,\ldots,m)$ are block
diagonal, $Z = \adj{A}y - C$ and its eigenvectors inherit
the same block diagonal structure, and it is important to
exploit this structure in all stages of the algorithm. To
the see the inadequacy of the standard Lanczos method, let
$s$ be the number of diagonal blocks, and denote the
$\supth{k}$ block $(k=1,\ldots,s)$ of these matrices by
$C(k), A_i(k)\;(i=1,\ldots,m)$ and $Z(k)$ respectively.
Suppose $p$ maximal eigenpairs of $Z$ are desired.  On the
one hand, the na\"{\i}ve approach (call this Method 1) of
computing the $p$ largest eigenvalues of each diagonal
block $Z(k)$ retains block structure, but is obviously
inefficient, since $ps$ eigenpairs are computed.  On the
other hand, a Lanczos process that treats the entire block
diagonal matrix $Z$ as a single linear operator on $\reals^{n_1 +
\ldots + n_s}$ can produce exactly $p$ maximal eigenpairs
(call this Method 2), but will likely converge to Ritz
vectors (and consequently, eigenvectors) without block
structure when the associated maximal eigenvalue is
multiple, and split across different blocks of $Z$.  For
example, if $Z$ has a multiple maximum eigenvalue that
occurs in every block, then the Lanczos method will likely
produce an associated eigenvector that is fully dense.
Although such eigenvectors, whether possessing block
structure or not, yield legitimate
$\epsilon$-subgradients, the loss of block structure
degrades the overall efficiency of the algorithm.

Fortunately, this can be rectified by incorporating simple
modifications into the Lanczos method; for ease of
reference, let us call the modified algorithm the Block
Structured Lanczos Method.\footnote{Not to be confused
with the block Lanczos method.}  Here we perform a
\emph{synchronized} set of $s$ Lanczos processes, one for
each block, maintaining Lanczos vectors independently for
each block (akin to Method 1).  However, none of these
processes is allowed to completely resolve $p$ maximal
eigenpairs of any block of $Z$.  Instead, after $p+d$
steps, the restarting process uses Ritz values from
\emph{all} blocks to determine unwanted eigenvalues, and
hence also implicit shifts for the entire matrix $Z$ (akin
to Method 2).  Thus the proposed method may be viewed as a
hybridization of Methods 1 and 2, with the restarting
procedure synchronizing the independent Lanczos processes.
Equivalently, it is a method that incorporates block
structure in the Lanczos basis.  This method requires
storage for $p+d$ Lanczos vectors of length $n_1 + \ldots
+ n_s$ (plus a bit more for some auxiliary arrays) and,
upon convergence, produces $p$ (approximate) maximal
eigenpairs of $Z$, not counting multiplicities, \ie{} the
method may resolve only one copy of a multiple eigenvalue
(since it is not a block Lanczos method), but a single
eigenvector corresponding to the maximum eigenvalue is all
that is needed for convergence of the bundle method.

\subsection{Active Block Strategy}
\label{subsec:abs}
It is obvious that blocks of $Z$ for which the Lanczos
process has converged (for instance, small blocks whose
sizes are smaller than the size chosen for the Lanczos
basis) may be excluded from further computation.  By
monitoring the Ritz values and their corresponding error
bounds~\cite{Par98} in every remaining block prior to each
implicit restart, we can ignore those blocks which are no
longer candidates for producing one of the largest $p$
eigenvalues; such blocks are deemed \emph{inactive}.  The
remaining blocks are called \emph{active}, and subsequent
computation is restricted to these blocks only.  This can
produce significant savings (as demonstrated
in~\secref{sec:numerical-results}) in problems with many
blocks, but where the largest $p$ eigenvalues are
concentrated in a small number of blocks in most of the
iterations.

\subsection{Null Steps}
\label{sec:null-steps}
We include a time-saving inexact evaluation procedure for
null steps suggested by Helmberg~\cite{Hel00}.  Since the
Ritz values provide progressively better lower bounds on
the maximal eigenvalues as the Lanczos method progresses,
the error bounds on the Ritz values may be used to
determine if a null step is inevitable.  In this case,
further resolution of eigenpairs is prematurely
terminated, and the algorithm performs a null step;
see~\cite{Hel00} for details.

\subsection{Parallelism}
\label{sec:parallelism}
Parallelism is easily exploited in Step~14 in the
matrix-vector products $w = Zv^{j+1}$.  Recall that $Z =
\sum_{i=1}^m y_i A_i - C$ is not available explicitly, as
the $C,\, A_i\;(i=1,\ldots,m)$ matrices reside on
different processors.  If $\mathcal{A} \subseteq
\setof{1,\ldots,s}$ denotes the index set of active blocks
at a particular iteration, then each processor $q$, including
the root processor $0$, uses the locally available
$A_i\;\;(i\in M_q)$ to compute the partial sum
\begin{equation}
\label{eq:mat-vec1}
w_q(k) = \sum_{i\in M_q} y_i A_i(k) v^{j+1}(k)\qquad
(k\in\mathcal{A}).
\end{equation}
These partial sums are subsequently collected on the root
processor which computes the desired matrix-vector product
\begin{equation}
\label{eq:mat-vec2}
w(k) = \sum_q w_q(k) - C(k) v^{j+1}(k)  
\qquad (k\in \mathcal{A}),
\end{equation}
without explicitly forming the blocks $Z(k)$.  Implicit
restarts and complete reorthogonalization are executed in
serial by the root processor.  These operations depend on
the number of Lanczos vectors stored (which is small in
comparison to the block sizes $n_1,\ldots,n_s$) and the
number of primal constraints $m$.  For large-scale problems,
these serial operations are dominated by the cost of the
matrix-vector products.

We now combine all these features to give a full
description of the eigenvector computation procedure
in~\algoref{algo:BSLM-ABS}.  Steps 4 -- 20 implement the
basic Lanczos process, while Steps 28 -- 40 deal with the
implicit restarts, for which \LAMBDA uses some
ARPACK~\cite{ARPACK} routines.  The portion in Steps 21 --
27 keeps track of the desired eigenvalues, computes the
implicit shifts for the unwanted ones, and implements the
active block strategy.  More precisely, the error bounds
for each eigenvalue produced by the Lanczos processes are
used to infer which eigenvalues could never possibly be
one of the $p$ largest ones (not counting multiplicities).
These unwanted eigenvalues then determine the implicit
shifts and the inactive blocks.

As a general-purpose SDP solver, \LAMBDA assumes no particular
prior information about the maximal eigenvector, and hence
initializes $v^1$ to be a random vector normalized within
each block.

\begin{algorithm}
\caption{Block Structured Lanczos Method with Active
  Block Strategy}
\label{algo:BSLM-ABS}
\begin{algorithmic}[1]

  \STATE \textbf{(Input)} A vector $y$, a blockwise
  normalized initial vector $v^1$, a convergence tolerance
  \texttt{tol}, the number $p$ of maximal eigenpairs
  desired, the number $d$ of additional Lanczos vectors
  allowed.

  \STATE \textbf{(Output)} $p$ triples $(\theta,
  x, k)$ such that $\theta$ is an (approximate)
  maximal eigenvalue of $Z$, with $(\theta, x)$
  being an (approximate) eigenpair of $Z(k)$.

  \STATE Initialize the active blocks index set
  $\mathcal{A} = \setof{1,\ldots,s}$, the inactive blocks
  index set $\mathcal{I} = \emptyset$, and the converged
  blocks index set $\mathcal{C} = \emptyset$.

  \STATE Compute $w = Z v^1 = (\adj{A}y - C)v^1$ in
  parallel using~\eqref{eq:mat-vec1}
  and~\eqref{eq:mat-vec2}.

  \STATE $\alpha^1 = \adj{(v^1)} w; \quad r^1 = w -
  \alpha^1 v^1; \quad V^1 = v^1; \quad T^1 =
  \alpha^1; \quad \beta^1 = \norm{r^1}{}$.

  \STATE $j = 1$.

  \FOR{$\mathrm{rst}=1,2,\ldots$ until convergence}

    \STATE \COMMENT{---- begin restart loop ----}

    \WHILE{$j < p+d$} 

      \STATE \COMMENT{---- begin Lanczos factorization ----}
 
      \FOR{$k \in \mathcal{A}$} 

        \STATE \COMMENT{---- augment Lanczos factorization
        for each block ----}
        
	\STATE $v^{j+1}(k) = r^j(k) / \beta^{j}(k); \quad
	V^{j+1}(k) = [V^j(k)\;\;v^{j+1}(k)]; \quad
	\widehat{T}^j(k) = \left[T^j(k);\; \beta^j(k) 
	\adj{(e_j(k))}\right]$.

        \STATE Compute $w(k) = Z(k) v^{j+1}(k)$ in
        parallel using~\eqref{eq:mat-vec1}
        and~\eqref{eq:mat-vec2}.

	\STATE $t(k) = \adj{V^{j+1}(k)} w(k); \quad
	r^{j+1}(k) = w(k) - V^{j+1}(k)\, t(k); \quad
	\beta^{j+1}(k) = \norm{r^{j+1}(k)}{}$.

	\STATE $T^{j+1}(k) = [\widehat{T}^j(k)\;\;t(k)]$.

	\ENDFOR \COMMENT{---- $k$ loop ----}

	\STATE $j = j+1$.

     \ENDWHILE
     
     \STATE Compute the spectral factorization $T^{p+d}(k) = S(k)
     \Theta(k) \adj{S(k)}$ for $k\in \mathcal{A}$.

     \STATE For each $(i,k) \in \setof{1,\ldots,p+d}
     \times \setof{1,\ldots,s}$, test the convergence
     condition $\beta^{p+d}(k) \abs{s_{p+d,i}(k)} <
     \mathtt{tol}$ to find indices $(i,k)$ corresponding
     to converged eigenpairs.  

     \STATE If all eigenpairs in any block $k \in
     \mathcal{A}$ have converged, set $\mathcal{C} =
     \mathcal{C}\cup\setof{k}$.

     \STATE Compute upper bounds
     \[
       \bar{\lambda}(k) =
       \max\limits_{i=1\,\ldots,p+d;\;k\in \mathcal{A}
       \cup \mathcal{C}} \theta_i(k) + \beta^{p+d}(k)
       \abs{s_{p+d,i}(k)}
     \]
     Set $\mathcal{W}$ to be the index pairs $(i,k)$ for
     the $p$ largest values among the upper bounds
     $\theta_i(k) +
     \beta^{p+d}(k)\abs{s_{p+d,i}(k)},\;\;(i=1,\,\ldots,p+d\; ;
     k \in \mathcal{A}\cup \mathcal{C})$.

     \STATE Define $\omega$ to be the $\supth{p}$ largest
     of the $\theta_i(k) -
     \beta^{p+d}(k)\abs{s_{p+d,i}}\;\; (i=1,\ldots,p+d; \;
     k\in\mathcal{A})$.  Set $\mathcal{I} = \mathcal{I}
     \cup \setof{k \in \mathcal{A}\;:\;\bar{\lambda}(k) <
     \omega}$.

     \STATE Update $\mathcal{A} = \mathcal{A} \setminus
     \setof{\mathcal{C} \cup
     \mathcal{I}}$. \COMMENT{---- active block strategy ----}

     \STATE Compute lower bound on maximal eigenvalue
     \[
       \underline{\lambda} =
       \max\limits_{i=1,\ldots,p+d;\;k\in\mathcal{A}}
       \theta_i(k) - \beta^{p+d}(k)\abs{s_{p+d,i}(k)}.
     \]
     If $\underline{\lambda}$ is large enough to
     violate~\eqref{eq:serious-step} by the margin
     required for a null step, break out of the
     $\mathrm{rst}$ loop.

     \STATE If the Ritz values $\theta_i(k)$ for all
     $(i,k)\in \mathcal{W}$ have converged to the
     prescribed tolerance $\mathtt{tol}$, break out of the
     $\mathrm{rst}$ loop.

     \COMMENT{---- synchronize Lanczos processes via
     implicit restarts ----}

     \STATE For each $k \in \mathcal{A}$, set
     $\widehat{Q}(k) = I_{p+d}$, and select the
     indices of the 'unwanted' eigenvalues 
     \begin{equation}
     \label{eq:implicit-shifts}
     \mathcal{U} = \setof{(i,k) \; : \; \theta_i(k) +
     \beta^{p+d}(k)\abs{s_{p+d,i}} < \omega}.
     \end{equation}
  
     \FOR{$(i,k)\in\mathcal{U}$} 

       \STATE $[Q(k),R] = \mathrm{qr}(T^{p+d}(k) -
       \theta_i(k)\, I_{p+d})$. \COMMENT{---- QR
       factorization with implicit shifts ----}
  
       \STATE $T^{p+d}(k) = \adj{Q(k)}T^{p+d}(k)Q(k)$.

       \STATE $\widehat{Q}(k) = \widehat{Q}(k)\,Q(k)$.
  
     \ENDFOR

     \FOR{$k \in \mathcal{A}$}

       \STATE $\beta^p_+(k) = T^{p+d}_{p+1,p}(k); \quad \sigma^p =
       \widehat{Q}_{p+d,p}(k)$.

       \STATE $r^p(k) = \beta^p_+(k) v^{p+1}(k) + \sigma^p r^p(k)$.

       \STATE $V^p(k) = V^{p+d}\widehat{Q}_{:,1:p}(k); \quad T^p(k) =
       T^{p+d}_{1:p,1:p}(k)$.

       \STATE Reorthogonalize $V^p(k)$.

    \ENDFOR \COMMENT{---- $k$ loop ----}

    \STATE $j = p$.

\ENDFOR \COMMENT{---- rst loop ----}

\STATE For each $k$ such that $(i,k) \in \mathcal{W}$, reorthogonalize
$V^{p+d}(k)$.

\STATE For each of the $p$ index pairs $(i,k) \in
\mathcal{W}$, compute an approximate eigenvector $x_i (k)= V^{p+d}(k)s_{i}(k)$.
Return the $p$ triples $(\theta_i, x_i(k), k)$ for each $(i,k)
\in \mathcal{W}$.

\end{algorithmic}
\end{algorithm}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Subdifferential Model}
\label{sec:subdifferential-model}

\LAMBDA employs a mixed polyhedral-semidefinite
$\epsilon$-subdifferential model.  At a given iteration,
the bundle consists of a polyhedral part $G =
[g_1,\ldots,g_l] \in \reals^{m\times l}$ and a
semidefinite part $P \in \reals^{n \times r}$, resulting
in a subdifferential model of the form
\[
\setof{\sum_{i=1}^l \alpha_i g_i + A(PV\adj{P}) -
  b\;:\;\alpha_i \geq 0, \; V\succeq 0,\;
  \alpha_1+\ldots+\alpha_l + \tr{V} = 1},
\]
with some columns in $G$ and in $P$ possibly retained from
earlier iterations, \ie each $g_i = A(W_i)$ for some
subgradient $W_i \in
\epssubdiff{\lambda_{\max}}{Z}{\epsilon}$, each column
$p_i$ of $P$ is a maximal eigenvector of $Z$, and $Z =
\adj{A}y - C$ is derived from the current iterate $y =
y^k$ or an earlier iterate $y = y^j \; (j\in J^k)$.

To enrich the model (in the Update step
in~\algoref{algo:PBM}), maximal eigenvectors to $Z =
\adj{A}\bar{y} - C$ are orthogonalized against, and
subsequently added to, existing columns in $P$.  This
model may be viewed as a hybrid combination of that used
by traditional bundle methods (obtained by restricting $V$
to be diagonal) and the one used by the spectral bundle
method (obtained by restricting $l = 1$, a minimal
requirement for aggregation).

Following~\cite{Nay99}, we can write the model
minimization subproblem~\eqref{eq:subproblem} as
\begin{align}
\min_{(\alpha; V)\in \reals^l\times\sym{r}} \quad & \beta \nonumber \\
\st & \half \ip{[\alpha; \svec(V)]}{Q \,[\alpha; \svec(V)]} +
\ip{\left[u;\;v\right]}{[\alpha; \svec(V)]} - \beta \leq 0
\label{eq:fos-cop} \\
    & \sum_{i=1}^l \alpha_i + \tr{V} = 1 \nonumber \\
    & \alpha \geq 0; \;\; V \succeq 0, \nonumber
\end{align}
where the matrix $Q$ and the vector $[u;\; v]$ are given by
\begin{align}
& u_i = \ip{y - \frac{1}{\rho}b}{A(W_i)} -\ip{C}{W_i}
\;\;(i=1,\ldots,l), \label{eq:u-def}\\
& v = \svec(\adj{P}(\adj{A}(y - \frac{1}{\rho}b) -
C) P)],\;\;\textrm{and} \label{eq:u-tilde-def} \\
& Q = \left[ 
\begin{array}{cc}
Q_{11} & Q_{12} \\
Q_{12}^T & Q_{22}
\end{array} \right], \;\; \textrm{with} \label{eq:Q-def} \\
& Q_{11} = \frac{1}{\rho} \adj{G}G, \label{eq:Q11-def} \\
& Q_{12} = \frac{1}{\rho} \left[
  \begin{array}{c}
    \adj{(\svec(\adj{P}\,(\adj{A}g_1)\,P))} \\
    \vdots\\
    \adj{(\svec(\adj{P}\,(\adj{A}g_l)\,P))}
  \end{array} \right], \;\;\textrm{and} \label{eq:Q12-def} \\
& Q_{22} = \frac{1}{\rho} \sum_{i=1}^m \svec(\adj{P}\,A_i\,P)\,
           \adj{(\svec(\adj{P}\,A_i\,P))}. \label{eq:Q22-def}
\end{align}
It is obvious from our data distribution scheme and the
formulas above that the evaluations of $A(\cdot)$ and
$\adj{A}(\cdot)$ in~\eqref{eq:u-def},
\eqref{eq:u-tilde-def} and~\eqref{eq:Q12-def}, as well as
the orthogonal conjugations in~\eqref{eq:Q12-def}
and~\eqref{eq:Q22-def} may be performed in parallel.
Further, observe that the polyhedral part of the bundle
can be distributed in exactly the same manner as the $C,\,
A_i \; (i=1,\ldots,m)$, \ie processor $k$ holds those rows
of $G$ whose indices $i \in M_k$.  Hence the Gram matrix
in~\eqref{eq:Q11-def} may also be computed in parallel.
Note that each processor needs the matrix $P$ to compute
the orthogonal conjugations defining $Q_{22}$, hence the
semidefinite portion of the bundle must be stored on all
processors, if the costs of communicating it are to be
avoided.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Subproblem}
\label{sec:subproblem}

The data defining the subproblem ($Q,u,v$) thus computed
in parallel are assembled on the root processor, which
solves the subproblem --- a conic program --- in serial
using SeQuL, a primal-dual interior-point code written in
the C language by J.-P.~Haeberly.  Conversion
of~\eqref{eq:fos-cop} to standard, dual form requires a
Cholesky factorization of $Q$.  Pivoting is employed to
deal with (nearly) linearly dependent subgradients that
cause rank deficiency in $Q$.  Finally, since SeQuL
accepts only inequalities in its dual formulation, we
incorporated modifications to the algorithm to accommodate
the equality constraint in~\eqref{eq:fos-cop}.

Whereas in a serial code, solving the subproblem takes a
negligible fraction of the time spent in any iteration,
this may not be the case in a parallel context.  Dominant
computational costs (subgradient computation) could be
sped up so much with multiple processors (as we will
demonstrate in~\secref{sec:numerical-results}) that, at
least for some range of problem sizes and number of
available processors, the time spent in solving the
subproblem and other serial portions of the code is no
longer negligible, and will eventually limit scalability.
The mixed polyhedral-semidefinite subdifferential model
helps in this regard by allowing a small semidefinite
model to be compensated, to some extent, by a large
polyhedral model, thus allowing the subproblem to be
solved quickly.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Update}
\label{sec:update}

To keep storage bounded, subgradients in both the
polyhedral and semidefinite parts of the bundle that are
less relevant to the model need to be lumped together to
make room for the addition of new subgradients in
subsequent iterations.  This update of the bundle is
termed \emph{aggregation}~\cite{Kiw83}, and is implemented
as described in~\cite{Nay99}.

The final ingredient is a suitable update of $\rho^k$, the
penalty parameter.  Roughly speaking, one would expect to
increase $\rho^k$ after a null step and decrease it after
a serious step, but in any case, ensuring that it never
violates the lower bound $\rho_{\min}$ required for
convergence.  Based on Helmberg's numerical
experiments~\cite{Hel03}, \LAMBDA uses a starting value of
$\rho^0 = \norm{g^0}{} / \sqrt{m}$, but updates it using a
rather involved scheme developed by Kiwiel~\cite{Kiw90}.

The updates to the bundle and to the penalty parameter are
executed in serial by the root processor.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Numerical Results}
\label{sec:numerical-results}

Our description of the experimental setup (collection of
test problems, hardware details, default parameters) is
followed by numerical experiments evaluating convergence
behavior, the proposed block structured Lanczos method for
eigenvector computation, the effect of the polyhedral
component in the subdifferential model, and finally,
parallel scalability.

\subsection{Experimental Setup}
\label{subsec:experimental-setup}

\subsubsection{Test Problems}
\label{subsubsec:test-problems}
Our overall test set consists of 13 problems
(see~\tabref{tab:test-set}), including:
\begin{itemize}
\item a fully dense, randomly generated SDP,
\item four SDP relaxations (max-cut and Lov\'asz
  $\vartheta$-function) of graph problems\footnote{Graphs
  generated using the graph generator program
  \texttt{rudy}, written by G.~Rinaldi.}, and
\item eight SDP relaxations arising in quantum chemistry
  from \textit{ab initio} calculations of electronic
  structure.\footnote{Available from
  \texttt{http://www.is.titech.ac.jp/\~{}mituhiro/software.html}
  in the SDPA format.}  
\end{itemize}
The maxG11 problem is a max-cut relaxation on a graph with
800 nodes.  For the theta problems, the names indicate the
sizes of the graphs, \eg theta-5k-67k is for a graph with
5,000 nodes and about 67,000 edges.  The quantum chemistry
problems have 23 blocks each.  Using $a(b)$ to denote $b$
blocks each of size $a \times a$, the block structure for
these problems is [6(1) 10(4) 45(4) 100(4) 120(2) 200(1)
230(1) 450(4) 1450(2)] with $m = 7230$ constraints, except
in \BHplus where the block structure is [6(5) 15(4) 20(2)
36(4) 72(1) 90(4) 94(1) 306(2)] with $m = 948$
constraints.  The total number of nonzeros in the
$C,\,A_i\;(i=1,\ldots,m)$ for the quantum chemistry
problems is about 0.85 million, except for \BHplus, where
it is 0.067 million.  For the dense, randomly generated
problem rand-1k-8k, the total number of nonzeros is about
1.17 millon.

The problem \BHplus is small enough to be solved quickly by
interior-point methods even in serial.  (Indeed SDPA solves it to high
accuracy within about 30 minutes, when executed in serial on one of
the nodes of the {\sc Kali} cluster described below.)  Nevertheless,
its solution provides insight into the convergence behavior of the
bundle method.  In particular, the quantum chemistry problems, though
modest in size, have solutions where the maximal eigenvalue has high
multiplicity --- a feature that poses a significant challenge for
bundle methods.

\begin{table}[t]
\caption{Description of problems in test set.}
\centering 
\label{tab:test-set}
\begin{tabular}{|l|c|c|p{6cm}|}\hline
Prob          &  $n_1,...,n_s$    & $m$      & Description \\ \hline \hline
rand-1k-8k    & [250 250 250 250] & 8192     & Fully dense, random \\ \hline
maxG11        & 800               & 800      & Max cut \\ \hline
theta-1k-4k   & 1,001             & 5996     & Lov\'asz $\vartheta$-function, random graph \\  \hline
theta-5k-67k  & 5001              & 67489    & Lov\'asz $\vartheta$-function, random graph \\  \hline
theta-5k-100k & 5001              & 105,106  & Lov\'asz $\vartheta$-function, random graph \\  \hline
\BHplus       &     & 948      & Quantum chemistry \\ \cline{1-1}\cline{3-4}
\Btwo         &     & 7230     & Quantum chemistry \\ \cline{1-1}\cline{3-4}
\BeO          &     & 7230     & Quantum chemistry \\ \cline{1-1}\cline{3-4}
\Ctwo         & (see~\secref{subsubsec:test-problems})    & 7230     & Quantum chemistry \\ \cline{1-1}\cline{3-4}
\Ctwoplus     &     & 7230     & Quantum chemistry \\ \cline{1-1}\cline{3-4}
\Litwo        &	    & 7230     & Quantum chemistry \\ \cline{1-1}\cline{3-4}
\LiF          &     & 7230     & Quantum chemistry \\ \cline{1-1}\cline{3-4}
\NaH          &     & 7230     & Quantum chemistry \\ \hline \hline
\end{tabular}
\end{table}

\subsubsection{Hardware}
\label{subsubsec:hardware}
All runs of our code were conducted on {\sc
Kali},\footnote{See \texttt{http://kali.math.umbc.edu} for
details.} a 64-cpu Beowulf cluster consisting of 32 dual
Intel Xeon nodes operating at 2~GHz.  Each processor has a
512~KB L2 cache.  Each node has 1~GB RAM (shared by two
processors), with the entire cluster interconnected by a
Myrinet network (2 Gbit/sec peak bandwidth).  One of these
32 nodes also acts as a \emph{storage node} with access to
a 0.5~TB RAID array.

\subsubsection{Default Parameter Values}
\label{subsubsec:default-parameters}
The default tolerance for eigenvalue computations is
$10^{-12}$.  The parameter $m_L$ for determining a serious
step is set to $0.1$.  The parameter $\rho_{\min} =
10^{-6}$, and Kiwiel's update~\cite{Kiw83} of $\rho$
requires another parameter $m_R \in \left(m_L,1\right)$,
which we set to $0.5$.  The maximum number of iterations
allowed is $10,000$.  The default accuracy level is
$\delta = 0.01$.  In every iteration, $p=5$ new
eigenvectors are computed using $p+d = 50$ Lanczos
vectors.  The subdifferential model contains $l=10$
subgradients in the polyhedral part, and $r=15$
eigenvectors in the semidefinite part.  The starting point
is always $y^0 = 0$.  Exceptions to these default values
are mentioned when appropriate.

In all tables, the following naming conventions are
employed: 'Prob' denotes problem name, 'Ser' and 'Tot'
denote the number of serious steps and the total number of
iterations, 'err' is the relative error in the objective
value, and 'Time' is computational time. Timings are
reported entirely in seconds or as hh:mm:ss
(hours:minutes:seconds).  All timings are elapsed
wall-clock time, and exclude the time taken to read
problem data.


\subsection{Convergence}
\label{subsec:convergence}
In~\figref{fig:convergence}, we plot the objective
function values for four of the smaller problems, with
each plotted point denoting a serious step.  The circles
on the plot indicate the number of iterations needed for
stopping tolerances of $\delta = 0.1, 0.01, 0.001$.  For
all problems, the plots exhibit the characteristic tailing
off effect of first order bundle methods, \ie rapid
initial progress, with the $\delta = 0.1$, and $\delta =
0.01$ accuracy levels attained within a few iterations,
and the $\delta = 0.001$ accuracy level attained between
$100$ and $1000$ iterations.

\begin{figure}[hbtp!]
\centering
\includegraphics[height=0.6\textheight,width=\textwidth]{conv}
\caption{Typical convergence behavior (serious steps only
  shown in plot) on four problems.}
\label{fig:convergence}
\end{figure}

A notable exception is \BHplus which takes an order of
magnitude more iterations to attain the $\delta = 0.001$
accuracy level.  This problem is block structured with 22
blocks making up a matrix of total dimension $1406$.  Yet
the multiplicity of the maximal eigenvalue of $Z^{\ast} =
\adj{A}y^{\ast} - C$ at any optimal solution $y^{\ast}$ is
at least $368$, as confirmed by a solution of high
accuracy computed by SDPA~\cite{SDPA}.  The slow
convergence witnessed here is an inevitable outcome of the
fact that the dimension of our subdifferential model
($l=10$ vectors in the polyhedral part $G$, and $r=15$
eigenvectors in the semidefinite part $P$) falls severely
short of the dimension of the subdifferential at the
optimal solution.  Using a computed solution $\tilde{y}$
and the true optimal value $f^{\ast}$, we calculate the
relative error as
\[
\text{err} = \frac{\abs{f(\tilde{y}) - f^{\ast}}}{1 +
\abs{f^{\ast}}}.
\]
For the three problems other than \BHplus, this quantity
is less than 1\% for $\delta = 0.01$, and less than 0.1\%
for $\delta = 0.001$.  For \BHplus, this error is 1.69\%
with $\delta = 0.01$.  With $\delta = 0.001$, this
relative error decreases to 0.7\% (at the cost of about
7900 extra iterations), but a solution satisfying the
$\delta = 0.0001$ accuracy level was unattainable within
10,000 iterations.  Relative errors for the other
instances in this problem class are similar and are given
in~\tabref{tab:rel-err}; for $f^{\ast}$, we used the
highly accurate optimal values computed in~\cite{ZBF04+}.

\begin{table}[t]
\caption{Relative errors in the optimal value for the
  quantum chemistry problems with $\delta = 0.01$.}
  \centering
\label{tab:rel-err}
\begin{tabular}{|l|l|} \hline
Problem     & err \\ \hline \hline
\BHplus     & 0.0169  \\ \hline
\Btwo	    & 0.0179  \\ \hline
\BeO	    & 0.0183  \\ \hline
\Ctwo	    & 0.0187  \\ \hline
\Ctwoplus   & 0.0195  \\ \hline
\Litwo	    & 0.0128  \\ \hline
\LiF	    & 0.0153  \\ \hline
\NaH        & 0.0255  \\ \hline
\end{tabular}
\end{table}

\subsection{Active Block Strategy}
\label{subsec:abs-results}
In~\tabref{tab:abs}, we show the effect of the active
block strategy (described in\secref{subsec:abs}) on the
time required for eigenvector calculation and overall run
time.  Each problem has two rows: (i) overall solution
time, and (ii) time for eigenvalue / eigenvector
computation.  The active block strategy produces savings
on every problem, and sometimes by significant amounts.
This is remarkable for the quantum chemistry problems,
since the maximal eigenvalue at the optimal solution not
only has high multiplicity but also occurs in most of the
blocks.  For example, this maximal eigenvalue for \BHplus,
earlier noted to have multiplicity at least $368$, occurs
in $18$ of the $22$ blocks.  Consequently, close to the
optimal solution, only $4$ blocks may be eliminated
prematurely from the matrix-vector products in the Lanczos
method.  Yet the savings on computation and communication
realized in the early iterations of the algorithm (when
the maximal eigenvalue is simple or occurs only in a few
blocks) are enough to speed up eigenvector computation and
total solution time by more than $3.5$ times.  Savings on
other problems range from barely noticeable (\Litwo,
\NaH), through significant (\Ctwoplus, \verb+rand-1k-8k+),
to dramatic (\Btwo, \BeO, \Ctwo, \LiF).

\begin{table}[t]
\caption{Effect of the active block strategy on run time.
  All problems were solved with $\delta=0.01$, $l=15$ and
  $r=10$ on 4 processors.}
\centering
\label{tab:abs}
\begin{tabular}{|c|c|c|} \hline
Name          & Time (without ABS) & Time (with ABS)\\ \hline \hline
\BHplus       & 00:25:09       & 00:07:40 \\ 
              & 00:25:00       & 00:07:00 \\ \hline
\Btwo         & 06:14:06       & 01:33:26 \\ 
              & 06:06:40       & 01:28:28 \\ \hline
\BeO          & 13:25:22       & 03:20:50 \\ 
              & 13:20:01       & 03:20:08 \\ \hline
\Ctwo         & 06:54:15       & 01:41:43 \\ 
              & 06:51:23       & 01:36:40 \\ \hline
\Ctwoplus     & 01:48:39       & 01:31:39 \\ 
              & 01:43:20       & 01:28:21 \\ \hline
\Litwo        & 02:31:03       & 02:28:28 \\ 
              & 02:26:00       & 02:23:20 \\ \hline
\LiF          & 01:25:13       & 00:50:03 \\ 
              & 01:21:40       & 00:48:21 \\ \hline
\NaH          & 01:28:21       & 01:23:44 \\ 
              & 01:25:01       & 01:20:01 \\ \hline
rand-1k-8k    & 00:20:36       & 00:15:36 \\ 
              & 00:19:51       & 00:14:42 \\ \hline
\end{tabular}
\end{table}

\subsection{Subdifferential Model}
\label{subsec:subdifferential-model-results}
Next, we show the effect of the mixed
polyhedral-semidefinite subdifferential model on parallel
run time using two sets of experiments on \BHplus.

In the first set (\tabref{tab:set1}), we use a coarse
accuracy level of $\delta = 0.01$.  The columns titled
'Total' and 'Sub' denote total solution time and time
occupied by the subproblem (in seconds); the last column 
expresses the subproblem solution time as a percentage 
of total solution time.  The first row is a serial run
taking $310$ seconds, of which subproblem solution
occupies under 12\%.  In a parallel context, solving the
subproblem could become a serial
bottleneck. \tabref{tab:set1} shows a series of runs, all
conducted on 8 processors, but with varying sizes of
polyhedral and semidefinite components in the
subdifferential model.  All runs require more or less the
same number of serious steps to achieve the specified
accuracy, but the richer models (larger $r$, smaller $l$)
do so in fewer total iterations, as expected.
Nevertheless, the simpler models (smaller $r$, larger $l$)
solve the problem faster.

\begin{table}[t]
\caption{Parallel runs on \BHplus for $\delta = 0.01$
  showing the effect of the subdifferential model on
  overall solution time.  The first row alone is a serial
  run; the remaining ones used 8 processors.}
\centering
\label{tab:set1}
\begin{tabular}{|l|l|l|l|l|l|l|} \hline
Ser / Tot & $l$ & $r$ & Total & Sub   & \%
\\ \hline \hline
24 / 282  & 10  & 25  & 310   & 37    & 11.94 \\ \hline \hline
23 / 283  & 10  & 25  & 130   & 38    & 29.23 \\ \hline
25 / 335  & 15  & 20  & 95    & 22    & 23.16 \\ \hline
25 / 351  & 20  & 15  & 84    & 11    & 13.10 \\ \hline
22 / 315  & 25  & 10  & 72    & 6.3   & 8.75 \\ \hline
24 / 402  & 30  & 5   & 65    & 2.8   & 4.31 \\ \hline
\end{tabular}
\end{table}

In the second set, we repeat a similar sequence of runs,
but using a finer accuracy level of $\delta =
0.001$. Since the optimal multiplicity is high, a large
subdifferential model will be required, so we choose $l =
10$ and $r = 50$ for the serial run shown in the first row
of~\tabref{tab:set2}.  Here solving the subproblem already
occupies a substantial $45\%$ of the total solution time,
since the blocks are small. Similar remarks as in the
first set apply to iteration counts and solution times,
but note that the last row (with $l = 50$ and $r = 10$)
shows an increase in the solution time achieved by the
penultimate row (with $l = 40$ and $r = 20$).  Thus there
is a tradeoff between saving computational time per
iteration by using a simpler model and incurring an
increase in iteration count, with a commensurate increase
in overall solution time.

\begin{table}[hbtp!]
\caption{Parallel runs on \BHplus for $\delta = 0.001$
  showing the effect of the subdifferential model on
  overall solution time.  The first row alone is a serial
  run; the remaining ones used 8 processors.}
\centering
\label{tab:set2}
\begin{tabular}{|l|l|l|l|l|l|l|} \hline
Ser / Tot & $l$ & $r$ & Total & Sub  & \%
\\ \hline \hline
44 / 1156 & 10  & 50  & 3085  & 1400 & 45.38 \\ \hline \hline
44 / 1236 & 10  & 50  & 1900  & 1400 & 73.68 \\ \hline
48 / 1554 & 20  & 40  & 1600  & 970  & 60.62 \\ \hline
45 / 1688 & 30  & 30  & 920   & 370  & 40.22 \\ \hline
42 / 1635 & 40  & 20  & 650   & 130  & 20.00 \\ \hline
44 / 2449 & 50  & 10  & 820   & 58   & 7.07  \\ \hline \hline
\end{tabular}
\end{table}

The speed-up factors of $4$ -- $5$ times obtained by a
suitable choice of the subdifferential model are due to
the small size of this problem; the subproblem solution
occupies a signicant fraction of overall solution time in
the serial run because the blocks are small.  With
increasing problem size, the speed-up factors become less
dramatic, and the choice of the subdifferential model
becomes less crucial, as evidenced by~\tabref{tab:set3}.

\begin{table}[t]
\caption{Parallel runs on \LiF\ for $\delta = 0.01$
  showing the effect of the subdifferential model on
  overall solution time.  The first row alone is a serial
  run; the remaining ones used 16 processors.}
\centering
\label{tab:set3}
\begin{tabular}{|l|l|l|l|l|l|l|} \hline
Ser / Tot & $l$ & $r$ & Total & Sub  & \%
\\ \hline \hline
32 / 354  & 10  & 25  & 5685  & 49   & 0.86 \\ \hline \hline
33 / 351  & 10  & 25  & 550   & 45   & 8.18 \\ \hline
31 / 367  & 15  & 20  & 474   & 22   & 4.64 \\ \hline
32 / 402  & 20  & 15  & 609   & 12   & 1.97 \\ \hline
30 / 433  & 25  & 10  & 676   & 5.4  & 0.80 \\ \hline
33 / 386  & 30  & 5   & 457   & 2.2  & 0.48 \\ \hline \hline
\end{tabular}
\end{table}

\subsection{Scalability}
\label{subsec:scalability}
To study parallel scalability, we solved the 10 largest
problems in our test set on processors numbering $1$
through $64$ in powers of two.  The results are reported
in~\tabref{tab:speedup}, where the columns denote the
number of processors used.  Each row contains three lines
in the following order: (i) total solution time; (ii) time
used in calculating eigenvalues (objective function
values) and eigenvectors (subgradients); and (iii) time
spent in forming the data for the subproblem.  Since
eigenvector computation is the dominant expense, the
scalability of overall solution time essentially follows
that of eigenvector computation.

\begin{table}[t]
\caption{Scalability studies on ten problems on up to 64
  processors.}
\centering
\label{tab:speedup}
\begin{tabular}{|l|l|l|l|l|l|l|l|} \hline
Name           & $1$      & $2$      & $4$      & $8$      & $16$      &  $32$     & $64$ \\ \hline \hline 
\Btwo          & 05:50:01 & 03:03:23 & 01:33:26 & 01:01:47 & 00:43:22 & 00:55:51 & 00:33:49 \\
               & 05:16:42 & 02:46:40 & 01:28:28 & 01:00:06 & 00:41:41 & 00:55:00 & 00:33:20 \\
               & 00:01:06 & 00:00:51 & 00:00:59 & 00:00:28 & 00:00:05 & 00:00:04 & 00:00:03 \\ \hline
\BeO           & 09:26:39 & 04:26:48 & 03:20:50 & 02:21:49 & 01:28:22 & 01:28:25 & 03:03:29 \\
               & 08:36:41 & 04:10:07 & 03:20:08 & 02:18:23 & 01:26:42 & 01:26:46 & 03:02:24 \\
               & 00:01:50 & 00:01:12 & 00:01:39 & 00:00:56 & 00:00:08 & 00:00:05 & 00:00:04 \\ \hline
\Ctwo          & 04:43:21 & 02:45:30 & 01:41:43 & 00:55:00 & 00:43:47 & 00:35:58 & 00:45:42 \\
               & 04:10:05 & 02:38:22 & 01:36:40 & 00:53:20 & 00:43:20 & 00:35:00 & 00:43:21 \\
               & 00:00:55 & 00:00:42 & 00:00:49 & 00:00:29 & 00:00:04 & 00:00:03 & 00:00:02 \\ \hline
\Ctwoplus      & 05:16:48 & 02:26:44 & 01:31:39 & 00:58:18 & 00:41:56 & 00:48:29 & 00:56:47 \\ 
               & 05:00:03 & 02:20:00 & 01:28:21 & 00:56:43 & 00:40:02 & 00:46:39 & 00:55:01 \\
               & 00:00:58 & 00:00:52 & 00:00:47 & 00:00:35 & 00:00:05 & 00:00:04 & 00:00:04 \\ \hline
\Litwo         & 07:46:51 & 04:10:48 & 02:28:28 & 01:31:39 & 00:53:28 & 01:01:06 & 01:10:20 \\
               & 07:30:11 & 03:53:23 & 02:23:20 & 01:30:02 & 00:51:41 & 01:00:01 & 01:08:31 \\
               & 00:01:40 & 00:01:08 & 00:01:27 & 00:00:47 & 00:00:09 & 00:00:07 & 00:00:06 \\ \hline
\LiF           & 02:36:42 & 01:01:40 & 00:50:03 & 00:35:52 & 00:23:25 & 00:20:52 & 00:39:20 \\
               & 02:23:17 & 00:58:16 & 00:48:21 & 00:33:20 & 00:21:45 & 00:18:27 & 00:38:20 \\
               & 00:00:31 & 00:00:29 & 00:00:31 & 00:00:17 & 00:00:03 & 00:00:02 & 00:00:02 \\ \hline
\NaH           & 03:36:47 & 02:01:39 & 01:23:44 & 00:39:34 & 00:41:01 & 00:28:20 & 00:08:20 \\
               & 03:20:09 & 01:55:04 & 01:20:01 & 00:38:18 & 00:40:00 & 00:28:20 & 00:07:40 \\
               & 00:00:39 & 00:00:42 & 00:00:44 & 00:00:25 & 00:00:04 & 00:00:02 & 00:00:02 \\ \hline
rand-          & 00:48:21 & 00:27:42 & 00:15:36 & 00:08:38 & 00:06:29 & 00:04:30 & 00:03:32 \\
1k-8k          & 00:43:20 & 00:26:38 & 00:14:42 & 00:08:20 & 00:06:10 & 00:04:18 & 00:03:20 \\
               & 00:00:37 & 00:00:21 & 00:00:12 & 00:00:06 & 00:00:03 & 00:00:02 & 00:00:01 \\ \hline
theta-         & 01:50:11 & 00:50:10 & 00:23:28 & 00:09:42 & 00:08:06 & 00:04:48 & 00:03:01 \\
5k-67k         & 01:45:03 & 00:46:41 & 00:18:22 & 00:08:13 & 00:06:10 & 00:04:10 & 00:02:43 \\
               & 00:00:30 & 00:00:14 & 00:01:09 & 00:00:32 & 00:00:21 & 00:00:11 & 00:00:06 \\ \hline
theta-         & 07:13:28 & 03:53:19 & 02:01:08 & 00:58:21 & 00:30:02 & 00:15:50 & 00:09:24 \\
5k-100k        & 06:56:42 & 03:36:40 & 01:56:40 & 00:56:26 & 00:26:40 & 00:13:41 & 00:08:20 \\
               & 00:08:01 & 00:04:13 & 00:02:23 & 00:01:31 & 00:01:11 & 00:00:55 & 00:00:17 \\ \hline \hline
\end{tabular}
\end{table}

A curious phenomenon in these experiments is that the same
problem solved on different numbers of processors could
yield quite different iterate sequences, producing
seemingly anomalous entries in~\tabref{tab:speedup}.
Rounding errors in the Lanczos process sometimes produce
different eigenvectors (albeit from the same eigenspace),
and hence different subgradients, and consequently a
different solution to the subproblem, which may further
trigger the conversion of a null step into a serious step
or vice versa.  In general, there could be significant
differences in iteration sequences, number of serious and
null steps, and running time across the columns
of~\tabref{tab:speedup} for any given problem.  In
particular, the suprisingly fast execution time for \NaH\
on 64 processors compared to that on 32 processors is
mostly due to the fact that the 64-processor run performed
17\% fewer serious steps with the benefit of the early
termination heuristic.  A null step is typically several
times faster than a serious step, and this has resulted in
a fortuitous speed-up.  In contrast, the 64-processor runs
of \BeO\ and \LiF\ consumed approximately the same number
of total iterations and serious steps as their respective
32-processor runs, yet the solution times have nearly
doubled.  This increase in run time as the number of
processors is increased --- typical in parallel computing
applications --- is purely due to communication losses.

The speed-up plots in~\figref{fig:speedup} for overall
solution time provide a quick picture of scalability.  In
these plots, the speed-up factor $S(p)$ for $p$ processors
is calculated as
\[
S(p) = \frac{\text{Time taken on } 1 \text{ processor
}}{\text{Time taken on } p \text{ processors }}.
\]
In summary, for the quantum chemistry problems, the speed-up
is acceptable for up to 8 processors.  (We include only
one speed-up plot from this class of problems, as the
plots for the remaining problems are more or less
similar.)  Improved solution times are obtained with up to
16 processors (and on some problems, up to 32 processors),
beyond which communication costs inevitably degrade run
time.  But we hasten to point out two facts:  (i) these
problems have relatively small blocks (with only two
moderately large blocks of size 1450 each) and a modest
value of $m = 7230$, and hence are too small to be fully
scalable up to the 64 processors available in our cluster;
and (ii) although these problems were the largest
instances available within their class at the time of
this work, they constitute toy problems in the realm of
quantum chemistry.  Realistic atomic systems would yield
much larger problems more amenable to scalability up to 64
processors.

The situation improves slightly with \verb+rand-1k-8k+.
Although comparable in size with the quantum chemistry
problems, this problem is fully dense, hence there is a
higher computation-to-communication ratio in the
matrix-vector products within the Lanczos method.
Acceptable scalability is observed for up to 8 processors,
with solution times improving up to 32 processors.

With increasing problem size, as in the two Lov\'asz
$\vartheta$-function SDP relaxations, the scalability
improves noticeably, with solution times improving
consistently all the way up to 64 processors.

\begin{figure}[hbtp!]
\caption{Speed-up plots for scalability studies. }
\centering
\label{fig:speedup}
\includegraphics[height=0.6\textheight,width=\textwidth]{speedup}
\end{figure}

\subsection{Fastest Solution Times}
\label{subsec:fastest}
Finally, we report the fastest run times that could be
obtained in solving the large quantum chemistry instances
to low accuracy with $\delta = 0.01$; this corresponds to
relative errors given in~\tabref{tab:rel-err}.  To obtain
the fastest solution times, we compute only $p=2$
eigenvectors in each iteration using $p+d=20$ Lanczos
vectors.  In~\tabref{tab:fastest}, we show
\begin{itemize}
\item the number of iterations (serious steps,
  \textit{ser}; total iterations, \textit{tot}),

\item details pertaining to eigenvalue computation (number
of implicit restarts, \textit{rst}; number of
matrix-vector products counting each block separately,
\textit{ops}; time spent in eigenvalue computation,
\textit{time}),

\item details pertaining to the subproblem (number of
  subproblems solved, \#; the total number of
  interior-point iterations, \textit{ip.it}; time spent in
  solving subproblems, \textit{time}), and

\item in the last column, the total solution time for each
  problem instance.
\end{itemize}

\begin{table}[t]
\caption{Fastest solution times on the large quantum
  chemistry problems using $\delta = 0.01, l = 10, r = 20,
  p = 2, q = 20$ on 16 processors.  The \textit{rst} and
  \textit{ip.it} columns are in thousands, and the \textit{ops}
  column is in millions.}  \centering
\label{tab:fastest}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|} \hline
Name          & \multicolumn{2}{c|}{It}  & \multicolumn{3}{c|}{Eig} &  \multicolumn{3}{c|}{Sub}   &  Time \\ \cline{2-9}
              & \textit{ser} & \textit{tot}   & \textit{rst} & \textit{ops} & \textit{time}     & \#   & \textit{ip.it} & \textit{time}     &         \\ \hline \hline
\Btwo         & 35  & 584   & 20.6 & 1.16 & 00:14:10 & 596  & 20.3 & 00:00:47 & 00:16:12\\ \hline
\BeO	      & 32  & 580   & 18.7 & 0.97 & 00:12:00 & 609  & 21.9 & 00:00:51 & 00:13:55 \\ \hline
\Ctwo	      & 34  & 547   & 18.2 & 0.94 & 00:12:30 & 564  & 19.1 & 00:00:41 & 00:14:16 \\ \hline
\Ctwoplus     & 35  & 1655  & 53.6 & 2.33 & 00:30:02 & 1701 & 59.1 & 00:02:00 & 00:36:33 \\ \hline
\Litwo	      & 30  & 761   & 33.1 & 1.60 & 00:32:03 & 791  & 26.7 & 00:00:50 & 00:34:21 \\ \hline
\LiF	      & 31  & 367   & 12.0 & 0.55 & 00:06:50 & 393  & 10.2 & 00:00:22 & 00:08:01 \\ \hline
\NaH          & 35  & 453   & 13.2 & 0.55 & 00:10:08 & 481  & 14.9 & 00:00:32 & 00:11:54 \\ \hline
\end{tabular}
\end{table}

As a rough comparison with interior-point methods, we
mention that these problems (in addition to others) were
originally solved in~\cite{ZBF04+} to high accuracy (6 or
7 digits) using the parallel interior-point code
SDPARA~\cite{SDPARA} on the NERSC SMP (shared memory
multiprocessor) cluster {\sc Seaborg}.\footnote{See
\texttt{http://www.nersc.gov/nusers/resources/SP} for
the hardware details of this cluster.}  Each problem required
about 14 hours of solution time and 27~GB of memory on a
16-cpu node with 375~MHz processors (slower than {\sc
Kali}'s 2~GHz processors), but 8~MB of L2 cache and 64~GB
of shared memory (much larger than {\sc Kali}'s 512~KB L2
cache per processor and 1~GB of memory per dual node).
Shared memory eliminates expensive communication over a
network.  A recent update~\cite{FBN05+} to these results,
using an improved version of SDPARA called SDPARA-SMP, has
reduced solution time to about 13 hours and memory usage
to about 5.7~GB.  With this improved code, the authors
were able to solve slightly larger atomic systems (but
significantly larger SDP's with $m \approx 20,000$ and a
largest block of size around $3200$) in about 5.3 days
using 73.9~GB of memory on a 64-cpu subcluster of {\sc
Seaborg}.  Recently, Mazziotti~\cite{Maz04a,Maz04b} has
reported results on even larger atomic systems (although
using a weaker SDP relaxation than employed
in~\cite{FBN05+}) based on a different type of first order
method for semidefinite programming.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Perspectives}
\label{sec:perspectives}

Subgradient methods, though used for eigenvalue
optimization over 30 years ago by Cullum, Donath and
Wolfe~\cite{CDW75}, had not seen widespread use in
semidefinite programming until revived by the spectral
bundle method of Helmberg and Rendl~\cite{HR99} and
improved variants~\cite{HK02} thereof.  Our implementation
draws much from Helmberg's serial code
\texttt{SBmethod}~\cite{SBmethod}, which has been
remarkably successful in solving large-scale problems,
albeit restricted to those arising in graph applications.

Our efforts to extend the applicability of the bundle
methodology center around efficiently handling block
diagonal structure, while resorting to parallelism to
solve very large-scale problems.  The proposed data
distribution scheme allows efficient storage of problem
data and subgradients together with ease of
implementation, without sacrificing performance on
algorithmic components (block structured Lanczos, the
active block strategy, implicit restarting, and a choice
of subdifferential models).  However, this scheme is not
without limitations.  In problems where sparsity levels
are vastly different among the data matrices
$C,\,A_i\;(i=1,\ldots,m)$, this scheme results in load
imbalance.  Thus some important problem classes (notably
maximum cut and graph bisection) cannot be effectively
parallelized.  No single scheme can work equally well for
all types of problems, and these problem classes are best
handled in a problem-dependent way.  A second limitation
in the present implementation is its inability to
effectively handle linear inequalities and bound
constraints on the variables. Such constraints serve to
tighten SDP relaxations of combinatorial optimization
problems, and \texttt{SBmethod} handles them by
incorporating a second (inner) iterative procedure within
the proximal bundle method.  Thus \LAMBDA is presently not
fully tailored for applications in combinatorial optimization.
Finally, \texttt{SBmethod} uses Lanczos vectors, even
those that have not converged to eigenvectors, to
construct cutting planes.  We believe this is not
essential in \LAMBDA, since it already has the option of
retaining many additional subgradients in the polyhedral
part of the bundle.

The Block Structured Lanczos Method offers an 
advantage\footnote{We thank an anonymous referee for this
observation.} even for problems without block structure,
when combined with a preprocessing technique introduced
in~\cite{FKMN00,NFFKM03}.  This technique, originally
proposed to exploit sparsity, uses positive semidefinite
matrix completion to convert a semidefinite program with a
single large block into one with several smaller blocks.
Thus subgradient computation with the Block Structured
Lanczos Method is likely to be more efficient on 
the modified problem with multiple blocks than on the 
original problem with a single large block.

We now comment on performance.  Throughout, we have used
basic blocking communication, \ie a processor executing a
'Send' ('Receive') will wait until the receiving (sending) 
processor executes a matching 'Receive' ('Send').  Using more
sophisticated MPI communication modes, it is possible to
overlap computation with communication in some parts of
the algorithm.  Further performance improvement could come
from selective orthogonalization and by using a block
Lanczos method.  The block Lanczos method is better at
resolving multiple eigenvalues.  Also, a single
communication of a block of vectors incurs a smaller
network latency cost than multiple communications
each involving a single vector.  However, the biggest
gains are to be realized by fully exploiting problem
structure, which is unfortunately completely lost by
encoding problems in the SDPA data format.  Thus, in our
experiments, the data matrices $C,\,A_i\;(i=1,\ldots,m)$
from all the quantum chemistry problems were treated as
unstructured, sparse matrices.  On the whole, given that
the Lanczos method is an intrinsically serial process (one
cannot compute $Z^2 v$ before computing $Zv$), and that
the parallelization is fine-grained at the linear algebra
level, the observed speed-up factors and solution times
are better than anticipated.

Finally, we fully acknowledge that the low accuracy
solutions computed for the quantum chemistry problems do
not meet the 6 or 7 digits of accuracy required for ground
state energy calculations.  We neither claim to have
'solved' these problems, nor do we suggest that first
order bundle methods are a viable solution methodology for
them.  We have merely used them as realistic, challenging
instances which allow all aspects of this general-purpose
code to be fully tested, and as such, the numerical
results for this problem set are to be viewed in that
light.  Attaining high accuracy by incorporating some
limited form of second order information, ideally in an
parallelizable way, remains a topic worthy of future
investigation.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section*{Acknowledgments}
We are grateful to Jean-Pierre Haeberly for allowing use
of the SeQuL code; to Mituhiro Fukuda for helpful
conversations about the quantum chemistry problems; and
especially to Jorge Mor\'e for access to the Chiba City
cluster at Argonne National Lab during the early stages of
this work.  We also thank the anonymous referees for their
constructive comments.

This work was done while the author was in the Department of
Mathematics \& Statistics at the University of Maryland, Baltimore
County, and was supported in part by NSF grants DMS-0238008 and
DMS-0215373.  The computational results presented here would not have
been possible without the high-performance cluster {\sc Kali}, which
was funded in part by the latter grant.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibliographystyle{spmpsci}
\bibliography{paper}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\end{document}
%% end of file `paper.tex'
