spacing=nonfrench

# Performance enhancements for a generic
conic interior point algorithm

Chris Coey  Operations Research Center, MIT, Cambridge, MA  Lea Kapelevich
Operations Research Center, MIT, Cambridge, MA  Juan Pablo Vielma  Google
Research and MIT Sloan School of Management, Cambridge, MA

###### Abstract

In recent work, we provide computational arguments for expanding the class of
proper cones recognized by conic optimization solvers, to permit simpler,
smaller, more natural conic formulations. We define an exotic cone as a proper
cone for which we can implement a small set of tractable (i.e. fast,
numerically stable, analytic) oracles for a logarithmically homogeneous self-
concordant barrier for the cone or for its dual cone. Our extensible, open
source conic interior point solver, Hypatia, allows modeling and solving any
conic optimization problem over a Cartesian product of exotic cones. In this
paper, we introduce Hypatia’s interior point algorithm. Our algorithm is based
on that of Skajaa and Ye (2015), which we generalize by handling exotic cones
without tractable primal oracles. With the goal of improving iteration count
and solve time in practice, we propose a sequence of four enhancements to the
interior point stepping procedure of Skajaa and Ye (2015), which alternates
between prediction and centering steps: (1) we use a less restrictive central
path proximity condition, (2) we adjust the prediction and centering
directions using a new third order directional derivative barrier oracle, (3)
we use a single backtracking search on a quadratic curve instead of two line
searches, and (4) we use a combined prediction and centering step. We
implement 23 useful exotic cones in Hypatia. We summarize the complexity of
computing oracles for these cones, showing that our new third order oracle is
not a bottleneck, and we derive efficient and numerically stable oracle
implementations for several cones. We generate a diverse benchmark set of 379
conic problems from 37 different applied examples. Our computational testing
shows that each stepping enhancement improves Hypatia’s iteration count and
solve time. Altogether, the enhancements reduce the shifted geometric means of
iteration count and solve time by over 80% and 70% respectively.

##### Funding

This work has been partially funded by the National Science Foundation under
grant OAC-1835443 and the Office of Naval Research under grant
N00014-18-1-2079.

###### Contents

  1. 1 Introduction
    1. 1.1 Conic optimization with Hypatia solver
    2. 1.2 The Skajaa-Ye algorithm
    3. 1.3 Practical algorithmic developments
    4. 1.4 Benchmark instances and computational testing
    5. 1.5 Overview
  2. 2 Notation
  3. 3 Exotic cones and oracles
  4. 4 General conic form and certificates
    1. 4.1 General conic form
    2. 4.2 Conic certificates
    3. 4.3 Homogeneous self-dual embedding
  5. 5 Central path following algorithm
    1. 5.1 Central path of the HSDE
    2. 5.2 Central path proximity
    3. 5.3 High level algorithm
    4. 5.4 Search directions
      1. 5.4.1 Centering
      2. 5.4.2 Prediction
    5. 5.5 Stepping procedures
      1. 5.5.1 Basic stepping procedure
      2. 5.5.2 Less restrictive proximity
      3. 5.5.3 Third order adjustments
      4. 5.5.4 Curve search
      5. 5.5.5 Combined directions
  6. 6 Oracles for predefined exotic cones
  7. 7 Computational testing
    1. 7.1 Exotic conic benchmark set
    2. 7.2 Methodology
    3. 7.3 Results
      1. 7.3.1 Less restrictive proximity
      2. 7.3.2 Third order adjustments
      3. 7.3.3 Curve search
      4. 7.3.4 Combined directions
  8. A Preprocessing and solving for search directions
  9. B Efficient proximity checks
  10. C Computing the TOO for some exotic cones
    1. C.1 Intersections of slices of the PSD cone
    2. C.2 LMI cone
    3. C.3 Matrix and scalar WSOS dual cones
    4. C.4 Sparse PSD cone
    5. C.5 Euclidean norm cone and Euclidean norm square cone

##  1 Introduction

Any convex optimization problem may be represented as a conic problem that
minimizes a linear function over the intersection of an affine subspace with a
Cartesian product of primitive proper cones (i.e. irreducible, closed, convex,
pointed, and full-dimensional conic sets). Under certain conditions, a conic
problem has a simple and easily checkable certificate of optimality, primal
infeasibility, or dual infeasibility (Permenter et al., 2017). Most commercial
and open-source conic solvers (such as CSDP (Borchers, 1999), CVXOPT (Andersen
et al., 2011), ECOS (Serrano, 2015), MOSEK (MOSEK ApS, 2020), SDPA (Yamashita
et al., 2003), Alfonso (Papp and Yıldız, 2017)) implement primal-dual interior
point methods (PDIPMs) based on the theory of logarithmically homogeneous
self-concordant barrier (LHSCB) functions. Compared to first order conic
methods (see O’Donoghue et al. (2016) on SCS solver), idealized PDIPMs
typically exhibit higher per-iteration cost, but have a lower iteration
complexity of $\mathcal{O}(\sqrt{\nu}\log(1/\varepsilon))$ iterations to
converge to $\varepsilon$ tolerance, where $\nu$ is the barrier parameter of
the LHSCB. We limit the scope of this paper to conic PDIPMs, but note that
there are other notions of duality and PDIPMs for convex problems outside of
the conic realm (see e.g. Karimi and Tunçel (2018)).

###  1.1 Conic optimization with Hypatia solver

Hypatia (Coey et al., 2020) is an open-source, extensible conic primal-dual
interior point solver.111 Hypatia is available at
github.com/chriscoey/Hypatia.jl under the MIT license. See Coey et al. (2021b)
for documentation, examples, and instructions for using Hypatia.  Hypatia is
written in the Julia language (Bezanson et al., 2017) and is accessible
through a flexible, low-level native interface or the modeling tool JuMP
(Dunning et al., 2017). A key feature of Hypatia is a generic cone interface
that allows users to define new _exotic cones_. An exotic cone is a proper
cone for which we can implement a small set of tractable LHSCB oracles (listed
in Section 1.2) for either the cone or its dual cone. Defining a new cone
through Hypatia’s cone interface makes both the cone and its dual available
for use in conic formulations. We have already predefined 23 useful exotic
cone types (some with multiple variants) in Hypatia. Several cones are new and
required the development of LHSCBs and efficient procedures for oracle
evaluations (see Coey et al. (2021c); Kapelevich et al. (2021)).

Advanced conic solvers such as MOSEK 9 currently recognize at most only a
handful of standard cones: nonnegative, (rotated) second order, positive
semidefinite (PSD), and three-dimensional exponential and power cones. In Coey
et al. (2020), we show across seven applied examples that modeling with the
much larger class of exotic cones often permits simpler, smaller, more natural
conic formulations. Our computational experiments with these examples
demonstrate the potential advantages, especially in terms of solve time and
memory usage, of solving natural formulations with Hypatia compared to solving
standard conic extended formulations with either Hypatia or MOSEK 9. However,
a description of Hypatia’s algorithms was outside the scope of Coey et al.
(2020). In contrast, the main purpose of this paper is to introduce some key
features of our exotic conic PDIPM that enable it to be both general and
performant.

###  1.2 The Skajaa-Ye algorithm

Most conic PDIPM solvers use efficient algorithms specialized for symmetric
cones, in particular, the nonnegative, (rotated) second order, and PSD cones.
Although non-symmetric conic PDIPMs proposed by Nesterov et al. (1996);
Nesterov (2012) can handle a broader class of cones, they have several
disadvantages compared to specialized symmetric methods (e.g. requiring larger
linear systems, strict feasibility of the initial iterate, or conjugate LHSCB
oracles). The algorithm by Skajaa and Ye (2015), henceforth referred to as
_SY_ , addresses these issues by approximately tracing the central path of the
homogeneous self-dual embedding (HSDE) (Andersen et al., 2003; Xu et al.,
1996) to an approximate solution for the HSDE. This final iterate provides an
approximate conic certificate for the conic problem, if a conic certificate
exists. The SY algorithm relies on an idea by Nesterov (2012) that a high
quality prediction direction (enabling a long step and rapid progress towards
a solution) can be computed if the current iterate is in close proximity to
the central path (i.e. it is an approximate scaling point). To restore
centrality after each prediction step, SY performs a series of centering
steps.

By using a different definition of central path proximity to Nesterov (2012),
SY avoids needing conjugate LHSCB oracles.222 Some proposed techniques such as
the Hessian scaling updates and central path proximity definitions of
Myklebust and Tunçel (2014); Dahl and Andersen (2021) require conjugate LHSCB
oracles.  Indeed, a major advantage of SY is that it only requires access to a
few tractable oracles for the primal cone: an initial interior point,
feasibility check, and gradient and Hessian evaluations for the LHSCB. In our
experience, for a large class of proper cones, these oracles can be evaluated
analytically, i.e. without requiring the implementation of iterative numerical
procedures (such as optimization) that can be expensive and may need numerical
tuning. Conjugate LHSCB oracles in general require optimization, and compared
to the analytic oracles, they are often significantly less efficient and more
prone to numerical instability.

###  1.3 Practical algorithmic developments

For many proper cones of interest, including most of Hypatia’s non-symmetric
cones, we are aware of LHSCBs with tractable oracles for either the cone or
its dual cone but not both. Suppose a problem involves a Cartesian product of
exotic cones, some with primal oracles implemented and some with dual oracles
implemented (as in several example formulations described in Coey et al.
(2020)). In this case, SY can solve neither the primal problem nor its conic
dual, as SY requires primal oracles. Our algorithm generalizes SY to allow a
conic formulation over any Cartesian product of exotic cones.

The focus of Skajaa and Ye (2015) is demonstrating that SY has the best known
iteration complexity for conic PDIPMs. This complexity analysis was corrected
by Papp and Yıldız (2017), who implemented SY in their recent MATLAB solver
Alfonso (Papp and Yıldız, 2020, 2021). It is well known that performant PDIPM
implementations tend to violate assumptions used in iteration complexity
analysis, so in this paper we are not concerned with iteration complexity. Our
goal is to reduce iteration counts and solve times in practice, by enhancing
the performance of the interior point stepping procedure proposed by SY and
implemented by Alfonso.

The basic SY-like stepping procedure computes a prediction or centering
direction by solving a large structured linear system, performs a backtracking
line search in the direction, and steps as far as possible given a restrictive
central path proximity condition. We propose a sequence of four practical
performance enhancements.

Less restrictive proximity.

We use a relaxed central path proximity condition, allowing longer prediction
steps and fewer centering steps.

Third order adjustments.

After computing the prediction or centering direction, we compute a _third
order adjustment_ (TOA) direction using a new _third order oracle_ (TOO) for
exotic cones. We use a line search in the unadjusted direction to determine
how to combine it with the TOA direction, before performing a second line
search and stepping in the new adjusted direction.

Curve search.

Due to the central path proximity checks, each backtracking line search can be
quite expensive. Instead of performing two line searches, we use a single
backtracking search along a particular quadratic curve of combinations of the
unadjusted and TOA directions.

Combined directions.

Unlike SY, most conic PDIPMs do not use separate prediction and centering
phases. We compute the prediction and centering directions and their
associated TOA directions, then perform a backtracking search along a
quadratic curve of combinations of all four directions.

Our TOA approach is distinct from the techniques by Mehrotra (1992); Dahl and
Andersen (2021) that also use higher order LHSCB information.333 To avoid
confusion, we do not use the term ‘corrector’ in this paper. In the
terminology of Mehrotra (1992); Dahl and Andersen (2021) our TOA approach is a
type of ‘higher order corrector’ technique, but also our unadjusted centering
direction is referred to by Skajaa and Ye (2015); Papp and Yıldız (2017) as
the ‘corrector’ direction.  Unlike these techniques, we derive adjustments
(using the TOO) for both the prediction and centering directions. Our TOO has
a simpler and more symmetric structure than the third order term used by Dahl
and Andersen (2021), and we leverage this for fast and numerically stable
evaluations. Whereas the method by Mehrotra (1992) only applies to symmetric
cones, and Dahl and Andersen (2021) test their technique only for the standard
exponential cone, we implement and test our TOO for all of Hypatia’s 23
predefined cones. In our experience, requiring a tractable TOO is only as
restrictive as requiring tractable gradient and Hessian oracles. We show that
the time complexity of the TOO is no higher than that of the other required
oracles for each of our cones. To illustrate, we describe efficient and
numerically stable TOO procedures for several cones that can be characterized
as intersections of slices of the PSD cone.

Although this paper is mainly concerned with the stepping procedures, we also
outline our implementations of other key algorithmic components. These include
preprocessing of problem data, finding an initial iterate, the solution of
structured linear systems for search directions, and efficient backtracking
searches with central path proximity checks. We note that Hypatia has a
variety of algorithmic options for these components; these different options
can have a dramatic impact on overall solve time and memory usage, but in most
cases they have minimal effect on the iteration count. For the purposes of
this paper, we only describe and test one set of (default) options for these
components.

###  1.4 Benchmark instances and computational testing

We implement and briefly describe 37 applied examples (available in Hypatia’s
examples folder), each of which has options for creating formulations of
different types and sizes. From these examples, we generate 379 problem
instances of a wide range of sizes. Since there is currently no conic
benchmark storage format that recognizes more than a handful of cone types, we
generate all instances on the fly using JuMP or Hypatia’s native interface.
All of Hypatia’s predefined cones are represented in these instances, so we
believe this is the most diverse conic benchmark set available.

On this benchmark set, we run five different stepping procedures: the basic
SY-like procedure (similar to Alfonso) and the sequence of four cumulative
enhancements to this procedure. Our results show that each enhancement tends
to improve Hypatia’s iteration count and solve time, with minimal impact on
the number of instances solved. We do not enforce time or iteration limits,
but we note that under strict limits the enhancements would greatly improve
the number of instances solved. The TOA enhancement alone leads to a
particularly consistent improvement of around 45% for iteration counts.
Overall, the enhancements together reduce the iterations and solve time by
more than 80% and 70% respectively. For instances that take more iterations or
solve time, the enhancements tend to yield greater relative improvements in
these measures.

###  1.5 Overview

In Section 2, we define our mathematical notation. In Section 3, we define
exotic cones, LHSCBs, and our required cone oracles (including the TOO). In
Section 4, we describe Hypatia’s general primal-dual conic form, associated
conic certificates, and the HSDE. In Section 5, we define the central path of
the HSDE and central path proximity measures, and we outline Hypatia’s high
level algorithm. We also derive the prediction and centering directions and
our new TOA directions, and we describe the SY-like stepping procedure and our
series of four enhancements to this procedure. In Appendices A to B, we
discuss advanced procedures for preprocessing and initial point finding,
solving structured linear systems for directions, and performing efficient
backtracking searches and proximity checks. In Section 6, we briefly introduce
Hypatia’s predefined exotic cones and show that our TOO is relatively cheap to
compute, and in Appendix C we describe some TOO implementations. In Section 7,
we summarize our applied examples and exotic conic benchmark instances, and
finally we present our computational results demonstrating the practical
efficacy of our stepping enhancements.

##  2 Notation

For a natural number $d$, we define the index set $\llbracket
d\rrbracket\coloneqq\\{1,2,\ldots,d\\}$. Often we construct vectors with round
parentheses, e.g. $(a,b,c)$, and matrices with square brackets, e.g.
$\begin{bmatrix}a&b\\\ c&d\end{bmatrix}$. For a set $\mathcal{C}$,
$\operatorname{cl}(\mathcal{C})$ and $\operatorname{int}(\mathcal{C})$ denote
the closure and interior of $\mathcal{C}$, respectively.

$\mathbb{R}$ denotes the space of reals, and $\mathbb{R}_{\geq}$,
$\mathbb{R}_{>}$, $\mathbb{R}_{\leq}$, $\mathbb{R}_{<}$ denote the
nonnegative, positive, nonpositive, and negative reals. $\mathbb{R}^{d}$ is
the space of $d$-dimensional real vectors, and $\mathbb{R}^{d_{1}\times
d_{2}}$ is the $d_{1}$-by-$d_{2}$-dimensional real matrices. The vectorization
operator $\operatorname{vec}:\mathbb{R}^{d_{1}\times
d_{2}}\to\mathbb{R}^{d_{1}d_{2}}$ maps matrices to vectors by stacking
columns, and its inverse operator is
$\operatorname{mat}_{d_{1},d_{2}}:\mathbb{R}^{d_{1}d_{2}}\to\mathbb{R}^{d_{1}\times
d_{2}}$.

$\mathbb{S}^{d}$ is the space of symmetric matrices with side dimension $d$,
and $\mathbb{S}^{d}_{\succeq}$ and $\mathbb{S}^{d}_{\succ}$ denote the
positive semidefinite and positive definite symmetric matrices. The inequality
$S\succeq Z$ is equivalent to $S-Z\in\mathbb{S}^{d}_{\succeq}$ (and similarly
for the strict inequality $\succ$ and $\mathbb{S}^{d}_{\succ}$). We let
$\operatorname{sd}(d)\coloneqq d(d+1)/2$ be the dimension of the vectorized
upper triangle of $\mathbb{S}^{d}$. We overload the vectorization operator
$\operatorname{vec}:\mathbb{S}^{d}\to\mathbb{R}^{\operatorname{sd}(d)}$ to
perform an _svec_ transformation, which rescales off-diagonal elements by
$\sqrt{2}$ and stacks columns of the upper triangle (or equivalently, rows of
the lower triangle). For example, for $S\in\mathbb{S}^{3}$ we have
$\operatorname{sd}(3)=6$ and
$\operatorname{vec}(S)=(S_{1,1},\sqrt{2}S_{1,2},S_{2,2},\sqrt{2}S_{1,3},\sqrt{2}S_{2,3},S_{3,3})\in\mathbb{R}^{\operatorname{sd}(3)}$.
The inverse mapping
$\operatorname{mat}:\mathbb{R}^{\operatorname{sd}(d)}\to\mathbb{S}^{d}$ is
well-defined.

For a vector or matrix $A$, the transpose is $A^{\prime}$ and the trace is
$\operatorname{tr}(A)$. We use the standard inner product on $\mathbb{R}^{d}$,
i.e. $s^{\prime}z={\textstyle\sum_{i\in\llbracket d\rrbracket}}s_{i}z_{i}$ for
$s,z\in\mathbb{R}^{d}$, which equips $\mathbb{R}^{d}$ with the standard norm
$\lVert s\rVert=(s^{\prime}s)^{1/2}$. The linear operators
$\operatorname{vec}$ and $\operatorname{mat}$ preserve inner products, e.g.
$\operatorname{vec}(S)^{\prime}\operatorname{vec}(Z)=\operatorname{tr}(S^{\prime}Z)$
for $S,Z\in\mathbb{R}^{d_{1}\times d_{2}}$ or $S,Z\in\mathbb{S}^{d}$.
$\operatorname{Diag}:\mathbb{R}^{d}\to\mathbb{S}^{d}$ is the diagonal matrix
of a given vector, and $\operatorname{diag}:\mathbb{S}^{d}\to\mathbb{R}^{d}$
is the vector of the diagonal of a given matrix. For dimensions implied by
context, $e$ is a vector of $1$s, $e_{i}$ is the $i$th unit vector, and $0$ is
a vector or matrix of $0$s.

$\lvert x\rvert$ is the absolute value of $x\in\mathbb{R}$ and $\log(x)$ is
the natural logarithm of $x>0$. $\det(X)$ is the determinant of
$X\in\mathbb{S}^{d}$, and $\operatorname{logdet}(X)$ is the log-determinant of
$X\succ 0$. For a vector $x\in\mathbb{R}^{d}$, $\lVert
x\rVert_{\infty}=\max_{i\in\llbracket d\rrbracket}\lvert x_{i}\rvert$ is the
$\ell_{\infty}$ norm and $\lVert x\rVert_{1}=\sum_{i\in\llbracket
d\rrbracket}\lvert x_{i}\rvert$ is the $\ell_{1}$ norm.

Suppose the function $f:\operatorname{int}(\mathcal{C})\rightarrow\mathbb{R}$
is strictly convex and three times continuously differentiable on the interior
of a set $\mathcal{C}\subset\mathbb{R}^{d}$. For a point
$p\in\operatorname{int}(\mathcal{C})$, we denote the gradient and Hessian of
$f$ at $p$ as $\nabla f(p)\in\mathbb{R}^{d}$ and
$\nabla^{2}f(p)\in\mathbb{S}^{d}_{\succ}$. Given an $h\in\mathbb{R}^{d}$, the
first, second, and third order directional derivatives of $f$ at $p$ in
direction $h$ are $\nabla f(p)[h]\in\mathbb{R}$,
$\nabla^{2}f(p)[h,h]\in\mathbb{R}_{\geq}$, and
$\nabla^{3}f(p)[h,h,h]\in\mathbb{R}$.

##  3 Exotic cones and oracles

Let $\mathcal{K}$ be a proper cone in $\mathbb{R}^{q}$, i.e. a conic subset of
$\mathbb{R}^{q}$ that is closed, convex, pointed, and full-dimensional (see
Skajaa and Ye (2015)). Note that requiring $\mathcal{K}$ to be a subset of
$\mathbb{R}^{q}$ simplifies our notation but is not restrictive, e.g. for the
PSD cone, we use the standard svec vectorization (see Section 2). The dual
cone of $\mathcal{K}$ is $\mathcal{K}^{\ast}$, which is also a proper cone in
$\mathbb{R}^{q}$:

$\mathcal{K}^{\ast}\coloneqq\\{z\in\mathbb{R}^{q}:s^{\prime}z\geq 0,\forall
s\in\mathcal{K}\\}.$ (1)

Following Nesterov and Nemirovskii (1994, Sections 2.3.1 and 2.3.3),
$f:\operatorname{int}(\mathcal{K})\to\mathbb{R}$ is a $\nu$-LHSCB for
$\mathcal{K}$, where $\nu\geq 1$ is the _LHSCB parameter_ , if it is three
times continuously differentiable, strictly convex, satisfies
$f(s_{i})\to\infty$ along every sequence
$s_{i}\in\operatorname{int}(\mathcal{K})$ converging to the boundary of
$\mathcal{K}$, and:

$\displaystyle\big{\lvert}\nabla^{3}f(s)[h,h,h]\big{\rvert}$
$\displaystyle\leq 2\bigl{(}\nabla^{2}f(s)[h,h]\bigr{)}^{3/2}$
$\displaystyle\forall s\in\mathcal{K},h\in\mathbb{R}^{q},$ (2a) $\displaystyle
f(\theta s)$ $\displaystyle=f(s)-\nu\log(\theta)$ $\displaystyle\forall
s\in\mathcal{K},\theta\in\mathbb{R}.$ (2b)

Following Renegar (2001, Section 3.3), we define the _conjugate_ of $f$,
$f^{\ast}:\operatorname{int}(\mathcal{K}^{\ast})\to\mathbb{R}$, as:

$f^{\ast}(z)\coloneqq-\textstyle\inf_{s\in\operatorname{int}(\mathcal{K})}\\{s^{\prime}z+f(s)\\},$
(3)

which is a $\nu$-LHSCB for $\mathcal{K}^{\ast}$.

A Cartesian product
$\mathcal{K}=\mathcal{K}_{1}\times\cdots\times\mathcal{K}_{K}$ of $K$ proper
cones is a proper cone, and its dual cone is
$\mathcal{K}^{\ast}=\mathcal{K}_{1}^{\ast}\times\cdots\times\mathcal{K}_{K}^{\ast}$.
In this case, if $f_{k}$ is a $\nu_{k}$-LHSCB for $\mathcal{K}_{k}$, then
$\sum_{k\in\llbracket K\rrbracket}f_{k}$ is an LHSCB for $\mathcal{K}$ with
parameter $\sum_{k\in\llbracket K\rrbracket}\nu_{k}$ (Nesterov and
Nemirovskii, 1994, Proposition 2.3.3). We call $\mathcal{K}$ a primitive cone
if it cannot be written as a Cartesian product of two or more lower-
dimensional cones (i.e. $K$ must equal $1$). Note $\mathcal{K}^{\ast}$ is
primitive if and only if $\mathcal{K}$ is primitive. Primitive proper cones
are the fundamental building blocks of conic formulations.

We call a proper cone $\mathcal{K}$ an exotic cone if we can implement a
particular set of tractable oracles for either $\mathcal{K}$ or
$\mathcal{K}^{\ast}$. Suppose we have tractable oracles for
$\mathcal{K}\subset\mathbb{R}^{q}$ and let
$f:\operatorname{int}(\mathcal{K})\to\mathbb{R}$ denote the $\nu$-LHSCB for
$\mathcal{K}$. The oracles for $\mathcal{K}$ that we require in this paper are
as follows.

Feasibility check.

The strict feasibility oracle checks whether a given point
$s\in\mathbb{R}^{q}$ satisfies $s\in\operatorname{int}(\mathcal{K})$.

Gradient and Hessian evaluations.

Given a point $s\in\operatorname{int}(\mathcal{K})$, the gradient oracle $g$
and Hessian oracle $H$ evaluated at $s$ are:

$\displaystyle g(s)$ $\displaystyle\coloneqq\nabla f(s)\in\mathbb{R}^{q},$
(4a) $\displaystyle H(s)$
$\displaystyle\coloneqq\nabla^{2}f(s)\in\mathbb{S}^{q}_{\succ}.$ (4b)

Third order directional derivative.

Given a point $s\in\operatorname{int}(\mathcal{K})$ and a direction
$\delta_{s}\in\mathbb{R}^{q}$, our new _third order oracle_ (TOO), denoted
$\mathrm{T}$, is a rescaled third order directional derivative vector:

$\mathrm{T}(s,\delta_{s})\coloneqq-\tfrac{1}{2}\nabla^{3}f(s)[\delta_{s},\delta_{s}]\in\mathbb{R}^{q}.$
(5)

Initial interior point.

The initial interior point $t\in\operatorname{int}(\mathcal{K})$ is an
arbitrary point in the interior of $\mathcal{K}$ (which is nonempty since
$\mathcal{K}$ is proper).

In Section 6, we introduce Hypatia’s predefined cones and discuss the time
complexity of computing the feasibility check, gradient, Hessian, and TOO
oracles. In Appendix C, we describe efficient and numerically stable
techniques for computing these oracles for a handful of our cones. Although
Hypatia’s generic cone interface allows specifying additional oracles that can
improve speed and numerical performance (e.g. a dual cone feasibility check,
Hessian product, and inverse Hessian product), these optional oracles are
outside the scope of this paper.

For the initial interior point (which Hypatia only calls once, when finding an
initial iterate), we prefer to use the _central point_ of $\mathcal{K}$. This
is the unique point satisfying
$t\in\operatorname{int}(\mathcal{K})\cap\operatorname{int}(\mathcal{K}^{\ast})$
and $t=-g(t)$ (Dahl and Andersen, 2021). It can also be characterized as the
solution to the following strictly convex problem:

$\textstyle\operatorname*{arg\,min}_{s\in\operatorname{int}(\mathcal{K})}\bigl{(}f(s)+\tfrac{1}{2}\lVert
s\rVert^{2}\bigr{)}.$ (6)

For the nonnegative cone $\mathcal{K}=\mathbb{R}_{\geq}$, $f(s)=-\log(s)$ is
an LHSCB with $\nu=1$, and we have $g(s)=-s^{-1}$ and the central point
$t=1=-g(1)$. For some of Hypatia’s cones, we are not aware of a simple
analytic expression for the central point, in which case we typically use a
non-central interior point.

##  4 General conic form and certificates

In Sections 4.1 and 4.2, we describe our general conic primal-dual form and
the associated conic certificates. In Section 4.3, we introduce the
homogeneous self-dual embedding (HSDE) conic feasibility problem, a solution
of which may provide a conic certificate.

###  4.1 General conic form

Hypatia uses the following primal conic form over variable
$x\in\mathbb{R}^{n}$:

$\displaystyle\textstyle\inf_{x}\quad c^{\prime}x$ $\displaystyle:$ (7a)
$\displaystyle b-Ax$ $\displaystyle=0,$ (7b) $\displaystyle h-Gx$
$\displaystyle\in\mathcal{K},$ (7c)

where $c\in\mathbb{R}^{n}$, $b\in\mathbb{R}^{p}$, and $h\in\mathbb{R}^{q}$ are
vectors, $A:\mathbb{R}^{n}\to\mathbb{R}^{p}$ and
$G:\mathbb{R}^{n}\to\mathbb{R}^{q}$ are linear maps, and
$\mathcal{K}\subset\mathbb{R}^{q}$ is a Cartesian product
$\mathcal{K}=\mathcal{K}_{1}\times\cdots\times\mathcal{K}_{K}$ of exotic
cones. For $k\in\llbracket K\rrbracket$, we let $q_{k}=\dim(\mathcal{K}_{k})$,
so $\sum_{k}q_{k}=q=\dim(\mathcal{K})$. Henceforth we use $n,p,q$ to denote
respectively the variable, equality, and conic constraint dimensions of a
conic problem.

Once a proper cone $\mathcal{K}_{k}$ is defined through Hypatia’s generic cone
interface, both $\mathcal{K}_{k}$ and $\mathcal{K}_{k}^{\ast}$ may be used in
any combination with other cones recognized by Hypatia to construct the
Cartesian product cone $\mathcal{K}$ in 7c. The primal form 7 matches CVXOPT’s
form, however CVXOPT only recognizes symmetric cones (Vandenberghe, 2010).
Unlike the conic form used by Skajaa and Ye (2015); Papp and Yıldız (2021),
which recognizes conic constraints of the form $x\in\mathcal{K}$, our form
does not require introducing slack variables to represent a more general
constraint $h-Gx\in\mathcal{K}$.

The conic dual problem of 7, over variables $y\in\mathbb{R}^{p}$ and
$z\in\mathbb{R}^{q}$ associated with 7b and 7c, is:

$\displaystyle\textstyle\sup_{y,z}\quad-b^{\prime}y-h^{\prime}z$
$\displaystyle:$ (8a) $\displaystyle c+A^{\prime}y+G^{\prime}z$
$\displaystyle=0,$ (8b) $\displaystyle z$
$\displaystyle\in\mathcal{K}^{\ast},$ (8c)

where 8b is associated with the primal variable $x\in\mathbb{R}^{n}$.

###  4.2 Conic certificates

Under certain conditions, there exists a simple conic certificate providing an
easily verifiable proof of infeasibility of the primal 7 or dual 8 problem
(via the conic generalization of Farkas’ lemma) or optimality of a given
primal-dual solution.

A primal improving ray

$x$ is a feasible direction for the primal along which the objective improves:

$\displaystyle c^{\prime}x$ $\displaystyle<0,$ (9a) $\displaystyle-Ax$
$\displaystyle=0,$ (9b) $\displaystyle-Gx$ $\displaystyle\in\mathcal{K},$ (9c)

and hence it certifies dual infeasibility.

A dual improving ray

$(y,z)$ is a feasible direction for the dual along which the objective
improves:

$\displaystyle-b^{\prime}y-h^{\prime}z$ $\displaystyle>0,$ (10a)
$\displaystyle A^{\prime}y+G^{\prime}z$ $\displaystyle=0,$ (10b)
$\displaystyle z$ $\displaystyle\in\mathcal{K}^{\ast},$ (10c)

and hence it certifies primal infeasibility.

A complementary solution

$(x,y,z)$ satisfies the primal-dual feasibility conditions 7b, 7c, 8b and 8c,
and has equal and attained primal and dual objective values:

$c^{\prime}x=-b^{\prime}y-h^{\prime}z,$ (11)

and hence certifies optimality of $(x,y,z)$ via conic weak duality.

One of these certificates exists if neither the primal nor the dual is _ill-
posed_. Intuitively, according to MOSEK ApS (2020, Section 7.2), a conic
problem is ill-posed if a small perturbation of the problem data can change
the feasibility status of the problem or cause arbitrarily large perturbations
to the optimal solution (see Permenter et al. (2017) for more details).

###  4.3 Homogeneous self-dual embedding

The HSDE is a self-dual conic feasibility problem in variables
$x\in\mathbb{R}^{n},y\in\mathbb{R}^{p},z\in\mathbb{R}^{q},\tau\in\mathbb{R},s\in\mathbb{R}^{q},\kappa\in\mathbb{R}$
(see Vandenberghe (2010, Section 6)), derived from a homogenization of the
primal-dual optimality conditions 7b, 7c, 8b, 8c and 11:

$\displaystyle\begin{bmatrix}0\\\ 0\\\ s\\\ \kappa\end{bmatrix}$
$\displaystyle=\begin{bmatrix}0&A^{\prime}&G^{\prime}&c\\\ -A&0&0&b\\\
-G&0&0&h\\\
-c^{\prime}&-b^{\prime}&-h^{\prime}&0\end{bmatrix}\begin{bmatrix}x\\\ y\\\
z\\\ \tau\end{bmatrix},$ (12a) $\displaystyle(z,\tau,s,\kappa)$
$\displaystyle\in\bigl{(}\mathcal{K}^{\ast}\times\mathbb{R}_{\geq}\times\mathcal{K}\times\mathbb{R}_{\geq}\bigr{)}.$
(12b)

For convenience we let
$\omega\coloneqq(x,y,z,\tau,s,\kappa)\in\mathbb{R}^{n+p+2q+2}$ represent a
point. We define the structured $4\times 6$ block matrix
$E\in\mathbb{R}^{(n+p+q+1)\times\dim(\omega)}$ such that 12a is equivalent to:

$E\omega=0.$ (13)

Here we assume $E$ has full row rank; in Appendix A we discuss preprocessing
techniques that handle linearly dependent rows. Note that $\omega=0$ satisfies
12, so the HSDE is always feasible. A point $\omega$ is an interior point if
it is strictly feasible for the conic constraints 12b, i.e. $\omega$ satisfies
$(z,\tau,s,\kappa)\in\operatorname{int}\bigl{(}\mathcal{K}^{\ast}\times\mathbb{R}_{\geq}\times\mathcal{K}\times\mathbb{R}_{\geq}\bigr{)}$.

Suppose a point $\omega$ is feasible for the HSDE 12. From skew symmetry of
the square $4\times 4$ block matrix in 12a, we have
$s^{\prime}z+\kappa\tau=0$. From the conic constraints 12b and the dual cone
inequality 1 we have $s^{\prime}z\geq 0$ and $\kappa\tau\geq 0$. Hence
$s^{\prime}z=\kappa\tau=0$. We consider an exhaustive list of cases below.

Optimality.

If $\tau>0,\kappa=0$, then $(x,y,z)/\tau$ is a complementary solution
satisfying the primal-dual optimality conditions 7b, 7c, 8b, 8c and 11.

Infeasibility.

If $\tau=0,\kappa>0$, then $c^{\prime}x+b^{\prime}y+h^{\prime}z<0$ and we
consider two sub-cases.

Of primal.

If $b^{\prime}y+h^{\prime}z<0$, then $(y,z)$ is a primal infeasibility
certificate satisfying 10.

Of dual.

If $c^{\prime}x<0$, then $x$ is a dual infeasibility certificate satisfying 9.

No information.

If $\tau=\kappa=0$, then $\omega$ provides no information about the
feasibility or optimal values of the primal or dual.

Thus an HSDE solution $\omega$ satisfying $\kappa+\tau>0$ provides an
optimality or infeasibility certificate (see Skajaa and Ye (2015, Lemma 1) and
Vandenberghe (2010, Section 6.1)).

According to Skajaa and Ye (2015, Section 2), if the primal and dual problems
are both feasible and have zero duality gap, SY (their algorithm) finds an
HSDE solution with $\tau>0$ (yielding a complementary solution), and if the
primal or dual (possibly both) is infeasible, SY finds an HSDE solution with
$\kappa>0$ (yielding an infeasibility certificate). This implies that if SY
finds a solution with $\kappa=\tau=0$, then $\kappa=\tau=0$ for all solutions
to the HSDE; in this case, no complementary solution or improving ray exists,
and the primal or dual (possibly both) is ill-posed (Permenter et al., 2017).
The algorithm we describe in Section 5 is an extension of SY that inherits
these properties.

##  5 Central path following algorithm

In Section 5.1, we describe the central path of the HSDE, and in Section 5.2
we define central path proximity measures. In Section 5.3, we outline a high
level PDIPM that maintains iterates close to the central path, and we give
numerical convergence criteria for detecting approximate conic certificates.
In Section 5.4, we derive prediction and centering directions and our
corresponding TOA directions using the TOO. Finally in Section 5.5, we
summarize an SY-like stepping procedure and describe our sequence of four
enhancements to this procedure.

###  5.1 Central path of the HSDE

We define the HSDE in 12. Recall that $\mathcal{K}$ in our primal conic form 7
is a Cartesian product
$\mathcal{K}=\mathcal{K}_{1}\times\cdots\times\mathcal{K}_{K}$ of $K$ exotic
cones. We partition the exotic cone indices $\llbracket K\rrbracket$ into two
sets: $K_{\text{pr}}$ for cones with primal oracles (i.e. for
$\mathcal{K}_{k}$) and $K_{\text{du}}$ for cones with dual oracles (i.e. for
$\mathcal{K}_{k}^{\ast}$). For convenience, we append the $\tau$ and $\kappa$
variables onto the $s$ and $z$ variables. Letting $\bar{K}=K+1$, we define for
$k\in\llbracket\bar{K}\rrbracket$:

$\displaystyle\bar{\mathcal{K}}_{k}$
$\displaystyle\coloneqq\begin{cases}\mathcal{K}_{k}&k\in K_{\text{pr}},\\\
\mathcal{K}_{k}^{\ast}&k\in K_{\text{du}},\\\
\mathbb{R}_{\geq}&k=\bar{K},\end{cases}$ (14a)
$\displaystyle(\bar{z}_{k},\bar{s}_{k})$
$\displaystyle\coloneqq\begin{cases}(z_{k},s_{k})&k\in K_{\text{pr}},\\\
(s_{k},z_{k})&k\in K_{\text{du}},\\\ (\kappa,\tau)&k=\bar{K}.\end{cases}$
(14b)

For a given initial interior point
$\omega^{0}=(x^{0},y^{0},z^{0},\tau^{0},s^{0},\kappa^{0})$, the central path
of the HSDE is the trajectory of solutions
$\omega_{\mu}=(x_{\mu},y_{\mu},z_{\mu},\tau_{\mu},s_{\mu},\kappa_{\mu})$,
parameterized by $\mu>0$, satisfying:

$\displaystyle E\omega_{\mu}$ $\displaystyle=\mu E\omega^{0},$ (15a)
$\displaystyle\bar{z}_{\mu,k}+\mu g_{k}(\bar{s}_{\mu,k})$
$\displaystyle=0\quad\forall k\in\llbracket\bar{K}\rrbracket,$ (15b)
$\displaystyle(\bar{z}_{\mu},\bar{s}_{\mu})$
$\displaystyle\in\operatorname{int}(\bar{\mathcal{K}}^{\ast}\times\bar{\mathcal{K}}).$
(15c)

When all exotic cones have primal oracles (i.e. $K_{\text{du}}$ is empty), our
definition 15 exactly matches the central path defined in Vandenberghe (2010,
Equation 32), and only differs from the definition in Skajaa and Ye (2015,
Equations 7-8) in the affine form (i.e. the variable names and affine
constraint structure). Unlike SY, our central path condition 15b allows cones
with dual oracles ($K_{\text{du}}$ may be nonempty).

To obtain an initial point $\omega^{0}$, we first let:

$\bigl{(}\bar{z}^{0}_{k},\bar{s}^{0}_{k}\bigr{)}=(-g_{k}(t_{k}),t_{k})\quad\forall
k\in\llbracket\bar{K}\rrbracket,$ (16)

where $t_{k}\in\operatorname{int}\bigl{(}\bar{\mathcal{K}}_{k}\bigr{)}$ is the
initial interior point oracle (note that $\tau^{0}=\kappa^{0}=1$). Although
$x^{0}$ and $y^{0}$ can be chosen arbitrarily, we let $x^{0}$ be the solution
of:

$\displaystyle\textstyle\min_{x\in\mathbb{R}^{n}}\quad\lVert x\rVert$
$\displaystyle:$ (17a) $\displaystyle-Ax+b\tau^{0}$ $\displaystyle=0,$ (17b)
$\displaystyle-Gx+h\tau^{0}-s^{0}$ $\displaystyle=0,$ (17c)

and we let $y^{0}$ be the solution of:

$\displaystyle\textstyle\min_{y\in\mathbb{R}^{p}}\quad\lVert y\rVert$
$\displaystyle:$ (18a) $\displaystyle A^{\prime}y+G^{\prime}z^{0}+c\tau^{0}$
$\displaystyle=0.$ (18b)

In Appendix A, we outline a QR-factorization-based procedure for preprocessing
the affine data of the conic model and solving for $\omega^{0}$.

Like Skajaa and Ye (2015, Section 4.1), we define the _complementarity gap_
function:

$\mu(\omega)\coloneqq\bar{s}^{\prime}\bar{z}/{\textstyle\sum_{k\in\llbracket\bar{K}\rrbracket}}\nu_{k},$
(19)

where $\nu_{k}$ is the LHSCB parameter of the LHSCB $f_{k}$ for
$\bar{\mathcal{K}}_{k}$ (see 2b). Note that $\mu(\omega)>0$ if
$(\bar{z},\bar{s})\in\bar{\mathcal{K}}^{\ast}\times\bar{\mathcal{K}}$, by a
strict version of the dual cone inequality 1. From 16, $\mu(\omega^{0})=1$,
since in 19 we have
$(\bar{s}^{0})^{\prime}\bar{z}^{0}={\textstyle\sum_{k\in\llbracket\bar{K}\rrbracket}}t_{k}^{\prime}(-g_{k}(t_{k}))$,
and $t_{k}^{\prime}(-g_{k}(t_{k}))=\nu_{k}$ by logarithmic homogeneity of
$f_{k}$ (Nesterov and Nemirovskii, 1994, Proposition 2.3.4). Hence
$\omega^{0}$ satisfies the central path conditions 15 for parameter value
$\mu=1$. The central path is therefore a trajectory that starts at
$\omega^{0}$ with complementarity gap $\mu=1$ and approaches a solution for
the HSDE as $\mu$ decreases to zero.

###  5.2 Central path proximity

Given a point $\omega$, we define the central path _proximity_ $\pi_{k}$ for
exotic cone $k\in\llbracket\bar{K}\rrbracket$ as:

$\pi_{k}(\omega)\coloneqq\begin{cases}\big{\lVert}(H_{k}(\bar{s}_{k}))^{-1/2}(\bar{z}_{k}/\mu(\omega)+g_{k}(\bar{s}_{k}))\big{\rVert}&\text{if
}\mu(\omega)>0,\bar{s}_{k}\in\operatorname{int}\bigl{(}\bar{\mathcal{K}}_{k}\bigr{)},\\\
\infty&\text{otherwise}.\end{cases}$ (20)

Hence $\pi_{k}$ is a measure of the distance from $\bar{s}_{k}$ and
$\bar{z}_{k}$ to the surface defined by the central path condition 15b
(compare to Skajaa and Ye (2015, Equation 9) and Nesterov and Todd (1998,
Section 4)).

In Lemma 5.1, we show that for exotic cone $k\in\llbracket\bar{K}\rrbracket$,
if $\pi_{k}(\omega)<1$, then
$\bar{s}_{k}\in\operatorname{int}\bigl{(}\bar{\mathcal{K}}_{k}\bigr{)}$ and
$\bar{z}_{k}\in\operatorname{int}\bigl{(}\bar{\mathcal{K}}^{\ast}_{k}\bigr{)}$.
This condition is sufficient but not necessary for strict cone feasibility. If
it holds for all $k\in\llbracket\bar{K}\rrbracket$, then $\omega$ is an
interior point (by definition) and 15c is satisfied. From 20,
$\pi_{k}(\omega)$ can be computed by evaluating the feasibility check,
gradient, and Hessian oracles for $\bar{\mathcal{K}}_{k}$ at $\bar{s}_{k}$.

######  Lemma 5.1.

Given a point $\omega$, for each $k\in\llbracket\bar{K}\rrbracket$,
$\pi_{k}(\omega)<1$ implies
$\bar{s}_{k}\in\operatorname{int}\bigl{(}\bar{\mathcal{K}}_{k}\bigr{)}$ and
$\bar{z}_{k}\in\operatorname{int}\bigl{(}\bar{\mathcal{K}}^{\ast}_{k}\bigr{)}$.

###### Proof.

We adapt Papp and Yıldız (2017, Lemma 15). Fix $\mu=\mu(\omega)$ for
convenience, and suppose $\pi_{k}(\omega)<1$ for exotic cone
$k\in\llbracket\bar{K}\rrbracket$. Then by 20, $\mu>0$ and
$\bar{s}_{k}\in\operatorname{int}\bigl{(}\bar{\mathcal{K}}_{k}\bigr{)}$. By
Papp and Yıldız (2017, Theorem 8),
$\bar{s}_{k}\in\operatorname{int}\bigl{(}\bar{\mathcal{K}}_{k}\bigr{)}$
implies
$-g_{k}(\bar{s}_{k})\in\operatorname{int}\bigl{(}\bar{\mathcal{K}}^{\ast}_{k}\bigr{)}$.
Let $f_{k}$ be the LHSCB for $\bar{\mathcal{K}}_{k}$, and let
$H^{\ast}_{k}\coloneqq\nabla^{2}f^{\ast}_{k}$ denote the Hessian operator for
the conjugate $f^{\ast}_{k}$ (see 3) of $f_{k}$. By Papp and Yıldız (2017,
Equation 13), $H^{\ast}_{k}(-g_{k}(\bar{s}_{k}))=(H_{k}(\bar{s}_{k}))^{-1}$,
so:

$\displaystyle\big{\lVert}(H^{\ast}_{k}(-g_{k}(\bar{s}_{k})))^{1/2}(\bar{z}_{k}/\mu+g_{k}(\bar{s}_{k}))\big{\rVert}$
(21a)
$\displaystyle=\big{\lVert}(H_{k}(\bar{s}_{k}))^{-1/2}(\bar{z}_{k}/\mu+g_{k}(\bar{s}_{k}))\big{\rVert}$
(21b) $\displaystyle=\pi_{k}(\omega)<1.$ (21c)

So by Papp and Yıldız (2017, Definition 1),
$\bar{z}_{k}/\mu\in\operatorname{int}\bigl{(}\bar{\mathcal{K}}^{\ast}_{k}\bigr{)}$,
hence
$\bar{z}_{k}\in\operatorname{int}\bigl{(}\bar{\mathcal{K}}^{\ast}_{k}\bigr{)}$.
∎

We now define a proximity function that aggregates the exotic cone central
path proximity values $\pi_{k}(\omega)\geq 0,\forall
k\in\llbracket\bar{K}\rrbracket$. SY aggregates by taking the $\ell_{2}$ norm:

$\pi_{\ell_{2}}(\omega)\coloneqq\big{\lVert}(\pi_{k}(\omega))_{k\in\llbracket\bar{K}\rrbracket}\big{\rVert}.$
(22)

An alternative aggregated proximity uses the $\ell_{\infty}$ norm (maximum):

$\pi_{\ell_{\infty}}(\omega)\coloneqq\big{\lVert}(\pi_{k}(\omega))_{k\in\llbracket\bar{K}\rrbracket}\big{\rVert}_{\infty}.$
(23)

Clearly,
$0\leq\pi_{k}(\omega)\leq\pi_{\ell_{\infty}}(\omega)\leq\pi_{\ell_{2}}(\omega),\forall
k\in\llbracket\bar{K}\rrbracket$. Both conditions $\pi_{\ell_{2}}(\omega)<1$
and $\pi_{\ell_{\infty}}(\omega)<1$ guarantee by Lemma 5.1 that $\omega$ is an
interior point, however using $\pi_{\ell_{2}}$ leads to a more restrictive
condition on $\omega$.

###  5.3 High level algorithm

We describe a high level algorithm for approximately solving the HSDE. The
method starts at the initial interior point $\omega^{0}$ with complementarity
gap $\mu(\omega^{0})=1$ and approximately tracks the central path trajectory
15 through a series of iterations. It maintains feasibility for the linear
equality conditions 15a and strict cone feasibility conditions 15c, but allows
violation of the nonlinear equality conditions 15b. On the $i$th iteration,
the current interior point is $\omega^{i-1}$ satisfying
$\pi_{k}(\omega^{i-1})<1,\forall k\in\llbracket\bar{K}\rrbracket$, and the
complementarity gap is $\mu(\omega^{i-1})$. The method searches for a new
point $\omega^{i}$ that maintains the proximity condition
$\pi_{k}(\omega^{i})<1,\forall k\in\llbracket\bar{K}\rrbracket$ (and hence is
an interior point) and either has a smaller complementarity gap
$\mu(\omega^{i})<\mu(\omega^{i-1})$ or a smaller aggregate proximity value
$\pi(\omega^{i})<\pi(\omega^{i-1})$ (where $\pi$ is $\pi_{\ell_{2}}$ or
$\pi_{\ell_{\infty}}$), or both. As the complementarity gap decreases towards
zero, the RHS of 15a approaches the origin, so the iterates approach a
solution of the HSDE 12.

To detect an approximate conic certificate and terminate the iterations, we
check if the current iterate $\omega$ satisfies any of the following numerical
convergence criteria. These conditions use positive tolerance values for
feasibility $\varepsilon_{f}$, infeasibility $\varepsilon_{i}$, absolute gap
$\varepsilon_{a}$, relative gap $\varepsilon_{r}$, and ill-posedness
$\varepsilon_{p}$ (see Section 7.2 for the tolerance values we use in
computational testing).

Optimality.

We terminate with a complementary solution $(x,y,z)/\tau$ approximately
satisfying the primal-dual optimality conditions 7b, 7c, 8b, 8c and 11 if:

$\max\biggl{(}\frac{\lVert
A^{\prime}y+G^{\prime}z+c\tau\rVert_{\infty}}{1+\lVert
c\rVert_{\infty}},\frac{\lVert-Ax+b\tau\rVert_{\infty}}{1+\lVert
b\rVert_{\infty}},\frac{\lVert-Gx+h\tau-s\rVert_{\infty}}{1+\lVert
h\rVert_{\infty}}\biggr{)}\leq\varepsilon_{f}\tau,$ (24a) and at least one of
the following two conditions holds: $\displaystyle s^{\prime}z$
$\displaystyle\leq\varepsilon_{a},$ (24b)
$\displaystyle\min(s^{\prime}z/\tau,\lvert
c^{\prime}x+b^{\prime}y+h^{\prime}z\rvert)$
$\displaystyle\leq\varepsilon_{r}\max(\tau,\min(\lvert
c^{\prime}x\rvert,\lvert b^{\prime}y+h^{\prime}z\rvert)).$ (24c)

Note that 24b and 24c are absolute and relative optimality gap conditions
respectively.

Primal infeasibility.

We terminate with a dual improving ray $(y,z)$ approximately satisfying 10 if:

$b^{\prime}y+h^{\prime}z<0,\qquad\lVert
A^{\prime}y+G^{\prime}z\rVert_{\infty}\leq-\varepsilon_{i}(b^{\prime}y+h^{\prime}z).$
(25)

Dual infeasibility.

We terminate with a primal improving ray $x$ approximately satisfying 9 if:

$c^{\prime}x<0,\qquad\max(\lVert Ax\rVert_{\infty},\lVert
Gx+s\rVert_{\infty})\leq-\varepsilon_{i}c^{\prime}x.$ (26)

Ill-posed primal or dual.

If $\tau$ and $\kappa$ are approximately $0$, the primal and dual problem
statuses cannot be determined (see Section 4.3). We terminate with an ill-
posed status if:

$\mu(\omega)\leq\varepsilon_{p},\qquad\tau\leq\varepsilon_{p}\min(1,\kappa).$
(27)

The high level path following algorithm below computes an approximate solution
to the HSDE. In Section 5.5, we describe specific stepping procedures for 5.

1:procedure SolveHSDE

2: compute initial interior point $\omega^{0}$

3: $i\leftarrow 1$

4: while $\omega^{i-1}$ does not satisfy any of the convergence conditions 24
to 27 do

5: $\omega^{i}\leftarrow$ Step($\omega^{i-1}$)

6: $i\leftarrow i+1$

7: end while

8: return $\omega^{i}$

9:end procedure

###  5.4 Search directions

At a given iteration of the path following method, let $\omega$ be the current
interior point and fix $\mu=\mu(\omega)$ for convenience. The stepping
procedures we describe in Section 5.5 first compute one or more search
directions, which depend on $\omega$. We derive the _centering_ direction in
Section 5.4.1 and the _prediction_ direction in Section 5.4.2. The goal of
centering is to step to a point with a smaller aggregate central path
proximity than the current point, i.e. to step towards the central path. The
goal of prediction is to step to a point with a smaller complementarity gap,
i.e. to step closer to a solution of the HSDE. The centering and prediction
directions match those used by SY. We associate with each of these directions
a new _third order adjustment_ (TOA) direction, which depends on the TOO and
helps to correct the corresponding unadjusted direction (which must be
computed before the TOA direction). Hence we derive four types of directions
here.

Each direction is computed as the solution to a linear system with a
structured square $6\times 6$ block matrix left hand side (LHS) and a
particular right hand side (RHS) vector. The LHS, which depends only on
$\omega$ and the problem data, is the same for all four directions at a given
iteration. We let
$r\coloneqq(r_{E},r_{1},\ldots,r_{\bar{K}})\in\mathbb{R}^{\dim(\omega)}$
represent an RHS, where $r_{E}\in\mathbb{R}^{n+p+q+1}$ corresponds to the
linear equalities 15a and $r_{k}\in\mathbb{R}^{q_{k}},\forall
k\in\llbracket\bar{K}\rrbracket$ corresponds to the nonlinear equalities 15b.
The direction
$\delta\coloneqq(\delta_{x},\delta_{y},\delta_{z},\delta_{\tau},\delta_{s},\delta_{\kappa})\in\mathbb{R}^{\dim(\omega)}$
corresponding to $r$ is the solution to:

$\displaystyle E\delta$ $\displaystyle=r_{E},$ (28a)
$\displaystyle\delta_{\bar{z},k}+\mu H_{k}(\bar{s}_{k})\delta_{\bar{s},k}$
$\displaystyle=r_{k}\quad\forall k\in\llbracket\bar{K}\rrbracket.$ (28b)

Since $E$ is assumed to have full row rank and each $H_{k}$ is positive
definite, this square system is nonsingular and hence has a unique solution.
In Appendix A, we describe a particular method for solving 28.

####  5.4.1 Centering

The centering direction $\delta^{c}$ is analogous to the definition of Skajaa
and Ye (2015, Section 3.2). It reduces the violation on the central path
nonlinear equality condition 15b (and can be interpreted as a Newton step),
while keeping the complementarity gap $\mu$ (approximately) constant. We
denote the centering TOA direction $\delta^{ct}$. To maintain feasibility for
the linear equality condition 15a, we ensure $E\delta^{c}=E\delta^{ct}=0$ in
28a.

Dropping the index $k\in\llbracket\bar{K}\rrbracket$ for conciseness, recall
that 15b expresses $\bar{z}+\mu g(\bar{s})=0$. A first order approximation of
this condition gives:

$\displaystyle\bar{z}+\delta_{\bar{z}}+\mu(g(\bar{s})+H(\bar{s})\delta_{\bar{s}})$
$\displaystyle=0$ (29a) $\displaystyle\Rightarrow\quad\delta_{\bar{z}}+\mu
H(\bar{s})\delta_{\bar{s}}$ $\displaystyle=-\bar{z}-\mu g(\bar{s}),$ (29b)

which matches the form of 28b. Hence we let the centering direction
$\delta^{c}$ be the solution to:

$\displaystyle E\delta$ $\displaystyle=0,$ (30a)
$\displaystyle\delta_{\bar{z},k}+\mu H_{k}(\bar{s}_{k})\delta_{\bar{s},k}$
$\displaystyle=-\bar{z}_{k}-\mu g_{k}(\bar{s}_{k})\quad\forall
k\in\llbracket\bar{K}\rrbracket.$ (30b)

Similarly, a second order approximation of $\bar{z}+\mu g(\bar{s})=0$ gives:

$\displaystyle\bar{z}+\delta_{\bar{z}}+\mu\bigl{(}g(\bar{s})+H(\bar{s})\delta_{\bar{s}}+\tfrac{1}{2}\nabla^{3}f(\bar{s})[\delta_{\bar{s}},\delta_{\bar{s}}]\bigr{)}$
$\displaystyle=0$ (31a) $\displaystyle\Rightarrow\quad\delta_{\bar{z}}+\mu
H(\bar{s})\delta_{\bar{s}}$ $\displaystyle=-\bar{z}-\mu
g(\bar{s})+\mu\mathrm{T}(\bar{s},\delta_{\bar{s}}),$ (31b)

where 31b uses the definition of the TOO in 5. Note that the RHSs of 29b and
31b differ only by $\mu\mathrm{T}(\bar{s},\delta_{\bar{s}})$, which depends on
$\delta_{\bar{s}}$. To remove this dependency, we substitute the centering
direction $\delta^{c}$, which we assume is already computed, into the RHS of
31b. Hence we let the centering TOA direction $\delta^{ct}$, which adjusts the
centering direction, be the solution to:

$\displaystyle E\delta$ $\displaystyle=0,$ (32a)
$\displaystyle\delta_{\bar{z},k}+\mu H_{k}(\bar{s}_{k})\delta_{\bar{s},k}$
$\displaystyle=\mu\mathrm{T}_{k}\bigl{(}\bar{s}_{k},\delta^{c}_{\bar{s},k}\bigr{)}\quad\forall
k\in\llbracket\bar{K}\rrbracket.$ (32b)

We note that for a rescaling factor $\alpha\in(0,1)$, the TOA direction
corresponding to $\alpha\delta^{c}$ (a rescaling of the centering direction)
is $\alpha^{2}\delta^{ct}$ (a rescaling of the centering TOA direction).

####  5.4.2 Prediction

The prediction direction $\delta^{p}$ reduces the complementarity gap and is
analogous to the definition of Skajaa and Ye (2015, Section 3.1). We derive
$\delta^{p}$ and its corresponding TOA direction $\delta^{pt}$ by considering
the central path conditions 15 as a dynamical system parametrized by $\mu>0$,
and differentiating the linear and nonlinear equalities 15a and 15b.

Differentiating 15a once gives:

$E\dot{\omega}_{\mu}=E\omega^{0}.$ (33)

Rescaling 33 by $-\mu$ and substituting 15a gives:

$E(-\mu\dot{\omega}_{\mu})=-\mu E\omega^{0}=-E\omega_{\mu}.$ (34)

Dropping the index $k\in\llbracket\bar{K}\rrbracket$ for conciseness, we
differentiate $\bar{z}_{\mu}+\mu g(\bar{s}_{\mu})=0$ from 15b once to get:

$\dot{\bar{z}}_{\mu}+g(\bar{s}_{\mu})+\mu
H(\bar{s}_{\mu})\dot{\bar{s}}_{\mu}=0.$ (35)

Rescaling 35 by $-\mu$ and substituting $\bar{z}_{\mu}=-\mu g(\bar{s}_{\mu})$
from 15b gives:

${-\mu}\dot{\bar{z}}_{\mu}+\mu
H(\bar{s}_{\mu})(-\mu\dot{\bar{s}}_{\mu})=-\bar{z}_{\mu}.$ (36)

The direction $\dot{\omega}_{\mu}$ is tangent to the central path. Like SY, we
interpret the prediction direction as $\delta^{p}=-\mu\dot{\omega}_{\mu}$, so
34 and 36 become:

$\displaystyle E\delta^{p}$ $\displaystyle=-E\omega_{\mu},$ (37a)
$\displaystyle\delta^{p}_{\bar{z}}+\mu H(\bar{s}_{\mu})\delta^{p}_{\bar{s}}$
$\displaystyle=-\bar{z}_{\mu},$ (37b)

which matches the form 28. So we let $\delta^{p}$ be the solution to:

$\displaystyle E\delta$ $\displaystyle=-E\omega,$ (38a)
$\displaystyle\delta_{\bar{z},k}+\mu H_{k}(\bar{s}_{k})\delta_{\bar{s},k}$
$\displaystyle=-\bar{z}_{k}\quad\forall k\in\llbracket\bar{K}\rrbracket.$
(38b)

Differentiating 15a twice and rescaling by $\frac{1}{2}\mu^{2}$ gives:

$E\bigl{(}\tfrac{1}{2}\mu^{2}\ddot{\omega}_{\mu}\bigr{)}=0.$ (39)

Differentiating $\bar{z}_{\mu}+\mu g(\bar{s}_{\mu})=0$ twice gives:

$\ddot{\bar{z}}_{\mu}+2H(\bar{s}_{\mu})\dot{\bar{s}}_{\mu}+\mu\nabla^{3}f(\bar{s}_{\mu})[\dot{\bar{s}}_{\mu},\dot{\bar{s}}_{\mu}]+\mu
H(\bar{s}_{\mu})\ddot{\bar{s}}_{\mu}=0.$ (40)

Rescaling 40 by $\frac{1}{2}\mu^{2}$ and substituting the TOO definition 5, we
have:

$\displaystyle\tfrac{1}{2}\mu^{2}\ddot{\bar{z}}_{\mu}+\mu
H(\bar{s}_{\mu})\bigl{(}\tfrac{1}{2}\mu^{2}\ddot{\bar{s}}_{\mu}\bigr{)}$
$\displaystyle=\mu
H(\bar{s}_{\mu})(-\mu\dot{\bar{s}}_{\mu})-\tfrac{1}{2}\mu\nabla^{3}f(\bar{s}_{\mu})[-\mu\dot{\bar{s}}_{\mu},-\mu\dot{\bar{s}}_{\mu}]$
(41a) $\displaystyle=\mu
H(\bar{s}_{\mu})(-\mu\dot{\bar{s}}_{\mu})+\mu\mathrm{T}(\bar{s}_{\mu},-\mu\dot{\bar{s}}_{\mu}).$
(41b)

We interpret the prediction TOA direction, which adjusts the prediction
direction, as $\delta^{pt}=\frac{1}{2}\mu^{2}\ddot{\omega}$. The RHS of 41b
depends on $\dot{\bar{s}}_{\mu}$, so we remove this dependency by substituting
the prediction direction $\delta^{p}=-\mu\dot{\omega}_{\mu}$, which we assume
is already computed. Hence using 39 and 41b, we let $\delta^{pt}$ be the
solution to:

$\displaystyle E\delta$ $\displaystyle=0,$ (42a)
$\displaystyle\delta_{\bar{z},k}+\mu H_{k}(\bar{s}_{k})\delta_{\bar{s},k}$
$\displaystyle=\mu
H_{k}(\bar{s}_{k})\delta^{p}_{\bar{s},k}+\mu\mathrm{T}_{k}\bigl{(}\bar{s}_{k},\delta^{p}_{\bar{s},k}\bigr{)}\quad\forall
k\in\llbracket\bar{K}\rrbracket.$ (42b)

We note that the RHS in 42b differs from the ‘higher order corrector’ RHS
proposed by Dahl and Andersen (2021, Equation 16), which has the form
$\frac{1}{2}\nabla^{3}f_{k}\bigl{[}\delta_{\bar{s},k}^{p},(H_{k}(\bar{s}_{k}))^{-1}\delta_{\bar{z},k}^{p}\bigr{]}$.
For example, our form does not satisfy all of the properties in Dahl and
Andersen (2021, Lemmas 3 and 4).

###  5.5 Stepping procedures

A stepping procedure computes one or more directions from Section 5.4 and uses
the directions to search for a new interior point. Recall from 5 of the high
level PDIPM in Section 5.3 that on iteration $i$ with current iterate
$\omega^{i-1}$, Step computes $\omega^{i}$ satisfying $\pi(\omega^{i})<1$ and
either $\mu(\omega^{i})<\mu(\omega^{i-1})$ (prediction) or
$\pi(\omega^{i})<\pi(\omega^{i-1})$ (centering) or both. In Section 5.5.1, we
describe a stepping procedure similar to that of Alfonso (Papp and Yıldız,
2021), which is a practical implementation of SY. This procedure, which we
call _basic_ , alternates between prediction and centering steps and does not
use the TOA directions. In Sections 5.5.2 to 5.5.5, we describe a sequence of
four cumulative enhancements to the _basic_ procedure, with the goal of
improving iteration counts and per-iteration computational efficiency in
practice. The main purpose of our computational testing in Section 7 is to
assess the value of these enhancements on a diverse set of benchmark
instances.

####  5.5.1 Basic stepping procedure

First, we decide whether to perform a centering step or a prediction step. If
the current iterate $\omega^{i-1}$ (at the $i$th iteration) is very close to
the central path, i.e. if the sum proximity 22 does not exceed $\eta=0.0332$
(from Alfonso (Papp and Yıldız, 2020)), or if the most recent $N=4$ steps have
all been centering steps, then we compute the prediction direction
$\delta^{p}$ from 38. Otherwise, we compute the centering direction
$\delta^{c}$ from 30. Letting $j$ be the number of consecutive centering steps
taken immediately before the current $i$th iteration, the search direction is:

$\delta\coloneqq\begin{cases}\delta^{p}&\text{if
$\pi_{\ell_{2}}(\omega^{i-1})\leq\eta$ or $j\geq N$},\\\
\delta^{c}&\text{otherwise}.\end{cases}$ (43)

Next, we perform a backtracking line search in the direction $\delta$. The
search finds a step length $\hat{\alpha}\in(0,1)$ from a fixed schedule of
decreasing values $\mathcal{A}=\\{\alpha_{l}\\}_{l\in\llbracket L\rrbracket}$,
where $L=18$, $\alpha_{1}=0.9999$, and $\alpha_{L}=0.0005$. The next iterate
$\omega^{i}=\omega^{i-1}+\hat{\alpha}\delta$ becomes the first point in the
backtracking line search that satisfies
$\pi_{\ell_{2}}(\omega^{i})\leq\beta_{1}$ for $\beta_{1}=0.2844$ (from Alfonso
(Papp and Yıldız, 2020)), which guarantees interiority by Lemma 5.1. If the
backtracking search terminates without a step length satisfying the proximity
condition (i.e. $\alpha_{L}$ is too large), the PDIPM algorithm terminates
without a solution. In Appendix B we discuss our implementation of the
proximity check that we run for each candidate point in the backtracking
search.

The _basic_ stepping procedure is summarized as follows. Note the centering
step count $j$ is initialized to zero before the first iteration $i=1$. Since
$\omega^{0}$ is exactly on the central path (i.e. the proximity is zero), the
first iteration uses a prediction step.

1:procedure BasicStep($\omega^{i-1}$, $j$)

2: if $\pi_{\ell_{2}}(\omega^{i-1})\leq\eta$ or $j\geq N$ then
$\triangleright$ choose predict or center

3: $\delta\leftarrow\delta^{p}$ from 38 $\triangleright$ compute prediction
direction

4: $j\leftarrow 0$

5: else

6: $\delta\leftarrow\delta^{c}$ from 30 $\triangleright$ compute centering
direction

7: $j\leftarrow j+1$

8: end if

9:
$\hat{\alpha}\leftarrow\max\\{\alpha\in\mathcal{A}:\pi_{\ell_{2}}(\omega^{i-1}+\alpha\delta)\leq\beta_{1}\bigr{\\}}$
$\triangleright$ compute step length by backtracking search

10: $\omega^{i}\leftarrow\omega^{i-1}+\hat{\alpha}\delta$ $\triangleright$
update current iterate

11:end procedure

####  5.5.2 Less restrictive proximity

The _basic_ stepping procedure in Section 5.5.1 requires iterates to remain in
close proximity to the central path and usually only takes prediction steps
from iterates that are very close to the central path. Although conservative
proximity conditions are used to prove polynomial iteration complexity in Papp
and Yıldız (2017), they may be too restrictive from the perspective of
practical performance. To allow prediction steps from a larger neighborhood of
the central path, we use the $\pi_{\ell_{\infty}}$ proximity measure from 23
instead of $\pi_{\ell_{2}}$ to compute the proximity of $\omega^{i-1}$, though
we do not change the proximity bound $\eta$. To allow longer step lengths, we
also use $\pi_{\ell_{\infty}}$ instead of $\pi_{\ell_{2}}$ for the
backtracking search proximity checks, and we increase this proximity bound to
$\beta_{2}=0.99$ (by Lemma 5.1, $\beta_{2}<1$ guarantees interiority).

The _prox_ stepping procedure, which enhances the _basic_ stepping procedure
by relaxing the proximity conditions somewhat, is summarized as follows.

1:procedure ProxStep($\omega^{i-1}$, $j$)

2: if $\pi_{\ell_{\infty}}(\omega^{i-1})\leq\eta$ or $j\geq N$ then
$\triangleright$ use less restrictive proximity measure $\pi_{\ell_{\infty}}$

3: $\delta\leftarrow\delta^{p}$ from 38

4: $j\leftarrow 0$

5: else

6: $\delta\leftarrow\delta^{c}$ from 30

7: $j\leftarrow j+1$

8: end if

9:
$\hat{\alpha}\leftarrow\max\\{\alpha\in\mathcal{A}:\pi_{\ell_{\infty}}(\omega^{i-1}+\alpha\delta)\leq\beta_{2}\\}$
$\triangleright$ use $\pi_{\ell_{\infty}}$ and larger proximity bound
$\beta_{2}$

10: $\omega^{i}\leftarrow\omega^{i-1}+\hat{\alpha}\delta$

11:end procedure

####  5.5.3 Third order adjustments

We modify the _prox_ stepping procedure in Section 5.5.2 to incorporate the
new TOA directions associated with the prediction and centering directions.
After deciding whether to predict or center (using the same criteria as
_prox_), we compute the unadjusted direction $\delta^{u}$ (i.e. $\delta^{p}$
or $\delta^{c}$) and its associated TOA direction $\delta^{t}$ (i.e.
$\delta^{pt}$ or $\delta^{ct}$). We perform a backtracking line search in
direction $\delta^{u}$, just like _prox_ , and we use this step length
$\hat{\alpha}^{u}\in(0,1)$ to scale down the TOA direction. We let the final
direction be $\delta^{u}+\hat{\alpha}^{u}\delta^{t}$. The rescaling of
$\delta^{t}$ helps to prevent over-adjustment. Finally, we perform a second
backtracking line search, using the same techniques and proximity condition as
the first line search.

The _TOA_ stepping procedure, which enhances the _prox_ stepping procedure by
incorporating the TOA directions, is summarized as follows.

1:procedure TOAStep($\omega^{i-1}$, $j$)

2: if $\pi_{\ell_{\infty}}(\omega^{i-1})\leq\eta$ or $j\geq N$ then

3: $\delta^{u}\leftarrow\delta^{p}$ from 38

4: $\delta^{t}\leftarrow\delta^{pt}$ from 42 $\triangleright$ compute
prediction TOA direction

5: $j\leftarrow 0$

6: else

7: $\delta^{u}\leftarrow\delta^{c}$ from 30

8: $\delta^{t}\leftarrow\delta^{ct}$ from 32 $\triangleright$ compute
centering TOA direction

9: $j\leftarrow j+1$

10: end if

11:
$\hat{\alpha}^{u}\leftarrow\max\\{\alpha\in\mathcal{A}:\pi_{\ell_{\infty}}(\omega^{i-1}+\alpha\delta^{u})\leq\beta_{2}\\}$
$\triangleright$ perform line search for unadjusted direction

12: $\delta\leftarrow\delta^{u}+\hat{\alpha}^{u}\delta^{t}$ $\triangleright$
compute final direction

13:
$\hat{\alpha}\leftarrow\max\\{\alpha\in\mathcal{A}:\pi_{\ell_{\infty}}(\omega^{i-1}+\alpha\delta)\leq\beta_{2}\\}$

14: $\omega^{i}\leftarrow\omega^{i-1}+\hat{\alpha}\delta$

15:end procedure

####  5.5.4 Curve search

The _TOA_ stepping procedure in Section 5.5.3 performs two backtracking line
searches, which can be quite expensive. We propose using a single backtracking
search along a curve that is quadratic in the step parameter $\alpha$ and
linear in the unadjusted and TOA directions. Recall from 12 of the _TOA_
procedure that we compute a direction $\delta$ as a linear function of the
step parameter from the first line search. Substituting this $\delta$ function
into the usual linear trajectory gives the curved trajectory
$\omega^{i-1}+\alpha(\delta^{u}+\alpha\delta^{t})$ for $\alpha\in(0,1)$, where
$\delta^{u}$ and $\delta^{t}$ are the unadjusted and TOA directions (as in the
_TOA_ procedure). Intuitively, a backtracking search along this curve achieves
a more dynamic rescaling of the TOA direction.

The _curve_ stepping procedure, which enhances the _TOA_ stepping procedure by
using a search on a curve instead of two line searches, is summarized as
follows.

1:procedure CurveStep($\omega^{i-1}$, $j$)

2: if $\pi_{\ell_{\infty}}(\omega^{i-1})\leq\eta$ or $j\geq N$ then

3: $\delta^{u}\leftarrow\delta^{p}$ from 38

4: $\delta^{t}\leftarrow\delta^{pt}$ from 42

5: $j\leftarrow 0$

6: else

7: $\delta^{u}\leftarrow\delta^{c}$ from 30

8: $\delta^{t}\leftarrow\delta^{ct}$ from 32

9: $j\leftarrow j+1$

10: end if

11: let
$\hat{\omega}(\alpha)\coloneqq\omega^{i-1}+\alpha(\delta^{u}+\alpha\delta^{t})$
$\triangleright$ use curved trajectory

12:
$\hat{\alpha}\leftarrow\max\\{\alpha\in\mathcal{A}:\pi_{\ell_{\infty}}(\hat{\omega}(\alpha))\leq\beta_{2}\\}$

13: $\omega^{i}\leftarrow\hat{\omega}(\hat{\alpha})$

14:end procedure

####  5.5.5 Combined directions

Unlike Skajaa and Ye (2015); Papp and Yıldız (2021), most conic PDIPMs combine
the prediction and centering phases (e.g. Vandenberghe (2010); Dahl and
Andersen (2021)). We propose using a single search on a curve that is
quadratic in the step parameter $\alpha$ and linear in all four directions
$\delta^{c},\delta^{ct},\delta^{p},\delta^{pt}$ from Section 5.5.3.
Intuitively, we can step further in a convex combination of the prediction and
centering directions than we can in just the prediction direction. In
practice, a step length of one is usually ideal for the centering phase, so we
can imagine performing a backtracking search from the point obtained from a
pure prediction step (with step length one) towards the point obtained from a
pure centering step, terminating when we are close enough to the centering
point to satisfy the proximity condition. This approach fundamentally differs
from the previous procedures we have described because the search trajectory
does not finish at the current iterate $\omega^{i-1}$. If
$\hat{\omega}^{p}(\alpha)$ and $\hat{\omega}^{c}(\alpha)$ are the prediction
and centering curve search trajectories from 11 of the _curve_ procedure, then
we define the combined trajectory as
$\hat{\omega}(\alpha)=\hat{\omega}^{p}(\alpha)+\hat{\omega}^{c}(1-\alpha)$.
Note that $\alpha=1$ corresponds to a full step in the adjusted prediction
direction $\delta^{p}+\delta^{pt}$, and $\alpha=0$ corresponds to a full step
in the adjusted centering direction $\delta^{c}+\delta^{ct}$.

The _comb_ stepping procedure, which enhances the _curve_ stepping procedure
by combining the prediction and centering phases, is summarized as follows.
Note that unlike the previous procedures, there is no parameter $j$ counting
consecutive centering steps. Occasionally in practice, the backtracking search
on 4 below fails to find a positive step value, in which case we perform a
centering step according to 11 to 13 of the _curve_ procedure.

1:procedure CombStep($\omega^{i-1}$)

2: compute $\delta^{c},\delta^{ct},\delta^{p},\delta^{pt}$ from 30, 32, 38 and
42 $\triangleright$ use four directions instead of two

3: let
$\hat{\omega}(\alpha)\coloneqq\omega^{i-1}+\alpha(\delta^{p}+\alpha\delta^{pt})+(1-\alpha)(\delta^{c}+(1-\alpha)\delta^{ct})$
$\triangleright$ use combined trajectory

4:
$\hat{\alpha}\leftarrow\max\\{\alpha\in\mathcal{A}:\pi_{\ell_{\infty}}(\hat{\omega}(\alpha))\leq\beta_{2}\\}$

5: $\omega^{i}\leftarrow\hat{\omega}(\hat{\alpha})$

6:end procedure

##  6 Oracles for predefined exotic cones

Below we list 23 exotic cone types that we have predefined through Hypatia’s
generic cone interface (see Section 3). Each of these cones is represented in
the benchmark set of conic instances that we introduce in Section 7.1. Recall
that we write any exotic cone $\mathcal{K}$ in vectorized form, i.e. as a
subset of $\mathbb{R}^{q}$, where $q=\dim(\mathcal{K})\geq 1$ is the cone
dimension. For cones typically defined using symmetric matrices, we use the
standard svec vectorization (see Section 2) to ensure the vectorized cone is
proper, to preserve inner products, and to simplify the dual cone definition.
Each cone is parametrized by at least one dimension and several cones have
additional parameters such as numerical data. For convenience, we drop these
parameters from the symbols we use to represent cone types. For several cones,
we have implemented additional variants over complex numbers (for example, a
Hermitian PSD cone), but we omit these definitions here for simplicity. We
defer a more complete description of Hypatia’s exotic cones and LHSCBs to Coey
et al. (2021a, c); Kapelevich et al. (2021).

Nonnegative cone.

$\mathcal{K}_{\geq}\coloneqq\mathbb{R}_{\geq}^{d}$ is the (self-dual)
nonnegative real vectors (note for $d>1$, $\mathcal{K}_{\geq}$ is not a
primitive cone).

PSD cone.

$\mathcal{K}_{\succeq}\coloneqq\bigl{\\{}w\in\mathbb{R}^{\operatorname{sd}(d)}:\operatorname{mat}(w)\in\mathbb{S}^{d}_{\succeq}\bigr{\\}}$
is the (self-dual) PSD matrices of side dimension $d$.

Doubly nonnegative cone.

$\mathcal{K}_{\operatorname{DNN}}\coloneqq\mathcal{K}_{\geq}\cap\mathcal{K}_{\succeq}$
is the PSD matrices with all nonnegative entries of side dimension $d$.

Sparse PSD cone.

$\mathcal{K}_{\operatorname{sPSD}}$ is the PSD matrices of side dimension $s$
with a fixed sparsity pattern $\mathcal{S}$ containing $d\geq s$ nonzeros
(including all diagonal elements); see Section C.4. The dual cone
$\mathcal{K}_{\operatorname{sPSD}}^{\ast}$ is the symmetric matrices with
pattern $\mathcal{S}$ for which there exists a PSD completion, i.e. an
assignment of the elements not in $\mathcal{S}$ such that the full matrix is
PSD. For simplicity, the complexity estimates in Table 1 assume the nonzeros
are grouped under $J\geq 1$ supernodes, each containing at most $l$ nodes, and
the monotone degree of each node is no greater than a constant $D$ (Andersen
et al., 2013).

Linear matrix inequality cone.

$\mathcal{K}_{\operatorname{LMI}}\coloneqq\bigl{\\{}w\in\mathbb{R}^{d}:{\textstyle\sum_{i\in\llbracket
d\rrbracket}}w_{i}P_{i}\in\mathbb{S}^{s}_{\succeq}\bigr{\\}}$ are the vectors
for which the matrix pencil of $d$ matrices $P_{i}\in\mathbb{S}^{s},\forall
i\in\llbracket d\rrbracket$ is PSD. We assume $P_{1}\succ 0$ so that we can
use the initial interior point $e_{1}$.

Infinity norm cone.

$\mathcal{K}_{\operatorname{\ell_{\infty}}}\coloneqq\\{(u,w)\in\mathbb{R}_{\geq}\times\mathbb{R}^{d}:u\geq\lVert
w\rVert_{\infty}\\}$ is the epigraph of the $\ell_{\infty}$ norm on
$\mathbb{R}^{d}$. Similarly, the dual cone
$\mathcal{K}_{\operatorname{\ell_{\infty}}}^{\ast}$ is the epigraph of the
$\ell_{1}$ norm.

Euclidean norm cone.

$\mathcal{K}_{\operatorname{\ell_{2}}}\coloneqq\\{(u,w)\in\mathbb{R}_{\geq}\times\mathbb{R}^{d}:u\geq\lVert
w\rVert\\}$ is the (self-dual) epigraph of the $\ell_{2}$ norm on
$\mathbb{R}^{d}$ (AKA second-order cone).

Euclidean norm square cone.

$\mathcal{K}_{\operatorname{sqr}}\coloneqq\\{(u,v,w)\in\mathbb{R}_{\geq}\times\mathbb{R}_{\geq}\times\mathbb{R}^{d}:2uv\geq\lVert
w\rVert^{2}\\}$ is the (self-dual) epigraph of the perspective of the square
of the $\ell_{2}$ norm on $\mathbb{R}^{d}$ (AKA rotated second-order cone).

Spectral norm cone.

$\mathcal{K}_{\operatorname{\ell_{\operatorname{spec}}}}\coloneqq\\{(u,w)\in\mathbb{R}_{\geq}\times\mathbb{R}^{rs}:u\geq\sigma_{1}(\operatorname{mat}(w))\\}$,
where $\sigma_{1}$ is the largest singular value function, is the epigraph of
the spectral norm on $\mathbb{R}^{r\times s}$, assuming $r\leq s$ without loss
of generality. Similarly,
$\mathcal{K}_{\operatorname{\ell_{\operatorname{spec}}}}^{\ast}$ is the
epigraph of the matrix nuclear norm (i.e. the sum of singular values).

Matrix square cone.

$\mathcal{K}_{\operatorname{matsqr}}\coloneqq\bigl{\\{}(u,v,w)\in\mathbb{R}^{\operatorname{sd}(r)}\times\mathbb{R}_{\geq}\times\mathbb{R}^{rs}:U\in\mathbb{S}^{r}_{\succeq},2Uv\succeq
WW^{\prime}\bigr{\\}}$, where $U\coloneqq\operatorname{mat}(u)$ and
$W\coloneqq\operatorname{mat}(w)\in\mathbb{R}^{r\times s}$, is the homogenized
symmetric matrix epigraph of the symmetric outer product, assuming $r\leq s$
without loss of generality (Güler and Tunçel, 1998).

Generalized power cone.

$\mathcal{K}_{\operatorname{gpow}}\coloneqq\bigl{\\{}(u,w)\in\mathbb{R}^{r}_{\geq}\times\mathbb{R}^{s}:{\textstyle\prod_{i\in\llbracket
r\rrbracket}}u_{i}^{\alpha_{i}}\geq\lVert w\rVert\bigr{\\}}$, parametrized by
$\alpha\in\mathbb{R}^{r}_{>}$ with $e^{\prime}\alpha=1$, is the generalized
power cone (Chares, 2009, Section 3.1.2).

Power mean cone.

$\mathcal{K}_{\operatorname{pow}}\coloneqq\bigl{\\{}(u,w)\in\mathbb{R}\times\mathbb{R}^{d}_{\geq}:u\leq{\textstyle\prod_{i\in\llbracket
d\rrbracket}}w_{i}^{\alpha_{i}}\bigr{\\}}$, parametrized by exponents
$\alpha\in\mathbb{R}^{d}_{>}$ with $e^{\prime}\alpha=1$, is the hypograph of
the power mean on $\mathbb{R}_{\geq}^{d}$.

Geometric mean cone.

$\mathcal{K}_{\operatorname{geo}}$ is the hypograph of the geometric mean on
$\mathbb{R}_{\geq}^{d}$, a special case of $\mathcal{K}_{\operatorname{pow}}$
with equal exponents.

Root-determinant cone.

$\mathcal{K}_{\operatorname{rtdet}}\coloneqq\bigl{\\{}(u,w)\in\mathbb{R}\times\mathbb{R}^{\operatorname{sd}(d)}:W\in\mathbb{S}^{d}_{\succeq},u\leq(\det(W))^{1/d}\bigr{\\}}$,
where $W\coloneqq\operatorname{mat}(w)$, is the hypograph of the $d$th-root-
determinant on $\mathbb{S}_{\succeq}^{d}$.

Logarithm cone.

$\mathcal{K}_{\log}\coloneqq\operatorname{cl}\bigl{\\{}(u,v,w)\in\mathbb{R}\times\mathbb{R}_{>}\times\mathbb{R}^{d}_{>}:u\leq{\textstyle\sum_{i\in\llbracket
d\rrbracket}}v\log(w_{i}/v)\bigr{\\}}$ is the hypograph of the perspective of
the sum of logarithms on $\mathbb{R}_{>}^{d}$.

Log-determinant cone.

$\mathcal{K}_{\operatorname{logdet}}\coloneqq\operatorname{cl}\bigl{\\{}(u,v,w)\in\mathbb{R}\times\mathbb{R}_{>}\times\mathbb{R}^{\operatorname{sd}(d)}:W\in\mathbb{S}^{d}_{\succ},u\leq
v\operatorname{logdet}(W/v)\bigr{\\}}$, where
$W\coloneqq\operatorname{mat}(w)$, is the hypograph of the perspective of the
log-determinant on $\mathbb{S}_{\succ}^{d}$.

Separable spectral function cone.

$\mathcal{K}_{\operatorname{sepspec}}\coloneqq\operatorname{cl}\\{(u,v,w)\in\mathbb{R}\times\mathbb{R}_{>}\times\operatorname{int}(\mathcal{Q}):u\geq
v\varphi(w/v)\\}$, where $\mathcal{Q}$ is $\mathcal{K}_{\geq}$ or
$\mathcal{K}_{\succeq}$ (a cone of squares of a Jordan algebra), is the
epigraph of the perspective of a convex separable spectral function
$\varphi:\operatorname{int}(\mathcal{Q})\to\mathbb{R}$, such as the sum or
trace of the negative logarithm, negative entropy, or power in $(1,2]$ (see
Coey et al. (2021c) for more details). The complexity estimates in Table 1
depend on whether $\mathcal{Q}$ is $\mathcal{K}_{\geq}$ or
$\mathcal{K}_{\succeq}$.

Relative entropy cone.

$\mathcal{K}_{\operatorname{relent}}\coloneqq\operatorname{cl}\big{\\{}(u,v,w)\in\mathbb{R}\times\mathbb{R}^{d}_{>}\times\mathbb{R}^{d}_{>}:\allowbreak
u\geq{\textstyle\sum_{i\in\llbracket
d\rrbracket}}w_{i}\log(w_{i}/v_{i})\big{\\}}$ is the epigraph of vector
relative entropy.

Matrix relative entropy cone.

$\mathcal{K}_{\operatorname{matrelent}}\coloneqq\operatorname{cl}\bigl{\\{}(u,v,w)\in\mathbb{R}\times\mathbb{R}^{\operatorname{sd}(d)}\times\mathbb{R}^{\operatorname{sd}(d)}:V\in\mathbb{S}^{d}_{\succ},W\in\mathbb{S}^{d}_{\succ},u\geq\operatorname{tr}(W(\log(W)-\log(V)))\bigr{\\}}$,
where $V\coloneqq\operatorname{mat}(v)$ and $W\coloneqq\operatorname{mat}(w)$,
is the epigraph of matrix relative entropy.444 The logarithmically homogeneous
barrier for $\mathcal{K}_{\operatorname{matrelent}}$ that Hypatia uses is
conjectured by Karimi and Tunçel (2020) to be self-concordant.

Weighted sum-of-squares (WSOS) cones.

An interpolant basis represents a polynomial implicitly by its evaluations at
a fixed set of $d$ points. Given a basic semialgebraic domain defined by $r$
polynomial inequalities, the four WSOS cones below are parameterized by
matrices $P_{l}\in\mathbb{R}^{d\times s_{l}}$ for $l\in\llbracket
r\rrbracket$. Each $P_{l}$ is constructed by evaluating $s_{l}$ independent
polynomials (columns) at the $d$ points (rows), following Papp and Yildiz
(2019). For simplicity, the complexity estimates in Table 1 assume
$s_{l}=s,\forall l\in\llbracket r\rrbracket$. Note that $s<d\leq s^{2}$. We
define $\mathcal{K}_{\operatorname{SOS}}$ and
$\mathcal{K}_{\operatorname{matSOS}}$ in Coey et al. (2020), and
$\mathcal{K}_{\ell_{1}\\!\operatorname{SOS}}$ and
$\mathcal{K}_{\ell_{2}\\!\operatorname{SOS}}$ in Kapelevich et al. (2021,
Equations 2.7 and 4.10).

Scalar WSOS cone.

$\mathcal{K}_{\operatorname{SOS}}$ is a cone of polynomials that are
guaranteed to be nonnegative pointwise on the domain.

Symmetric matrix WSOS cone.

$\mathcal{K}_{\operatorname{matSOS}}$ is a cone of polynomial symmetric
matrices (in an svec-like format) of side dimension $t$ that are guaranteed to
belong to $\mathbb{S}_{\succeq}$ pointwise on the domain. We let $m\coloneqq
st+d$ in Table 1 for succinctness.

$\ell_{1}$ epigraph WSOS cone.

$\mathcal{K}_{\ell_{1}\\!\operatorname{SOS}}$ is a cone of polynomial vectors
of length $1+t$ that are guaranteed to belong to
$\mathcal{K}_{\operatorname{\ell_{\infty}}}^{\ast}$ pointwise on the domain.

$\ell_{2}$ epigraph WSOS cone.

$\mathcal{K}_{\ell_{2}\\!\operatorname{SOS}}$ is a cone of polynomial vectors
of length $1+t$ that are guaranteed to belong to
$\mathcal{K}_{\operatorname{\ell_{2}}}$ pointwise on the domain.

For each cone, we have an analytic form for the feasibility check, gradient,
Hessian, and TOO oracles defined in Section 3. That is, we always avoid
iterative numerical procedures such as optimization, which are typically slow,
numerically unstable, and require tuning. Hypatia’s algorithm always evaluates
the feasibility check before the gradient, Hessian, and TOO (which are only
defined at strictly feasible points), and the gradient is evaluated before the
Hessian and TOO. For most of these cones, the feasibility check and gradient
oracles compute values and factorizations that are also useful for computing
the Hessian and TOO, so this data is cached in the cone data structures and
re-used where possible. In Table 1, we estimate the time complexities
(ignoring constants) of these four oracles for each cone, counting the cost of
cached values and factorizations only once (for the oracle that actually
computes them). Table 1 shows that the TOO is never more expensive than the
feasibility check, gradient, and Hessian oracles (i.e. the oracles needed by
SY). Indeed, our computational results in Section 7.3 demonstrate that the TOO
is very rarely an algorithmic bottleneck in practice.

cone | $\dim(\mathcal{K})$ | $\nu$ | feasibility | gradient | Hessian | TOO
---|---|---|---|---|---|---
$\mathcal{K}_{\geq}$ | $d$ | $d$ | $d$ | $d$ | $d$ | $d$
$\mathcal{K}_{\succeq}$ | $\operatorname{sd}(d)$ | $d$ | $d^{3}$ | $d^{3}$ | $d^{4}$ | $d^{3}$
$\mathcal{K}_{\operatorname{DNN}}$ | $\operatorname{sd}(d)$ | $\operatorname{sd}(d)$ | $d^{3}$ | $d^{3}$ | $d^{4}$ | $d^{3}$
$\mathcal{K}_{\operatorname{sPSD}}$ | $d$ | $s$ | $JD^{2}l$ | $JD^{2}l$ | $dJD^{2}l$ | $JD^{2}l$
$\mathcal{K}_{\operatorname{LMI}}$ | $d$ | $s$ | $ds^{2}+s^{3}$ | $ds^{3}$ | $d^{2}s^{2}$ | $ds^{2}+s^{3}$
$\mathcal{K}_{\operatorname{\ell_{\infty}}}$ | $1+d$ | $1+d$ | $d$ | $d$ | $d$ | $d$
$\mathcal{K}_{\operatorname{\ell_{2}}}$, $\mathcal{K}_{\operatorname{sqr}}$ | $1+d$ | $2$ | $d$ | $d$ | $d^{2}$ | $d$
$\mathcal{K}_{\operatorname{\ell_{\operatorname{spec}}}}$ | $1+rs$ | $1+r$ | $r^{2}s+r^{3}$ | $r^{2}s+r^{3}$ | $r^{2}s^{2}$ | $rs^{2}$
$\mathcal{K}_{\operatorname{matsqr}}$ | $\operatorname{sd}(r)+1+rs$ | $1+r$ | $r^{2}s+r^{3}$ | $r^{2}s+r^{3}$ | $r^{2}s^{2}$ | $rs^{2}$
$\mathcal{K}_{\operatorname{gpow}}$ | $r+s$ | $1+r$ | $r+s$ | $r+s$ | $r^{2}+s^{2}$ | $r+s$
$\mathcal{K}_{\operatorname{pow}}$, $\mathcal{K}_{\operatorname{geo}}$ | $1+d$ | $1+d$ | $d$ | $d$ | $d^{2}$ | $d$
$\mathcal{K}_{\operatorname{rtdet}}$ | $1+\operatorname{sd}(d)$ | $1+d$ | $d^{3}$ | $d^{3}$ | $d^{4}$ | $d^{3}$
$\mathcal{K}_{\log}$ | $2+d$ | $2+d$ | $d$ | $d$ | $d^{2}$ | $d$
$\mathcal{K}_{\operatorname{logdet}}$ | $2+\operatorname{sd}(d)$ | $2+d$ | $d^{3}$ | $d^{3}$ | $d^{4}$ | $d^{3}$
$\mathcal{K}_{\operatorname{sepspec}}\mathchar 45\relax\mathcal{K}_{\geq}$ | $2+d$ | $2+d$ | $d$ | $d$ | $d^{2}$ | $d$
$\mathcal{K}_{\operatorname{sepspec}}\mathchar 45\relax\mathcal{K}_{\succeq}$ | $2+\operatorname{sd}(d)$ | $2+d$ | $d^{3}$ | $d^{3}$ | $d^{5}$ | $d^{3}$
$\mathcal{K}_{\operatorname{relent}}$ | $1+2d$ | $1+2d$ | $d$ | $d$ | $d^{2}$ | $d$
$\mathcal{K}_{\operatorname{matrelent}}$ | $1+2\operatorname{sd}(d)$ | $1+2d$ | $d^{3}$ | $d^{3}$ | $d^{5}$ | $d^{4}$
$\mathcal{K}_{\operatorname{SOS}}$ | $d$ | $sr$ | $ds^{2}r$ | $ds^{2}r$ | $d^{2}sr$ | $ds^{2}r$
$\mathcal{K}_{\operatorname{matSOS}}$ | $d\operatorname{sd}(t)$ | $str$ | $ms^{2}t^{2}r$ | $ds^{2}t^{2}r$ | $d^{2}st^{3}r$ | $ms^{2}t^{2}r$
$\mathcal{K}_{\ell_{1}\\!\operatorname{SOS}}$ | $d(1+t)$ | $str$ | $ds^{2}tr$ | $ds^{2}tr$ | $d^{2}str$ | $ds^{2}tr$
$\mathcal{K}_{\ell_{2}\\!\operatorname{SOS}}$ | $d(1+t)$ | $2sr$ | $ds^{2}tr$ | $ds^{2}tr$ | $d^{2}st^{2}r$ | $ds^{2}t^{2}r$
Table 1:  Cone dimension $\dim(\mathcal{K})$, LHSCB parameter $\nu$, and time
complexity estimates (ignoring constants) for our feasibility check, gradient,
Hessian, and TOO implementations, for the exotic cones defined in Section 6.

Our TOO in 5 is distinct from the ‘higher order corrector’ terms proposed by
Mehrotra (1992); Dahl and Andersen (2021). The method by Mehrotra (1992) only
applies to symmetric cones, and Dahl and Andersen (2021) test their technique
only for the standard exponential cone. Compared to the third order term
proposed by Dahl and Andersen (2021), our TOO has a simpler and more symmetric
structure, as it relies on only one direction $\delta_{\bar{s}}$ rather than
two. Like the gradient and Hessian oracles, our TOO is additive for sums of
LHSCBs, which can be useful for cones (such as
$\mathcal{K}_{\operatorname{DNN}}$ and $\mathcal{K}_{\operatorname{SOS}}$)
that are defined as intersections of other cones. We leverage these properties
to obtain fast and numerically stable TOO implementations.

To illustrate, in Section C.1 we define LHSCBs and derive efficient TOO
procedures for a class of cones that can be characterized as intersections of
slices of the PSD cone $\mathcal{K}_{\succeq}$. We consider
$\mathcal{K}_{\operatorname{LMI}}$ in Section C.2 and
$\mathcal{K}_{\operatorname{SOS}}^{\ast}$ and
$\mathcal{K}_{\operatorname{matSOS}}^{\ast}$ in Section C.3. In Section C.4,
we handle $\mathcal{K}_{\operatorname{sPSD}}$ by differentiating a procedure
by Andersen et al. (2013) for computing Hessian products. In Section C.5 we
also show how to compute the TOO for $\mathcal{K}_{\operatorname{\ell_{2}}}$
and $\mathcal{K}_{\operatorname{sqr}}$. In Coey et al. (2021c), we derive
efficient TOO procedures for a class of spectral function cones on positive
domains ($\mathcal{K}_{\operatorname{sepspec}}$, $\mathcal{K}_{\log}$,
$\mathcal{K}_{\operatorname{logdet}}$, $\mathcal{K}_{\operatorname{geo}}$,
$\mathcal{K}_{\operatorname{rtdet}}$).

##  7 Computational testing

In Section 7.1, we introduce a diverse set of exotic conic benchmark instances
generated from a variety of applied examples. In Section 7.2, we describe our
methodology for comparing the stepping procedures from Section 5.5, and in
Section 7.3 we examine our computational results.

###  7.1 Exotic conic benchmark set

We generate 379 instances (in our primal general form 7) from 37 applied
examples in Hypatia’s examples folder. All instances are primal-dual feasible
except for 12 that are primal infeasible and one that is dual infeasible. For
most examples, we construct multiple formulations using different predefined
exotic cones from the list in Section 6. Each cone from this list appears in
at least one instance, so we consider our benchmark set to be the most diverse
collection of conic instances available.

We generate most instances using JuMP, but for some we use Hypatia’s native
model interface. Due to the size of some instances and the lack of a standard
instance storage format recognizing our cone types, we generate all instances
on the fly in Julia. For instances that use random data, we set random seeds
to ensure reproducibility. Figure 1 shows the distributions of instance
dimensions and exotic cone counts. All instances have at least one cone (note
any $\mathcal{K}_{\geq}$ cones are concatenated together, so
$\mathcal{K}_{\geq}$ is counted at most once) and take at least one iteration
to solve with Hypatia.

Below we briefly introduce each example. In Table 2, we summarize for each
example the number of corresponding instances and the cone types represented
in at least one of the instances. We do not distinguish dual cones and primal
cones in this summary (for example, instances that use
$\mathcal{K}_{\operatorname{\ell_{\infty}}}^{\ast}$ are only listed as using
$\mathcal{K}_{\operatorname{\ell_{\infty}}}$). For some examples, we describe
a subset of formulations in Coey et al. (2020); Kapelevich et al. (2021). Our
benchmark set includes ten instances from CBLIB (a conic benchmark instance
library, see Friberg (2016)). We chose to avoid running a larger sample of
instances from CBLIB so that the relatively few cone types supported by CBLIB
version 3 are not over-represented in our benchmark set.

Central polynomial matrix.

Minimize a spectral function of a gram matrix of a polynomial.

Classical-quantum capacity.

Compute the capacity of a classical-to-quantum channel (adapted from Fawzi and
Fawzi (2018, Section 3.1)).

Condition number.

Minimize the condition number of a matrix pencil subject to a linear matrix
inequality (adapted from Boyd et al. (1994, Section 3.2)).

Contraction analysis.

Find a contraction metric that guarantees global stability of a dynamical
system (adapted from Aylward et al. (2008, Section 5.3)). Six instances are
primal infeasible.

Convexity parameter.

Find the strong convexity parameter of a polynomial function over a domain.

Covariance estimation.

Estimate a covariance matrix that satisfies some given prior information and
minimizes a given convex spectral function.

Density estimation.

Find a valid polynomial density function maximizing the likelihood of a set of
observations (compare to (Papp and Alizadeh, 2014, Section 4.3), see Coey et
al. (2020, Section 5.6)).

Discrete maximum likelihood.

Maximize the likelihood of some observations at discrete points, subject to
the probability vector being close to a uniform prior.

D-optimal design.

Solve a D-optimal experiment design problem, i.e. maximize the determinant of
the information matrix subject to side constraints (adapted from Boyd and
Vandenberghe (2004, Section 7.5); see Coey et al. (2020, Section 5.4)).

Entanglement-assisted capacity.

Compute the entanglement-assisted classical capacity of a quantum channel
(adapted from Fawzi and Fawzi (2018, Section 3.2)).

Experiment design.

Solve a general experiment design problem that minimizes a given convex
spectral function of the information matrix subject to side constraints
(adapted from Boyd and Vandenberghe (2004, Section 7.5)).

Linear program.

Solve a simple linear program.

Lotka-Volterra.

Find an optimal controller for a Lotka-Volterra model of population dynamics
(adapted from Korda et al. (2016, Section 7.2)).

Lyapunov stability.

Minimize an upper bound on the root mean square gain of a dynamical system
(adapted from Boyd et al. (1994, Section 6.3.2) and Boyd (2009, Page 6)).

Matrix completion.

Complete a rectangular matrix by minimizing the nuclear norm and constraining
the missing entries (compare to Agrawal et al. (2019, Equation 8); see Coey et
al. (2020, Section 5.2)).

Matrix quadratic.

Find a rectangular matrix that minimizes a linear function and satisfies a
constraint on the outer product of the matrix.

Matrix regression.

Solve a multiple-output (or matrix) regression problem with regularization
terms, such as $\ell_{1}$, $\ell_{2}$, or nuclear norm (see Coey et al. (2020,
Section 5.3)).

Maximum volume hypercube.

Find a maximum volume hypercube (with edges parallel to the axes) inside a
given polyhedron or ellipsoid (adapted from MOSEK ApS (2020, Section 4.3.2)).

Nearest correlation matrix.

Compute the nearest correlation matrix in the quantum relative entropy sense
(adapted from Fawzi et al. (2019)).

Nearest polynomial matrix.

Given a symmetric matrix of polynomials $H$, find a polynomial matrix $Q$ that
minimizes the sum of the integrals of its elements over the unit box and
guarantees $Q-H$ is pointwise PSD on the unit box.

Nearest PSD matrix.

Find a sparse PSD matrix or a PSD-completable matrix (with a given sparsity
pattern) with constant trace that maximizes a linear function (adapted from
Sun and Vandenberghe (2015)).

Nonparametric distribution.

Given a random variable taking values in a finite set, compute the
distribution minimizing a given convex spectral function over all
distributions satisfying some prior information.

Norm cone polynomial.

Given a vector of polynomials, check a sufficient condition for pointwise
membership in $\mathcal{K}_{\operatorname{\ell_{2}}}$ or
$\mathcal{K}_{\operatorname{\ell_{\infty}}}^{\ast}$. Four instances are primal
infeasible.

Polynomial envelope.

Find a polynomial that closely approximates, over the unit box, the lower
envelope of a given list of polynomials (see Papp and Yildiz (2019, Section
7.2.1)).

Polynomial minimization.

Compute a lower bound for a given polynomial over a given semialgebraic set
(see Papp and Yildiz (2019, Section 7.3.1) and Coey et al. (2020, Section
5.5)). Some instances use polynomials with known optimal values from Burkardt
(2016).

Polynomial norm.

Find a polynomial that, over the unit box, has minimal integral and belongs
pointwise to the epigraph of the $\ell_{1}$ or $\ell_{2}$ norm of other given
polynomials (see Kapelevich et al. (2021, Section 6)).

Portfolio.

Maximize the expected returns of a stock portfolio and satisfy various risk
constraints (see Coey et al. (2020, Section 5.1)).

Region of attraction.

Find the region of attraction of a polynomial control system (see Henrion and
Korda (2013, Section 9.1)).

Relative entropy of entanglement.

Compute a lower bound on relative entropy of entanglement with a positive
partial transpose relaxation (adapted from Fawzi and Fawzi (2018, Section 4)).

Robust geometric programming.

Bound the worst-case optimal value of an uncertain signomial function with a
given coefficient uncertainty set (adapted from Chandrasekaran and Shah (2017,
Equation 39)).

Semidefinite polynomial matrix.

Check a sufficient condition for global convexity of a given polynomial. Two
instances are primal infeasible and one is dual infeasible.

Shape constrained regression.

Given a dataset, fit a polynomial function that satisfies shape constraints
such as monotonicity or convexity over a domain (see Coey et al. (2020,
Section 5.7)). Several instances use real datasets from Mazumder et al.
(2019).

Signomial minimization.

Compute a global lower bound for a given signomial function (see Murray et al.
(2020)). Several instances use signomials with known optimal values from
Murray et al. (2020); Chandrasekaran and Shah (2016).

Sparse LMI.

Optimize over a simple linear matrix inequality with sparse data.

Sparse principal components.

Solve a convex relaxation of the problem of approximating a symmetric matrix
by a rank-one matrix with a cardinality-constrained eigenvector (adapted from
d’Aspremont et al. (2007, Section 2)).

Stability number.

Given a graph, solve for a particular strengthening of the theta function
towards the stability number (adapted from Laurent and Piovesan (2015,
Equation 2.4)).

$10^{0}$$10^{1}$$10^{2}$$10^{3}$$10^{4}$$10^{5}$$0$$10$$20$$30$$40$$50$instance size ($n+p+q$)$10^{0}$$10^{1}$$10^{2}$$10^{3}$$0$$20$$40$$60$$80$$100$$120$exotic cone count ($K$) Figure 1:  Histograms summarizing the benchmark instances in the primal conic form 7. Instance size (log scale) is the sum of the primal variable, equality, and conic constraint dimensions. Exotic cone count (log scale) is the number of exotic cones comprising the Cartesian product cone.  example | # | cones in at least one instance
---|---|---
CBLIB | 10 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{\ell_{2}}}\;\mathcal{K}_{\operatorname{sqr}}\;\mathcal{K}_{\log}\;\mathcal{K}_{\operatorname{gpow}}$
central polynomial matrix | 24 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{sqr}}\;\mathcal{K}_{\operatorname{gpow}}\;\mathcal{K}_{\operatorname{rtdet}}\;\mathcal{K}_{\log}\;\mathcal{K}_{\operatorname{sepspec}}$
classical-quantum capacity | 9 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\succeq}\;\mathcal{K}_{\log}\;\mathcal{K}_{\operatorname{sepspec}}$
condition number | 6 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{LMI}}$
contraction analysis | 8 | $\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{matSOS}}$
convexity parameter | 7 | $\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{matSOS}}$
covariance estimation | 13 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{sqr}}\;\mathcal{K}_{\operatorname{gpow}}\;\mathcal{K}_{\operatorname{rtdet}}\;\mathcal{K}_{\log}\;\mathcal{K}_{\operatorname{sepspec}}$
density estimation | 16 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{sqr}}\;\mathcal{K}_{\operatorname{geo}}\;\mathcal{K}_{\log}\;\mathcal{K}_{\operatorname{SOS}}$
discrete maximum likelihood | 7 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\operatorname{pow}}\;\mathcal{K}_{\log}\;\mathcal{K}_{\operatorname{sepspec}}$
D-optimal design | 16 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{\ell_{\infty}}}\;\mathcal{K}_{\operatorname{\ell_{2}}}\;\mathcal{K}_{\operatorname{sqr}}\;\mathcal{K}_{\operatorname{geo}}\;\mathcal{K}_{\operatorname{rtdet}}\;\mathcal{K}_{\log}\;\mathcal{K}_{\operatorname{logdet}}$
entanglement-assisted capacity | 3 | $\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{sepspec}}\;\mathcal{K}_{\operatorname{matrelent}}$
experiment design | 13 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{sqr}}\;\mathcal{K}_{\operatorname{gpow}}\;\mathcal{K}_{\operatorname{rtdet}}\;\mathcal{K}_{\log}\;\mathcal{K}_{\operatorname{sepspec}}$
linear program | 3 | $\mathcal{K}_{\geq}$
Lotka-Volterra | 3 | $\mathcal{K}_{\succeq}$
Lyapunov stability | 10 | $\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{matsqr}}$
matrix completion | 11 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{sqr}}\;\mathcal{K}_{\operatorname{\ell_{\operatorname{spec}}}}\;\mathcal{K}_{\operatorname{gpow}}\;\mathcal{K}_{\operatorname{geo}}\;\mathcal{K}_{\log}$
matrix quadratic | 8 | $\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{matsqr}}$
matrix regression | 11 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{\ell_{\infty}}}\;\mathcal{K}_{\operatorname{\ell_{2}}}\;\mathcal{K}_{\operatorname{sqr}}\;\mathcal{K}_{\operatorname{\ell_{\operatorname{spec}}}}$
maximum volume hypercube | 15 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\operatorname{\ell_{\infty}}}\;\mathcal{K}_{\operatorname{\ell_{2}}}\;\mathcal{K}_{\operatorname{sqr}}\;\mathcal{K}_{\operatorname{geo}}$
nearest correlation matrix | 3 | $\mathcal{K}_{\operatorname{matrelent}}$
nearest polynomial matrix | 8 | $\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{SOS}}\;\mathcal{K}_{\operatorname{matSOS}}$
nearest PSD matrix | 28 | $\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{sPSD}}$
nonparametric distribution | 10 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\operatorname{sqr}}\;\mathcal{K}_{\operatorname{geo}}\;\mathcal{K}_{\log}\;\mathcal{K}_{\operatorname{sepspec}}$
norm cone polynomial | 10 | $\mathcal{K}_{\ell_{1}\\!\operatorname{SOS}}\;\mathcal{K}_{\ell_{2}\\!\operatorname{SOS}}$
polynomial envelope | 7 | $\mathcal{K}_{\operatorname{SOS}}$
polynomial minimization | 15 | $\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{SOS}}$
polynomial norm | 10 | $\mathcal{K}_{\operatorname{SOS}}\;\mathcal{K}_{\operatorname{matSOS}}\;\mathcal{K}_{\ell_{1}\\!\operatorname{SOS}}\;\mathcal{K}_{\ell_{2}\\!\operatorname{SOS}}$
portfolio | 9 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\operatorname{\ell_{\infty}}}\;\mathcal{K}_{\operatorname{\ell_{2}}}$
region of attraction | 6 | $\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{SOS}}$
relative entropy of entanglement | 6 | $\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{matrelent}}$
robust geometric programming | 6 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\operatorname{\ell_{\infty}}}\;\mathcal{K}_{\log}\;\mathcal{K}_{\operatorname{relent}}$
semidefinite polynomial matrix | 18 | $\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{\ell_{2}}}\;\mathcal{K}_{\operatorname{matSOS}}$
shape constrained regression | 11 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{\ell_{\infty}}}\;\mathcal{K}_{\operatorname{\ell_{2}}}\;\mathcal{K}_{\operatorname{SOS}}\;\mathcal{K}_{\operatorname{matSOS}}$
signomial minimization | 13 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\log}\;\mathcal{K}_{\operatorname{relent}}$
sparse LMI | 15 | $\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{sPSD}}\;\mathcal{K}_{\operatorname{LMI}}$
sparse principal components | 6 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{\ell_{\infty}}}$
stability number | 6 | $\mathcal{K}_{\geq}\;\mathcal{K}_{\succeq}\;\mathcal{K}_{\operatorname{DNN}}$
Table 2:  For each example, the count of instances and list of exotic cones
(defined in Section 6) used in at least one instance.

###  7.2 Methodology

We can assess the practical performance of a stepping procedure on a given
benchmark instance according to several metrics: whether the correct conic
certificate (satisfying our numerical tolerances, discussed below) is found,
and if so, the PDIPM iteration count and solve time. Across the benchmark set,
we compare performance between consecutive pairs of the five stepping
procedures outlined in Section 5.5.

basic.

The basic prediction or centering stepping procedure without any enhancements;
described in Section 5.5.1, this is similar to the method in Alfonso solver
(Papp and Yıldız, 2021), which is a practical implementation of the algorithm
by Skajaa and Ye (2015); Papp and Yıldız (2017).

prox.

The _basic_ procedure modified to use a less restrictive central path
proximity condition; described in Section 5.5.2.

TOA.

The _prox_ procedure with the TOA enhancement to incorporate third order LHSCB
information; described in Section 5.5.3.

curve.

The _TOA_ procedure adapted for a single backtracking search on a curve
instead of two backtracking line searches; described in Section 5.5.4.

comb.

The _curve_ procedure modified to search along a curve of combinations of both
the prediction and centering directions and their corresponding adjustment
directions; described in Section 5.5.5.

We perform all instance generation, computational experiments, and results
analysis using double precision floating point format, with Ubuntu 21.04,
Julia 1.7, and Hypatia 0.5.1 (with default options), on dedicated hardware
with an AMD Ryzen 9 3950X 16-core processor (32 threads) and 128GB of RAM. In
Appendix A, we outline the default procedures Hypatia uses for preprocessing,
initial point finding, and linear system solving for search directions. Simple
scripts and instructions for reproducing all results are available in
Hypatia’s benchmarks/stepper folder. The benchmark script runs all solves
twice and uses results from the second run, to exclude Julia compilation
overhead. A CSV file containing raw results is available at the Hypatia wiki
page.

When Hypatia converges for an instance, i.e. claims it has found a certificate
of optimality, primal infeasibility, or dual infeasibility, our scripts verify
that this is the correct type of certificate for that instance. For some
instances, our scripts also check additional conditions, for example that the
objective value of an optimality certificate approximately equals the known
true optimal value. We do not set restrictive time or iteration limits. All
failures to converge are caused by Hypatia ‘stalling’ during the stepping
iterations: either the backtracking search cannot step a distance of at least
the minimal value in the $\alpha$ schedule, or across several prediction steps
or combined directions steps, Hypatia fails to make sufficient progress
towards meeting the convergence conditions in Section 5.3.

Since some instances are more numerically challenging than others, we set the
termination tolerances (described in Section 5.3) separately for each
instance. Let $\epsilon\approx 2.22\times 10^{-16}$ be the machine epsilon.
For most instances, we use
$\varepsilon_{f}=\varepsilon_{r}=10\epsilon^{1/2}\approx 1.49\times 10^{-7}$
for the feasibility and relative gap tolerances,
$\varepsilon_{i}=\varepsilon_{a}=10\epsilon^{3/4}\approx 1.82\times 10^{-11}$
for the infeasibility and absolute gap tolerances, and
$\varepsilon_{p}=0.1\epsilon^{3/4}\approx 1.82\times 10^{-13}$ for the ill-
posedness tolerance. For 50 instances that are particularly numerically
challenging, we loosen all of these tolerances by a factor of either 10 or
100, and for two challenging primal infeasible instances of the _contraction
analysis_ example, we set $\varepsilon_{i}=10^{-9}$. This ensures that for
every benchmark instance, at least one of the five stepping procedures
converges.

Following Fleming and Wallace (1986), we define the _shifted geometric mean_
with shift $s\geq 0$, for $d$ values $v\in\mathbb{R}_{>}^{d}$, as:

$M(v,s)\coloneqq\textstyle\prod_{i\in\llbracket
d\rrbracket}(v_{i}+s)^{1/d}-s.$ (44)

We always apply a shift of one for iteration counts. Since different stepping
procedures converge on different subsets of instances, in tables we show three
types of shifted geometric means, each computed from a vector of values ($v$
in 44) obtained using one of the following approaches.

every.

Values for the 353 instances on which every stepping procedure converged.

this.

Values for instances on which this stepping procedure (corresponding to the
row of the table) converged.

all.

Values for all instances, but for any instances for which this stepping
procedure (corresponding to the row of the table) failed to converge, the
value is replaced with two times the maximum value for that instance across
the stepping procedures that converged.

The shifted geometric means for the _every_ approach are the most directly
comparable because they are computed on a fixed subset of instances, so we
usually quote the _every_ results in our discussion in Section 7.3.

Table 3 shows counts of converged instances and shifted geometric means of
iteration count and total solve time (in milliseconds), for the five stepping
procedures. We use a shift of one millisecond for the solve times in Table 3,
as some instances solve very quickly (see Figure 2).

Table 4 shows shifted geometric means of the time (in milliseconds) Hypatia
spends performing each of the following key algorithmic components, for the
five stepping procedures.

init.

Performed once during an entire solve run, independently of the stepping
iterations. Includes rescaling and preprocessing of model data, initial
interior point finding, and linear system solver setup (see Appendix A).

LHS.

Performed at the start of each iteration. Includes updating data that the
linear system solver (which has a fixed LHS in each iteration) uses to
efficiently compute at least one direction (such as updating and factorizing
the positive definite matrix in Appendix A).

RHS.

Performed between one and four times per iteration, depending on the stepping
procedure. Includes updating an RHS vector (see 28) for the linear system for
search directions. Note that the TOO is only evaluated while computing the
centering TOA RHS 32b and the prediction TOA RHS 42b.

direc.

Performed for each RHS vector. Includes solving the linear system for a search
direction (see 28) using the data computed during _LHS_ and a single RHS
vector computed during _RHS_ , and performing iterative refinement on the
direction (see Appendix A).

search.

Performed once or twice per iteration (occasionally more if the step length is
near zero), depending on the stepping procedure. Includes searching using
backtracking along a line or curve to find an interior point satisfying the
proximity conditions (see Appendix B).

For some instances that solve extremely quickly, these subtimings sum to only
around half of the total solve time due to extraneous overhead. However for
slower instances, these components account for almost the entire solve time.
In Table 4, _total_ is the time over all iterations, and _per iteration_ is
the average time per iteration (the arithmetic means are computed before the
shifted geometric mean). We use a shift of 0.1 milliseconds for the _init_ and
_total_ subtimings (left columns) and a shift of 0.01 milliseconds for the
_per iteration_ subtimings (right columns).

Finally, in Figures 3 and 6 we use _performance profiles_ (Dolan and Moré,
2002; Gould and Scott, 2016) to compare iteration counts and solve times
between pairs of stepping procedures. These should be interpreted as follows.
The _performance ratio_ for procedure $i$ and instance $j$ is the value
(iterations or solve time) attained by procedure $i$ on instance $j$ divided
by the better/smaller value attained by the two procedures on instance $j$.
Hence a performance ratio is at least one, and smaller values indicate better
relative performance. For a point $(x,y)$ on a performance profile curve for a
particular procedure, $x$ is the logarithm (base $2$) of performance ratio and
$y$ is the proportion of instances for which the procedure attains that
performance ratio or better/smaller. For example, a curve crosses the vertical
axis at the proportion of instances on which the corresponding procedure
performed at least as well as the alternative procedure. We use the Julia
package BenchmarkProfiles.jl (Orban, 2019) to compute coordinates for the
performance profile curves.

###  7.3 Results

Tables 3 and 6 demonstrate that each of the four cumulative stepping
enhancements tends to improve Hypatia’s iteration count and solve time. The
enhancements do not have a significant impact on the number of instances
Hypatia converges on. However, if we had enforced time or iteration limits,
the enhancements would have also improved the number of instances solved. This
is clear from Figure 2, which shows the distributions of iteration counts and
solve times for the _basic_ and _comb_ stepping procedures. We note that
Figure 4 (left) supports the intuition that formulation size is strongly
positively correlated with solve time for _comb_.

Overall, Table 3 shows that on the subset of instances solved by every
stepping procedure (_every_), the enhancements together reduce the shifted
geometric means of iterations and solve time by more than 80% and 70%
respectively (i.e. comparing _comb_ to _basic_). Figure 3 shows that the
iteration count and solve time improve on nearly every instance solved by both
_basic_ and _comb_ , and the horizontal axis scale shows that the magnitude of
these improvements is large on most instances. Figure 5 shows that for
instances that take more iterations or solve time, the enhancements tend to
yield a greater improvement in these measures. On every instance, the
enhancements improve the iteration count by at least 33%. The few instances
for which solve time regressed with the enhancements all solve relatively
quickly.

$10^{0}$$10^{1}$$10^{2}$$10^{3}$$0$$20$$40$$60$$80$$100$iteration countbasiccomb$10^{-4}$$10^{-2}$$10^{0}$$10^{2}$$10^{4}$$0$$10$$20$$30$$40$solve time (seconds)basiccomb Figure 2:  Overlayed histograms of iteration count (left, log scale) and solve time (right, log scale, in seconds) for the _basic_ and _comb_ stepping procedures, excluding instances that fail to converge.  |  | iterations | solve time
---|---|---|---
step | conv | every | this | all | every | this | all
basic | 371 | $101.34$ | $100.93$ | $102.39$ | $2130.98$ | $2207.10$ | $2282.19$
prox | 369 | $64.73$ | $65.28$ | $67.23$ | $1316.50$ | $1390.14$ | $1451.24$
TOA | 374 | $34.98$ | $35.31$ | $36.08$ | $1014.26$ | $1062.66$ | $1103.01$
curve | 372 | $29.67$ | $30.01$ | $30.99$ | $742.49$ | $780.95$ | $820.16$
comb | 367 | $18.31$ | $18.55$ | $20.02$ | $623.82$ | $655.71$ | $706.49$
Table 3:  For each stepping procedure, the number of converged instances and shifted geometric means of iterations and solve times (in milliseconds).  $0$$1$$2$$3$$4$$5$$0$$0.2$$0.4$$0.6$$0.8$$1$iteration countbasiccomb$0$$1$$2$$3$$4$$0$$0.2$$0.4$$0.6$$0.8$$1$solve timebasiccomb Figure 3:  Performance profiles (see Section 7.2) of iteration count (left) and solve time (right) for the four stepping enhancements overall.  $10^{1}$$10^{2}$$10^{3}$$10^{4}$$10^{5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$10^{1}$$10^{2}$$10^{3}$instance size ($n+p+q$)solve time (seconds)$0$$0.1$$0.2$$0.3$$0.4$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$10^{1}$$10^{2}$$10^{3}$proportion of time in _RHS_ Figure 4:  Solve time (log scale, in seconds) for the _comb_ stepping procedure against (left) instance size (log scale) and (right) the proportion of solve time spent in _RHS_ , excluding instances that fail to converge.  |  |  | total | per iteration
---|---|---|---|---
set | step | init | LHS | RHS | direc | search | LHS | RHS | direc | search
every | basic | $29.49$ | $741.38$ | $1.45$ | $75.59$ | $125.64$ | $7.73$ | $0.02$ | $0.81$ | $1.29$
prox | $29.45$ | $486.34$ | $1.11$ | $50.33$ | $67.74$ | $7.88$ | $0.02$ | $0.84$ | $1.10$
TOA | $29.42$ | $284.68$ | $10.96$ | $52.21$ | $73.33$ | $8.28$ | $0.32$ | $1.53$ | $2.14$
curve | $29.55$ | $244.38$ | $9.24$ | $44.60$ | $33.43$ | $8.33$ | $0.32$ | $1.53$ | $1.15$
comb | $29.31$ | $159.94$ | $10.51$ | $57.57$ | $35.23$ | $8.74$ | $0.58$ | $3.14$ | $1.94$
this | basic | $30.33$ | $784.33$ | $1.48$ | $78.82$ | $131.78$ | $8.20$ | $0.02$ | $0.85$ | $1.36$
prox | $30.07$ | $519.48$ | $1.12$ | $53.43$ | $72.60$ | $8.35$ | $0.02$ | $0.88$ | $1.16$
TOA | $30.16$ | $301.54$ | $11.97$ | $55.44$ | $78.34$ | $8.70$ | $0.35$ | $1.61$ | $2.26$
curve | $30.52$ | $260.67$ | $9.99$ | $47.30$ | $35.10$ | $8.80$ | $0.34$ | $1.60$ | $1.20$
comb | $30.46$ | $171.04$ | $11.03$ | $60.72$ | $36.40$ | $9.23$ | $0.60$ | $3.27$ | $1.98$
all | basic | $31.13$ | $814.44$ | $1.62$ | $82.65$ | $134.65$ | $8.52$ | $0.02$ | $0.91$ | $1.40$
prox | $31.29$ | $549.29$ | $1.25$ | $56.22$ | $75.09$ | $8.74$ | $0.02$ | $0.94$ | $1.20$
TOA | $31.16$ | $316.91$ | $12.23$ | $57.52$ | $79.71$ | $9.04$ | $0.36$ | $1.66$ | $2.28$
curve | $31.38$ | $275.95$ | $10.40$ | $49.63$ | $37.34$ | $9.17$ | $0.36$ | $1.68$ | $1.26$
comb | $31.36$ | $188.07$ | $11.88$ | $64.18$ | $40.02$ | $9.66$ | $0.63$ | $3.35$ | $2.10$
Table 4:  For each stepping procedure, the shifted geometric means of
subtimings (in milliseconds) for the key algorithmic components.
$0$$20$$40$$60$$80$$0.4$$0.6$$0.8$$1$iteration countrelative
improvement$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$10^{1}$$10^{2}$$10^{3}$$-1$$-0.5$$0$$0.5$$1$solve
time (seconds) Figure 5:  Relative improvement, from _basic_ to _comb_ , in
iteration count (left) or solve time (right) against iteration count or solve
time (in seconds) respectively for _comb_ , over the 356 instances on which
both _basic_ and _comb_ converge.
$0$$0.5$$1$$1.5$$2$$2.5$$0$$0.2$$0.4$$0.6$$0.8$$1$iteration
countbasicprox$0$$0.5$$1$$1.5$$2$$2.5$$0$$0.2$$0.4$$0.6$$0.8$$1$solve
timebasicprox$0$$0.5$$1$$1.5$$2$$0$$0.2$$0.4$$0.6$$0.8$$1$proxTOA$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$1.4$$0$$0.2$$0.4$$0.6$$0.8$$1$proxTOA$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$TOAcurve$0$$0.5$$1$$1.5$$2$$0$$0.2$$0.4$$0.6$$0.8$$1$TOAcurve$0$$0.5$$1$$1.5$$2$$0$$0.2$$0.4$$0.6$$0.8$$1$curvecomb$0$$0.5$$1$$1.5$$2$$0$$0.2$$0.4$$0.6$$0.8$$1$curvecomb
Figure 6:  Performance profiles (see Section 7.2) of iteration count (left
column) and solve time (right column) for the four stepping enhancements
(rows).

Each enhancement, by design, changes one modular component or aspect of the
stepping procedure. Below, we examine the impact of our algorithmic choices by
discussing pairwise comparisons of consecutive stepping procedures.

####  7.3.1 Less restrictive proximity

We compare _basic_ and _prox_ to evaluate the central path proximity
enhancement introduced in Section 5.5.2. Figure 6 (first row) shows that the
iteration count and solve time improve for nearly all instances. From Table 3,
the shifted geometric means of iteration count and solve time improve by over
35%.

The similarity between the iteration count and solve time performance profiles
in Figure 6 and also between the per iteration subtimings in Table 4 suggests
that the solve time improvement is driven mainly by the reduction in iteration
count. The per iteration _search_ time decreases slightly, since on average
fewer backtracking search steps are needed per iteration for _prox_ (because
it tends to step further in the prediction directions, as evidenced by the
smaller iteration counts). These results suggest that the central path
proximity restrictions in the algorithms by Skajaa and Ye (2015); Papp and
Yıldız (2021) are too conservative from the perspective of practical
performance, and that we need not restrict iterates to a very small
neighborhood of the central path in order to obtain high quality prediction
directions in practice.

####  7.3.2 Third order adjustments

We compare _prox_ and _TOA_ to evaluate the TOA enhancement introduced in
Section 5.5.3. Figure 6 (second row) shows that the iteration count improves
for all instances and by a fairly consistent magnitude, and the solve time
improves for nearly 80% of instances. From Table 3, the shifted geometric
means of iteration count and solve time improve by over 45% and over 20%
respectively.

Since _TOA_ computes an additional direction and performs an additional
backtracking search every iteration, the per iteration times for _direc_ and
_search_ in Table 4 nearly double. The _RHS_ time increases substantially,
because the TOO is evaluated for the second RHS vector (used to compute the
TOA direction), but _RHS_ is still much faster than the other components. Per
iteration, _direc_ and _search_ also remain fast compared to _LHS_. We see an
overall solve time improvement because the reduction in iteration count
usually outweighs the additional cost at each iteration. This suggests that
the TOO is generally relatively cheap to compute, and our TOA approach very
reliably improves the quality of the search directions.

####  7.3.3 Curve search

We compare _TOA_ and _curve_ to evaluate the curve search enhancement
introduced in Section 5.5.4. Figure 6 (third row) shows that the iteration
count and solve time improve for most instances, with larger and more
consistent improvements for the solve time. From Table 3, the shifted
geometric means of iteration count and solve time improve by over 15% and over
25% respectively.

Since _curve_ performs one backtracking search along a curve instead of the
two backtracking line searches needed by _TOA_ , the per iteration _search_
time in Table 4 nearly halves. The other subtimings are unaffected, so _curve_
improves the speed of each iteration. The improvement in iteration count may
stem from the more dynamic nature of the curve search compared to _TOA_ ’s
approach of computing a fixed combination of the unadjusted and TOA directions
as a function of the step distance in the unadjusted direction.

####  7.3.4 Combined directions

Finally, we compare _curve_ and _comb_ to evaluate the combined directions
enhancement introduced in Section 5.5.5. Figure 6 (fourth row) shows that the
iteration count and solve time improve on around 90% and 70% of instances
respectively. From Table 3, the shifted geometric means of iteration count and
solve time improve by nearly 40% and over 15% respectively.

Since _comb_ computes four directions per iteration (unadjusted and TOA
directions for both prediction and centering) instead of two, the per
iteration times for _RHS_ and _direc_ approximately double in Table 4. The
_search_ time increases because on average more backtracking curve search
steps are needed per iteration (for _curve_ , the centering phase typically
does not require multiple backtracking steps). Per iteration, _LHS_ remains
slower than the other components combined. Hence combining the prediction and
centering phases generally improves practical performance, and should be more
helpful when _LHS_ is particularly expensive (such as when $n-p$, the side
dimension of the PSD matrix we factorize during _LHS_ , is large; see Appendix
A). Furthermore, Figure 4 (right) shows that for most instances, _RHS_
accounts for a small proportion of the overall solve time for _comb_ ,
especially for instances that take longer to solve. This suggests that the TOO
is rarely a bottleneck for our _comb_ stepping procedure.

##  Appendix A Preprocessing and solving for search directions

We discuss preprocessing and initial point finding procedures and solving
structured linear systems for directions. Although Hypatia has various
alternative options for these procedures, we only describe the set of options
we fix in our computational experiments in Section 7, to give context for
these results. These techniques are likely to be useful for other conic PDIPM
implementations.

Given a conic model specified in the general primal conic form 7, we first
rescale the primal and dual equality constraints 7b and 8b to improve the
conditioning of the affine data. Next, we perform a QR factorization of
$A^{\prime}$ and check whether any primal equalities are inconsistent
(terminating if so). We use this factorization to modify $c,G,h$ and eliminate
all $p$ primal equalities (removing dual variable $y$), reducing the dimension
of the primal variable $x$ from $n$ to $n-p$. Next, we perform a QR
factorization of the modified $G$. We use this factorization to check whether
any dual equalities are inconsistent (terminating if so) and to remove any
redundant dual equalities, further reducing the dimension of $x$. This
factorization also allows us to cheaply compute an initial $x^{0}$ satisfying
17c. Since $y$ is eliminated, we do not need to solve 18b for $y^{0}$.

Starting from the initial interior point $\omega^{0}$ defined in Section 5.1,
we perform PDIPM iterations until the convergence conditions in Section 5.3
(in the preprocessed space) are met. Finally, we re-use the two QR
factorizations to lift the approximate certificate for the preprocessed model
to one for the original model. The residual norms for the lifted certificate
could violate the convergence tolerances, but we have not found such
violations to be significant on our benchmark instances.

During each PDIPM iteration, we solve the linear system 28 for a single LHS
matrix and between one and four RHS vectors, to obtain directions vectors
needed for one of the stepping procedures described in Section 5.5. Instead of
factorizing the large square nonsymmetric block-sparse LHS matrix, we utilize
its structure to reduce the size of the factorization needed. Some of these
techniques are adapted from methods in CVXOPT (see Vandenberghe (2010, Section
10.3)).

First we eliminate $s$ and $\kappa$, yielding a square nonsymmetric system,
then we eliminate $\tau$ to get a symmetric indefinite system in $x$ and $z$.
Most interior point solvers use a sparse LDL factorization (with precomputed
symbolic factorization) to solve this system. Although Hypatia can optionally
do the same, we see improved performance on our benchmark instances by further
reducing the system. After eliminating $z$, we have a (generally dense)
positive definite system, which we solve via a dense Cholesky factorization.
In terms of the original dimensions of the model before preprocessing
(assuming no redundant equalities), the side dimension of this system is
$n-p$. Finally, after finding a solution to 28, we apply several rounds of
iterative refinement in working precision to improve the solution quality.

We note that this Cholesky-based system solver method does not require
explicit Hessian oracles, only oracles for left-multiplication by the Hessian
or inverse Hessian. As we discuss in Appendix C and Coey et al. (2021c), these
optional oracles can be more efficient and numerically stable to compute for
many exotic cones. For cones without these oracles, Hypatia calls the explicit
Hessian matrix oracle, performing a Cholesky factorization of the Hessian if
necessary. A deeper discussion of Hypatia’s linear system solving techniques
and optional cone oracles is outside the scope of this paper.

##  Appendix B Efficient proximity checks

Recall that each stepping procedure in Section 5.5 uses at least one
backtracking search (on a line or a curve) to find a point $\omega$ satisfying
an aggregate proximity condition: $\pi_{\ell_{2}}(\omega)\leq\beta_{1}$ for
the _basic_ procedure in Section 5.5.1 or
$\pi_{\ell_{\infty}}(\omega)\leq\beta_{2}$ for the procedures in Sections
5.5.2 to 5.5.5. In Section 5.2, we define $\pi_{\ell_{2}}$ and
$\pi_{\ell_{\infty}}$ in 22 and 23. For each primitive cone
$k\in\llbracket\bar{K}\rrbracket$,
$0\leq\pi_{k}(\omega)\leq\pi_{\ell_{\infty}}(\omega)\leq\pi_{\ell_{2}}(\omega)$,
and by Lemma 5.1, $\pi_{k}(\omega)<1$ implies
$\bar{s}_{k}\in\operatorname{int}\bigl{(}\bar{\mathcal{K}}_{k}\bigr{)}$ and
$\bar{z}_{k}\in\operatorname{int}\bigl{(}\bar{\mathcal{K}}^{\ast}_{k}\bigr{)}$.
We use a schedule of decreasing trial values for the step parameter $\alpha$
and accept the first value that yields a candidate point satisfying the
aggregate proximity condition.

Suppose at a particular iteration of the backtracking search, we have the
candidate point $\omega$. We check a sequence of increasingly expensive
conditions that are necessary for the proximity condition to hold for
$\omega$. First, we verify that $\bar{s}_{k}^{\prime}\bar{z}_{k}>0,\forall
k\in\llbracket\bar{K}\rrbracket$, which is necessary for interiority (by a
strict version of the dual cone inequality 1). Note that this condition
implies $\mu(\omega)>0$. Next, we verify that $\rho_{k}(\omega)<\beta,\forall
k\in\llbracket\bar{K}\rrbracket$, where $\rho_{k}(\omega)$ is:

$\rho_{k}(\omega)\coloneqq\nu_{k}^{-1/2}\lvert\bar{s}_{k}^{\prime}\bar{z}_{k}/\mu-\nu_{k}\rvert\geq
0.$ (45)

In Lemma B.1 below, we show that $\rho_{k}(\omega)$ is a lower bound on
$\pi_{k}(\omega)$, so if $\rho_{k}(\omega)>\beta$ then then
$\pi_{k}(\omega)>\beta$. Computing $\rho_{k}$ is much cheaper than computing
$\pi_{k}(\omega)$ as it does not require evaluating any cone oracles.

Next, we iterate over $k\in\llbracket\bar{K}\rrbracket$ to check first the
primal feasibility oracle, then the optional dual feasibility oracle if
implemented, and finally the proximity condition $\pi_{k}(\omega)<\beta$.
Before computing $\pi_{k}(\omega)$, we check that the gradient and Hessian
oracle evaluations approximately satisfy two logarithmic homogeneity
conditions (Nesterov and Nemirovskii, 1994, Proposition 2.3.4):

$(g_{k}(\bar{s}_{k}))^{\prime}(H_{k}(\bar{s}_{k}))^{-1}g_{k}(\bar{s}_{k})=-\bar{s}_{k}^{\prime}g_{k}(\bar{s}_{k})=\nu_{k}.$
(46)

This allows us to reject $\omega$ if the cone oracles and the proximity value
$\pi_{k}(\omega)$ are likely to be numerically inaccurate.

######  Lemma B.1.

Given a point $\omega$ for which $\mu(\omega)>0$, for each
$k\in\llbracket\bar{K}\rrbracket$, $0\leq\rho_{k}(\omega)\leq\pi_{k}(\omega)$.

###### Proof.

We fix $\mu=\mu(\omega)>0$ for convenience. Let $f_{k}$ be the $\nu_{k}$-LHSCB
for $\bar{\mathcal{K}}_{k}$, and let the conjugate of $f_{k}$ be
$f^{\ast}_{k}$ (see 3), which is a $\nu_{k}$-LHSCB for
$\bar{\mathcal{K}}_{k}^{\ast}$. Let $g^{\ast}_{k}\coloneqq\nabla f^{\ast}_{k}$
and $H^{\ast}_{k}\coloneqq\nabla^{2}f^{\ast}_{k}$ denote the gradient and
Hessian operators for $f^{\ast}_{k}$. Using the logarithmic homogeneity
properties from Nesterov and Nemirovskii (1994, Proposition 2.3.4), and from
the definition of $\pi_{k}(\omega)$ in 20, we have:

$\displaystyle(\pi_{k}(\omega))^{2}$
$\displaystyle=(\bar{z}_{k}/\mu+g_{k}(\bar{s}_{k}))^{\prime}(H_{k}(\bar{s}_{k}))^{-1}(\bar{z}_{k}/\mu+g_{k}(\bar{s}_{k}))$
(47a)
$\displaystyle=\mu^{-2}\bar{z}_{k}^{\prime}(H_{k}(\bar{s}_{k}))^{-1}\bar{z}_{k}+2\mu^{-1}\bar{z}_{k}^{\prime}(H_{k}(\bar{s}_{k}))^{-1}g_{k}(\bar{s}_{k})+(g_{k}(\bar{s}_{k}))^{\prime}(H_{k}(\bar{s}_{k}))^{-1}g_{k}(\bar{s}_{k})$
(47b)
$\displaystyle=\mu^{-2}\bar{z}_{k}^{\prime}(H_{k}(\bar{s}_{k}))^{-1}\bar{z}_{k}-2\mu^{-1}\bar{z}_{k}^{\prime}\bar{s}_{k}+\nu_{k}.$
(47c)

By Papp and Yıldız (2017, Equation 13),
$(H_{k}(\bar{s}_{k}))^{-1}=H^{\ast}_{k}(-g_{k}(\bar{s}_{k}))$. Since
$f^{\ast}_{k}$ is a self-concordant barrier with parameter $\nu_{k}$, by
Nesterov et al. (2018, Equation 5.3.6) we have:
$(\bar{z}_{k}^{\prime}g^{\ast}_{k}(-g_{k}(\bar{s}_{k})))^{2}\leq\nu_{k}\bar{z}_{k}^{\prime}H^{\ast}_{k}(-g_{k}(\bar{s}_{k}))\bar{z}_{k}$.
Furthermore, $g^{\ast}_{k}(-g_{k}(\bar{s}_{k}))=\bar{s}_{k}$. Using these
facts, from 47 we have $\rho_{k}(\omega)\geq 0$ and:

$\displaystyle(\pi_{k}(\omega))^{2}$
$\displaystyle=\mu^{-2}\bar{z}_{k}^{\prime}H^{\ast}_{k}(-g_{k}(\bar{s}_{k}))\bar{z}_{k}-2\mu^{-1}\bar{z}_{k}^{\prime}\bar{s}_{k}+\nu_{k}$
(48a)
$\displaystyle\geq\nu_{k}^{-1}\mu^{-2}(\bar{z}_{k}^{\prime}\bar{s}_{k})^{2}-2\mu^{-1}\bar{z}_{k}^{\prime}\bar{s}_{k}+\nu_{k}$
(48b)
$\displaystyle=\nu_{k}^{-1}(\bar{s}_{k}^{\prime}\bar{z}_{k}/\mu-\nu_{k})$
(48c) $\displaystyle=(\rho_{k}(\omega))^{2}.$ (48d)

Therefore, $\pi_{k}(\omega)\geq\rho_{k}(\omega)\geq 0$ for all
$k\in\llbracket\bar{K}\rrbracket$. ∎

As an aside, we can use similar arguments to Lemma B.1 to show that
$\rho_{k}(\omega)$ also symmetrically bounds a conjugate proximity measure
$\pi_{k}^{\ast}(\omega)$, which we define as:

$\pi_{k}^{\ast}(\omega)\coloneqq\big{\lVert}(H_{k}^{\ast}(\bar{z}_{k}))^{-1/2}(\bar{s}_{k}/\mu+g_{k}^{\ast}(\bar{z}_{k}))\big{\rVert}\geq\nu_{k}^{-1/2}\lvert\bar{z}_{k}^{\prime}(\bar{s}_{k}/\mu+g_{k}^{\ast}(\bar{z}_{k}))\rvert=\rho_{k}(\omega).$
(49)

In general, we cannot check whether $\pi_{k}^{\ast}(\omega)<\beta$ because as
we discuss in Section 1.3 we do not have access to fast and numerically stable
conjugate barrier oracles ($g_{k}^{\ast}$ and $H_{k}^{\ast}$).

##  Appendix C Computing the TOO for some exotic cones

###  C.1 Intersections of slices of the PSD cone

First, we consider a proper cone $\mathcal{K}\subset\mathbb{R}^{q}$ that is an
inverse linear image (or slice) of the PSD cone
$\mathbb{S}^{\jmath}_{\succeq}$ of side dimension $\jmath$. Suppose:

$\mathcal{K}\coloneqq\\{s\in\mathbb{R}^{q}:\Lambda(s)\succeq 0\\},$ (50)

where $\Lambda:\mathbb{R}^{d}\to\mathbb{S}^{\jmath}$ is a linear operator,
with adjoint linear operator
$\Lambda^{\ast}:\mathbb{S}^{\jmath}\to\mathbb{R}^{d}$. Then the dual cone can
be characterized as:

$\mathcal{K}^{\ast}\coloneqq\\{s\in\mathbb{R}^{q}:\exists S\succeq
0,s=\Lambda^{\ast}(S)\\}.$ (51)

We note that for $\mathcal{K}_{\succeq}$ (the self-dual vectorized PSD cone),
we can let $q=d=\operatorname{sd}(\jmath)$,
$\Lambda(s)=\operatorname{mat}(s)$, and
$\Lambda^{\ast}(S)=\operatorname{vec}(S)$. Given a point $s\in\mathbb{R}^{q}$,
strict feasibility for $\mathcal{K}$ can be checked, for example, by
attempting a Cholesky factorization $\Lambda(s)=LL^{\prime}$, where $L$ is
lower triangular.

For $\mathcal{K}$ we have the LHSCB $f(s)=-\operatorname{logdet}(\Lambda(s))$
with parameter $\nu=\jmath$. Given a point
$s\in\operatorname{int}(\mathcal{K})$, we have
$\Lambda(s)\in\mathbb{S}^{\jmath}_{\succ}$ and its inverse
$\Lambda^{-1}(s)\in\mathbb{S}^{\jmath}_{\succ}$. For a direction
$\delta\in\mathbb{R}^{q}$, for $f$ at $s$ we can write the gradient, and the
Hessian and TOO applied to $\delta$, as (compare to Papp and Yildiz (2019,
Section 3)):

$\displaystyle g(s)$ $\displaystyle=-\Lambda^{\ast}(\Lambda^{-1}(s)),$ (52a)
$\displaystyle H(s)\delta$
$\displaystyle=\Lambda^{\ast}(\Lambda^{-1}(s)\Lambda(\delta)\Lambda^{-1}(s)),$
(52b) $\displaystyle\mathrm{T}(s,\delta)$
$\displaystyle=\Lambda^{\ast}(\Lambda^{-1}(s)\Lambda(\delta)\Lambda^{-1}(s)\Lambda(\delta)\Lambda^{-1}(s)).$
(52c)

If we have, for example, a Cholesky factorization $\Lambda(s)=LL^{\prime}$
(computed during the feasibility check), then the oracles in 52 are easy to
compute if $\Lambda$ and $\Lambda^{\ast}$ are easy to apply. We can compute
the TOO 52c using the following steps:

$\displaystyle Y$ $\displaystyle\coloneqq