# A Primal-Dual Framework for Symmetric Cone Programming

Jiaqi Zheng Department of Computer Science, National University of Singapore
Antonios Varvitsiotis Engineering Systems and Design Pillar, Singapore
University of Technology and Design Tiow-Seng Tan11footnotemark: 1 Wayne
Lin22footnotemark: 2

###### Abstract

In this paper, we introduce a primal-dual algorithmic framework for solving
Symmetric Cone Programs (SCPs), a versatile optimization model that unifies
and extends Linear, Second-Order Cone (SOCP), and Semidefinite Programming
(SDP). Our work generalizes the primal-dual framework for SDPs introduced by
Arora and Kale [6], leveraging a recent extension of the Multiplicative
Weights Update method (MWU) to symmetric cones. Going beyond existing works,
our framework can handle SOCPs and mixed SCPs, exhibits nearly linear time
complexity, and can be effectively parallelized. To illustrate the efficacy of
our framework, we employ it to develop approximation algorithms for two
geometric optimization problems: the Smallest Enclosing Sphere problem and the
Support Vector Machine problem. Our theoretical analyses demonstrate that the
two algorithms compute approximate solutions in nearly linear running time and
with parallel depth scaling polylogarithmically with the input size. We
compare our algorithms against CGAL as well as interior point solvers applied
to these problems. Experiments show that our algorithms are highly efficient
when implemented on a CPU and achieve substantial speedups when parallelized
on a GPU, allowing us to solve large-scale instances of these problems.

Keywords: Symmetric Cone Programming, Approximation Algorithms, Geometric
Optimization

##  Introduction

A linear conic program (LCP) is an optimization problem that minimizes a
linear objective function over the intersection of a convex cone $\mathcal{K}$
and a finite number of halfspaces:

$\begin{aligned} \max\quad&{\bm{c}}\bullet{\bm{x}}\\\
\operatorname*{s.t.}\quad&{\bm{a}}_{j}\bullet{\bm{x}}\leq b_{j}\quad\forall
j\in[m]\\\ &{\bm{x}}\succeq_{\mathcal{K}}0\end{aligned}\hskip
40.00006pt\begin{aligned} \min\quad&{\bm{b}}^{T}{\bm{y}}\\\
\operatorname*{s.t.}\quad&\textstyle\sum_{j=1}^{m}{\bm{a}}_{j}y_{j}\succeq_{\mathcal{K}^{*}}{\bm{c}}\\\
&{\bm{y}}\in\mathbb{R}^{m}_{+}\end{aligned}$ (LCP)

Here $\mathcal{K}^{*}$ denotes the dual cone of $\mathcal{K}$ and
$\succeq_{\mathcal{K}}$ denotes the generalized inequality induced by the cone
$\mathcal{K}$, i.e. ${\bm{x}}\succeq_{\mathcal{K}}{\bm{y}}$ iff
${\bm{x}}-{\bm{y}}\in\mathcal{K}$. Some of the most important convex
optimization models are instances of LCPs with respect to appropriate convex
cones. For example, LCPs over the $d$-dimensional non-negative orthant
$\mathbb{R}^{d}_{+}$ correspond to linear programs (LPs), and LCPs over the
cone of $d\times d$ positive semidefinite (PSD) matrices $\mathcal{S}^{d}_{+}$
correspond to semidefinite programs (SDPs). Additionally, LCPs related to
second-order cone $\mathcal{Q}^{d+1}$ correspond to second-order cone programs
(SOCPs). These optimization models – LPs, SDPs, and SOCPs – hold a prominent
place in the realm of convex optimization, finding applications in a wide
array of domains.

A symmetric cone is a closed convex cone that is self-dual and homogeneous.
The non-negative orthant, the cone of (real or Hermitian) PSD matrices, the
second-order cone, and Cartesian products thereof are all examples of
symmetric cones, and a symmetric cone program (SCP) is simply an LCP over a
symmetric cone. Consequently, SCPs provide a unifying framework for studying
LPs, SDPs, and SOCPs. From an algorithmic standpoint, SCPs admit self-
concordant barrier functions, which facilitate the development and analysis of
interior point methods (IPMs) [20, 35].

Many approximation algorithms have been introduced in the literature for LPs
and SDPs based on the Multiplicative Weight Update (MWU) method over the
simplex, e.g. [39, 4], and its matrix variant over the set of density
matrices, e.g. [6]. In this work, we focus on algorithms based on the MWU
framework that are primal-dual, in the sense that, at the end of execution, we
obtain a pair of primal/dual feasible solutions with an optimality gap bounded
by the desired error tolerance. An important instance of this is the primal-
dual SDP framework by Arora and Kale for SDPs [6]. MWU-based algorithms have a
running time that is nearly linear in the input size and are amenable to
parallelization due to the nature of the multiplicative update.

Motivated by the success of such methods for LPs and SDPs, in this work we try
to understand:

Can we develop primal-dual MWU-based algorithms for symmetric cone programs
that have nearly linear time complexity and are effectively parallelizable?

In this work, we introduce the first primal-dual approximation framework for
SCPs, which is detailed in Section 3. Our framework is motivated by, and
generalizes, the primal-dual framework for SDPs developed by Arora and Kale
[6]. To devise our framework, we rely on the recently introduced symmetric
cone multiplicative weight update method (SCMWU) [13] that is designed for
online optimization over symmetric cones.

In contrast to existing MWU-based algorithms designed solely for LPs or SDPs,
our new framework can handle SOCPs and, additionally, SCPs with mixed conic
constraints, e.g., problems featuring both linear constraints and positive-
semidefinite constraints. Compared to IPMs, which require solving systems of
linear equations at each iteration, our framework has significantly lower
computational costs and can be effectively parallelized. Additionally, we
offer geometric insights into our primal-dual framework that are also
applicable to the previous MWU-based frameworks for LPs [39] and SDPs [6].

As a practical demonstration of our framework we develop approximation
algorithms for two computational geometry problems. The first one is the
smallest enclosing sphere (SES) problem of finding the sphere of minimum
radius that encloses an input set of spheres, and the second one is the
support vector machine (SVM) problem of finding a hyperplane that separates
two sets of data with maximum margin. Our framework allows various degrees of
freedom that necessitate further specification for each individual problem,
which we provide in Sections 4.1 and 4.2 respectively. When applied to the SES
problem, our framework yields a parallelizable approach (dual-only algorithm
in this context) that computes a $(1+\varepsilon)$-approximate solution with a
running time of $O(\frac{nd\log n}{\varepsilon^{2}})$, where $n$ represents
the number of inputs and $d$ is the dimensionality. When applied to the SVM
problem, our framework yields a parallelizable primal-dual algorithm that
computes a $(1-\varepsilon)$-approximate solution for SVM and a
$(1+\varepsilon)$-approximate solution for its dual counterpart in
$O(\frac{End\log n}{\varepsilon^{2}})$ time, where $n$ is the number of data
points, $d$ is the dimensionality, and $E$ is a instance-specific parameter
that measures its difficulty.

In addition to the theoretical analyses, our algorithms for SES and SVM have
been implemented in both sequential (CPU) and parallel (GPU) settings.
Extensive experiments have been conducted to compare them with CGAL and IPM
solvers in Section 5. The results show that our algorithms are GPU-friendly
and possess fairly large parallelism. For large-scale input instances, our GPU
implementations outperform the commercial solvers.

Related work.  Plotkin, Shmoys, and Tardos [39] introduced a framework based
on the classical MWU method for computing approximate solutions to packing and
covering LPs with a running time that is nearly linear in the input size. The
framework of Plotkin et al. [39] was later employed by Klein and Lu [28] to
devise fast approximation algorithms for the Max Cut problem, which was
further extended by Arora et al. [4] for solving general SDPs. These methods
can be categorized as primal-only methods. They can only find approximately
feasible primal solutions that satisfy every constraint up to an additive
error $\varepsilon$, and the running time is proportional to
$1/\varepsilon^{2}$, which poses challenges when a strictly feasible solution
is required or when $\varepsilon$ is very small. The problem is addressed when
Arora and Kale [6] introduce their primal-dual SDP framework. Their framework
is based on the matrix multiplicative weight update method (MMWU) that has
been discovered in different fields including optimization [8], machine
learning [44], and online learning [47]. It produces a pair of solutions to
both primal and dual SDPs, which are strictly feasible and have multiplicative
approximation ratios. Together with specific combinatorial techniques, Arora
and Kale also developed efficient approximation algorithms for NP-hard graph
problems such as Sparsest Cut and Balanced Separator.

The MWU-based frameworks (for both LPs and SDPs) rely on the “width” of the
problem, a parameter related to the largest entry or eigenvalue of the
constraints. As a closely related area of study, the literature also contains
a range of sequential and parallel algorithms that are designed for specific
families of LPs and SDPs, known as positive LPs (see, e.g. [30, 50, 2]) and
positive SDPs [26, 1]. Unlike the MWU-based algorithms, the running time of
these algorithms is typically independent of the width. However, these width-
independent algorithms often necessitate more sophisticated analysis, and most
of the techniques developed in this line of research are not applicable to the
generic problems.

##  Preliminaries

###  Symmetric cones and Euclidean Jordan algebras

A _symmetric cone_ is a closed convex cone $\mathcal{K}$ in a finite-
dimensional inner product space $\mathcal{J}$ that is self-dual (i.e.,
$\mathcal{K}^{*}\triangleq\\{{\bm{y}}\in\mathcal{J}:{\bm{y}}\bullet{\bm{x}}\geq
0,\forall{\bm{x}}\in\mathcal{K}\\}=\mathcal{K}$) and homogeneous (i.e., for
any ${\bm{u}},{\bm{v}}\in\mathbf{int}(\mathcal{K})$ there exists an invertible
linear transformation $\mathcal{T}:\mathcal{J}\to\mathcal{J}$ such that
$\mathcal{T}({\bm{u}})={\bm{v}}$ and $\mathcal{T}(\mathcal{K})=\mathcal{K}$).
Importantly, a symmetric cone can also be characterized as the cone of squares
of an _Euclidean Jordan algebra_ (EJA, e.g. see Proposition 2.5.8 in [46]):

$\mathcal{K}=\\{{\bm{x}}^{2}\triangleq{\bm{x}}\circ{\bm{x}}:\forall{\bm{x}}\in\mathcal{J}\\}.$

Here, we recall that an EJA is a finite-dimensional vector space $\mathcal{J}$
equipped with a bilinear product
$\circ:\mathcal{J}\times\mathcal{J}\to\mathcal{J}$ that satisfies
${\bm{x}}\circ{\bm{y}}={\bm{y}}\circ{\bm{x}}$ and
${\bm{x}}^{2}\circ({\bm{x}}\circ{\bm{y}})={\bm{x}}\circ({\bm{x}}^{2}\circ{\bm{y}})$
for all ${\bm{x}},{\bm{y}}\in\mathcal{J}$, an inner product $\bullet$ that
satisfies
$({\bm{x}}\circ{\bm{y}})\bullet{\bm{z}}={\bm{x}}\bullet({\bm{y}}\circ{\bm{z}})$
for all ${\bm{x}},{\bm{y}},{\bm{z}}\in\mathcal{J}$, and an identity element
${\bm{e}}$ satisfying ${\bm{e}}\circ{\bm{x}}={\bm{x}}\circ{\bm{e}}={\bm{x}}$
for all ${\bm{x}}\in\mathcal{J}$. Every EJA has a finite _rank_ , which is
defined as the largest degree of the minimal polynomial of any of its
elements.

Every element ${\bm{x}}$ of an EJA has a spectral decomposition, in the sense
that there exists unique (up to ordering and multiplicities) real numbers
$\lambda_{1},\cdots,\lambda_{r}$ and a Jordan frame
$\\{{\bm{q}}_{1},\ldots,{\bm{q}}_{r}\\}$ such that
${\bm{x}}=\sum_{i=1}^{r}\lambda_{i}{\bm{q}}_{i}$, where a Jordan frame is a
collection of primitive idempotents ${\bm{q}}_{1},\ldots,{\bm{q}}_{r}$ that
satisfy ${\bm{q}}_{i}\circ{\bm{q}}_{j}={\bm{0}}\ \forall i\neq j$ and
$\sum_{i=1}^{r}{\bm{q}}_{i}={\bm{e}},$ and $r$ is the rank of the EJA. An
element ${\bm{q}}\in\mathcal{J}$ is an idempotent if ${\bm{q}}^{2}={\bm{q}}$
and is said to be primitive if it is nonzero and cannot be written as the sum
of two nonzero idempotents. The _trace_ of an EJA element with spectral
decomposition ${\bm{x}}=\sum_{i}^{r}\lambda_{i}{\bm{q}}_{i}$ is given by the
sum of its eigenvalues $\mathbf{Tr}({\bm{x}})=\sum_{i}^{r}\lambda_{i}$. In the
following, the inner product of the EJA is denoted by
${\bm{x}}\bullet{\bm{y}}=\mathbf{Tr}({\bm{x}}\circ{\bm{y}})$. The cone of
squares of an EJA is then the set of elements with non-negative eigenvalues.

The characterization of symmetric cones as the cone of squares of EJAs, taken
together with the classification result for finite-dimensional EJAs (e.g. see
[19]), implies that all symmetric cones over the real field are isomorphic to
Cartesian products of PSD cones and second-order cones. We now present some
examples of symmetric cones, together with their algebraic properties which
shall be important in our applications:

  * •

The non-negative orthant $\mathbbm{R}^{d}_{+}$. The EJA product associated
with $\mathbbm{R}^{d}_{+}$ is the component-wise product
${\bm{x}}\circ{\bm{y}}=\mathbf{diag}({\bm{x}})\ {\bm{y}}$, and the associated
identity element is ${\bm{e}}={\bm{1}}_{d}$. The rank of $\mathbbm{R}^{d}_{+}$
is $d$, and the spectral decomposition of an element ${\bm{x}}$ is
$\sum_{i=1}^{d}x_{i}{\bm{e}}_{i}$, where ${\bm{e}}_{i}$ is the $i$-th standard
basis vector. The trace is then $\mathbf{Tr}({\bm{x}})=\sum_{i=1}^{d}x_{i}$.

  * •

The second-order cone $\mathcal{Q}^{d+1}$. We consider $\mathbbm{R}^{d+1}$ as
an EJA with Jordan product
$({\bm{u}};u_{0})\circ({\bm{v}};v_{0})=\frac{1}{\sqrt{2}}(v_{0}{\bm{u}}+u_{0}{\bm{v}};{\bm{u}}^{T}{\bm{v}}+u_{0}v_{0})$.
The corresponding cone of squares is the second-order cone
$\mathcal{Q}^{d+1}\triangleq\\{({\bm{u}};u_{0})\in\mathbbm{R}^{d+1}:\|{\bm{u}}\|_{2}\leq
u_{0}\\}.$ The rank of this EJA is 2 and its identity element is
${\bm{e}}=({\bm{0}};\sqrt{2})$. For an element $({\bm{u}};u_{0})$ of this EJA,
its eigenvalues are
$\lambda_{1,2}=\frac{1}{\sqrt{2}}(u_{0}\pm\|{\bm{u}}\|_{2})$ and the
corresponding idempotents are
${\bm{q}}_{1,2}=\frac{1}{\sqrt{2}}(\pm{\bm{v}};1)$, where
${\bm{v}}=\frac{{\bm{u}}}{\|{\bm{u}}\|}$ if ${\bm{u}}\neq{\bm{0}}$ and any
vector with unit $\ell_{2}$-norm otherwise. The trace of an element
$({\bm{u}};u_{0})$ in this EJA is given by
$\mathbf{Tr}({\bm{u}};u_{0})=\sqrt{2}u_{0}$.

  * •

The PSD cone $\mathcal{S}^{d}_{+}$. The underlying EJA of the cone of $d\times
d$ Hermitian PSD matrices is the vector space $\mathbbm{R}^{d(d+1)/2}$ endowed
with the Jordan product ${\bm{x}}\circ{\bm{y}}={1\over 2}{\rm
vec}({\bm{X}}{\bm{Y}}+{\bm{Y}}{\bm{X}})$, where ${\bm{X}}={\rm mat}({\bm{x}})$
and ${\bm{Y}}={\rm mat}({\bm{y}})$. The rank of the EJA is $d$ and its
identity element ${\bm{e}}={\rm vec}(I)$, where $I$ is the identity matrix.
For an element ${\bm{x}}$ in the EJA, its spectral decomposition is
${\bm{x}}=\sum_{i=1}^{d}\lambda_{i}{\rm vec}({\bm{v}}_{i}{\bm{v}}_{i}^{T})$,
where $\lambda_{i}$ and ${\bm{v}}_{i}$ are from the eigendecomposition of
${\rm mat}({\bm{x}})$. The trace of ${\bm{x}}$ is simply $\mathbf{Tr}({\rm
mat}({\bm{x}}))$.

  * •

Product cones. The Cartesian products of symmetric cones are also symmetric
cones. We denote the Cartesian product of $n$ symmetric cones
$\mathcal{K}_{1},\dots,\mathcal{K}_{n}$ as
$\mathcal{K}=\operatorname*{\operatorname*{\oplus}}_{i=1}^{n}\mathcal{K}_{i}$.
(We will also use $\operatorname*{\operatorname*{\oplus}}$ to denote vector
concatenations.) The rank/trace of a product cone is defined as the sum of the
ranks/traces of its components. The spectral decomposition of the product cone
follows from the decomposition of each component.

The spectral decomposition allows us to define the _Löwner extension_ of the
exponentiation function $\exp:\mathbb{R}\to\mathbb{R}$, which gives a function
that maps elements of the EJA with spectrum in $\mathbb{R}$ to elements of the
EJA as follows:

$\mathbf{exp}:\quad\sum_{i=1}^{r}\lambda_{i}{\bm{q}}_{i}\mapsto\sum_{i=1}^{r}\exp(\lambda_{i}){\bm{q}}_{i}\
\in\mathcal{J}.$

Let $(\mathcal{J},\circ)$ be an EJA and ${\bm{x}},{\bm{y}}\in\mathcal{J}$, the
generalized Golden-Thompson inequality holds [42]:

$\mathbf{Tr}(\mathbf{exp}({\bm{x}}+{\bm{y}}))\leq\mathbf{Tr}(\mathbf{exp}({\bm{x}})\circ\mathbf{exp}({\bm{y}})).$
(1)

Moreover, for any element ${\bm{x}}\in\mathcal{J}$, let
$\sum_{i}\lambda_{i}{\bm{q}}_{i}$ be its spectral decomposition. We define the
infinity norm of ${\bm{x}}$ as the maximum magnitude of its eigenvalues,
namely $\|{\bm{x}}\|_{\infty}\triangleq\max_{i\in[r]}|\lambda_{i}|$.

Additional details on symmetric cones and EJAs can be found in [19], [46], and
[45].

###  Symmetric cone programming

A Symmetric Cone Program (SCP) corresponds to linear optimization over the
intersection of a symmetric cone with an affine subspace. A pair of
primal/dual SCPs over a symmetric cone $\mathcal{K}$ (in standard form) is
given by:

$\begin{aligned} \max\quad&{\bm{c}}\bullet{\bm{x}}\\\
\operatorname*{s.t.}\quad&{\bm{a}}_{j}\bullet{\bm{x}}=b_{j}\quad\forall
j\in[m]\\\ &{\bm{x}}\succeq_{\mathcal{K}}0\end{aligned}\hskip
40.00006pt\begin{aligned} \min\quad&{\bm{b}}^{T}{\bm{y}}\\\
\operatorname*{s.t.}\quad&\sum_{j=1}^{m}{\bm{a}}_{j}y_{j}\succeq_{\mathcal{K}}{\bm{c}}\end{aligned}$
(2)

As a consequence of the classification of EJAs, symmetric cones provide a
unified framework to study Linear Programs (LPs), Semidefinite Programs
(SDPs), and Second-Order Cone Programs (SOCPs). Specifically, when
$\mathcal{K}=\mathbb{R}^{d}_{+}$, SCPs correspond to LPs; when $\mathcal{K}$
is a Cartesian product of second-order cones, they correspond to SOCPs; and
when $\mathcal{K}$ is the cone of Hermitian PSD matrices
$\mathcal{S}^{d}_{+}$, they correspond to SDPs.

###  Symmetric cone multiplicative weights update (SCMWU)

SCMWU is a recently introduced algorithm for online linear optimization (OLO)
over the trace-one slice of a symmetric cone [13]. It extends the seminal MWU
method for OLO over the simplex [5] and the MMWU method for OLO over the set
of density matrices [27, 44].

We now describe the online linear optimization framework where the SCMWU finds
its application. At each decision epoch $t$, a decision maker needs to choose
a trace-one element ${{\bm{p}}}^{(t)}$ within a symmetric cone $\mathcal{K}$.
Once this choice is made, a linear loss function
${\ell}^{(t)}({\bm{p}})={{\bm{m}}}^{(t)}\bullet{\bm{p}}$, where
${{\bm{m}}}^{(t)}$ is in the underlying EJA space of $\mathcal{K}$, is
revealed. Employing SCMWU, the next iterate is determined by

${{\bm{p}}}^{(t+1)}=\frac{\mathbf{exp}(-\eta\sum_{\tau=1}^{t}{{\bm{m}}}^{(\tau)})}{\mathbf{Tr}(\mathbf{exp}(-\eta\sum_{\tau=1}^{t}{{\bm{m}}}^{(\tau)}))},$
(SCMWU)

and ${{\bm{p}}}^{(1)}$ is initialized as $\frac{{\bm{e}}}{r}$, where $r$
represents the rank of the underlying EJA.

In the asymptotic limit $t\to+\infty$, the cumulative losses incurred by SCMWU
become comparable to the losses of the best fixed action that the decision
maker can make in hindsight. This is commonly referred to as the “no-regret”
property. Here, we provide the regret bound for SCMWU presented in [13] for
the sake of comprehensiveness.

######  Theorem 1.

Let $(\mathcal{J},\circ)$ be an EJA of rank $r$, $\mathcal{K}$ be its cone of
squares, and ${\bm{x}}\bullet{\bm{y}}=\mathbf{Tr}({\bm{x}}\circ{\bm{y}})$ be
the EJA inner product. For any $\eta\in(0,1]$ and any sequence of loss vectors
${{\bm{m}}}^{(t)}$ satisfy $\|{{\bm{m}}}^{(t)}\|_{\infty}\leq 1$, the iterates
${{\bm{p}}}^{(t)}$ generated from (SCMWU) satisfy

$\sum_{t=1}^{T}{{\bm{m}}}^{(t)}\bullet{{\bm{p}}}^{(t)}\leq\lambda_{\min}({\sum_{t=1}^{T}{{\bm{m}}}^{(t)}})+\eta
T+\frac{\ln r}{\eta},$

where $\lambda_{\min}({\cdot})$ is the minimum eigenvalue of an EJA vector.

###### Proof.

Let
${{\bm{w}}}^{(t)}=\mathbf{exp}(-\eta\sum_{\tau=1}^{t-1}{{\bm{m}}}^{(\tau)})$,
${{\bm{w}}}^{(1)}={\bm{e}}$ and
${{\bm{p}}}^{(t)}=\frac{{{\bm{w}}}^{(t)}}{\mathbf{Tr}({{\bm{w}}}^{(t)})}$. We
can establish an upper bound on $\mathbf{Tr}({{\bm{w}}}^{(T+1)})$ via the
following steps:

$\displaystyle\mathbf{Tr}({{\bm{w}}}^{(T+1)})$
$\displaystyle=\mathbf{Tr}(\mathbf{exp}(-\eta\sum_{t=1}^{T}{{\bm{m}}}^{(t)}))$
(3)
$\displaystyle\overset{(a)}{\leq}\mathbf{Tr}({{\bm{w}}}^{(T)}\circ\mathbf{exp}(-\eta{{\bm{m}}}^{(T)}))$
$\displaystyle={{\bm{w}}}^{(T)}\bullet\mathbf{exp}(-\eta{{\bm{m}}}^{(T)})$
$\displaystyle\overset{(b)}{\leq}{{\bm{w}}}^{(T)}\bullet({\bm{e}}-\eta{{\bm{m}}}^{(T)}+\eta^{2}({{\bm{m}}}^{(T)})^{2})$
$\displaystyle=\mathbf{Tr}({{\bm{w}}}^{(T)})-\eta{{\bm{m}}}^{(T)}\bullet{{\bm{w}}}^{(T)}+\eta^{2}({{\bm{m}}}^{(T)})^{2}\bullet{{\bm{w}}}^{(T)}$
$\displaystyle=\mathbf{Tr}({{\bm{w}}}^{(T)})\cdot(1-\eta{{\bm{m}}}^{(T)}\bullet{{\bm{p}}}^{(T)}+\eta^{2}({{\bm{m}}}^{(T)})^{2}\bullet{{\bm{p}}}^{(T)})$
$\displaystyle\overset{(c)}{\leq}\mathbf{Tr}({{\bm{w}}}^{(T)})\cdot\exp(-\eta{{\bm{m}}}^{(T)}\bullet{{\bm{p}}}^{(T)}+\eta^{2}({{\bm{m}}}^{(T)})^{2}\bullet{{\bm{p}}}^{(T)})$
$\displaystyle\overset{(d)}{\leq}r\cdot\exp(-\eta\sum_{t=1}^{T}{{\bm{m}}}^{(t)}\bullet{{\bm{p}}}^{(t)}+\eta^{2}\sum_{t=1}^{T}({{\bm{m}}}^{(t)})^{2}\bullet{{\bm{p}}}^{(t)}).$

Here $(a)$ is due to the generalized Golden-Thompson inequality (see (1)),
$(b)$ is due to the self-duality of the symmetric cone and the fact that
$e^{x}\leq 1+x+x^{2}\ \forall x[-1,1]$, $(c)$ is because $1+x\leq e^{x}$, and
$(d)$ is by induction. On the other hand, a lower bound on
$\mathbf{Tr}({{\bm{w}}}^{(T+1)})$ can also be established:

$\displaystyle\mathbf{Tr}({{\bm{w}}}^{(T+1)})=\mathbf{Tr}(\mathbf{exp}(-\eta\sum_{t=1}^{T}{{\bm{m}}}^{(t)}))=\sum_{k=1}^{r}\lambda_{k}({\mathbf{exp}(-\eta\sum_{t=1}^{T}{{\bm{m}}}^{(t)})})$
(4)
$\displaystyle=\sum_{k=1}^{r}\exp(\lambda_{k}({-\eta\sum_{t=1}^{T}{{\bm{m}}}^{(t)}}))\geq\exp(\lambda_{\max}({-\eta\sum_{t=1}^{T}{{\bm{m}}}^{(t)}}))=\exp(-\eta\lambda_{\min}({\sum_{t=1}^{T}{{\bm{m}}}^{(t)}})),$

where $\lambda_{k}({\cdot})$ represents the $k$-th eigenvalue. Combining (3)
and (4) gives

$\displaystyle\exp(-\eta\lambda_{\min}({\sum_{t=1}^{T}{{\bm{m}}}^{(t)}}))$
$\displaystyle\leq
r\cdot\exp(-\eta\sum_{t=1}^{T}{{\bm{m}}}^{(t)}\bullet{{\bm{p}}}^{(t)}+\eta^{2}\sum_{t=1}^{T}({{\bm{m}}}^{(t)})^{2}\bullet{{\bm{p}}}^{(t)})$
(5) $\displaystyle\leq
r\cdot\exp(-\eta\sum_{t=1}^{T}{{\bm{m}}}^{(t)}\bullet{{\bm{p}}}^{(t)}+\eta^{2}T).$

The last inequality is because
$({{\bm{m}}}^{(t)})^{2}\bullet{{\bm{p}}}^{(t)}\leq 1$, since
$-{\bm{e}}\preceq_{\mathcal{K}}{{\bm{m}}}^{(t)}\preceq_{\mathcal{K}}{\bm{e}}$
and $\mathbf{Tr}({{\bm{p}}}^{(t)})=1$. Taking logarithms on both sides of (5),
scaling by $\frac{1}{\eta}$, and rearranging gives the desired result. ∎

###  Parallel computation

The work-span (or work-depth) model is a framework for conceptualizing and
describing parallel algorithms. In this model, a computation can be viewed as
a directed acyclic graph (DAG). The work $W$ of a parallel algorithm is the
total number of operations, and the span $S$ is the length of the longest path
in the DAG. Given $P$ processors with shared memory, the computation can be
executed in $W/P+O(S)$ time with high probability [3, 12]. Since $P$ is often
asymptotically much smaller than $W$, the running time is often dominated by
the total work of the algorithm. On the other hand, the span measures the
ability of a parallel algorithm to effectively utilize more processors to
improve the performance [11].

In this work, we will use the work-span model to analyze our algorithms under
parallel settings. A primitive problem that frequently arises as a subproblem
when designing our parallel algorithms is the information aggregation problem.
Given an array $(x_{1},\dots,x_{n})$ of length $n$ and an associative operator
$\diamond$, we want to compute $x_{1}\diamond\dots\diamond x_{n}$. Let
$\mathcal{I}$ be the set of odd indices in the array. For each
$i\in\mathcal{I}$, we compute $y_{\lceil i/2\rceil}=x_{i}\diamond x_{i+1}$ (if
$n$ is odd, we let $y_{\lceil n/2\rceil}=x_{n}$). Then we get a new array $y$
with length $\lceil n/2\rceil$ and the problem becomes computing
$y_{1}\diamond\dots\diamond y_{\lceil n/2\rceil}$. The problem can be solved
recursively using the same approach until the length of the array is reduced
to 1. The associated computational DAG of this process form a binary tree,
where the total number of operations is $O(n)$ and the longest path in the DAG
(the height of the binary tree) is $O(\log n)$. Therefore, under the setting
of the work-span model, the information aggregation problem can be solved in
$O(n)$ work and $O(\log n)$ span.

The information aggregation problem covers many basic operations such as
summation and extremum, which are building blocks for our algorithms. Here we
elucidate the process of parallelizing the SCMWU method within the work-span
model, given its frequent status as the most computationally intensive
component in our algorithms. Let $\mathcal{K}$ be the Cartesian product of
several primitive symmetric cones and let $\mathcal{J}$ be the underlying EJA.
The computation of the exponential of an element ${\bm{x}}\in\mathcal{J}$ can
be decomposed into problems of computing the exponentials of each of its
components. Let $N$ be the number of non-zero entries in
$\mathbf{exp}({\bm{x}})$ and let $r$ be the rank of $\mathcal{J}$. Given the
eigenvalues (or the trace) of each component of $\mathbf{exp}({\bm{x}})$, the
computation of $\mathbf{Tr}(\mathbf{exp}({\bm{x}}))$ is simply a summation
problem that is solvable in $O(r)$ work and $O(\log r)$ span, and consequently
the normalization of $\mathbf{exp}({\bm{x}})$ can be performed in $O(N)$ work
and $O(1)$ span. By the classification theorem of symmetric cones, any
symmetric cone over the real field is isomorphic to the Cartesian product of
PSD cones and second-order cones. Thus to analyze the computation of
$\mathbf{exp}({\bm{x}})$, it only remains analyzing the exponentials of
elements in the spaces of PSD cones and second-order cones. Let
$\mathcal{Q}^{d+1}$ be a second-order cone and let $({\bm{u}};u_{0})$ be an
element in the underlying space $\mathbbm{R}^{d+1}$. The spectral
decomposition of $({\bm{u}};u_{0})$ is given by
$\lambda_{1,2}={1\over\sqrt{2}}(u_{0}\pm\|{\bm{u}}\|_{2})$ and
${\bm{q}}_{1,2}={1\over\sqrt{2}}(\pm{\bm{v}};1)$, where
${\bm{v}}=\frac{{\bm{u}}}{\|{\bm{u}}\|_{2}}$ or any unit vector if
${\bm{u}}=0$. If ${\bm{u}}\neq 0$, To compute the eigenvalues and
eigenvectors, the most intensive part is to compute $\|{\bm{u}}\|_{2}$, which
is again a summation problem with some simple additional arithmetic and can be
computed in $O(d)$ work and $O(\log d)$ span. Once the spectral decomposition
is obtained, the exponential can be directly computed in $O(d)$ work and
$O(1)$ span. In summary, an SCMWU step for $\mathcal{Q}^{d+1}$ can be
performed with $O(d)$ work and $O(\log d)$ span. In [6, Sec. 7], Arora and
Kale proposed an approximation algorithm for computing the exponential of
matrices that lie in the space of the PSD cone. Their algorithm utilizes the
power-series expansion of the exponential and the Johnson-Lindenstrauss lemma
and only involves matrix multiplications, which can be reduced to a series of
inner products and summations and can be directly parallelized within the
work-span model. Given a $d\times d$ matrix, Arora and Kale’s algorithm can be
parallelized in ${\sf poly}(d)$ work and ${\sf polylog}(d)$ span.

##  Primal-dual meta algorithm for SCPs

In this section, and throughout this work, we consider a pair of primal/dual
SCPs of the form

$\begin{aligned}
\max\quad&{\bm{c}}\bullet{\bm{x}}+{\bm{c}}^{\prime}\bullet{\bm{x}}^{\prime}\\\
\operatorname*{s.t.}\quad&{\bm{a}}_{j}\bullet{\bm{x}}+{\bm{a}}_{j}^{\prime}\bullet{\bm{x}}^{\prime}=b_{j}\quad\forall
j\in[m]\\\ &{\bm{x}}\in\mathcal{K},\
{\bm{x}}^{\prime}\in\mathcal{K}^{\prime}\end{aligned}\hskip
30.00005pt\begin{aligned} \min\quad&{\bm{b}}^{T}{\bm{y}}\\\
\operatorname*{s.t.}\quad&\textstyle\sum_{j=1}^{m}{\bm{a}}_{j}y_{j}\succeq_{\mathcal{K}}{\bm{c}}\\\
&\textstyle\sum_{j=1}^{m}{\bm{a}}^{\prime}_{j}y_{j}\succeq_{\mathcal{K}^{\prime}}{\bm{c}}^{\prime}\end{aligned}$
(SCP)

where $\mathcal{K},\mathcal{K}^{\prime}$ are symmetric cones,
${\bm{c}},{\bm{a}}_{j}\in\mathcal{J}$,
${\bm{c}}^{\prime},{\bm{a}}_{j}^{\prime}\in\mathcal{J}^{\prime}$, and
$\mathcal{J},\mathcal{J}^{\prime}$ are the EJAs associated with the cones
$\mathcal{K},\mathcal{K}^{\prime}$ respectively. We assume throughout that
(SCP) exhibit strong duality (e.g. see [31]), and use ${\sf OPT}$ to denote
the shared optimal value of the primal and dual problem.

Notice that the dual feasible region is defined by two generalized inequality
constraints,

$\mathcal{F}\ {\triangleq}\
\\{\mathbf{y}\in\mathbb{R}^{m}:\sum_{j=1}^{m}\mathbf{a}_{j}y_{j}\succeq_{\mathcal{K}}\mathbf{c}\\}\
\text{ and }\ \mathcal{C}\ {\triangleq}\
\\{\mathbf{y}\in\mathbb{R}^{m}:\sum_{j=1}^{m}\mathbf{a}^{\prime}_{j}y_{j}\succeq_{\mathcal{K}^{\prime}}\mathbf{c}^{\prime}\\},$
(6)

where $\mathcal{F}$ is a constraint set that is “hard to deal with” while
$\mathcal{C}$ is a constraint set that is “easy to satisfy”. In our algorithm,
the hard constraint set $\mathcal{F}$ is relaxed to a halfspace while for the
easy part we will need to to perform linear optimization over the intersection
of the set $\mathcal{C}$ and a given halfspace.

As a concrete example, consider the case where
$\mathcal{K}=\mathbbm{R}^{n}_{+}$, $\mathcal{K}^{\prime}=\mathbb{R}^{m}_{+}$,
${\bm{a}}^{\prime}_{j}\in\mathbb{R}^{m}$ is the $j$-th unit vector, and
${\bm{c}}^{\prime}={\bm{0}}$. Thehe pair of primal/dual SCPs is then just a
pair of primal/dual LPs:

$\max\\{{\bm{c}}^{T}{\bm{x}}:{\bm{A}}{\bm{x}}\leq{\bm{b}},\ {\bm{x}}\geq
0\\}\text{ and }\min\\{{\bm{b}}^{T}{\bm{y}}:{\bm{A}}^{T}{\bm{y}}-{\bm{c}}\geq
0,\ {\bm{y}}\geq 0\\}.$

Here, the dual feasible region is the intersection of a polyhedron
$\mathcal{F}=\\{{\bm{y}}:{\bm{A}}^{T}{\bm{y}}-{\bm{c}}\geq 0\\}$ corresponding
to the hard set and the non-negative orthant $\mathcal{C}$, which corresponds
to the easy set. Note that optimizing a linear function over the intersection
of the non-negative orthant and a halfspace can be accomplished easily.

As a second example, we consider the case where
$\mathcal{K}=\mathcal{S}^{n}_{+}$, $\mathcal{K}^{\prime}=\mathbb{R}^{m}_{+}$,
${\bm{a}}^{\prime}_{j}\in\mathbb{R}^{m}$ is the $j$-th unit vector, and
${\bm{c}}^{\prime}={\bm{0}}$. In this setting, the pair of primal/dual SCPs
defined above correspond to a pair of primal/dual SDPs:

$\max\\{{\bm{c}}\bullet{\bm{x}}:{\bm{a}}_{j}\bullet{\bm{x}}\leq b_{j}\ \forall
j\in[m],\ {\bm{x}}\succeq 0\\}\text{ and
}\min\\{{\bm{b}}^{T}{\bm{y}}:\textstyle\sum_{j=1}^{m}{\bm{a}}_{j}y_{j}\succeq{\bm{c}},\
{\bm{y}}\geq 0\\}.$

In this setting the hard part of the dual feasible region is carved out by a
linear matrix inequality, i.e.
$\mathcal{F}=\\{{\bm{y}}:\sum_{j}{\bm{a}}_{j}y_{j}\succeq{\bm{c}}\\}$, while
the easy part is again the non-negative orthant. This is exactly the setting
that Arora and Kale consider in their primal-dual framework for SDPs [6].

Finally, our analysis necessitates that the SCPs adhere to the following
assumption:

$\exists j\in[m]\text{ such that }{\bm{a}}_{j}={\bm{e}}\text{ and
}{\bm{a}}^{\prime}_{j}\in\mathcal{K}^{\prime}.$

This property allows us to “project” an approximately feasible dual solution
to a strict one. For simplicity, we assume that this assumption is true for
the last primal constraint, i.e.

${\bm{a}}_{m}={\bm{e}}\ \text{ and }\
{\bm{a}}^{\prime}_{m}\in\mathcal{K}^{\prime}.$ (7)

As a canonical example, this assumption holds when we have a bound $R$ on the
trace of primal solutions. Indeed, this case is captured by setting
${\bm{a}}_{m}={\bm{e}}$, ${\bm{a}}^{\prime}_{m}\in\mathcal{K}^{\prime}$, and
$b_{m}=R$.

In order to solve the SCP problem we reduce it to a feasibility problem, i.e.
we consider the following problem: can the dual problem in (SCP) achieve a
specified objective value $\alpha$? More specifically, does
$(1+\varepsilon)\alpha\geq{\sf OPT}$ or $\alpha<{\sf OPT}$? In addition to the
answer, for the former case we also want to find a dual feasible solution to
(SCP) with objective value at most $(1+\varepsilon)\alpha$, and for the latter
case a primal feasible solution with objective value at least $\alpha$. In the
following, the decision problem, together with the problem of finding the
primal/dual solutions, is referred to as the $\alpha$-feasibility test
problem, $\alpha$-FTP. Given an algorithm that solves $\alpha$-FTP, the SCP
problem can then be solved by employing a binary search (or other searching
methods) to narrow the range of $\alpha$. Therefore, in this section, we
discuss mainly the algorithm for solving $\alpha$-FTP.

Define
$\mathcal{S}_{\alpha}=\\{{\bm{y}}\in\mathbbm{R}^{m}:{\bm{b}}^{T}{\bm{y}}\leq\alpha\\}$
to be the $\alpha$-sublevel set of the dual objective function. Recalling the
definitions of $\mathcal{C}$ and $\mathcal{F}$ (see (6)) we have that

$\mathcal{F}\cap\mathcal{C}\cap\mathcal{S}_{\alpha}\neq\varnothing\iff{\sf
OPT}\leq\alpha.$

Indeed, if $\alpha\geq{\sf OPT}$, the dual optimal solution lies in
$\mathcal{F}\cap\mathcal{C}\cap\mathcal{S}_{\alpha}$, so it is nonempty.
Conversely, if $\alpha<{\sf OPT}$, all feasible solutions have value
$>\alpha$, so we have
$\mathcal{F}\cap\mathcal{C}\cap\mathcal{S}_{\alpha}=\varnothing$. Therefore,
$\alpha$-FTP is equivalent to determine whether the intersection between
$\mathcal{F}$ and $\mathcal{C}\cap\mathcal{S}_{\alpha}$ is empty or not. Since
$\mathcal{F}$, $\mathcal{C}$ and $\mathcal{S}_{\alpha}$ are convex sets,
$\mathcal{F}\cap\mathcal{C}\cap\mathcal{S}_{\alpha}=\varnothing$ is equivalent
to saying that the hard constraint set $\mathcal{F}$ can be separated from
$\mathcal{C}\cap\mathcal{S}_{\alpha}$, which can be certified by a separating
hyperplane. Our meta-algorithm is designed to discover a separating
hyperplane, thereby certifying that $\alpha<{\sf OPT}$. To this end, we relax
the hard constraint set $\mathcal{F}$ to a halfspace $\mathcal{H}_{+}$ and
test whether the intersection
$\mathcal{H}_{+}\cap\mathcal{C}\cap\mathcal{S}_{\alpha}$ is empty. If this
intersection is empty, then $\mathcal{H}$ serves as a certificate that
$\alpha<{\sf OPT}$. As long as the intersection is non-empty, we adjst the
halfspace to ensure it contains the hard constraint set $\mathcal{F}$ while
also endeavoring to separate it from the average of all previous iterates. If
we succeed to do this for enough times, we will find an approximate feasible
solution. The high-level outline of the algorithm is presented in the boxed
environment below.

Primal-dual meta-algorithm for $\alpha$-FTP Let ${\mathcal{H}}^{(1)}$ be a
hyperplane such that $\mathcal{F}$ is contained in the closed positive
halfspace, ${\mathcal{H}}^{(1)}_{+}$. For each iteration $t$, we test whether
the intersection of ${\mathcal{H}}^{(t)}_{+}$ with
$\mathcal{C}\cap\mathcal{S}_{\alpha}$ is empty: • If yes, conclude that
$\alpha<{\sf OPT}$, find a primal feasible solution to (SCP) with objective
value at least $\alpha$, and terminate. • If no, find a point
${{\bm{y}}}^{(t)}\in{\mathcal{H}}^{(t)}_{+}\cap\mathcal{C}\cap\mathcal{S}_{\alpha}$.
Let ${\bar{{\bm{y}}}}^{(t)}=\frac{1}{t}\sum_{\tau=1}^{t}{{\bm{y}}}^{(\tau)}$,
which is the mean of the past iterates. Adjust ${\mathcal{H}}^{(t)}$ to get a
new hyperplane ${\mathcal{H}}^{(t+1)}$ (maintaining
$\mathcal{F}\subseteq{\mathcal{H}}^{(t+1)}_{+}$), with the goal of separating
${\bar{{\bm{y}}}}^{(t)}$ from $\mathcal{F}$. If the algorithm does not
terminate within a sufficiently large number of iterations, output a dual
feasible solution to (SCP) with objective value at most
$(1+\varepsilon)\alpha$.

Figure 1: Illustration of the process of the meta-algorithm, where
$\alpha_{1}>{\sf OPT}>\alpha_{2}$.

Figure 1 illustrates three iterations (in two settings) of our meta-algorithm
in both non-separable and separable cases. In the former case (the top row),
where $\mathcal{F}\cap\mathcal{C}\cap\mathcal{S}_{\alpha_{1}}\not=\varnothing$
and $\alpha_{1}>{\sf OPT}$, our algorithm would not find a separating
hyperplane, and would thus continue till $t$ reaches some large number and
then output a dual feasible solution. In the latter case, indicated in the
bottom row, the algorithm identifies a separating hyperplane, concluding that
$\alpha_{2}<{\sf OPT}$.

We designate this approach as a meta-algorithm because it contains several
degrees of freedom that must be specified to develop a concrete algorithm for
the $\alpha$-feasibility test problem. We will now discuss how each of the
steps within the meta-algorithmic framework can be implemented.

Parameterizing the hyperplanes ${\mathcal{H}}^{(t)}$. For each iteration $t$,
we need to find a halfspace ${\mathcal{H}}^{(t)}_{+}$ satisfying
$\mathcal{F}\subseteq{\mathcal{H}}^{(t)}_{+}$. Recall that $\mathcal{F}$ is
defined by the generalized inequality
$\sum_{j}{\bm{a}}_{j}y_{j}-{\bm{c}}\succeq_{\mathcal{K}}{\bm{0}}$. A
generalized inequality is equivalent to an infinite number of halfspace
constraints, i.e.

$\sum_{j}{\bm{a}}_{j}y_{j}-{\bm{c}}\succeq_{\mathcal{K}}{\bm{0}}\iff(\sum_{j}{\bm{a}}_{j}y_{j}-{\bm{c}})\bullet{\bm{p}}\geq
0,\ \forall{\bm{p}}\in\mathcal{K}\text{ with }\mathbf{Tr}({\bm{p}})=1,$

by the self-duality and homogeneous properties of symmetric cones. Thus, it is
a natural choice to consider hyperplanes ${\mathcal{H}}^{(t)}$ that are
parametrized by trace-one vectors ${{\bm{p}}}^{(t)}\in\mathcal{K}$. That is,
we define

${\mathcal{H}}^{(t)}=\\{{\bm{y}}\in\mathbbm{R}^{m}:(\textstyle\sum_{j}{\bm{a}}_{j}y_{j}-{\bm{c}})\bullet{{\bm{p}}}^{(t)}=0\\},$

and $\mathcal{F}$ indeed lies in the closed positive halfspace
${\mathcal{H}}^{(t)}_{+}=\\{{\bm{y}}\in\mathbbm{R}^{m}:(\sum_{j}{\bm{a}}_{j}y_{j}-{\bm{c}})\bullet{{\bm{p}}}^{(t)}\geq
0\\}$.

Testing of separation. Given a halfspace ${\mathcal{H}}^{(t)}_{+}$ satisfying
$\mathcal{F}\subseteq{\mathcal{H}}^{(t)}_{+}$, we try to determine whether the
intersection of ${\mathcal{H}}^{(t)}_{+}$ and
$\mathcal{C}\cap\mathcal{S}_{\alpha}$ is empty. If yes, then $\mathcal{F}$,
which is contained in ${\mathcal{H}}^{(t)}_{+}$, also has no intersection with
$\mathcal{C}\cap\mathcal{S}_{\alpha}$ and we can conclude that $\alpha<{\sf
OPT}$. On the other hand, if
${\mathcal{H}}^{(t)}_{+}\cap\mathcal{C}\cap\mathcal{S}_{\alpha}\neq\varnothing$,
we also want to output a point ${{\bm{y}}}^{(t)}$ inside the region as a
witness to this fact. As this decision problem is the core of the meta-
algorithm, we abstract it as an ORACLE:

Given ${{\bm{p}}}^{(t)}\in\mathcal{K}$ with $\mathbf{Tr}({{\bm{p}}}^{(t)})=1$,
determine if there exists ${\bm{y}}\in\mathbbm{R}^{m}$ such that
$(\textstyle\sum_{j}{\bm{a}}_{j}y_{j}-{\bm{c}})\bullet{{\bm{p}}}^{(t)}\geq 0,\
{\bm{b}}^{T}{\bm{y}}\leq\alpha\text{ and }{\bm{y}}\in\mathcal{C}.$

The ORACLE can be implemented in a variety of ways. One possible
implementation of the ORACLE is to solve the following convex problem (with
${{\bm{p}}}^{(t)}$ as the input at iteration $t$):

$\displaystyle\max_{\bm{y}}$
$\displaystyle\textstyle\sum_{j}({\bm{a}}_{j}\bullet{{\bm{p}}}^{(t)})y_{j}-{\bm{c}}\bullet{{\bm{p}}}^{(t)}$
(8) $\displaystyle\operatorname*{s.t.}$
$\displaystyle{\bm{b}}^{T}{\bm{y}}\leq\alpha,\quad{\bm{y}}\in\mathcal{C}.$

If the optimal value is negative, this means we have
$(\sum_{j}{\bm{a}}_{j}{y}_{j}-{\bm{c}})\bullet{{\bm{p}}}^{(t)}<0$ for all the
points ${\bm{y}}\in\mathcal{C}\cap\mathcal{S}_{\alpha}$, in which case the
entire region of $\mathcal{C}\cap\mathcal{S}_{\alpha}$ is in the open negative
halfspace ${\mathcal{H}}^{(t)}_{-}$ and separated from $\mathcal{F}$ in
${\mathcal{H}}^{(t)}_{+}$. Thus, the ORACLE can output affirmatively that
${\mathcal{H}}^{(t)}_{+}\cap\mathcal{C}\cap\mathcal{S}_{\alpha}$ is empty. On
the other hand, if the optimal value of (8) is non-negative and
${{\bm{y}}}^{(t)}$ is an optimal solution, we clearly have that
${{\bm{y}}}^{(t)}\in{\mathcal{H}}^{(t)}_{+}\cap\mathcal{C}\cap\mathcal{S}_{\alpha}$.

The complexity of each iteration of our algorithm depends on the complexity of
implementing the ORACLE, which in turn depends on the “complexity” of the set
$\mathcal{C}$. Additionally, let
$\mathcal{D}=\\{{\bm{p}}\in\mathcal{K}:\mathbf{Tr}({\bm{p}})=1\\}$ be the set
of all possible inputs to the ORACLE. We assume that for any input
${\bm{p}}\in\mathcal{D}$, the output ${\bm{y}}$ of the ORACLE satisfies
$\|\sum_{j}{\bm{a}}_{j}y_{j}-{\bm{c}}\|_{\infty}\leq\rho$. The quantity $\rho$
is designated as the width of the ORACLE. The number of iterations of our
algorithm will depend on the width of the ORACLE as will be revealed in what
follows. These two considerations on the complexity of implementing the ORACLE
and its width offer insights into the properties that the constraint set
$\mathcal{C}$ should possess to be considered easy.

Finding a primal solution. Consider the pair of primal/dual optimization
problems:

$\begin{aligned} \min_{\bm{y}}\quad&{\bm{b}}^{T}{\bm{y}}\\\
\operatorname*{s.t.}\quad&\textstyle\sum_{j}({\bm{a}}_{j}\bullet{{\bm{p}}}^{(t)})y_{j}\geq{\bm{c}}\bullet{{\bm{p}}}^{(t)}\\\
&{\bm{y}}\in\mathcal{C}\end{aligned}\qquad\begin{aligned}
\max_{z,{\bm{x}}^{\prime}}\quad&({\bm{c}}\bullet{{\bm{p}}}^{(t)})z+{\bm{c}}^{\prime}\bullet{\bm{x}}^{\prime}\\\
\operatorname*{s.t.}\quad&({\bm{a}}_{j}\bullet{{\bm{p}}}^{(t)})z+{\bm{a}}_{j}^{\prime}\bullet{\bm{x}}^{\prime}=b_{j}\quad\forall
j\in[m]\\\ &z\geq 0,\ {\bm{x}}^{\prime}\in\mathcal{K}^{\prime}\end{aligned}$
(9)

denoted P and D, respectively. Clearly, if the optimal value of (8) is
negative (i.e. when the ORACLE outputs affirmatively), the value of P is
greater than $\alpha$. Then, by strong duality, the value of D is also greater
than $\alpha$. Let $(\tilde{z};\tilde{{\bm{x}}}^{\prime})$ be the optimal
solution of D. Since $\tilde{z}\geq 0$ and ${{\bm{p}}}^{(t)}\in\mathcal{K}$,
$(\tilde{z}{{\bm{p}}}^{(t)};\tilde{{\bm{x}}}^{\prime})$ is a primal feasible
solution to (SCP) with value at least $\alpha$.

Adjusting the hyperplane ${\mathcal{H}}^{(t)}$. Consider the case when
${\mathcal{H}}^{(t)}_{+}$ intersects $\mathcal{C}\cap\mathcal{S}_{\alpha}$ and
the ORACLE returns ${{\bm{y}}}^{(t)}\in\mathcal{C}\cap\mathcal{S}_{\alpha}$.
Then, we need to adjust ${\mathcal{H}}^{(t)}$ to a new hyperplane
${\mathcal{H}}^{(t+1)}$. This is done by using the symmetric cone
multiplicative weights update (SCMWU) method introduced in Section 2.3 to
update the vector ${{\bm{p}}}^{(t)}\in\mathcal{K}$ that parameterizes
${\mathcal{H}}^{(t)}$. Specifically, we define the loss vectors
${{\bm{m}}}^{(\tau)}=\frac{1}{\rho}(\sum_{j}{\bm{a}}_{j}{y}^{(\tau)}_{j}-{\bm{c}})$
for $\tau=1,\cdots,t$, where ${{\bm{y}}}^{(\tau)}$ is the output of the ORACLE
at the $\tau$-th iteration and $\rho$ is the width of the ORACLE. The next
vector ${{\bm{p}}}^{(t+1)}$ is then obtained from (SCMWU). Note that since
$\sum_{\tau}^{t}{{\bm{m}}}^{(\tau)}$ equals
$\frac{t}{\rho}(\sum_{j}{\bm{a}}_{j}{\bar{y}}^{(t)}_{j}-{\bm{c}})$, where
${\bar{{\bm{y}}}}^{(t)}$ is the mean of the past points, the next vector
${{\bm{p}}}^{(t+1)}$ can be computed by:

${{\bm{p}}}^{(t+1)}=\frac{1}{C}\cdot\mathbf{exp}\Big{(}{-\eta
t\over\rho}(\sum_{j}{\bm{a}}_{j}{\bar{y}}^{(t)}_{j}-{\bm{c}})\Big{)},$ (10)

where $C$ normalizes ${{\bm{p}}}^{(t+1)}$ to have trace one.

We have two notes to make here. First, for any ${\bm{y}}\in\mathcal{F}$, we
have $(\sum_{j}{\bm{a}}_{j}{\bm{y}}-{\bm{c}})\bullet{{\bm{p}}}^{(t+1)}\geq 0$
and thus $\mathcal{F}$ is still in ${\mathcal{H}}^{(t+1)}_{+}$. Second, here
is the intuition on how the new hyperplane ${\mathcal{H}}^{(t+1)}$ endeavors
to separate ${\bar{{\bm{y}}}}^{(t)}$ and $\mathcal{F}$: since the vectors
$\sum_{j}{\bm{a}}_{j}{\bar{y}}^{(t)}_{j}-{\bm{c}}$ and ${{\bm{p}}}^{(t+1)}$
commute, let $\lambda_{k}$ be the eigenvalues of
$\sum_{j}{\bm{a}}_{j}{\bar{y}}^{(t)}_{j}-{\bm{c}}$, we have:

$(\sum_{j}{\bm{a}}_{j}{\bar{y}}^{(t)}_{j}-{\bm{c}})\bullet{{\bm{p}}}^{(t+1)}=\frac{1}{C}\sum_{k}\lambda_{k}\exp({-\eta
t\lambda_{k}\over\rho}).$

Hence, on the right-hand side above, significantly negative eigenvalues are
given large positive weights in an effort to make the above summation
negative. From a geometric perspective, the adjustment endeavors to place
${\bar{{\bm{y}}}}^{(t)}$ in the open negative halfspace
${\mathcal{H}}^{(t+1)}_{-}\ \triangleq\
\\{{\bm{y}}\in\mathbbm{R}^{m}:(\sum_{j}{\bm{a}}_{j}y_{j}-{\bm{c}})\bullet{{\bm{p}}}^{(t+1)}<0\\}$
so that it is separated from $\mathcal{F}$ in ${\mathcal{H}}^{(t+1)}_{+}$.

Having outlined all the steps in the meta-algorithm, we are now ready to
present our central theorem, specifying the required number of steps to reach
a conclusion regarding whether OPT surpasses or falls short of our current
estimate $\alpha$.

######  Theorem 2.

Consider a pair of primal/dual SCPs over the symmetric cone
$\mathcal{K}\operatorname*{\operatorname*{\oplus}}\mathcal{K}^{\prime}$, where
the rank of $\mathcal{K}$ is $r$, and assumption (7) holds. Fix an error
parameter $\varepsilon>0$ and a guess $\alpha>0$ for the optimal value OPT.
Execute the primal-dual meta-algorithm for
$T=\left\lceil\frac{4R^{2}\rho^{2}\ln
r}{\varepsilon^{2}\alpha^{2}}\right\rceil$ iterations, adjusting the
hyperplane using SCMWU (as defined in (10)) with step size
$\eta=\sqrt{\frac{\ln r}{T}}$ and an ORACLE with width $\rho.$ Then, we either
correctly conclude that ${\sf OPT}>\alpha$, or, if the algorithm does not
terminate within $T$ iterations we conclude that ${\sf
OPT}\leq(1+\varepsilon)\alpha$ and the vector

$\tilde{{\bm{y}}}=\left({\bar{y}}^{(T)}_{1},\ldots,{\bar{y}}^{(T)}_{m-1},{\bar{y}}^{(T)}_{m}+\frac{\varepsilon\alpha}{R}\right)$

is a dual feasible solution to SCP with objective value
${\bm{b}}^{T}\tilde{{\bm{y}}}\leq(1+\varepsilon)\alpha$.

###### Proof.

Suppose that the algorithm does not terminate after $t=1,\ldots,T$ iterations.
By our choice of ${{\bm{y}}}^{(t)}$ and
${{\bm{m}}}^{(t)}=\frac{1}{\rho}(\sum_{j}{\bm{a}}_{j}{y}^{(t)}_{j}-{\bm{c}})$,
the loss ${{\bm{m}}}^{(t)}\bullet{{\bm{p}}}^{(t)}$ for SCMWU is always non-
negative. By the regret bound of SCMWU (Theorem 1), we have:

$0\leq\sum_{t=1}^{T}\frac{1}{\rho}(\sum_{j}{\bm{a}}_{j}{y}^{(t)}_{j}-{\bm{c}})\bullet{{\bm{p}}}^{(t)}\leq\lambda_{\min}({\sum_{t=1}^{T}\frac{1}{\rho}(\sum_{j}{\bm{a}}_{j}{y}^{(t)}_{j}-{\bm{c}})})+\eta
T+\frac{\ln r}{\eta}.$

Then, scaling both sides of the inequality by $\frac{\rho}{T}$, we have

$0\leq\lambda_{\min}({\sum_{j}{\bm{a}}_{j}{\bar{y}}^{(T)}_{j}-{\bm{c}}})+\eta\rho+\frac{\rho\ln
r}{\eta T},$ (11)

and substituting $T=\left\lceil\frac{4R^{2}\rho^{2}\ln
r}{\varepsilon^{2}\alpha^{2}}\right\rceil$ and $\eta=\sqrt{\frac{\ln r}{T}}$
into the above gives us

$\lambda_{\min}({\sum_{j}{\bm{a}}_{j}{\bar{y}}^{(T)}_{j}-{\bm{c}}})+\frac{\varepsilon\alpha}{R}\geq
0,$

or equivalently

$\sum_{j}{\bm{a}}_{j}{\bar{y}}^{(T)}_{j}\succeq_{\mathcal{K}}{\bm{c}}-\frac{\varepsilon\alpha}{R}{\bm{e}}.$
(12)

We now show that the vector
$\tilde{{\bm{y}}}=({\bar{y}}^{(T)}_{1},\ldots,{\bar{y}}^{(T)}_{m-1},{\bar{y}}^{(T)}_{m}+\frac{\varepsilon\alpha}{R})$
is dual feasible. First, to show that $\tilde{{\bm{y}}}\in\mathcal{F}$, we use
(12) and recall that by assumption (7) we have that ${\bm{a}}_{m}={\bm{e}}$.
Then, (12) becomes:

$\sum_{j=1}^{m-1}{\bm{a}}_{j}{\bar{y}}^{(T)}_{j}+({\bar{y}}^{(T)}_{m}+\frac{\varepsilon\alpha}{R}){\bm{e}}=\sum_{j=1}^{m}{\bm{a}}_{j}\tilde{y}_{j}\succeq_{\mathcal{K}}{\bm{c}}.$

Next, we show that $\tilde{{\bm{y}}}\in\mathcal{C}$. For this, by definition
of the ORACLE we have that ${{{\bm{y}}}}^{(t)}\in\mathcal{C}$, for all
$t=1,\ldots,T$, and by averaging we get that
${\bar{{\bm{y}}}}^{(T)}\in\mathcal{C}$, i.e.

$\sum_{j}{\bm{a}}_{j}^{\prime}{\bar{y}}^{(T)}_{j}\succeq_{\mathcal{K}^{\prime}}{\bm{c}}^{\prime}.$
(13)

Then, by assumption (7) we have that
${\bm{a}}_{m}^{\prime}\in\mathcal{K}^{\prime}$, and thus (13) implies that

$\sum_{j=1}^{m}{\bm{a}}_{j}^{\prime}\tilde{y}_{j}=\sum_{j=1}^{m-1}{\bm{a}}_{j}^{\prime}{\bar{y}}^{(T)}_{j}+{\bm{a}}_{m}^{\prime}({\bar{y}}^{(T)}_{m}+{\epsilon\alpha\over
R})\succeq_{\mathcal{K}^{\prime}}\sum_{j=1}^{m-1}{\bm{a}}_{j}^{\prime}{\bar{y}}^{(T)}_{j}+{\bm{a}}_{m}^{\prime}{\bar{y}}^{(T)}_{m}\succeq_{\mathcal{K}^{\prime}}{\bm{c}}^{\prime}.$

Finally, by assumption (7) we have that $b_{m}=R$. Thus, the objective value
of $\tilde{{\bm{y}}}$ satisfies:

${\bm{b}}^{T}\tilde{{\bm{y}}}={\bm{b}}^{T}{\bar{{\bm{y}}}}^{(T)}+R\frac{\varepsilon\alpha}{R}\leq\alpha+\varepsilon\alpha=(1+\varepsilon)\alpha,$

which completes the proof. ∎

Remarks on $\mathcal{C}$. The iteration complexity is quadratic in the width
$\rho$ of the ORACLE used. Thus, it is desirable to design an ORACLE with a
small width. As $\mathcal{C}$ is a constraint of the convex problem solved by
the ORACLE, $\rho$ in some way depends on the choice of $\mathcal{C}$.
Therefore, to reduce $\rho$, we may shrink the space of $\mathcal{C}$ by
adding extra constraints to the problem. These constraints should be redundant
to the original problem (SCP), but strong enough to restrict the result of (8)
so that ORACLE has a smaller $\rho$, as we will see in an example in Section
4.2. Moreover, if it is not required to find a primal solution (by solving
(9)), the choice of $\mathcal{C}$ can be more aggressive: we can replace the
set $\mathcal{C}$ by any convex set $\mathcal{C^{\prime}}\subseteq\mathcal{C}$
that contains the dual optimal solution ${\bm{y}}^{*}$, since whenever the
ORACLE confirms that
${\mathcal{H}}^{(t)}_{+}\cap\mathcal{C}^{\prime}\cap\mathcal{S}_{\alpha}$ is
empty we will have the correct outcome that ${\sf OPT}>\alpha$ (since
$y^{*}\in\mathcal{F}\cap\mathcal{C}^{\prime}\subseteq{\mathcal{H}}^{(t)}_{+}\cap\mathcal{C}^{\prime}$)
and $\mathcal{F}\cap\mathcal{C}\cap\mathcal{S}_{\alpha}=\varnothing$.

From feasibility to optimization. Assume an upper/lower bound for the value of
${\sf OPT}$. (Otherwise, we can employ an exponential search or a constant-
approximation algorithm to obtain one.) With our meta-algorithm, a
$(1+\varepsilon)$-approximate solution to (SCP) can then be computed using a
binary search. However, using the desired accuracy $\varepsilon$ as the error
tolerance in each binary search step results in an iteration complexity of
$O(T\log\frac{1}{\varepsilon})$. Instead of doing that, we can choose the
error tolerance $\varepsilon_{i}$ adaptively for each search step by setting
$\varepsilon_{i}\alpha_{i}$ to be within a constant fraction of the size of
its searching range. The total number of iterations is then dominated by the
last step (with the smallest $\varepsilon_{i}\alpha_{i}=\Theta(\varepsilon{\sf
OPT})$), and is improved to $O(T)$. Assume we know that ${\sf OPT}>0$ and lies
in the range $[L_{0},U_{0}]$. In the first step, we run our meta-algorithm,
with $\alpha_{1}=L_{0}+{1\over 3}(U_{0}-L_{0})$ and $\varepsilon_{1}={1\over
3\alpha_{1}}(U_{0}-L_{0})$. According to Theorem 2, we either correctly
conclude that this problem is infeasible, in which case ${\sf OPT}>\alpha_{1}$
and the search range becomes:

$L_{0}+{1\over 3}(U_{0}-L_{0})\leq{\sf OPT}\leq U_{0},$

or we conclude that ${\sf OPT}\leq(1+\varepsilon_{1})\alpha_{1}=L_{0}+{2\over
3}(U_{0}-L_{0})$ (and find a solution with this value), and the search range
becomes:

$L_{0}\leq{\sf OPT}\leq L_{0}+{2\over 3}(U_{0}-L_{0}).$

Setting $L_{1}=L_{0}+{1\over 3}(U_{0}-L_{0}),\ U_{1}=U_{0}$ if the first case
materializes and $L_{1}=L_{0},\ U_{1}=L_{0}+{2\over 3}(U_{0}-L_{0})$ if the
second case materializes. Both cases can be summarized as follows:

$L_{1}\leq{\sf OPT}\leq U_{1},\ \text{ where }(U_{1}-L_{1})={2\over
3}(U_{0}-L_{0}).$

Adopting this iterative refinement approach, at the $i$-th step we set
$\alpha_{i}=L_{i-1}+{1\over 3}(U_{i-1}-L_{i-1})$ and $\varepsilon_{i}={1\over
3\alpha_{i}}(U_{i-1}-L_{i-1})$. For the resultant interval $[L_{i},U_{i}]$ we
have

$(U_{i}-L_{i})=\Big{(}{2\over 3}\Big{)}^{i}(U_{0}-L_{0}).$

Thus, in $O\big{(}\log(\frac{U_{0}-L_{0}}{\varepsilon{\sf OPT}})\big{)}$
iterations the range size $(U_{i}-L_{i})$ will be less than the desired
(additive) error $\varepsilon{\sf OPT}$. The very last dual solution will have
value at most $(1+\varepsilon){\sf OPT}$, and a scaled version of the
${{\bm{p}}}^{(t)}$ that parameterized the last separating hyperplane is a
primal solution with value at least $(1-\varepsilon){\sf OPT}$. At each step
of the binary search, we invoke Theorem 2 with
$\varepsilon_{i}\alpha_{i}={1\over 3}\left({2\over
3}\right)^{i-1}(U_{0}-L_{0})$ that decreases geometrically. The total number
of iterations needed to solve all $\alpha_{i}$-feasibility test during the
search is then upper bounded by the complexity needed for the very last
feasibility test, which is $O\big{(}\frac{R^{2}\rho^{2}\ln
r}{\varepsilon^{2}{\sf OPT}^{2}}\big{)}$.

Extension to minimization problems. With slight changes, our meta-algorithm
also works for minimization problems. Consider the following variant of (SCP)
where the primal problem is minimizing and the dual problem is maximizing:

$\begin{aligned}
\min\quad&{\bm{c}}\bullet{\bm{x}}+{\bm{c}}^{\prime}\bullet{\bm{x}}^{\prime}\\\
\operatorname*{s.t.}\quad&{\bm{a}}_{j}\bullet{\bm{x}}+{\bm{a}}_{j}^{\prime}\bullet{\bm{x}}^{\prime}=b_{j}\quad\forall
j\in[m]\\\ &{\bm{x}}\in\mathcal{K},\
{\bm{x}}^{\prime}\in\mathcal{K}^{\prime}\end{aligned}\hskip
30.00005pt\begin{aligned} \max\quad&{\bm{b}}^{T}{\bm{y}}\\\
\operatorname*{s.t.}\quad&\textstyle\sum_{j=1}^{m}{\bm{a}}_{j}y_{j}\preceq_{\mathcal{K}}{\bm{c}}\\\
&\textstyle\sum_{j=1}^{m}{\bm{a}}^{\prime}_{j}y_{j}\preceq_{\mathcal{K}^{\prime}}{\bm{c}}^{\prime}\end{aligned}$
(14)

This time, given a specified objective value $\alpha$, the
$\alpha$-feasibility test problem is to determine whether
$(1-\varepsilon)\alpha\leq{\sf OPT}$ or $\alpha>{\sf OPT}$. For the former
case, we provide a dual solution with objective value at least
$(1-\varepsilon)\alpha$, and for the latter case we want a primal solution
with value at most $\alpha$. Define the constraint sets

$\mathcal{F}\ {\triangleq}\
\\{\mathbf{y}\in\mathbb{R}^{m}:\sum_{j=1}^{m}\mathbf{a}_{j}y_{j}\preceq_{\mathcal{K}}\mathbf{c}\\}\
\text{ and }\ \mathcal{C}\ {\triangleq}\
\\{\mathbf{y}\in\mathbb{R}^{m}:\sum_{j=1}^{m}\mathbf{a}^{\prime}_{j}y_{j}\preceq_{\mathcal{K}^{\prime}}\mathbf{c}^{\prime}\\},$
(15)

and the $\alpha$-superlevel set of the dual objective function
$\mathcal{S}_{\alpha}=\\{{\bm{y}}\in\mathbb{R}^{m}:{\bm{b}}^{T}{\bm{y}}\geq\alpha\\}$.
Similar to the maximization problems, we have

$\mathcal{F}\cap\mathcal{C}\cap\mathcal{S}_{\alpha}\neq\varnothing\iff{\sf
OPT}\geq\alpha.$

Thus $\alpha$-FTP can be solved by checking the emptiness of the intersection
of $\mathcal{F},\mathcal{C}$ and $\mathcal{S}_{\alpha}$. Given a
${{\bm{p}}}^{(t)}\in\mathcal{K}$ with $\mathbf{Tr}({{\bm{p}}}^{(t)})=1$, we
parameterize the hyperplane ${\mathcal{H}}^{(t)}$ as

${\mathcal{H}}^{(t)}=\\{{\bm{y}}\in\mathbb{R}^{m}:({\bm{c}}-\textstyle\sum_{j}{\bm{a}}_{j}y_{j})\bullet{{\bm{p}}}^{(t)}=0\\},$

so that the region $\mathcal{F}$ is kept inside the closed positive halfspace
${\mathcal{H}}^{(t)}_{+}$. To test if the hyperplane ${\mathcal{H}}^{(t)}$
separates $\mathcal{F}$ from $\mathcal{C}\cap\mathcal{S}_{\alpha}$ (or
equivalently whether
${\mathcal{H}}^{(t)}_{+}\cap\mathcal{C}\cap\mathcal{S}_{\alpha}=\varnothing$),
the ORACLE for minimization problems becomes:

Given ${{\bm{p}}}^{(t)}\in\mathcal{K}$ with $\mathbf{Tr}({{\bm{p}}}^{(t)})=1$,
determine if there exists ${\bm{y}}\in\mathbbm{R}^{m}$ such that
$({\bm{c}}-\textstyle\sum_{j}{\bm{a}}_{j}y_{j})\bullet{{\bm{p}}}^{(t)}\geq 0,\
{\bm{b}}^{T}{\bm{y}}\geq\alpha\text{ and }{\bm{y}}\in\mathcal{C}.$

Similarly, the ORACLE can be implemented by solving the following convex
problem:

$\displaystyle\max_{\bm{y}}$
$\displaystyle{\bm{c}}\bullet{{\bm{p}}}^{(t)}-\textstyle\sum_{j}({\bm{a}}_{j}\bullet{{\bm{p}}}^{(t)})y_{j}$
(16) $\displaystyle\operatorname*{s.t.}$
$\displaystyle{\bm{b}}^{T}{\bm{y}}\geq\alpha,\quad{\bm{y}}\in\mathcal{C}.$

If the optimal value of the above problem is non-negative, and
${{\bm{y}}}^{(t)}$ is the solution. Let $\rho$ be the width of the above
ORACLE which is defined as an upperbound on
$\|{\bm{c}}-\sum_{j}{\bm{a}}_{j}y_{j}\|_{\infty}$. We define the loss vector
${{\bm{m}}}^{(t)}={1\over\rho}({\bm{c}}-\sum_{j}{\bm{a}}_{j}{{\bm{y}}}^{(t)})$
and update ${{\bm{p}}}^{(t)}$ using SCMWU:

${{\bm{p}}}^{(t+1)}=\frac{1}{C}\cdot\mathbf{exp}\Big{(}{-\eta
t\over\rho}({\bm{c}}-\sum_{j}{\bm{a}}_{j}{\bar{y}}^{(t)}_{j})\Big{)},$ (17)

where $C$ normalizes ${{\bm{p}}}^{(t+1)}$ to be trace-one. If the value of
problem (16) is negative, we know that $\alpha>{\sf OPT}$ and a scaled version
of ${{\bm{p}}}^{(t)}$ serves as a primal solution of (14) with objective value
at most $\alpha$. The running time of the meta-algorithm for solving
$\alpha$-FTP for the minimization problem is concluded as the following
theorem.

######  Theorem 3.

Consider problem (14) over the symmetric cone
$\mathcal{K}\operatorname*{\operatorname*{\oplus}}\mathcal{K}^{\prime}$, where
the rank of $\mathcal{K}$ is $r$, and assumption (7) holds. Fix an error
parameter $\varepsilon>0$ and a guess $\alpha>0$ for the optimal value OPT.
Execute the primal-dual meta-algorithm for
$T=\left\lceil\frac{4R^{2}\rho^{2}\ln
r}{\varepsilon^{2}\alpha^{2}}\right\rceil$ iterations, adjusting the
hyperplane using SCMWU (as defined in (17)) with step size
$\eta=\sqrt{\frac{\ln r}{T}}$ and an ORACLE with width $\rho$. Then, we either
correctly conclude that ${\sf OPT}<\alpha$, or, if the algorithm does not
terminate within $T$ iterations, we conclude that ${\sf
OPT}\geq(1-\varepsilon)\alpha$ and the vector

$\tilde{{\bm{y}}}=\left({\bar{y}}^{(T)}_{1},\ldots,{\bar{y}}^{(T)}_{m-1},{\bar{y}}^{(T)}_{m}-\frac{\varepsilon\alpha}{R}\right)$

is a dual feasible solution of (14) with objective value
${\bm{b}}^{T}\tilde{{\bm{y}}}\geq(1-\varepsilon)\alpha$.

##  Applications

We apply the meta-algorithm to derive novel algorithms for two geometry
applications: the smallest enclosing sphere problem (Section 4.1) and the
support vector machine problem (a.k.a. polytope distance, Section 4.2). By
meta-algorithm, we refer to the meta-algorithm for $\alpha$-FTP that
incorporates the necessary search of a good $\alpha$ to solve (SCP) or (14).
For each application, we also show the analysis of the algorithm in the
parallel setting.

###  Smallest enclosing sphere

In the Smallest Enclosing Sphere (SES) problem, the goal is to find the sphere
with the minimum radius that encloses all the spheres in the input set
$\mathcal{I}=\\{S({\bm{v}}_{1},\gamma_{1}),\dots,S({\bm{v}}_{n},\gamma_{n})\\}$,
where $S({\bm{v}}_{i},\gamma_{i})$ is a sphere centered at ${\bm{v}}_{i}$ with
radius $\gamma_{i}$. SES is an important problem arising in many branches of
computer science such as computer graphics, machine learning, and operations
research. In the past few decades, numerous methods have been proposed for
solving the SES problem. These methods can be grouped under two classes: exact
and approximation methods. The exact methods (for example, [32, 33, 18, 48,
22, 24]) mainly focus on solving the SES problem (for a set of points) in some
fixed (low) dimension where the running time is usually linear in the input
size $n$. However, when extending these exact methods to higher dimensions,
the dependence on dimension could be sub-exponential or worse. On the other
hand, the approximation methods are considerably faster in general (high)
dimensions, and can be readily adapted to diverse input types such as spheres
and ellipsoids. A notable set of approximation methods uses core-sets [7, 29,
14, 49, 36], to compute an $(1+\varepsilon)$-approximate solution using only a
subset of the input points of size $O(\frac{1}{\varepsilon})$ or
$O(\frac{1}{\varepsilon^{2}})$, and the running time is usually linear in the
input size, i.e. $O(\frac{nd}{\varepsilon})$ or
$O(\frac{nd}{\varepsilon^{2}})$, where $d$ is the dimensionality of the
problem. Note that the SES problem can also be formulated as a SOCP. However,
employing the interior point methods on SOCPs often requires solving a linear
system in each iteration and results in a super-quadratic or quadratic
dependence in the dimensionality: see, for example, [29, 51].

We now formulate the SES problem as an SCP and use our meta-algorithm to
develop an algorithm with a nearly linear running time. Clearly, SES can be
formulated as the following optimization problem:

$\displaystyle\min_{{\bm{u}},\gamma}$ $\displaystyle\gamma$ (SES)
$\displaystyle\operatorname*{s.t.}$
$\displaystyle\|{\bm{u}}-{\bm{v}}_{i}\|_{2}\leq\gamma-\gamma_{i}\quad\forall
i\in[n]$

Since the constraint $\|{\bm{u}}-{\bm{v}}_{i}\|_{2}\leq\gamma-\gamma_{i}$ can
be interpreted as a generalized inequality constraint w.r.t.
$\mathcal{Q}^{d+1}$, the problem can be viewed as a SCP where the cone is the
Cartesian product of $n$ second-order cones,
$\operatorname*{\operatorname*{\oplus}}_{i=1}^{n}\mathcal{Q}^{d+1}$. Let
$({\bm{u}}^{*},\gamma^{*})$ be the optimal solution to the above program, then
the SES of $\mathcal{I}$ is centered at ${\bm{u}}^{*}$ and the optimal radius
is $\gamma^{*}$.

In the algorithm, we use the constraint
$\|{\bm{u}}-{\bm{v}}_{1}\|_{2}\leq\gamma-\gamma_{1}$ to define the set
$\mathcal{C}$ and the remaining constraints to define the set $\mathcal{F}$.
Notice that the vector ${{\bm{p}}}^{(t)}$ that parameterizes the hyperplane
${\mathcal{H}}^{(t)}$ lies in the product cone
$\operatorname*{\operatorname*{\oplus}}_{i=2}^{n}\mathcal{Q}^{d-1}$. The
running time of our algorithm is stated below.

######  Theorem 4.

Applying the meta-algorithm to problem (SES) for $n$ spheres in a
$d$-dimensional space, we can compute an enclosing sphere with radius at most
$(1+\varepsilon)\gamma^{*}$ in $O\left(\frac{nd\log
n}{\varepsilon^{2}}\right)$ time.

###### Proof.

Noting that

$\|{\bm{u}}-{\bm{v}}_{i}\|_{2}\leq\gamma-\gamma_{i}\iff({\bm{u}};\gamma)-({\bm{v}}_{i};\gamma_{i})\in\mathcal{Q}^{d+1},$

we can write the (SES) problem as an SOCP:

${\min_{{\bm{u}},\gamma}}\left\\{\gamma:({\bm{u}};\gamma)\succeq_{\mathcal{Q}^{d+1}}({\bm{v}}_{i};\gamma_{i}),\
i=1,\ldots,n\right\\}.$

For the (SES) problem, distinguishing between easy and hard constraints is
very natural. The “easy” portion is derived by keeping only one of the second-
order cone constraints (for instance, the first one), while the “hard” part is
defined by the intersection of all other constraints, i.e.

$\displaystyle\mathcal{C}$
$\displaystyle=\\{({\bm{u}};\gamma)\in\mathbbm{R}^{d+1}:\|{\bm{u}}-{\bm{v}}_{1}\|_{2}\leq\gamma-\gamma_{1}\\},$
$\displaystyle\text{and}\quad\mathcal{F}$
$\displaystyle=\\{({\bm{u}};\gamma)\in\mathbbm{R}^{d+1}:\|{\bm{u}}-{\bm{v}}_{i}\|_{2}\leq\gamma-\gamma_{i},\
\text{for}\ i=2,\dots,n\\}.$

Following this splitting, we express the (SES) problem as a pair of
primal/dual SCPs in the form of (SCP). Setting
$\mathcal{K}=\operatorname*{\operatorname*{\oplus}}_{i=2}^{n}\mathcal{Q}^{d+1}$
and $\mathcal{K}^{\prime}=\mathcal{Q}^{d+1}$, the dual formulation is given
by:

$\displaystyle\min$ $\displaystyle{\underbrace{({\bf
0};1)}_{{\bm{b}}}}^{T}\underbrace{({\bm{u}};\gamma)}_{{\bm{y}}}$ (18)
$\displaystyle\operatorname*{s.t.}$
$\displaystyle\Big{(}\sum_{j=1}^{d}u_{j}\underbrace{\operatorname*{\operatorname*{\oplus}}_{i=2}^{n}({\bm{e}}_{j};0)}_{{\bm{a}}_{j}}\Big{)}+\gamma\underbrace{\operatorname*{\operatorname*{\oplus}}_{i=2}^{n}({\bf
0};1)}_{{\bm{a}}_{d+1}}\succeq_{\mathcal{K}}\underbrace{\operatorname*{\operatorname*{\oplus}}_{i=2}^{n}({\bm{v}}_{i};\gamma_{i})}_{{\bm{c}}}$
$\displaystyle\Big{(}\sum_{j=1}^{d}u_{j}\underbrace{({\bm{e}}_{j};0)}_{{\bm{a}}^{\prime}_{j}}\Big{)}+\gamma\underbrace{({\bf
0};1)}_{{\bm{a}}^{\prime}_{d+1}}\succeq_{\mathcal{K}^{\prime}}\underbrace{({\bm{v}}_{1};\gamma_{1})}_{{\bm{c}}^{\prime}}.$

On the other hand, the primal formulation is given by:

$\displaystyle\max$
$\displaystyle\sum_{i=1}^{n}({\bm{v}}_{i};\gamma_{i})\bullet{\bm{x}}_{i}$ (19)
$\displaystyle\operatorname*{s.t.}$
$\displaystyle\underbrace{\Big{(}\operatorname*{\operatorname*{\oplus}}_{i=2}^{n}({\bm{e}}_{j};0)\Big{)}}_{{\bm{a}}_{j}}\bullet{\bm{x}}+\underbrace{({\bm{e}}_{j};0)}_{{\bm{a}}^{\prime}_{j}}\bullet{\bm{x}}^{\prime}=0,\
j=1,\ldots,d$
$\displaystyle\underbrace{\Big{(}\operatorname*{\operatorname*{\oplus}}_{i=2}^{n}({\bf
0};1)\Big{)}}_{{\bm{a}}_{d+1}}\bullet{\bm{x}}+\underbrace{({\bf
0};1)}_{{\bm{a}}^{\prime}_{d+1}}\bullet{\bm{x}}^{\prime}=1$
$\displaystyle{\bm{x}}=({\bm{x}}_{1};\ldots;{\bm{x}}_{n-1})\in\mathcal{K},\
{\bm{x}}^{\prime}={\bm{x}}_{n}\in\mathcal{K}^{\prime}$

Examining (19) and (18) we see that assumption (7) is implicitly satisfied.
Specifically, if we scale the last constraint in (19) by $\sqrt{2}$, we get an
equivalent problem where the constraint becomes

$\underbrace{\Big{(}\operatorname*{\operatorname*{\oplus}}_{i=2}^{n}({\bf
0};\sqrt{2})\Big{)}}_{{\bm{a}}_{d+1}}\bullet{\bm{x}}+\underbrace{({\bf
0};\sqrt{2})}_{{\bm{a}}^{\prime}_{d+1}}\bullet{\bm{x}}^{\prime}=\sqrt{2},$

and we have $m=d+1$,
${\bm{a}}_{d+1}=\operatorname*{\operatorname*{\oplus}}_{i=2}^{n}({\bf
0};\sqrt{2})={\bm{e}}$ (the identity in $\mathcal{K}$) and
${\bm{a}}^{\prime}_{d+1}=({\bf 0};\sqrt{2})\in\mathcal{K}^{\prime}$. Moreover,
we see $b_{d+1}=R=\sqrt{2}$. Note that we don’t need to perform the scaling
explicitly, but the scaled constraint helps in deriving the final result.

Defining $D=\max_{i\geq
2}\|{\bm{v}}_{1}-{\bm{v}}_{i}\|_{2}+\gamma_{1}+\gamma_{i}$, then a sphere
centered at ${\bm{v}}_{1}$ with radius $D$ will enclose $\mathcal{I}$, so
$\gamma^{*}\leq D$. The distance of the farthest pair of points in
$\mathcal{I}$ is at least $D$, so $\gamma^{*}\geq\frac{D}{2}$. To solve the
SES problem, we can search for the optimal radius in the range
$[\frac{D}{2},D]$.

Fix an error parameter $\varepsilon>0$ and a guess $\alpha>0$ for the optimal
value OPT. Let ${{\bm{p}}}^{(t)}$ be the vector that parameterizes
${\mathcal{H}}^{(t)}$ in the $t$-th iteration: it lies in the product cone
$\mathcal{K}=\operatorname*{\operatorname*{\oplus}}_{i=2}^{n}\mathcal{Q}^{d+1}$
that corresponds to $\mathcal{F}$. Let
$({{\bm{w}}}^{(t)}_{i};{s}^{(t)}_{i})\in\mathcal{Q}^{d+1}$, where
${{\bm{w}}}^{(t)}_{i}\in\mathbbm{R}^{d}$ and ${s}^{(t)}_{i}\in\mathbbm{R}$, be
the $(i-1)^{th}$ component of ${{\bm{p}}}^{(t)}$. In other words,
${{\bm{p}}}^{(t)}=\operatorname*{\operatorname*{\oplus}}_{i=2}^{n}({{\bm{w}}}^{(t)}_{i};{s}^{(t)}_{i})$.
The convex program that the ORACLE needs to solve to check whether
${\mathcal{H}}^{(t)}_{+}\cap\mathcal{C}\cap\mathcal{S}_{\alpha}=\varnothing$
(recall (8)) is given by:

$\displaystyle\max_{{\bm{u}},\gamma}$
$\displaystyle\sum_{i=2}^{n}({\bm{u}};\gamma)\bullet({{\bm{w}}}^{(t)}_{i};{s}^{(t)}_{i})-\sum_{i=2}^{n}({\bm{v}}_{i};\gamma_{i})\bullet({{\bm{w}}}^{(t)}_{i};{s}^{(t)}_{i})$
$\displaystyle\operatorname*{s.t.}$
$\displaystyle\gamma\leq\alpha,\quad\|{\bm{u}}-{\bm{v}}_{1}\|_{2}\leq\gamma-\gamma_{1}.$

Since
$\sum_{i=2}^{n}{s}^{(t)}_{i}=\frac{1}{\sqrt{2}}\mathbf{Tr}({{\bm{p}}}^{(t)})=\frac{1}{\sqrt{2}}$
(see Section 2.1), the problem can be simplified to:

$\displaystyle\max_{{\bm{u}},\gamma}$
$\displaystyle\left(\sum_{i=2}^{n}{{\bm{w}}}^{(t)}_{i}\right)^{T}{\bm{u}}+\frac{\gamma}{\sqrt{2}}-\sum_{i=2}^{n}{\bm{v}}_{i}^{T}{{\bm{w}}}^{(t)}_{i}-\sum_{i=2}^{n}\gamma_{i}{s}^{(t)}_{i}$
(20) $\displaystyle\operatorname*{s.t.}$
$\displaystyle\gamma\leq\alpha,\quad\|{\bm{u}}-{\bm{v}}_{1}\|_{2}\leq\gamma-\gamma_{1}$

Notice that once $({{\bm{w}}}^{(t)}_{i};{s}^{(t)}_{i})$ are fixed, the last
two terms in the objective are constants, and increasing the value of $\gamma$
only relaxes the constraint for ${\bm{u}}$, resulting in a larger objective
value. Thus, the optimal value of (20) will be achieved at $\gamma=\alpha$.
Once $\gamma$ is fixed, the problem boils down to optimizing a linear function
over a ball, which admits a simple closed-form solution.

Next, we show that the width $\rho$ of the ORACLE is at most
$\frac{3D}{\sqrt{2}}$. As explained above, at each iteration the ORACLE
returns a point $({\bm{u}};\alpha)$ that satisfies

$\|{\bm{u}}-{\bm{v}}_{1}\|_{2}\leq\alpha-\gamma_{1}.$ (21)

The width then equals the maximum value of

$\Big{|}\lambda_{k}\big{(}\operatorname*{\operatorname*{\oplus}}_{i=2}^{n}({{\bm{u}}}-{\bm{v}}_{i};{\alpha}-\gamma_{i})\big{)}\Big{|},$

over all ${\bm{u}},k$ that satisfy
$\|{\bm{u}}-{\bm{v}}_{1}\|_{2}\leq\alpha-\gamma_{1}$ and $k\in[r]$, which is

$\max_{i\geq
2}\frac{1}{\sqrt{2}}(\alpha-\gamma_{i}+\|{{\bm{u}}}-{\bm{v}}_{i}\|_{2}).$

Finally, we have

$\displaystyle\max_{i\geq
2}\frac{1}{\sqrt{2}}(\alpha-\gamma_{i}+\|{\bm{u}}-{\bm{v}}_{i}\|_{2})\overset{(a)}{\leq}\frac{1}{\sqrt{2}}(\alpha+\|{\bm{u}}-{\bm{v}}_{1}\|_{2}+\max_{i\geq
2}\|{\bm{v}}_{1}-{\bm{v}}_{i}\|_{2})$
$\displaystyle\overset{(b)}{\leq}\frac{1}{\sqrt{2}}(\alpha+\alpha-\gamma_{1}+\max_{i\geq
2}\|{\bm{v}}_{1}-{\bm{v}}_{i}\|_{2})\overset{(c)}{\leq}\frac{1}{\sqrt{2}}(2\alpha+D)\overset{(d)}{\leq}\frac{3D}{\sqrt{2}},$

where $(a)$ is due to the triangle inequality, $(b)$ follows from (21), $(c)$
follows from the definition of $D$, and $(d)$ is because $\alpha\leq D$.

By Theorem 2 we know that if we execute the primal-dual meta-algorithm for
$T=\left\lceil\frac{4R^{2}\rho^{2}\ln
r}{\varepsilon^{2}\alpha^{2}}\right\rceil$ iterations, we either correctly
conclude infeasibility or we find a solution with value
$(1+\varepsilon)\alpha$. For the latter case, a sphere centered at
$\bar{{\bm{u}}}=\frac{1}{T}\sum_{t}{{\bm{u}}}^{(t)}$ with radius
$(1+\varepsilon)\alpha$ will enclose all spheres in $\mathcal{I}$. As
discussed above, for problem (SES), we have
$R=\sqrt{2},\rho\leq\frac{3D}{\sqrt{2}}$ and $r=2(n-1)$. The number of
iterations $T$ can be simplified to
$\left\lceil\frac{36D^{2}\ln(2n-2)}{\varepsilon^{2}\alpha^{2}}\right\rceil$,
which is $O\left(\frac{D^{2}\log n}{\varepsilon^{2}\alpha^{2}}\right)$.

We now analyze the overall running time of the algorithm. For the SES problem,
since we obtained a natural lower and upper bound on the value of $\alpha$,
where $L_{0}=\frac{D}{2}$ and $U_{0}=D$. Therefore, the technique of reducing
optimization problems to feasibility problems discussed in Section 3 applies.
The total number of iterations required for solving the optimization problem
is dominated by the number of iterations required for the last feasibility
problem, which is $O(\frac{D^{2}\log n}{\varepsilon^{2}{\sf OPT}^{2}})$ and
can be further simplified to $O(\frac{\log n}{\varepsilon^{2}})$ because ${\sf
OPT}\geq\frac{D}{2}$. The computation in each iteration of the algorithm
includes: $(a)$ Solving (20) for test of separation. This involves computing
the coefficients of the objective function and testing the intersection of a
halfspace and a hyperball, which requires $O(nd)$ time. $(b)$ Performing a
SCWMU step for adjusting the hyperplane. This involves computing the loss
vector and its spectral decomposition and performing the exponentiation and
normalization step, which requires $O(nd)$ time in total. Combining the
iteration complexity and the running time of each iteration, we conclude that
the overall running time of our algorithm is $O(\frac{nd\log
n}{\varepsilon^{2}})$. ∎

######  Theorem 5.

The algorithm described in Theorem 4 can be parallelized under the setting of
the work-span model. The parallel algorithm requires $O\Big{(}\frac{nd\log
n}{\varepsilon^{2}}\Big{)}$ work and has $O\Big{(}\frac{\log(nd)\log
n}{\varepsilon^{2}}\Big{)}$ span.

###### Proof.

To parallelize the algorithm for the SES problem, we parallelize the
computation during each iteration. The computation in each iteration of the
algorithm involves computing the analytical solution of (20) and updating an
element of the product cone
$\operatorname*{\operatorname*{\oplus}}_{i=2}^{n}\mathcal{Q}^{d+1}$ via SCMWU,
which only requires simple arithmetic and information aggregation operations
(summations or extremum queries) on vectors, and can be trivially
parallelized. The work of each iteration of the parallel algorithm is the same
as the time complexity of the sequential version, while the span (or parallel
depth) becomes $O\left(\log(nd)\right)$. ∎

###  Support vector machine

The Support Vector Machine (SVM) is one of the most successful machine
learning tools. Its efficacy is rooted in its inherent simplicity, coupled
with exceptional generalization properties applicable to both classification
and regression tasks. The vanilla SVM problem [15] is to find a hyperplane
that separates two sets of data points with maximum margin, which is closely
related to the geometric problem of finding the closest points in two convex
polytopes (a.k.a. polytope distance) [9]. The soft-margin variants of SVM,
such as $C$-SVM [15] and $\nu$-SVM [40], can also be interpreted as polytope
distance (PD) problems for reduced polytopes [9, 10, 17]. The algorithms
designed for SVM are usually based on optimization techniques. A milestone in
this direction is the proposal of the sequential minimal optimization (SMO)
method [38], which can be viewed as a block coordinate descent method and is
widely used in practice. Another line of research focuses on solving the PD
problem using geometric approaches: remarkable algorithms in this modality
include the Gilbert-Schlesinger-Kozinec (GSK) method [21] and the Mitchell-
Dem’yanov-Malozemov (MDM) method [34]. In [23], the authors also proved that
the size of core-set for the SVM/PD problem is
$\Theta(\frac{E}{\varepsilon})$, where $E$ is the excentricity of the problem
instance, which is a quantity that measures its difficulty.

In this section we study the SVM problem in the separable case. We formulate
the problem as a pair of primal/dual SCPs and apply the meta-algorithm to
develop an algorithm that produces solutions for both problems in nearly
linear time. Given two point sets
$\mathcal{P}=\\{{\bm{u}}_{1},\dots,{\bm{u}}_{n_{1}}\\}$ and
$\mathcal{Q}=\\{{\bm{v}}_{1},\dots,{\bm{v}}_{n_{2}}\\}$ in $\mathbbm{R}^{d}$,
we define the matrices ${\bm{P}}=({\bm{u}}_{1},\dots,{\bm{u}}_{n_{1}})$ and
${\bm{Q}}=({\bm{v}}_{1},\dots,{\bm{v}}_{n_{2}})$, where ${\bm{u}}_{i}$ is the
$i$-th column of ${\bm{P}}$ and, similarly, ${\bm{v}}_{i}$ the $i$-th column
of ${\bm{Q}}$. We can write SVM as the optimization problem (e.g. see [41]):

$\displaystyle\max_{{\bm{w}},s_{1},s_{2}}$ $\displaystyle s_{1}+s_{2}$ (SVM)
$\displaystyle\operatorname*{s.t.}$
$\displaystyle{\bm{P}}^{T}{\bm{w}}-s_{1}{\bm{1}}\geq{\bm{0}}$
$\displaystyle-{\bm{Q}}^{T}{\bm{w}}-s_{2}{\bm{1}}\geq{\bm{0}}$
$\displaystyle\|{\bm{w}}\|_{2}\leq 1$

If $({\bm{w}}^{*};s_{1}^{*};s_{2}^{*})$ is the optimal solution to the above
problem, then the maximum margin of SVM is given by $s_{1}^{*}+s_{2}^{*}$
while ${\bm{w}}^{*}$ is the normal vector of the separating hyperplane. The
counterpart of (SVM) is the PD problem, whose objective is to compute the
minimum distance between the convex polytopes corresponding to $\mathcal{P}$
and $\mathcal{Q}$:

$\displaystyle\min$
$\displaystyle\|{\bm{P}}{\bm{\mu}}-{\bm{Q}}{\bm{\gamma}}\|_{2}$ (PD)
$\displaystyle\operatorname*{s.t.}$
$\displaystyle{\bm{1}}^{T}{\bm{\mu}}={\bm{1}}^{T}{\bm{\gamma}}=1$
$\displaystyle{\bm{\mu}},{\bm{\gamma}}\geq{\bm{0}}$

In the following, we view (SVM) as a dual SCP in the form of (14). The
following lemma is useful for us to add additional easy constraints to the
problem that help to design the ORACLE without altering the solution to the
original problem.

######  Lemma 6.

Let $D$ be the maximum norm of the input points in $\mathcal{P}$ and
$\mathcal{Q}$. If $({\bm{w}};s_{1};s_{2})$ is a feasible solution of (SVM),
then $s_{1},s_{2}\leq D$.

###### Proof.

Let $({\bm{w}};s_{1};s_{2})$ be any feasible solution to (SVM). From the
constraint ${\bm{P}}^{T}{\bm{w}}-s_{1}{\bm{1}}\geq 0$, we see that for any
${\bm{u}}_{i}\in\mathcal{P}$:

$s_{1}\leq{\bm{u}}_{i}^{T}{\bm{w}}\leq\|{\bm{u}}_{i}\|_{2}.$

Since $\|{\bm{u}}_{i}\|_{2}\leq D\ \forall{\bm{u}}_{i}\in\mathcal{P}$, this
gives:

$s_{1}\leq\max_{{\bm{u}}_{i}\in\mathcal{P}}\|{\bm{u}}_{i}\|_{2}\leq D.$

In the same way, one can prove that $s_{2}\leq D$. ∎

######  Theorem 7.

Let $E=\frac{D^{2}}{{\sf OPT}^{2}}$ be an instance-specified parameter.
Applying the primal-dual meta-algorithm to problem (SVM), we can compute a
$(1-\varepsilon)$-approximate solution to (SVM) and a
$(1+\varepsilon)$-approximate solution to (PD) in $O\Big{(}\frac{E\,nd\log
n}{\varepsilon^{2}}\Big{)}$ time, where $n=n_{1}+n_{2}$.

###### Proof.

By Lemma 6, problem (SVM) is equivalent to:

$\displaystyle\max_{{\bm{w}},s_{1},s_{2}}$ $\displaystyle s_{1}+s_{2}$ (22)
$\displaystyle\operatorname*{s.t.}$
$\displaystyle{\bm{P}}^{T}{\bm{w}}-s_{1}{\bm{1}}\geq{\bm{0}}$
$\displaystyle-{\bm{Q}}^{T}{\bm{w}}-s_{2}{\bm{1}}\geq{\bm{0}}$
$\displaystyle\|{\bm{w}}\|_{2}\leq 1,\ s_{1}\leq D,\ s_{2}\leq D.$

It is natural to split the constraints into “easy” and “hard” parts as
follows:

$\displaystyle\mathcal{C}$
$\displaystyle=\\{({\bm{w}};s_{1};s_{2})\in\mathbbm{R}^{d+2}:\|{\bm{w}}\|_{2}\leq
1,s_{1}\leq D,s_{2}\leq D\\},$ (23) $\displaystyle\text{and}\quad\mathcal{F}$
$\displaystyle=\\{({\bm{w}};s_{1};s_{2})\in\mathbbm{R}^{d+2}:{\bm{P}}^{T}{\bm{w}}-s_{1}{\bm{1}}\geq{\bm{0}},\
-{\bm{Q}}^{T}{\bm{w}}-s_{2}{\bm{1}}\geq{\bm{0}}\\}.$

Following this splitting, we express problem (22) as a dual SCP in the form of
(14). Let ${\bm{P}}_{j}$ be the $j$-th column of ${\bm{P}}$ and ${\bm{Q}}_{j}$
the $j$-th column of ${\bm{Q}}$. Setting $\mathcal{K}=\mathbb{R}^{n}_{+}$ and
$\mathcal{K}^{\prime}=\mathcal{Q}^{d+1}\operatorname*{\operatorname*{\oplus}}\mathbb{R}^{2}_{+}$,
the formulation is given by:

$\displaystyle\max$
$\displaystyle{\underbrace{({\bm{0}};1;1)}_{{\bm{b}}}}^{T}\underbrace{({\bm{w}};s_{1};s_{2})}_{{\bm{y}}}$
(24) $\displaystyle\operatorname*{s.t.}$
$\displaystyle\Big{(}\sum_{j=1}^{d}w_{j}\underbrace{(-{\bm{P}}^{T}_{j};{\bm{Q}}^{T}_{j})}_{{\bm{a}}_{j}}\Big{)}+s_{1}\underbrace{({\bm{1}}_{n_{1}};{\bm{0}})}_{{\bm{a}}_{d+1}}+s_{2}\underbrace{({\bm{0}};{\bm{1}}_{n_{2}})}_{{\bm{a}}_{d+2}}\preceq_{\mathcal{K}}\underbrace{({\bm{0}};{\bm{0}})}_{{\bm{c}}}$
$\displaystyle\Big{(}\sum_{j=1}^{d}w_{j}\underbrace{({\bm{e}}_{j};0)\operatorname*{\operatorname*{\oplus}}(0;0)}_{{\bm{a}}_{j}^{\prime}}\Big{)}+s_{1}\underbrace{({\bm{0}};0)\operatorname*{\operatorname*{\oplus}}(1;0)}_{{\bm{a}}_{d+1}^{\prime}}+s_{2}\underbrace{({\bm{0}};0)\operatorname*{\operatorname*{\oplus}}(0;1)}_{{\bm{a}}_{d+2}^{\prime}}\preceq_{\mathcal{K}^{\prime}}\underbrace{({\bm{0}};1)\operatorname*{\operatorname*{\oplus}}(D;D)}_{{\bm{c}}^{\prime}}$

On the other hand, the primal formulation is given by:

$\displaystyle\min$
$\displaystyle\big{(}({\bm{0}};1)\operatorname*{\operatorname*{\oplus}}(D;D)\big{)}\bullet{\bm{x}}^{\prime}$
(25) $\displaystyle\operatorname*{s.t.}$
$\displaystyle\underbrace{(-{\bm{P}}_{j}^{T};{\bm{Q}}_{j}^{T})}_{{\bm{a}}_{j}}\bullet{\bm{x}}+\underbrace{\big{(}({\bm{e}}_{j};0)\operatorname*{\operatorname*{\oplus}}(0;0)\big{)}}_{{\bm{a}}_{j}^{\prime}}\bullet{\bm{x}}^{\prime}=0,\
j=1,\ldots,d$
$\displaystyle\underbrace{({\bm{1}}_{n_{1}};{\bm{0}})}_{{\bm{a}}_{d+1}}\bullet{\bm{x}}+\underbrace{\big{(}({\bm{0}};0)\operatorname*{\operatorname*{\oplus}}(1;0)\big{)}}_{{\bm{a}}_{d+1}^{\prime}}\bullet{\bm{x}}^{\prime}=1$
$\displaystyle\underbrace{({\bm{0}};{\bm{1}}_{n_{2}})}_{{\bm{a}}_{d+2}}\bullet{\bm{x}}+\underbrace{\big{(}({\bm{0}};0)\operatorname*{\operatorname*{\oplus}}(0;1)\big{)}}_{{\bm{a}}_{d+2}^{\prime}}\bullet{\bm{x}}^{\prime}=1$
$\displaystyle{\bm{x}}\in\mathcal{K},\
{\bm{x}}^{\prime}\in\mathcal{K}^{\prime}$

Notice that by combining the last two constraints, we obtain an implicit
constraint

${\bm{1}}_{n}\bullet{\bm{x}}+\big{(}({\bm{0}};0)\operatorname*{\operatorname*{\oplus}}(1;1)\big{)}\bullet{\bm{x}}^{\prime}=2.$

Here ${\bm{1}}_{n}={\bm{e}}$ (the identity in $\mathcal{K}$) and
$({\bm{0}};0)\operatorname*{\operatorname*{\oplus}}(1;1)\in\mathcal{K}^{\prime}$.
Therefore, assumption (7) is implicitly satisfied by (24) and (25). Moreover,
we have $R=2$.

As we assume that the point sets are separable, 0 is a lower bound for the
optimal margin. From Lemma (6), we also obtain an upperbound
$s_{1}^{*}+s_{2}^{*}\leq 2D$. To solve the SVM problem, we search for the
optimal margin in the range $(0,2D]$.

Fix an error tolerance $\varepsilon>0$ and a guess $\alpha\in(0,2D]$ for the
optimal value ${\sf OPT}$. Let ${{\bm{p}}}^{(t)}$ be the vector that
parameterizes the hyperplane ${\mathcal{H}}^{(t)}$: it lies in the cone
$\mathcal{K}=\mathbbm{R}^{n}_{+}$. Let ${{\bm{\mu}}}^{(t)}$ and
${{\bm{\gamma}}}^{(t)}$ be the components of ${{\bm{p}}}^{(t)}$, where
${{\bm{\mu}}}^{(t)}\in\mathbbm{R}^{n_{1}}_{+}$ and
${{\bm{\gamma}}}^{(t)}\in\mathbbm{R}^{n_{2}}_{+}$. In other words,
${{\bm{p}}}^{(t)}=({{\bm{\mu}}}^{(t)};{{\bm{\gamma}}}^{(t)})$. The convex
program that the ORACLE needs to solve to check whether
${\mathcal{H}}^{(t)}_{+}\cap\mathcal{C}\cap\mathcal{S}_{\alpha}=\varnothing$
(recall (16)) is given by:

$\displaystyle\max_{{\bm{w}},s_{1},s_{2}}$
$\displaystyle({\bm{P}}^{T}{\bm{w}}-s_{1}{\bm{1}})\bullet{{\bm{\mu}}}^{(t)}+({-}{\bm{Q}}^{T}{\bm{w}}-s_{2}{\bm{1}})\bullet{{\bm{\gamma}}}^{(t)}$
$\displaystyle\operatorname*{s.t.}$ $\displaystyle
s_{1}+s_{2}\geq\alpha,\quad\|{\bm{w}}\|_{2}\leq 1,\quad s_{1},s_{2}\leq D,$

which can be re-organized as:

$\displaystyle\max_{{\bm{w}},s_{1},s_{2}}$
$\displaystyle({\bm{P}}{{\bm{\mu}}}^{(t)}-{\bm{Q}}{{\bm{\gamma}}}^{(t)})^{T}{\bm{w}}-\mathbf{Tr}({{\bm{\mu}}}^{(t)})s_{1}-\mathbf{Tr}({{\bm{\gamma}}}^{(t)})s_{2}$
(26) $\displaystyle\operatorname*{s.t.}$ $\displaystyle
s_{1}+s_{2}\geq\alpha,\quad\|{\bm{w}}\|_{2}\leq 1,\quad s_{1},s_{2}\leq D.$

This problem has the following analytical solution:

${{\bm{w}}}^{(t)}=\frac{{\bm{P}}{{\bm{\mu}}}^{(t)}-{\bm{Q}}{{\bm{\gamma}}}^{(t)}}{\|{\bm{P}}{{\bm{\mu}}}^{(t)}-{\bm{Q}}{{\bm{\gamma}}}^{(t)}\|_{2}},\quad{s}^{(t)}_{1}=\begin{cases}D&\text{if
$\mathbf{Tr}({{\bm{\mu}}}^{(t)})\leq\mathbf{Tr}({{\bm{\gamma}}}^{(t)})$},\\\
\alpha-D&\text{otherwise},\end{cases}\quad{s}^{(t)}_{2}=\alpha-{s}^{(t)}_{1}.$

######  Claim 8.

Let ${{\bm{p}}}^{(t)}$ be the vector that parameterizes the separating
hyperplane ${\mathcal{H}}^{(t)}$ in the $t$-th iteration, and let
${{\bm{\mu}}}^{(t)},{{\bm{\gamma}}}^{(t)}$ be its components. If the ORACLE
fails to return a point in
${\mathcal{H}}^{(t)}_{+}\cap\mathcal{C}\cap\mathcal{S}_{\alpha}$, then
$(\frac{{{\bm{\mu}}}^{(t)}}{\mathbf{Tr}({{\bm{\mu}}}^{(t)})};\frac{{{\bm{\gamma}}}^{(t)}}{\mathbf{Tr}({{\bm{\gamma}}}^{(t)})})$
is a solution of (PD) with objective value at most $\alpha$.

We give two distinct proofs of the above claim:

  * •

The ORACLE fails when the optimal value of (26) is negative, i.e.

$({\bm{P}}{{\bm{\mu}}}^{(t)})^{T}{\bm{w}}-\mathbf{Tr}({{\bm{\mu}}}^{(t)}){s}^{(t)}_{1}<({\bm{Q}}{{\bm{\gamma}}}^{(t)})^{T}{\bm{w}}+\mathbf{Tr}({{\bm{\gamma}}}^{(t)}){s}^{(t)}_{2}\quad\text{for
all }\|{\bm{w}}\|_{2}\leq 1.$ (27)

If $\mathbf{Tr}({{\bm{\mu}}}^{(t)})\leq\mathbf{Tr}({{\bm{\gamma}}}^{(t)})$, we
have ${s}^{(t)}_{1}=D$ and ${s}^{(t)}_{2}=\alpha-D$. For the left-hand side of
(27), we have:

$\max_{{\bm{w}}:\|{\bm{w}}\|_{2}\leq
1}({\bm{P}}{{\bm{\mu}}}^{(t)})^{T}{\bm{w}}-\mathbf{Tr}({{\bm{\mu}}}^{(t)}){s}^{(t)}_{1}\leq\max_{{\bm{u}}\in\mathcal{P}}\mathbf{Tr}({{\bm{\mu}}}^{(t)})\|{\bm{u}}\|_{2}-\mathbf{Tr}({{\bm{\mu}}}^{(t)})D\leq
0.$

Therefore, if we scale the left-hand and right-hand side of (27) by
$\frac{1}{\mathbf{Tr}({{\bm{\mu}}}^{(t)})}$ and
$\frac{1}{\mathbf{Tr}({{\bm{\gamma}}}^{(t)})}$ respectively the inequality
still holds, and we have:

$(\textstyle\frac{{\bm{P}}{{\bm{\mu}}}^{(t)}}{\mathbf{Tr}({{\bm{\mu}}}^{(t)})})^{T}{\bm{w}}-{s}^{(t)}_{1}<(\textstyle\frac{{\bm{Q}}{{\bm{\gamma}}}^{(t)}}{\mathbf{Tr}({{\bm{\gamma}}}^{(t)})})^{T}{\bm{w}}+{s}^{(t)}_{2}\quad\text{for
all }\|{\bm{w}}\|_{2}\leq 1.$

Re-arranging the inequality gives:

$\left\|\frac{{\bm{P}}{{\bm{\mu}}}^{(t)}}{\mathbf{Tr}({{\bm{\mu}}}^{(t)})}-\frac{{\bm{Q}}{{\bm{\gamma}}}^{(t)}}{\mathbf{Tr}({{\bm{\gamma}}}^{(t)})}\right\|_{2}<{s}^{(t)}_{1}+{s}^{(t)}_{2}=\alpha.$

Since
${\mathbf{Tr}(}\frac{{{\bm{\mu}}}^{(t)}}{\mathbf{Tr}({{\bm{\mu}}}^{(t)})}{)}={\mathbf{Tr}(}\frac{{{\bm{\gamma}}}^{(t)}}{\mathbf{Tr}({{\bm{\gamma}}}^{(t)})}{)}=1$,
we conclude that
$(\frac{{{\bm{\mu}}}^{(t)}}{\mathbf{Tr}({{\bm{\mu}}}^{(t)})};\frac{{{\bm{\gamma}}}^{(t)}}{\mathbf{Tr}({{\bm{\gamma}}}^{(t)})})$
is a feasible solution to (PD) with objective value at most $\alpha$.

One can get the same result using the same technique for the
$\mathbf{Tr}({{\bm{\mu}}}^{(t)})>\mathbf{Tr}({{\bm{\gamma}}}^{(t)})$ case.

  * •

Since the optimal value of (26) is negative, it is equivalent to saying that
the optimal value of the following problem is less than $\alpha$:

$\displaystyle\max_{{\bm{w}},s_{1},s_{2}}$ $\displaystyle s_{1}+s_{2}$ (28)
$\displaystyle\operatorname*{s.t.}$
$\displaystyle({\bm{P}}{{\bm{\mu}}}^{(t)}-{\bm{Q}}{{\bm{\gamma}}}^{(t)})^{T}{\bm{w}}-\mathbf{Tr}({{\bm{\mu}}}^{(t)})s_{1}-\mathbf{Tr}({{\bm{\gamma}}}^{(t)})s_{2}\geq
0$ $\displaystyle\|{\bm{w}}\|_{2}\leq 1,\quad s_{1},s_{2}\leq D.$

Observe that any point $({\bm{w}};s_{1};s_{2})$ that satisfies the constraints

$\displaystyle({\bm{P}}{{\bm{\mu}}}^{(t)})^{T}{\bm{w}}-\mathbf{Tr}({{\bm{\mu}}}^{(t)})s_{1}\geq
0,$ (29)
$\displaystyle-({\bm{Q}}{{\bm{\gamma}}}^{(t)})^{T}{\bm{w}}-\mathbf{Tr}({{\bm{\gamma}}}^{(t)})s_{2}\geq
0,$ $\displaystyle\|{\bm{w}}\|_{2}\leq 1,$

also satisfies the constraints in (28). Thus the constraints in (29) are more
restrictive than (28), and the optimal value of the following problem is also
less than $\alpha$:

$\displaystyle\max_{{\bm{w}},s_{1},s_{2}}$ $\displaystyle s_{1}+s_{2}$ (30)
$\displaystyle\operatorname*{s.t.}$
$\displaystyle({\bm{P}}{{\bm{\mu}}}^{(t)})^{T}{\bm{w}}-\mathbf{Tr}({{\bm{\mu}}}^{(t)})s_{1}\geq
0,$
$\displaystyle-({\bm{Q}}{{\bm{\gamma}}}^{(t)})^{T}{\bm{w}}-\mathbf{Tr}({{\bm{\gamma}}}^{(t)})s_{2}\geq
0,$ $\displaystyle\|{\bm{w}}\|_{2}\leq 1.$

The dual problem of (30) is:

$\displaystyle\min_{z_{1},z_{2},{\bm{\xi}},\xi_{0}}$ $\displaystyle\xi_{0}$
(31) $\displaystyle\operatorname*{s.t.}$
$\displaystyle{-}({\bm{P}}{{\bm{\mu}}}^{(t)})z_{1}{+}({\bm{Q}}{{\bm{\gamma}}}^{(t)})z_{2}-{\bm{\xi}}={\bm{0}},$
$\displaystyle\mathbf{Tr}({{\bm{\mu}}}^{(t)})z_{1}=\mathbf{Tr}({{\bm{\gamma}}}^{(t)})z_{2}=1,$
$\displaystyle{z_{1},z_{2}\geq 0,\quad}\|{\bm{\xi}}\|_{2}\leq\xi_{0}.$

Let $(\tilde{z}_{1};\tilde{z}_{2};\tilde{{\bm{\xi}}};\tilde{\xi}_{0})$ be the
optimal solution to the above problem. It is not difficult to see that
$\tilde{z}_{1}=\frac{1}{\mathbf{Tr}({{\bm{\mu}}}^{(t)})}$ and
$\tilde{z}_{2}=\frac{1}{\mathbf{Tr}({{\bm{\gamma}}}^{(t)})}$. Since the
optimal value is less than $\alpha$, we have:

$\|\tilde{{\bm{\xi}}}\|_{2}=\left\|{-}({\bm{P}}{{\bm{\mu}}}^{(t)})\tilde{z}_{1}{+}({\bm{Q}}{{\bm{\gamma}}}^{(t)})\tilde{z}_{2}\right\|_{2}\leq\tilde{\xi}_{0}<\alpha.$

As
$\mathbf{Tr}(\tilde{z}_{1}{{\bm{\mu}}}^{(t)})=\mathbf{Tr}(\tilde{z}_{2}{{\bm{\gamma}}}^{(t)})=1$,
the above result implies that
$(\tilde{z}_{1}{{\bm{\mu}}}^{(t)};\tilde{z}_{2}{{\bm{\gamma}}}^{(t)})$ is a
feasible solution to (PD) with objective value at most $\alpha$.

Next, we show that the width $\rho$ of the ORACLE is at most 2D. As explained
above, at each iteration the ORACLE returns a point $({\bm{w}};s_{1};s_{2})$
that satisfies

$\|{\bm{w}}\|_{2}\leq 1,\ s_{1}\leq D,\ s_{2}\leq D,\ s_{1}+s_{2}\geq\alpha.$
(32)

The last constraint, together with the upper bounds of $s_{1}$ and $s_{2}$
implies

$s_{1}\geq\alpha-s_{2}\geq-D\text{ and }s_{2}\geq\alpha-s_{1}\geq-D.$ (33)

The width then equals the maximum value of

$\Big{|}\lambda_{k}\big{(}{\bm{P}}^{T}{\bm{w}}-s_{1}{\bm{1}};-{\bm{Q}}^{T}{\bm{w}}-s_{2}{\bm{1}}\big{)}\Big{|},$

over all $k\in[r]$ and $({\bm{w}};s_{1};s_{2})$ that satisfies (32) and (33),
which is

$\max\Big{\\{}\max_{{\bm{u}}\in\mathcal{P}}\big{|}{\bm{u}}^{T}{\bm{w}}-s_{1}\big{|},\
\max_{{\bm{v}}\in\mathcal{Q}}\big{|}-{\bm{v}}^{T}{\bm{w}}-s_{2}\big{|}\Big{\\}}.$

Finally, we have

$\displaystyle\max_{{\bm{u}}\in\mathcal{P}}\big{|}{\bm{u}}^{T}{\bm{w}}-s_{1}\big{|}\overset{(a)}{\leq}\max_{{\bm{u}}\in\mathcal{P}}\big{|}{\bm{u}}^{T}{\bm{w}}\big{|}+|s_{1}|\overset{(b)}{\leq}D+|s_{1}|\overset{(c)}{\leq}2D,$
and
$\displaystyle\max_{{\bm{v}}\in\mathcal{Q}}\big{|}-{\bm{v}}^{T}{\bm{w}}-s_{2}\big{|}\overset{(a)}{\leq}\max_{{\bm{v}}\in\mathcal{Q}}\big{|}{\bm{v}}^{T}{\bm{w}}\big{|}+|s_{2}|\overset{(b)}{\leq}D+|s_{2}|\overset{(c)}{\leq}2D,$

where $(a)$ are due to the triangle inequality, $(b)$ follow from (32) and the
definition of $D$, and $(c)$ are due to the bounds on $s_{1}$ and $s_{2}$ in
(32) and (33).

By Theorem 3 we know that if we execute the primal-dual meta-algorithm for
$T=\left\lceil\frac{4R^{2}\rho^{2}\ln
r}{\varepsilon^{2}\alpha^{2}}\right\rceil$ iterations, we either correctly
conclude infeasibility or we find a solution with value
$(1-\varepsilon)\alpha$. For the latter case, a hyperplane with normal vector
$\bar{{\bm{w}}}={1\over T}\sum_{t}{{\bm{w}}}^{(t)}$ can separate the point
sets $\mathcal{P}$ and $\mathcal{Q}$ with a margin at least
$(1-\varepsilon)\alpha$. As discussed above, for problem (22), we have
$R=2,\rho\leq 2D$ and $r=n$. The number of iterations $T$ can be simplified to
$\left\lceil\frac{64D^{2}\ln n}{\varepsilon^{2}\alpha^{2}}\right\rceil$.

We now analyze the overall running time of the algorithm for solving the
optimization problem. Since we obtained lower and upper bounds for ${\sf
OPT}$, where $L_{0}=0$ and $U_{0}=2D$. The technique that reduces optimization
problems to $\alpha$-FTP introduced in Section 3 applies. The total number of
iterations required for solving the optimization problem is $O(\frac{D^{2}\log
n}{\varepsilon^{2}{\sf OPT}^{2}})$. Define $E=\frac{D^{2}}{{\sf OPT}^{2}}$ as
a instance-specific parameter. The iteration complexity can be simplified to
$O(\frac{E\log n}{\varepsilon^{2}})$. The computation in each iteration of our
algorithm includes: $(a)$ Solving (26) for test of separation. As in our
previous discussion, the optimal solution to (26) has an analytical solution.
Computing the optimal solution and its objective value takes $O(nd)$ time.
$(b)$ Performing a SCMWU step for adjusting the hyperplane. This involves
computing the loss vector, and performing the update and the normalization
step. Since the vector ${{\bm{p}}}^{(t)}$ lies in $\mathbbm{R}^{n}_{+}$, this
also requires $O(nd)$ time. Combining the iteration complexity and the running
time for each iteration, we conclude that the total time complexity of our
algorithm is $O(\frac{E\,nd\log n}{\varepsilon^{2}})$. ∎

######  Theorem 9.

The algorithm described in Theorem 7 can be parallelized under the work-span
model. The parallel algorithm takes $O\Big{(}\frac{End\log
n}{\varepsilon^{2}}\Big{)}$ work and has $O\Big{(}\frac{E\log(nd)\log
n}{\varepsilon^{2}}\Big{)}$ span.

###### Proof.

Similar to Theorem 5, we parallelize the computation during each iteration of
the algorithm for the SVM/PD problem. This involves computing the analytical
solution of (26) and adjusting the hyperplane using SCMWU for an element in
$\mathbbm{R}^{n}_{+}$. Under the work-span model, the work of the parallel
algorithm is the same as the sequential time complexity, while the span of
each iteration becomes $O(\log(nd))$. When the searching range is small and a
solution of PD is required, we only need to normalize the vectors
${{\bm{\mu}}}^{(t)}$ and ${{\bm{\gamma}}}^{(t)}$, which can be done in $O(n)$
work and $O(\log n)$ span. Therefore, the total work of the parallel algorithm
is $O\Big{(}\frac{End\log n}{\varepsilon^{2}}\Big{)}$ and the span is
$O\Big{(}\frac{E\log(nd)\log n}{\varepsilon^{2}}\Big{)}$. ∎

##  Experimental results

Implementation strategies. We implemented our algorithms for the SES and the
SVM/PD problem in both sequential (CPU) and parallel (GPU) settings. The CPU
versions are developed using C++ and denoted as PDSCP-ST, while the GPU
versions use NVIDIA’s CUDA programming framework [37] and are denoted as
PDSCP. For both CPU and GPU versions, we employ the same implementation
strategies. We examine the solution we obtained in each iteration. For the SES
problem, letting ${\bar{{\bm{u}}}}^{(t)}$ be the mean of the first $t$
iterates, we compute

${f}^{(t)}\triangleq\max_{i\in[n]}\|{\bar{{\bm{u}}}}^{(t)}-{\bm{v}}_{i}\|_{2}+\gamma_{i}.$

For the SVM problem, letting ${\bar{{\bm{w}}}}^{(t)}$ be the first $t$
iterates, we compute

${f}^{(t)}\triangleq\max\Big{\\{}\textstyle\frac{1}{\|{\bar{{\bm{w}}}}^{(t)}\|_{2}}\displaystyle\big{(}\min_{{\bm{v}}\in\mathcal{P}}{\bm{v}}^{T}{\bar{{\bm{w}}}}^{(t)}-\max_{{\bm{u}}\in\mathcal{Q}}{\bm{u}}^{T}{\bar{{\bm{w}}}}^{(t)}\big{)},\
0\Big{\\}}.$

If the current result is already better than our guess $\alpha$ (for SES, when
the current radius ${f}^{(t)}$ is smaller than $\alpha$; and for SVM, when the
current margin ${f}^{(t)}$ is larger than $\alpha$), we can terminate the
process and get an affirmative conclusion for the $\alpha$-feasibility test
problem. Moreover, instead of running the algorithms for a sufficiently large
number of iterations as stated in their theoretical analyses, we use an early-
stopping technique that is often used in practice, and terminate the process
when the result is stabilized. Specifically, in each iteration, we check
whether

$\frac{|{f}^{(t)}-{f}^{(t-1)}|}{{f}^{(t-1)}}<\delta,$

where $\delta$ is a small constant (in our experiments,
$\delta=10^{-4}$).111Note that when ${f}^{(t-1)}=0$, we don’t compute the
ratio and the algorithm continues. If the criterion has been satisfied
consecutively (in our experiments, 10 consecutive iterations), the algorithm
will be terminated.

Experimental setups. We compare the performance of our algorithms against the
two popular commercial softwares Gurobi Optimizer [25] and IBM ILOG Cplex
Optimizer [16], as well as the computational geometry algorithms library CGAL
[43]. The main results of Gurobi and Cplex are executed with the default
settings of using the barrier method on up to 8 threads. For simplicity, the
input instances of the SES problem are randomly generated point sets (i.e.,
$\gamma_{i}=0$), in which case the problem can also be formulated as the
following quadratic program [24]:

$\displaystyle\min_{{\bm{x}},{\bm{y}}}$
$\displaystyle{\bm{y}}^{T}{\bm{y}}-\sum_{i=1}^{n}{\bm{v}}_{i}^{T}{\bm{v}}_{i}x_{i}$
(34) $\displaystyle\operatorname*{s.t.}$
$\displaystyle{\bm{y}}=\sum_{i=1}^{n}{\bm{v}}_{i}x_{i},$
$\displaystyle{\bm{1}}^{T}{\bm{x}}=1,$ $\displaystyle{\bm{x}}\geq{\bm{0}}.$

As for the SVM/PD problem, the input instances are generated to have similar
parameters $E$. The problem also has a quadratic programming formulation:

$\displaystyle\min_{{\bm{\mu}},{\bm{\gamma}}}$
$\displaystyle{\bm{\mu}}^{T}{\bm{P}}^{T}{\bm{P}}{\bm{\mu}}+{\bm{\gamma}}^{T}{\bm{Q}}^{T}{\bm{Q}}{\bm{\gamma}}-2{\bm{\mu}}^{T}{\bm{P}}^{T}{\bm{Q}}{\bm{\gamma}}$
(35) $\displaystyle\operatorname*{s.t.}$
$\displaystyle{\bm{1}}^{T}{\bm{\mu}}=1,\ {\bm{1}}^{T}{\bm{\gamma}}=1,$
$\displaystyle{\bm{\mu}},{\bm{\gamma}}\geq{\bm{0}}.$

We test the performance of Gurobi and Cplex for solving the SES and SVM/PD
problems using both the SOCP formulations ((SES) and (SVM), marked with
“-SOCP” in the figures) and the quadratic program formulations ((34) and (35),
marked with “-QP” in the figures). To better understand the parallel speedup
they get from running on 8 threads, we also include results on them running on
just a single thread. The “-ST” suffix is used to denote the single-thread
variants of the solvers. The algorithms implemented in CGAL are sequential and
simplex-based methods for quadratic programming. On the other hand, we did not
manage to find any publicly available software adapting the core-set methods
for inclusion in our studies. For the SVM problem, almost all the other
existing solvers are designed for the soft-margin variants, which we do not
consider in this work. All the experiments are executed on a machine with an
Intel Core i7-9700K CPU with 32GB memory and 8 cores, and an NVIDIA RTX 2080Ti
GPU with 11GB memory and 4352 CUDA cores. Each resulting running time reported
here is averaged over 10 randomly generated instances. We terminate any
experiment that takes beyond 1000 seconds and do not include them in our
figures.

Different input sizes. The left chart in Figure 3 and Figure 3 shows the
running time of different solvers for solving the SES problem and the SVM
problem respectively, while the input size $n$ (for SVM, $n_{1}$ and $n_{2}$)
ranges from $2^{10}$ to $2^{20}$. In both charts, the running times are
illustrated in a logarithmic scale, and the number of dimensions $d$ is fixed
to $64$. The right chart in Figure 3 and Figure 3 shows the ratio of the
running time of the best competitive solver to that of PDSCP for the SES and
SVM problems. The results demonstrate that our GPU implementation PDSCP
significantly outperforms the multi-threaded SOCP solvers Gurobi-SOCP and
Cplex-SOCP, which use the same problem formulations. With a slower growth rate
to the problem size, PDSCP also outperforms the QP solvers Gurobi-QP and
Cplex-QP (which are better than their corresponding SOCP solvers) as the input
size increases. This is evidenced in the right chart of Figure 3 and Figure 3,
where PDSCP runs faster than the fastest QP solver once the input size is
larger than $2^{17}$ for the SES problem and $2^{14}$ for the SVM problem in
our experiments. The growth rate of the running time (w.r.t. the increase in
input size) of our sequential implementation PDSCP-ST is roughly the same as
the sequential interior point solvers, which is reasonable as in the
theoretical analyses our algorithms have a dependence of $O(n\log n)$ on the
input size and the interior point methods have a dependence of $O(n^{1.5})$.
The results of CGAL do not show up in Figure 3 since they exceeded the time
limit of 1000 seconds.

Figure 2: (Left) Running time of different solvers with inputs of different
sizes, and (Right) the ratio of the running time of the best competitive
solvers to that of PDSCP for the SES problem.

Figure 3:  (Left) Running time of different solvers with inputs of different
sizes, and (Right) the ratio of the running time of the best competitive
solvers to that of PDSCP for the SVM problem.

Different dimensionalities. The left charts in Figure 5 and Figure 5 show the
running time of the solvers for solving the SES problem and the SVM problem
respectively, with the dimensionality $d$ ranging from 2 to 512. In both
figures, the running time is demonstrated on a logarithmic scale. For the SES
problem, the number of input points $n$ is fixed to 100,000, and for the SVM
problem, the number of points in each set ($n_{1}$ and $n_{2}$) is fixed to
100,000. The right charts in Figure 5 and Figure 5 show the ratio of the
running time of the best competitive solver to that of PDSCP for the SES and
SVM problems. The results demonstrate that PDSCP outperforms the multi-
threaded SOCP solvers Gurobi-SOCP and Cplex-SOCP in most cases, which use the
same problem formulations. It also surpasses the QP solvers Gurobi-QP and
Cplex-QP as the dimensionality increases. This is evidenced in the right chart
of Figure 5 and Figure 5, where PDSCP runs faster than the best QP solver once
the dimensionality is larger than 128 for the SES problem and 16 for the SVM
problem. It is worth noting that PDSCP and PDSCP-ST have a significantly
slower growth rate concerning the increase in dimensionality as compared with
the other solvers (either the simplex-based solver of CGAL or interior point
solvers of Gurobi and Cplex), which reflects the results we obtained from the
theoretical analyses, where the running time of our algorithms have a $O(d)$
dependence on the number of dimensions and the dependences of simplex and
interior point methods are super-quadratic. It is reasonable to anticipate
that the advantage of PDCSP will persist in higher dimensions. The results for
Gurobi-SOCP are not available because it runs out of memory during execution.

Figure 4: (Left) Running time of the solvers under different dimensionalities,
and (Right) the ratio of the running time of the best competitive solvers to
that of PDSCP for the SES problem.

Figure 5: (Left) Running time of the solvers under different dimensionalities,
and (Right) the ratio of the running time of the best competitive solvers to
that of PDSCP for the SVM problem.

Parallel speedup. The parallelism of a parallel solver is to be understood
through its speedup, which is defined as the ratio of the running time of its
single-thread version to the running time of its parallel version. From the
left charts in Figure 3–5, it is evident that PDSCP demonstrates a
significantly superior capability, achieving speeds over an order of magnitude
faster than its sequential counterpart, PDSCP-ST. On the other hand, the
charts also depict that all the interior point solvers from Gurobi and Cplex
do not seem to gain much benefit from parallelization. The performances of
their multi-thread versions are sometimes even worse than the sequential
versions due to the extra overhead of multithreading. For example, the
performance of Cplex-QP v.s. Cplex-QP-ST in all the figures, and the
performance of Gurobi-SOCP v.s. Gurobi-SOCP-ST in Figure 3 and Figure 5.

Quality of results. To measure the quality of the results of our algorithms,
we compare them with the best result obtained from other solvers. Let
$\hat{f}$ be the result obtained from PDSCP (for SES, the smallest
${f}^{(t)}$; and for SVM, the largest ${f}^{(t)}$), and let $f^{*}$ be the
best result obtained from other solvers. For the SES problem, we define the
error as $(\hat{f}-f^{*})/f^{*}$; and for the SVM problem, we define the error
as $(f^{*}-\hat{f})/f^{*}$. Notice that the errors can have negative values,
in which cases our results are better than the best results of the other
solvers. The average errors of PDSCP for the SES problem are given in Table 1,
and the average errors for the SVM problem are shown in Table 2. As we can see
from the tables, the error exhibits a slow rate of growth as the input size
increases. The scenario concerning the increase of dimensionality differs. For
the SES problem, the error increases as the dimensionality grows, but for the
SVM problem, the error decreases. In general, PDSCP produces a desirable
result for most problem instances. For the SES problem, the 80th percentile of
the errors over all experimented cases is 0.0056; and for the SVM problem, the
80th percentile of the errors over all inputs is 0.0008.

Table 1: Average errors of PDSCP for the SES problem

Number of points | $2^{10}$ | $2^{11}$ | $2^{12}$ | $2^{13}$ | $2^{14}$ | $2^{15}$ | $2^{16}$ | $2^{17}$ | $2^{18}$ | $2^{19}$ | $2^{20}$
---|---|---|---|---|---|---|---|---|---|---|---
Average error | 0.0019 | 0.0021 | 0.0023 | 0.0024 | 0.0025 | 0.0029 | 0.0031 | 0.0042 | 0.0041 | 0.0044 | 0.0055

Number of dimensions | 2 | 4 | 8 | 16 | 32 | 64 | 128 | 256 | 512
---|---|---|---|---|---|---|---|---|---
Average error | 0.0003 | 0.0009 | 0.0008 | 0.0012 | 0.0024 | 0.0027 | 0.0058 | 0.0085 | 0.0098

Table 2: Average errors of PDSCP for the SVM problem

Number of points | $2^{10}$ | $2^{11}$ | $2^{12}$ | $2^{13}$ | $2^{14}$ | $2^{15}$ | $2^{16}$ | $2^{17}$ | $2^{18}$ | $2^{19}$ | $2^{20}$
---|---|---|---|---|---|---|---|---|---|---|---
Average error | 0.0004 | 0.0004 | 0.0005 | 0.0005 | 0.0006 | 0.0002 | 0.0006 | 0.0006 | 0.0007 | 0.0007 | 0.0007

Number of dimensions | 2 | 4 | 8 | 16 | 32 | 64 | 128 | 256 | 512
---|---|---|---|---|---|---|---|---|---
Average error | 0.0025 | 0.0023 | 0.0015 | 0.0011 | 0.0008 | 0.0006 | 0.0006 | 0.0006 | 0.0007

##  Concluding remark

In this work, we introduce a novel primal-dual framework for symmetric cone
programming (SCP) that utilizes a recent extension of the multiplicative
weights update method to symmetric cones. We applied this framework to devise
parallel algorithms for smallest enclosing sphere and support vector machine
(polytope distance). Our implementation of these two algorithms, tested in
both sequential and parallel settings, demonstrated excellent efficiency and
scalability with large-scale inputs. Looking forward, the versatility of SCP
allows our framework to have broad applications in both theoretical and
practical aspects. Besides the width-dependent method proposed in this paper,
width-independent algorithms for special classes of SCPs (such as positive
SCPs) would also be an interesting topic to study.

## References

  * [1] Zeyuan Allen-Zhu, Yin Tat Lee, and Lorenzo Orecchia.  Using optimization to obtain a width-independent, parallel, simpler, and faster positive SDP solver.  In Proceedings of the 27th annual ACM-SIAM symposium on Discrete algorithms, pages 1824–1831. SIAM, 2016.
  * [2] Zeyuan Allen-Zhu and Lorenzo Orecchia.  Using optimization to break the epsilon barrier: A faster and simpler width-independent algorithm for solving positive linear programs in parallel.  In Proceedings of the 26th annual ACM-SIAM symposium on Discrete algorithms, pages 1439–1456, 2014.
  * [3] Nimar S. Arora, Robert D. Blumofe, and C. Greg Plaxton.  Thread scheduling for multiprogrammed multiprocessors.  In Proceedings of the tenth annual ACM symposium on parallel algorithms and architectures, pages 119–129, 1998.
  * [4] Sanjeev Arora, Elad Hazan, and Satyen Kale.  Fast algorithms for approximate semidefinite programming using the multiplicative weights update method.  In Proceedings of IEEE Symposium on Foundations of Computer Science, 2005.
  * [5] Sanjeev Arora, Elad Hazan, and Satyen Kale.  The multiplicative weights update method: a meta-algorithm and applications.  Theory of Computing, 8:121–164, 2012.
  * [6] Sanjeev Arora and Satyen Kale.  A combinatorial, primal-dual approach to semidefinite programs.  Journal of the ACM, 63(2):1–35, 2016.
  * [7] Mihai Badoiu and Kenneth L. Clarkson.  Smaller core-sets for balls.  In SODA, volume 3, pages 801–802, 2003.
  * [8] Aharon Ben-Tal and Arkadi Nemirovski.  Non-Euclidean restricted memory level method for large-scale convex optimization.  Mathematical Programming, 102(3):407–456, 2005.
  * [9] Kristin P. Bennett and Erin J. Bredensteiner.  Duality and geometry in SVM classifiers.  In Proceedings of the Seventeenth International Conference on Machine Learning, pages 57–64, 2000.
  * [10] Jinbo Bi and Kristin P. Bennett.  A geometric approach to support vector regression.  Neurocomputing, 55(1-2):79–108, 2003.
  * [11] Guy E. Blelloch, Laxman Dhulipala, and Yihan Sun.  Introduction to parallel algorithms.  Lecture notes, Computer Science Department, Carnegie Mellon University, Spring, 2021.
  * [12] Guy E. Blelloch and Margaret Reid-Miller.  Fast set operations using treaps.  In Proceedings of the tenth annual ACM symposium on parallel algorithms and architectures, pages 16–26, 1998.
  * [13] Ilayda Canyakmaz, Wayne Lin, Georgios Piliouras, and Antonios Varvitsiotis.  Multiplicative updates for online convex optimization over symmetric cones.  arXiv:2307.03136, 2023.
  * [14] Kenneth L. Clarkson.  Coresets, sparse greedy approximation, and the Frank-Wolfe algorithm.  ACM Transactions on Algorithms (TALG), 6(4):1–30, 2010.
  * [15] Corinna Cortes and Vladimir Vapnik.  Support-vector networks.  Machine learning, 20:273–297, 1995.
  * [16] Cplex, IBM ILOG.  User’s Manual for CPLEX, 2023.  URL: https://www.ibm.com/support/pages/cplex-optimization-studio-reference-manuals.
  * [17] David Crisp and Christopher J. C. Burges.  A geometric interpretation of $\nu$-SVM classifiers.  In S. Solla, T. Leen, and K. Müller, editors, Advances in Neural Information Processing Systems, volume 12. MIT Press, 1999.
  * [18] Martin Dyer.  A class of convex programs with applications to computational geometry.  In Proceedings of the 8th annual symposium on Computational Geometry, pages 9–15, 1992.
  * [19] Jacques Faraut and Adam Korányi.  Analysis on Symmetric Cones.  Oxford mathematical monographs. Clarendon Press, 1994.
  * [20] Leonid Faybusovich.  Euclidean Jordan Algebras and interior-point algorithms.  Positivity, 1997.
  * [21] Vojtěch Franc and Václav Hlaváč.  An iterative algorithm learning the maximal margin classifier.  Pattern recognition, 36(9):1985–1996, 2003.
  * [22] Bernd Gärtner.  A subexponential algorithm for abstract optimization problems.  SIAM Journal on Computing, 24(5):1018–1035, 1995.
  * [23] Bernd Gärtner and Martin Jaggi.  Coresets for polytope distance.  In Proceedings of the 25th annual symposium on Computational Geometry, pages 33–42, 2009.
  * [24] Bernd Gärtner and Sven Schönherr.  An efficient, exact, and generic quadratic programming solver for geometric optimization.  In Proceedings of the sixteenth annual symposium on Computational geometry, pages 110–118, 2000.
  * [25] Gurobi Optimization, LLC.  Gurobi Optimizer Reference Manual, 2023.  URL: https://www.gurobi.com.
  * [26] Rahul Jain and Penghui Yao.  A parallel approximation algorithm for positive semidefinite programming.  In 2011 IEEE 52nd Annual Symposium on Foundations of Computer Science, pages 463–471. IEEE, 2011.
  * [27] Satyen Kale.  Efficient Algorithms Using The Multiplicative Weights Update Method.  PhD thesis, Princeton University, 2007.
  * [28] Philip Klein and Hsueh-I Lu.  Efficient approximation algorithms for semidefinite programs arising from MAX CUT and COLORING.  In STOC’96, 1996.
  * [29] Piyush Kumar, Joseph S.B. Mitchell, and Emre Alper Yildirim.  Computing core-sets and approximate smallest enclosing hyperspheres in high dimensions.  In Proceedings of the 5th Workshop on Algorithm Engineering and Experiments, pages 45–55, 2003.
  * [30] Michael Luby and Noam Nisan.  A parallel approximation algorithm for positive linear programming.  In Proceedings of the 25th annual ACM symposium on Theory of Computing, pages 448–457, 1993.
  * [31] David G. Luenberger and Yinyu Ye.  Linear and nonlinear programming, volume 2.  Springer, 1984.
  * [32] Nimrod Megiddo.  Linear-time algorithms for linear programming in $\mathbb{R}^{3}$ and related problems.  SIAM journal on computing, 12(4):759–776, 1983.
  * [33] Nimrod Megiddo.  Linear programming in linear time when the dimension is fixed.  Journal of the ACM, 31(1):114–127, 1984.
  * [34] B.F. Mitchell, Vladimir Fedorovich Dem’yanov, and V.N. Malozemov.  Finding the point of a polyhedron closest to the origin.  SIAM Journal on Control, 12(1):19–26, 1974.
  * [35] Arkadi S. Nemirovski and Michael J. Todd.  Self-scaled barriers and interior-point methods for convex programming.  Mathematics of Operations Research, 22(1):1–42, 1997.
  * [36] Frank Nielsen and Richard Nock.  Approximating smallest enclosing balls with applications to machine learning.  International Journal of Computational Geometry & Applications, 19(05):389–414, 2009.
  * [37] NVIDIA.  CUDA C++ Programming Guide, 2023.  URL: https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html.
  * [38] John Platt.  Sequential minimal optimization: A fast algorithm for training support vector machines.  Technical Report MSR-TR-98-14, Microsoft Research, 1998.
  * [39] Serge A. Plotkin, David B. Shmoys, and Éva Tardos.  Fast approximation algorithms for fractional packing and covering problems.  Mathematics of Operations Research, 20(2):257–301, 1995.
  * [40] Bernhard Schölkopf, Alex J. Smola, Robert C. Williamson, and Peter L Bartlett.  New support vector algorithms.  Neural computation, 12(5):1207–1245, 2000.
  * [41] Pannagadatta K. Shivaswamy, Chiranjib Bhattacharyya, and Alexander J. Smola.  Second order cone programming approaches for handling missing and uncertain data.  Journal of Machine Learning Research, pages 1283–1314, 2006.
  * [42] J. Tao, G. Q. Wang, and L. Kong.  The Araki-Lieb-Thirring inequality and the Golden-Thompson inequality in Euclidean Jordan algebras.  Linear and Multilinear Algebra, 70(19):4228–4243, 2022.
  * [43] The CGAL Project.  CGAL User and Reference Manual, 2023.  URL: https://doc.cgal.org/5.6/Manual/packages.html.
  * [44] Koji Tsuda, Gunnar Rätsch, and Manfred K. Warmuth.  Matrix exponentiated gradient updates for on-line learning and Bregman projection.  Journal of Machine Learning Research, 6:995–1018, 2005.
  * [45] Lieven Vandenberghe.  Symmetric cones.  Lecture notes, Electrical and Computer Engineering Department, UCLA, Spring, 2016.
  * [46] MVC Vieira.  Jordan algebraic approach to symmetric optimization.  PhD thesis. Electrical Engineering, Mathematics and Computer Science, Delft University of Technology, 2007.
  * [47] Manfred K Warmuth and Dima Kuzmin.  Online variance minimization.  In International conference on computational learning theory, pages 514–528. Springer, 2006.
  * [48] Emo Welzl.  Smallest enclosing disks (balls and ellipsoids).  In New Results and New Trends in Computer Science: Graz, Austria, June 20–21, 1991 Proceedings, pages 359–370. Springer, 2005.
  * [49] E. Alper Yildirim.  Two algorithms for the minimum enclosing ball problem.  SIAM Journal on Optimization, 19(3):1368–1391, 2008.
  * [50] Neal E. Young.  Sequential and parallel algorithms for mixed packing and covering.  In Proceedings of the 42nd IEEE symposium on Foundations of Computer Science, pages 538–546. IEEE, 2001.
  * [51] Guanglu Zhou, Kim-Chuan Toh, and Jie Sun.  Efficient algorithms for the smallest enclosing ball problem.  Computational Optimization and Applications, 30:147–160, 2005.