Conic Optimization

4 January 2026 · Download PDF

1. Introduction

This note introduces some central ideas behind the implementation of efficient interior-point methods for solving conic programs. A large portion of the theory section is based on the papers for the conic solvers Clarabel¹¹ You can play around with a few examples using a Wasm-compiled Clarabel here. [1] and CVXOPT [2]. What this note adds is some essential implementation details that were left out of the original two papers (e.g. iterative refinement, matrix equilibration, and handling of the structure of the scaling matrices).

2. Primal and dual problems

We will in this text consider optimization problems where the objective function is at most quadratic and the constraints are conic (including the zero cone in order to model equality constraints). The general form of such (primal) problems is the following

\begin{aligned} minimize & \frac{1}{2} 𝑥^{⊤} 𝑃 𝑥 + 𝑞^{⊤} 𝑥 \\ subject to & 𝐴 𝑥 + 𝑠 = 𝑏 \\ 𝑠 \in 𝒦︀, \end{aligned}

(1)

where $𝑥 \in ℝ^{𝑛}$ is the decision variable, $𝑠 \in ℝ^{𝑚}$ is the slack variable, and $𝑃 = 𝑃^{⊤} ⪰ 0 \in ℝ^{𝑛 \times 𝑛}$ , $𝑞 \in ℝ^{𝑛}$ , $𝐴 \in ℝ^{𝑚 \times 𝑛}$ , $𝑏 \in ℝ^{𝑚}$ . Furthermore, $𝒦︀$ is a composite i.e. the Cartesian product of $𝐾$ cones ( $𝒦︀ = {𝒦︀}_{1} \times \dots \times {𝒦︀}_{𝐾}$ ) with each of the cones ${𝒦︀}_{𝑘}$ of dimension $𝑙_{𝑘}$ being of the following types

The zero cone: ${𝒵︀}_{𝑙_{𝑘}} ≔ {0}^{𝑙_{𝑘}}$ (used to model equality constraints).
The nonnegative cone: $ℝ_{+}^{𝑙_{𝑘}} ≔ {𝑠 \in ℝ^{𝑙_{𝑘}} | 𝑠 \geq 0}$ .
The second-order cone: ${𝒬︀}_{𝑙_{𝑘}} ≔ {(𝑠_{1}, 𝑠_{2}) \in ℝ \times ℝ^{𝑙_{𝑘} - 1} | {‖ 𝑠_{2} ‖}_{2} \leq 𝑠_{1}}$ .
The exponential cone: ${𝒦︀}_{exp} ≔ {(𝑠_{1}, 𝑠_{2}, 𝑠_{3}) \in ℝ^{3} | 𝑠_{2} exp (\frac{𝑠_{1}}{𝑠_{2}}) \leq 𝑠_{3}, 𝑠_{2} > 0} \cup {(𝑥, 0, 𝑧) | 𝑥 \leq 0, 𝑧 \geq 0} .$
The power cone: ${𝒦︀}_{pow} ≔ {(𝑠_{1}, 𝑠_{2}, 𝑠_{3}) \in ℝ^{3} | 𝑠_{1}^{𝛼} 𝑠_{2}^{1 - 𝛼} \geq | 𝑠_{3} |, 𝑠_{1}, 𝑠_{2} \geq 0}$ with $𝛼 \in (0, 1)$ .
The positive semidefinite cone: ${𝒮︀}_{𝑙_{𝑘}} ≔ {𝑠 \in ℝ^{\frac{𝑙_{𝑘} (𝑙_{𝑘} + 1)}{2}} | smat (𝑠) ⪰ 0}$ .

where $smat (𝑥)$ creates a symmetric matrix from the vector $𝑥$ . This can be done in a plethora of ways, but we will use the definition from Clarabel, which fills the lower triangle row by row

smat (𝑥) = [\begin{matrix} 𝑥_{1} & 𝑥_{2} / \sqrt{2} & \dots & 𝑥_{𝑛 (𝑛 - 1) / 2 + 1} / \sqrt{2} \\ 𝑥_{2} / \sqrt{2} & 𝑥_{3} & \dots & 𝑥_{𝑛 (𝑛 - 1) / 2 + 2} / \sqrt{2} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 𝑥_{𝑛 (𝑛 - 1) / 2 + 1} / \sqrt{2} & 𝑥_{𝑛 (𝑛 - 1) / 2 + 2} / \sqrt{2} & \dots & 𝑥_{𝑛 (𝑛 + 1) / 2} \end{matrix}] .

(2)

where the scaling of $\sqrt{2}$ have been introduced on the off-diagonal elements in order to preserve the inner product, i.e. $𝑥^{⊤} 𝑦 = tr (smat (𝑥) smat (𝑦))$ . Furthermore, we denote the inverse operation by $svec (𝑋)$ , i.e. an operation that creates a vector from a symmetric matrix $𝑋$ . The nonnegative, second-order, and positive semidefinite cones are examples of symmetric cones which are a special class of cones that have a lot of nice properties. The exponential and power cones are examples of nonsymmetric cones which do not have the same nice properties as symmetric cones, but are still very useful in practice. The zero cone is a special case that is used to model equality constraints.

2.1. The dual problem

The primal optimization problem of interest can equivalently be written as

Primal Problem : {\begin{cases} minimize & \frac{1}{2} 𝑥^{⊤} 𝑃 𝑥 + 𝑞^{⊤} 𝑥, \\ subject to & 𝑠 = 𝑏 - 𝐴 𝑥 \in 𝒦︀ . \end{cases}

(3)

where the conic constraint $𝑏 - 𝐴 𝑥 \in 𝒦︀$ can be understood as a generalized inequality, i.e.

0 ⪯_{𝒦︀} 𝑠 = 𝑏 - 𝐴 𝑥 .

(4)

However, the standard form (similar to LPs) requires that the inequality is turned the other way. As such we negate the expression, resulting in

𝐴 𝑥 - 𝑏 ⪯_{𝒦︀} 0 .

(5)

We now write up the Lagrangian as

𝐿 (𝑥, 𝑧) = \frac{1}{2} 𝑥^{⊤} 𝑃 𝑥 + 𝑞^{⊤} 𝑥 + 𝑧^{⊤} (𝐴 𝑥 - 𝑏)

(6)

Since $- (𝐴 𝑥 - 𝑏) \in 𝒦︀$ we must have that $𝑧$ is in the dual cone ${𝒦︀}^{*}$ in order for the Lagrangian to be a lower bound for any feasible $𝑥$ . This comes directly from the definition of the dual cone, i.e.

{𝒦︀}^{*} = {𝑧 \in ℝ^{𝑚} | 𝑧^{⊤} 𝑠 \geq 0, \forall 𝑠 \in 𝒦︀},

(7)

so that, in particular, $𝑧^{⊤} (𝐴 𝑥 - 𝑏) \leq 0$ for every $𝑥$ with $𝑏 - 𝐴 𝑥 \in 𝒦︀$ .

Now since we want to find a stationary point we compute the gradient of the Lagrangian and set it to zero

\nabla_{𝑥} 𝐿 (𝑥, 𝑧) = 𝑃 𝑥 + 𝑞 + 𝐴^{⊤} 𝑧 = 0 .

(8)

Now notice that the Lagrangian can be written as

𝑥^{⊤} (\frac{1}{2} 𝑃 𝑥 + 𝑞 + 𝐴^{⊤} 𝑧) - 𝑏^{⊤} 𝑧

(9)

The thing inside the parentheses, is almost equal to the zero constraint on the gradient. By performing the mathematician’s favorite magic trick of adding zero followed by rearranging the expressions we end up with

\begin{aligned} 𝑥^{⊤} (\frac{1}{2} 𝑃 𝑥 + 𝑞 + 𝐴^{⊤} 𝑧) + (\underset{= 0}{\underset{⏟}{\frac{1}{2} 𝑥^{⊤} 𝑃 𝑥 - \frac{1}{2} 𝑥^{⊤} 𝑃 𝑥}}) - 𝑏^{⊤} 𝑧 \\ = & 𝑥^{⊤} \underset{\nabla_{𝑥} 𝐿 (𝑥, 𝑧) = 0}{\underset{⏟}{(𝑃 𝑥 + 𝑞 + 𝐴^{⊤} 𝑧)}} - 𝑏^{⊤} 𝑧 - \frac{1}{2} 𝑥^{⊤} 𝑃 𝑥 . \end{aligned}

(10)

As such the dual formulation of the problem becomes

Dual Problem: {\begin{cases} maximize & - \frac{1}{2} 𝑥^{⊤} 𝑃 𝑥 - 𝑏^{⊤} 𝑧, \\ subject to & 𝑃 𝑥 + 𝑞 + 𝐴^{⊤} 𝑧 = 0, \\ 𝑧 \in {𝒦︀}^{*} \end{cases}

(11)

2.2. Homogeneous embedding

The standard KKT conditions for the primal problem are

\begin{aligned} 𝐴 𝑥 + 𝑠 & = 𝑏 \\ 𝑃 𝑥 + 𝐴^{⊤} 𝑧 & = - 𝑞 \\ 𝑠^{⊤} 𝑧 & = 0 \\ (𝑠, 𝑧) & \in 𝒦︀ \times {𝒦︀}^{*} . \end{aligned}

(12)

However, we can take a different path in order to solve the same problem. First we define the duality gap $𝛾$ as the difference between the primal and dual objectives

\begin{aligned} 𝛾 : & = (\frac{1}{2} 𝑥^{⊤} 𝑃 𝑥 + 𝑞^{⊤} 𝑥) - (- \frac{1}{2} 𝑥^{⊤} 𝑃 𝑥 - 𝑏^{⊤} 𝑧) \\ = 𝑥^{⊤} 𝑃 𝑥 + 𝑞^{⊤} 𝑥 + 𝑏^{⊤} 𝑧 = 𝑠^{⊤} 𝑧 \end{aligned}

(13)

where we used that $𝑞 = - 𝑃 𝑥 - 𝐴^{⊤} 𝑧$ (from the gradient constraint) and $𝑠 = 𝑏 - 𝐴 𝑥$ (from the definition of the equality constraints). Because $(𝑠, 𝑧) \in (𝒦︀, {𝒦︀}^{*})$ the duality gap is always nonnegative for any feasible point, and is zero at the optimal point $(𝑥^{*}, 𝑠^{*}, 𝑧^{*})$ when strong duality holds. We can therefore come up with a different optimization problem by combining the constraints of the primal and dual problem and with an objective to minimize the duality gap

\begin{aligned} \underset{𝑥, 𝑠, 𝑧}{minimize} & 𝑠^{⊤} 𝑧, \\ subject to & 𝑃 𝑥 + 𝐴^{⊤} 𝑧 + 𝑞 = 0, \\ 𝐴 𝑥 + 𝑠 - 𝑏 = 0, \\ 𝑥^{⊤} 𝑃 𝑥 + 𝑞^{⊤} 𝑥 + 𝑏^{⊤} 𝑧 = 0, \\ (𝑠, 𝑧) \in 𝒦︀ \times {𝒦︀}^{*} . \end{aligned}

(14)

An issue with the above is that the duality gap is explicitly set to zero in the third equality constraint. This condition can be relaxed after a change of variables $𝑥 \to \frac{𝑥}{𝜏}, 𝑧 \to \frac{𝑧}{𝜏},$ and $𝑠 \to \frac{𝑠}{𝜏}$ and instead solve

Homogeneous embedding: {\begin{cases} \underset{𝑥, 𝑠, 𝑧, 𝜏, 𝜅}{minimize} & 𝑠^{⊤} 𝑧 + 𝜏 𝜅 \\ subject to & 𝑃 𝑥 + 𝐴^{⊤} 𝑧 + 𝑞 𝜏 = 0, \\ 𝐴 𝑥 + 𝑠 - 𝑏 𝜏 = 0, \\ 𝜏^{- 1} 𝑥^{⊤} 𝑃 𝑥 + 𝑞^{⊤} 𝑥 + 𝑏^{⊤} 𝑧 + 𝜅 = 0, \\ (𝑠, 𝑧, 𝜏, 𝜅) \in 𝒦︀ \times {𝒦︀}^{*} \times ℝ_{+} \times ℝ_{+} . \end{cases}

(15)

In [1] they go on to show the existence of solutions and infeasibility detection. This will be skipped here.

3. Barrier functions

Barrier functions are functions that are defined only inside of the feasible region (i.e. when $𝑠 = 𝑏 - 𝐴 𝑥 \in 𝒦︀$ ) and tend towards infinity when we move towards the boundary of the feasible set. They are useful as they can be used to approximate the constrained problem by an unconstrained problem

{\begin{cases} minimize & \frac{1}{2} 𝑥^{⊤} 𝑃 𝑥 + 𝑞^{⊤} 𝑥 \\ subject to & 𝐴 𝑥 + 𝑠 = 𝑏 \\ 𝑠 \in 𝒦︀, \end{cases} \approx minimize \frac{1}{2} 𝑥^{⊤} 𝑃 𝑥 + 𝑞^{⊤} 𝑥 + 𝜇 𝑓 (𝑏 - 𝐴 𝑥) .

(16)

As $𝜇 \to 0$ the solution to the unconstrained problem tends towards the solution of the constrained problem. A naive approach would simply be to pick $𝜇$ to be a small value in order to get a good approximation of the original problem. However, this is not a good idea as the resulting optimization problem is very ill-conditioned and hard to solve. Instead, interior-point methods start with a large value of $𝜇$ and then track the solution of the unconstrained problem as $𝜇$ decreases until we are reasonably close to the optimal solution of the original problem. Methods that rely on this strategy of solving a series of simpler problems are known as homotopy methods. As an example we will now describe the properties of the logarithmic barrier functions for the three symmetric cones²² That is the nonnegative, second-order, and positive semidefinite cone., but later we will also discuss the barrier functions for the exponential and power cones.

3.1. The logarithmic barrier function for the three symmetric cones

For a composite cone $𝒦︀ = {𝒦︀}_{1} \times {𝒦︀}_{2} \times \dots \times {𝒦︀}_{𝐾}$ where each component is one of the three symmetric cones we have that the logarithmic barrier function is given by

𝑓 (𝑠) = \sum_{𝑘 = 1}^{𝐾} 𝑓_{𝑘} (𝑠_{𝑘}), 𝑓_{𝑘} (𝑠) = {\begin{cases} - \sum_{𝑗 = 1}^{𝑛} log (𝑠_{𝑗}) & {𝒦︀}_{𝑘} = ℝ_{+}^{𝑙_{𝑘}} \\ - \frac{1}{2} log (𝑠^{⊤} 𝐽 𝑠) & {𝒦︀}_{𝑘} = {𝒬︀}_{𝑙_{𝑘}} \\ - log det smat (𝑠) & {𝒦︀}_{𝑘} = {𝒮︀}_{𝑙_{𝑘}} \end{cases},

(17)

where $𝑠 = {[\begin{matrix} 𝑠_{1} & 𝑠_{2} & \dots & 𝑠_{𝐾} \end{matrix}]}^{⊤}$ with $𝑠_{𝑘} \in {𝒦︀}_{𝑘}$ and

𝐽 = [\begin{matrix} 1 & 0 \\ 0 & - 𝐼_{𝑙_{𝑘} - 1} \end{matrix}],

(18)

meaning that $𝑠_{𝑘}^{⊤} 𝐽 𝑠_{𝑘} = 𝑠_{1}^{2} - {‖ 𝑠_{2} ‖}^{2} \geq 0$ ³³ The notation here is a little sloppy as we use that $𝑠_{𝑘} = (𝑠_{1}, 𝑠_{2}) \in ℝ \times ℝ^{𝑙_{𝑘} - 1}$ . for all $𝑠_{𝑘} \in {𝒬︀}_{𝑙_{𝑘}}$ . Furthermore, the logarithmic barrier function satisfies the relation $𝑓 (𝑡 𝑠) = 𝑓 (𝑠) - 𝑚 log (𝑡)$ for $𝑡 > 0$ where

𝑚 = 𝑚_{1} + \dots + 𝑚_{𝐾}, 𝑚_{𝑘} = {\begin{cases} 𝑙_{𝑘} & {𝒦︀}_{𝑘} = ℝ_{+}^{𝑙_{𝑘}} \\ 1 & {𝒦︀}_{𝑘} = {𝒬︀}_{𝑙_{𝑘}} \\ 𝑙_{𝑘} & {𝒦︀}_{𝑘} = {𝒮︀}_{𝑙_{𝑘}} \end{cases} .

(19)

We refer to $𝑚_{𝑘}$ as the degree of the $𝑘$ th cone, so that $𝑚$ is the degree of the composite cone $𝒦︀$ .

3.2. Gradients of the barrier function for the three symmetric cones

By linearity the gradient of the logarithmic barrier function is given by the sum of gradients of the individual barrier functions for each symmetric cone given by

\nabla 𝑓_{𝑘} (𝑠_{𝑘}) = {\begin{cases} - {diag (𝑠_{𝑘})}^{- 1} 1_{𝑙_{𝑘}} & {𝒦︀}_{𝑘} = ℝ_{+}^{𝑙_{𝑘}} \\ - {(𝑠_{𝑘}^{⊤} 𝐽 𝑠_{𝑘})}^{- 1} 𝐽 𝑠_{𝑘} & {𝒦︀}_{𝑘} = {𝒬︀}_{𝑙_{𝑘}} \\ - svec ({smat (𝑠_{𝑘})}^{- 1}) & {𝒦︀}_{𝑘} = {𝒮︀}_{𝑙_{𝑘}} \end{cases},

(20)

where $1_{𝑙}$ is a $𝑙$ -vector of ones. It can be shown that $- \nabla 𝑓 (𝑠) \in 𝒦︀$ and that $𝑠^{⊤} \nabla 𝑓 (𝑠) = - 𝑚$ for $𝑠 \in 𝒦︀$ .

3.3. Hessians of the barrier function for the three symmetric cones

We denote the Hessian of $𝑓$ and $𝑓_{𝑘}$ as $𝐻 (𝑠) = \nabla^{2} 𝑓 (𝑠)$ and $𝐻_{𝑘 (𝑠_{𝑘})} = \nabla^{2} 𝑓_{𝑘 (𝑠_{𝑘})}$ :

𝐻_{𝑘} (𝑠_{𝑘}) = {\begin{cases} {diag (𝑠_{𝑘})}^{- 2} & {𝒦︀}_{𝑘} = ℝ_{+}^{𝑙} \\ {(𝑠_{𝑘}^{⊤} 𝐽 𝑠_{𝑘})}^{- 2} (2 𝐽 𝑠_{𝑘} 𝑠_{𝑘}^{⊤} 𝐽 - (𝑠_{𝑘}^{⊤} 𝐽 𝑠_{𝑘}) 𝐽) & {𝒦︀}_{𝑘} = {𝒬︀}_{𝑙} \\ {smat (𝑠_{𝑘})}^{- 1} \otimes_{𝑠} {smat (𝑠_{𝑘})}^{- 1} & {𝒦︀}_{𝑘} = {𝒮︀}_{𝑙} \end{cases},

(21)

where $\otimes_{𝑠}$ is the symmetric Kronecker product, defined in the section on semidefinite cones at the end of the post. For the scaling matrices below, we also need $𝑅_{𝑘} = 𝐿_{1} 𝑉 Λ_{𝑘}^{- 1 / 2} = 𝐿_{2}^{- ⊤} 𝑈 Λ_{𝑘}^{1 / 2}$ , which depends on both $𝑠_{𝑘}$ and $𝑧_{𝑘}$ , with

smat (𝑠_{𝑘}) = 𝐿_{1} 𝐿_{1}^{⊤}, smat (𝑧_{𝑘}) = 𝐿_{2} 𝐿_{2}^{⊤}, 𝐿_{2}^{⊤} 𝐿_{1} = 𝑈 Λ_{𝑘} 𝑉^{⊤} .

(22)

In addition, we have that

𝐻_{𝑘} {(𝑠_{𝑘})}^{- 1} = {\begin{cases} {diag (𝑠_{𝑘})}^{2} & {𝒦︀}_{𝑘} = ℝ_{+}^{𝑙} \\ 2 𝑠_{𝑘} 𝑠_{𝑘}^{⊤} - (𝑠_{𝑘}^{⊤} 𝐽 𝑠_{𝑘}) 𝐽 & {𝒦︀}_{𝑘} = {𝒬︀}_{𝑙} \\ smat (𝑠_{𝑘}) \otimes_{𝑠} smat (𝑠_{𝑘}) & {𝒦︀}_{𝑘} = {𝒮︀}_{𝑙} \end{cases} .

(23)

Further properties of the Hessian of the logarithmic barrier function are

\begin{aligned} {𝐻 (𝑠)}^{- 1} = & 𝐻 (𝑠^{- 1}), {𝐻 (𝑠)}^{1 / 2} = 𝐻 (𝑠^{1 / 2}) \\ 𝐻 (𝑠) 𝑠 = 𝑠^{- 1}, & 𝐻 (𝑠^{1 / 2}) 𝑠 = 𝑒, {(𝐻 (𝑠) 𝑣)}^{- 1} = {𝐻 (𝑠)}^{- 1} 𝑣^{- 1} . \end{aligned}

(24)

3.4. The Jordan product for the three symmetric cones

The Jordan products for the three symmetric cones are given by

𝑠 ⚬ 𝑧 = {\begin{cases} {[\begin{matrix} 𝑠_{1} 𝑧_{1} & 𝑠_{2} 𝑧_{2} & \dots & 𝑠_{𝑙} 𝑧_{𝑙} \end{matrix}]}^{⊤} & {𝒦︀}_{𝑘} = ℝ_{+}^{𝑙_{𝑘}} \\ {[\begin{matrix} 𝑠^{⊤} 𝑧 & 𝑠_{1} 𝑧_{2}^{⊤} + 𝑠_{2}^{⊤} 𝑧_{1} \end{matrix}]}^{⊤} & {𝒦︀}_{𝑘} = {𝒬︀}_{𝑙_{𝑘}} \\ \frac{1}{2} svec (smat (𝑠) smat (𝑧) + smat (𝑧) smat (𝑠)) & {𝒦︀}_{𝑘} = {𝒮︀}_{𝑙_{𝑘}} \end{cases},

(25)

where $𝑒$ is the idempotent for the cones given by

𝑒_{𝑘} = {\begin{cases} {[\begin{matrix} 1 & 1 & \dots & 1 \end{matrix}]}^{⊤} & {𝒦︀}_{𝑘} = ℝ_{+}^{𝑙_{𝑘}} \\ {[\begin{matrix} 1 & 0 & \dots & 0 \end{matrix}]}^{⊤} & {𝒦︀}_{𝑘} = {𝒬︀}_{𝑙_{𝑘}} \\ svec (𝐼_{𝑙_{𝑘}}) & {𝒦︀}_{𝑘} = {𝒮︀}_{𝑙_{𝑘}} \end{cases} .

(26)

Note that $𝑒^{⊤} (𝑧 ⚬ 𝑠) = 𝑧^{⊤} 𝑠$ and $𝑒^{⊤} 𝑒 = 𝑚$ . In addition, the inverse, square, and square root are defined by the relations $𝑠^{- 1} ⚬ 𝑠 = 𝑒$ , $𝑠^{2} = 𝑠 ⚬ 𝑠$ , and $𝑠^{\frac{1}{2}} ⚬ 𝑠^{\frac{1}{2}} = 𝑠$ .

3.5. Logarithmically homogeneous self-concordant barrier

The ideas of the logarithmic barrier functions can be extended to the nonsymmetric cones as well. We often refer to the general class of logarithmically homogeneous self-concordant barrier functions, which are a special class of barrier functions that have nice properties that can be exploited in the design of interior-point methods as follows: A function $𝑓 : 𝒦︀ \to ℝ$ is called a logarithmically homogeneous self-concordant barrier with degree $𝑚 \geq 1$ ( $𝑚$ -LHSCB)⁴⁴ The logarithmic barrier functions from the previous section are all $𝑚$ -LHSCB functions. for the convex cone $𝒦︀$ if it satisfies the following properties

\begin{aligned} | \nabla^{3} 𝑓 (𝑥) [𝑟, 𝑟, 𝑟] | \leq 2 {(\nabla^{2} 𝑓 (𝑥) [𝑟, 𝑟])}^{3 / 2} & \forall 𝑥 \in int (𝒦︀), 𝑟 \in ℝ^{𝑑}, \\ 𝑓 (𝑡 𝑥) = 𝑓 (𝑥) - 𝑚 log (𝑡), & \forall 𝑥 \in 𝒦︀, 𝑡 > 0 . \end{aligned}

(27)

The conjugate barrier function $𝑓_{*} : {𝒦︀}^{*} \to ℝ$ is defined as ⁵⁵ This is not to be confused with the convex conjugate of a function, which is denoted by $𝑓^{*} (𝑦) = sup {𝑦^{⊤} 𝑥 - 𝑓 (𝑥), 𝑥 \in int (𝒦︀)}$ . The relation between the two conjugates are $𝑓_{*} (𝑦) = 𝑓^{*} (- 𝑦)$ and $\nabla 𝑓_{*} (𝑦) = - \nabla 𝑓^{*} (- 𝑦)$ .

𝑓_{*} (𝑦) : = sup_{𝑥 \in int (𝒦︀)} {- 𝑦^{⊤} 𝑥 - 𝑓 (𝑥)} .

(28)

The function $𝑓_{*}$ is $𝑚$ -LHSCB for ${𝒦︀}^{*}$ if the associated $𝑓$ is. The primal gradient $\nabla 𝑓$ satisfies

𝑥^{⊤} \nabla 𝑓 (𝑥) = - 𝑚, \forall 𝑥 \in int (𝒦︀) .

(29)

In cases where the conjugate barrier function is known only through the definition (i.e. rather than a closed-form representation) we can compute its gradient as the solution to

\nabla 𝑓_{*} (𝑦) : = - arg sup_{𝑥 \in int (𝒦︀)} {- 𝑦^{⊤} 𝑥 - 𝑓 (𝑥)} .

(30)

4. The central path

As described previously, the idea of interior-point methods is to track the solution of an approximate unconstrained problem as $𝜇$ decreases until we are reasonably close to the optimal solution of the original constrained problem. The central path denotes this path. In order to describe the central path mathematically we define the function $𝐺$ representing the equality constraints of the homogeneous embedding

𝐺 (𝑥, 𝑧, 𝑠, 𝜏, 𝜅) : = [\begin{matrix} 0 \\ 𝑠 \\ 𝜅 \end{matrix}] - [\begin{matrix} 𝑃 & 𝐴^{⊤} & 𝑞 \\ - 𝐴 & 0 & 𝑏 \\ - 𝑞^{⊤} & - 𝑏^{⊤} & 0 \end{matrix}] [\begin{matrix} 𝑥 \\ 𝑧 \\ 𝜏 \end{matrix}] + [\begin{matrix} 0 \\ 0 \\ 𝜏^{- 1} 𝑥^{⊤} 𝑃 𝑥 \end{matrix}] .

(31)

Now assume that $𝑓 : 𝒦︀ \to ℝ$ is an $𝑚$ -LHSCB function on $𝒦︀$ with conjugate barrier $𝑓_{*}$ for some $𝑚 \geq 1 .$ Given any initial $𝑣^{0}$ , we define the central path of $𝑣^{*} (𝜇)$ as the unique solution to

\begin{aligned} 𝐺 (𝑣) = & 𝜇 𝐺 (𝑣^{0}) \\ 𝑠 = - 𝜇 \nabla 𝑓_{*} (𝑧), & 𝑧 = - 𝜇 \nabla 𝑓 (𝑠), 𝜏 𝜅 = 𝜇 \end{aligned}

(32)

Now inserting the definition of $𝑧$ while using that $𝑠^{⊤} \nabla 𝑓 (𝑠) = - 𝑚$ it follows that

𝑠^{⊤} \frac{𝑧}{𝑚} = 𝑠^{⊤} \frac{- 𝜇 \nabla 𝑓 (𝑠)}{𝑚} = 𝜇 \frac{𝑚}{𝑚} = 𝜇 .

(33)

If the cone $𝒦︀$ is symmetric, then we can replace the condition with a simpler condition defined in terms of the Jordan product on $𝒦︀$ , i.e.

𝑠 ⚬ 𝑧 = 𝜇 𝑒,

(34)

where $𝑒$ is the idempotent for each of the symmetric cones.

Interactive figure 1: A central path in two dimensions

Interactive figure 1 (A central path in two dimensions) Consider $𝑥 \in ℝ^{2}$ , the objective $𝑞^{⊤} 𝑥$ , and the composite cone $𝒦︀ = ℝ_{+} \times {𝒬︀}_{3} \times {𝒦︀}_{exp}$ with

𝑠 = 𝑏 - 𝐴 𝑥 = [\begin{matrix} 1.2 \\ 2 \\ 0 \\ 0 \\ - 0.5 \\ 1 \\ 1.5 \end{matrix}] - [\begin{matrix} 0 & - 1 \\ 0 & 0 \\ - 1 & 0 \\ 0 & - 1 \\ - 1 & 0 \\ 0 & 0 \\ 0 & - 1 \end{matrix}] [\begin{matrix} 𝑥_{1} \\ 𝑥_{2} \end{matrix}] = [\begin{matrix} 𝑥_{2} + 1.2 \\ 2 \\ 𝑥_{1} \\ 𝑥_{2} \\ 𝑥_{1} - 0.5 \\ 1 \\ 𝑥_{2} + 1.5 \end{matrix}] \in 𝒦︀ .

(35)

The three cones give, respectively, the half-plane $𝑥_{2} \geq - 1.2$ , the disk ${‖ 𝑥 ‖}_{2} \leq 2$ , and the region $𝑥_{1} \leq 0.5 + log (𝑥_{2} + 1.5)$ , and the barrier is

\begin{aligned} 𝑓 (𝑠) & = - log (𝑠_{1}) - log (𝑠_{2}^{2} - 𝑠_{3}^{2} - 𝑠_{4}^{2}) \\ - log (𝑠_{6} log (𝑠_{7} / 𝑠_{6}) - 𝑠_{5}) - log (𝑠_{6}) - log (𝑠_{7}), \end{aligned}

(36)

of degree $𝑚 = 1 + 2 + 3 = 6$ (here the second-order cone uses $- log (𝑠^{⊤} 𝐽 𝑠)$ , twice the barrier of the previous section, so it contributes degree 2). The central path is the curve of minimizers $𝑥 (𝜇)$ of $𝑞^{⊤} 𝑥 + 𝜇 𝑓 (𝑏 - 𝐴 𝑥)$ and along it $𝑧 = - 𝜇 \nabla 𝑓 (𝑠)$ . For large $𝜇$ the barrier dominates, and $𝑥 (𝜇)$ lies close to the analytic center, well inside the feasible region. As $𝜇$ decreases the linear term takes over, and the level sets are pushed against the boundary, where they follow its shape.

- {log}_{10} 𝜇

-1

4.1. Scaling matrices

We now introduce the idea of a scaling matrix $𝐻$ . The choice of scaling matrix $𝐻$ depends on which way we linearize the central path. For symmetric cones the most common choice is the Nesterov-Todd scaling. For nonsymmetric cones the central path is defined by the set of points satisfying (32). For the zero cone we set $𝐻 = 0$ .

4.1.1. Symmetric cones

For symmetric cones the central path described with (34) is linearized. The Nesterov-Todd scaling exploits the self-scaled property of symmetric cones to compute a unique scaling point $𝑤 \in 𝒦︀$ satisfying

𝐻 (𝑤) 𝑠 = 𝑧 .

(37)

We can factorize $𝐻 (𝑤)$ as $𝐻^{- 1} (𝑤) = 𝑊^{⊤} 𝑊$ .

4.1.2. Nonsymmetric cones

As nonsymmetric cones are not self-scaled, we cannot use (34) to describe the central path and must instead linearize (32). A general primal-dual scaling strategy suitable for nonsymmetric cones was proposed in [3], and used later in [4], which relies on the satisfaction of the two secant equations

𝐻 𝑧 = 𝑠, 𝐻 \nabla 𝑓 (𝑠) = \nabla 𝑓_{*} (𝑧) .

(38)

We now define shadow iterates as

\tilde{𝑧} : = - \nabla 𝑓 (𝑠), \tilde{𝑠} : = - \nabla 𝑓_{*} (𝑧),

(39)

with

\tilde{𝜇} = {\tilde{𝑠}}^{⊤} \tilde{𝑧} / 𝑚

(40)

A scaling matrix $𝐻$ can be obtained from the rank-4 Broyden-Fletcher-Goldfarb-Shanno (BFGS) update, which is commonly used in quasi-Newton methods,

𝐻 : = 𝐻_{BFGS} : = {𝑆 (𝑆^{⊤} 𝑍)}^{- 1} 𝑆^{⊤} + 𝐻_{𝑎} - 𝐻_{𝑎} {𝑍 (𝑍^{⊤} 𝐻_{𝑎} 𝑍)}^{- 1} 𝑍^{⊤} 𝐻_{𝑎},

(41)

where $𝑍 : = [𝑧, \tilde{𝑧}]$ , $𝑆 : = [𝑠, \tilde{𝑠}]$ , and $𝐻_{𝑎} ≻ 0$ is an approximation of the Hessian. In our implementation we choose $𝐻_{𝑎} = 𝜇 \nabla^{2} 𝑓_{*} (𝑧)$ following [4]. Putting things together means that the computation of $𝐻_{BFGS}$ reduces to a rank-3-update as

\begin{aligned} 𝐻_{BFGS} = 𝜇 \nabla^{2} 𝑓_{*} (𝑧) + & {(2 𝜇 𝑚)}^{- 1} 𝛿_{𝑠} {(𝑠 + 𝜇 \tilde{𝑠} + {(𝜇 \tilde{𝜇} - 1)}^{- 1} 𝛿_{𝑠})}^{⊤} \\ + & {(2 𝜇 𝑚)}^{- 1} (𝑠 + 𝜇 \tilde{𝑠} + {(𝜇 \tilde{𝜇} - 1)}^{- 1} 𝛿_{𝑠}) 𝛿_{𝑠}^{⊤} \\ - & 𝜇 \frac{(\nabla^{2} 𝑓_{*} (𝑧) \tilde{𝑧} - \tilde{𝜇} \tilde{𝑠}) {(\nabla^{2} 𝑓_{*} (𝑧) \tilde{𝑧} - \tilde{𝜇} \tilde{𝑠})}^{⊤}}{⟨ \tilde{𝑧}, \nabla^{2} 𝑓_{*} (𝑧) \tilde{𝑧} ⟩ - 𝑚 {\tilde{𝜇}}^{2}}, \end{aligned}

(42)

where $𝜇 = 𝑠^{⊤} 𝑧 / 𝑚$ and $𝛿_{𝑠} ≔ 𝑠 - 𝜇 \tilde{𝑠}$ . A more detailed description can be found in [4].

The above is defined in terms of the conjugate barrier function $𝑓_{* (𝑧)}$ . This is in contrast to [4], which defines it in terms of the standard barrier function.

5. Computing step directions

The overall aim is to find the optimal solution to the homogeneous embedding. This amounts to solving the nonlinear system of equations $𝐺 (𝑣) = 0$ with the additional constraint $𝜅 𝜏 = 𝜇$ . The aim of most interior-point methods is to track the central path for $𝜇 \to 0$ until reasonably close to the optimum. In order to solve the non-linear system Newton’s method is applied. This requires a linearization of the central path which can be seen below,

\begin{aligned} [\begin{matrix} 0 \\ Δ 𝑠 \\ Δ 𝜅 \end{matrix}] - & [\begin{matrix} 𝑃 & 𝐴^{⊤} & 𝑞 \\ - 𝐴 & 0 & 𝑏 \\ - {(𝑞 + 2 𝑃 𝜉)}^{⊤} & - 𝑏^{⊤} & 𝜉^{⊤} 𝑃 𝜉 \end{matrix}] [\begin{matrix} Δ 𝑥 \\ Δ 𝑧 \\ Δ 𝜏 \end{matrix}] = - [\begin{matrix} 𝑑_{𝑥} \\ 𝑑_{𝑧} \\ 𝑑_{𝜏} \end{matrix}] \\ 𝐻 Δ 𝑧 + Δ 𝑠 = - 𝑑_{𝑠}, 𝜅 Δ 𝜏 + 𝜏 Δ 𝜅 = - 𝑑_{𝜅}, \end{aligned}

(43)

where $𝑑 = (𝑑_{𝑥}, 𝑑_{𝑧}, 𝑑_{𝜏}, 𝑑_{𝑠}, 𝑑_{𝜅})$ is the vector of residuals, while $𝜉 = \frac{𝑥}{𝜏}$ and $𝐻 \in ℝ^{𝑚 \times 𝑚}$ is a positive semidefinite block-diagonal matrix (zero on zero-cone blocks) with blocks equal to the scaling matrices ⁶⁶ Ignoring the sparse form of the scaling matrices that are usually required in practice. This will be explained in a later section.. Note that in order to allow for nonsymmetric cones we use $𝐻$ , which is in contrast to [2] that uses $𝜆 ⚬ (𝑊 Δ 𝑧 + 𝑊^{- ⊤} Δ 𝑠) = 𝑑_{𝑠}$ with $𝐻 = 𝑊^{⊤} 𝑊$ .

We now eliminate $(Δ 𝑠, Δ 𝜅)$ from the system resulting in

\begin{aligned} - & [\begin{matrix} 𝑃 & 𝐴^{⊤} & 𝑞 \\ - 𝐴 & 𝐻 & 𝑏 \\ - {(𝑞 + 2 𝑃 𝜉)}^{⊤} & - 𝑏^{⊤} & 𝜉^{⊤} 𝑃 𝜉 + 𝜅 𝜏^{- 1} \end{matrix}] [\begin{matrix} Δ 𝑥 \\ Δ 𝑧 \\ Δ 𝜏 \end{matrix}] = - [\begin{matrix} 𝑑_{𝑥} \\ 𝑑_{𝑧} - 𝑑_{𝑠} \\ 𝑑_{𝜏} - 𝑑_{𝜅} 𝜏^{- 1} \end{matrix}] \\ Δ 𝑠 = - 𝑑_{𝑠} - 𝐻 Δ 𝑧, Δ 𝜅 = - (𝑑_{𝜅} + 𝜅 Δ 𝜏) 𝜏^{- 1} . \end{aligned}

(44)

The above system can be solved by a pair of linear systems with a common coefficient matrix, meaning that it only has to be factorized once

[\begin{matrix} 𝑃 & 𝐴^{⊤} \\ 𝐴 & - 𝐻 \end{matrix}] [\begin{matrix} Δ 𝑥_{1} & Δ 𝑥_{2} \\ Δ 𝑧_{1} & Δ 𝑧_{2} \end{matrix}] = [\begin{matrix} 𝑑_{𝑥} & - 𝑞 \\ - (𝑑_{𝑧} - 𝑑_{𝑠}) & 𝑏 \end{matrix}]

(45)

From the above solution we can recover the search direction $(Δ 𝑥, Δ 𝑧, Δ 𝑠, Δ 𝜏, Δ 𝜅)$ as follows

\begin{aligned} Δ 𝑥 & = Δ 𝑥_{1} + Δ 𝜏 \cdot Δ 𝑥_{2}, \\ Δ 𝑧 & = Δ 𝑧_{1} + Δ 𝜏 \cdot Δ 𝑧_{2} \\ Δ 𝑠 & = - 𝑑_{𝑠} - 𝐻 Δ 𝑧, \\ Δ 𝜏 & = \frac{𝑑_{𝜏} - 𝑑_{𝜅} 𝜏^{- 1} + {(2 𝑃 𝜉 + 𝑞)}^{⊤} Δ 𝑥_{1} + 𝑏^{⊤} Δ 𝑧_{1}}{𝜅 𝜏^{- 1} + 𝜉^{⊤} 𝑃 𝜉 - {(2 𝑃 𝜉 + 𝑞)}^{⊤} Δ 𝑥_{2} - 𝑏^{⊤} Δ 𝑧_{2}} \\ = \frac{𝑑_{𝜏} - 𝑑_{𝜅} 𝜏^{- 1} + 𝑞^{⊤} Δ 𝑥_{1} + 𝑏^{⊤} Δ 𝑧_{1} + 2 𝜉^{⊤} 𝑃 Δ 𝑥_{1}}{𝜅 𝜏^{- 1} + ‖ Δ 𝑥_{2} - 𝜉 ‖_{𝑃}^{2} - ‖ Δ 𝑥_{2} ‖_{𝑃}^{2} - 𝑞^{⊤} Δ 𝑥_{2} - 𝑏^{⊤} Δ 𝑧_{2}} \\ Δ 𝜅 & = - (𝑑_{𝜅} + 𝜅 Δ 𝜏) 𝜏^{- 1} . \end{aligned}

(46)

The actual step is then obtained by computing the maximal step size $𝛼$ that ensures that the new update is still in the interior of the conic constraints.

Interactive figure 2: The largest step inside the cones

Interactive figure 2 (The largest step inside the cones) The feasible region of the central path example, with an iterate $𝑥$ and a search direction $Δ 𝑥$ , both of which can be dragged. Following $Δ 𝑥$ from $𝑥$ , the first cone boundary is reached at the step $𝛼_{max}$ (dashed). Clarabel moves only a fraction of the way, $𝑥 + 𝛼 Δ 𝑥$ with $𝛼 = 0.99 min (1, 𝛼_{max})$ , so that the iterate stays strictly inside the cones.

5.1. The affine and centering directions

At every interior-point iteration we solve the KKT system for two sets of right-hand sides each, corresponding to the so-called affine and centering directions. In short this means solving (45) twice, however since the system matrix stays the same and the right-hand side $(- 𝑞, 𝑏)$ stays the same the computation requires only a single numerical factorization ⁷⁷ As the sparsity pattern does not change the symbolic factorization can be reused and three solves.

The two sets of right-hand sides correspond to a predictor step and a corrector step. The steps are also known as the affine step and the centering step. In our case the affine step has the right-hand side of

𝑑_{𝑥} = 𝑟_{𝑥}, 𝑑_{𝑧} = 𝑟_{𝑧}, 𝑑_{𝜏} = 𝑟_{𝜏}, 𝑑_{𝑠} = 𝑠, 𝑑_{𝜅} = 𝜅 𝜏,

(47)

where $(𝑟_{𝑥}, 𝑟_{𝑧}, 𝑟_{𝜏})$ are the three components of the residual $𝐺 (𝑣)$ at the current iterate $𝑣$ , while the corrector step has the right-hand side of

\begin{aligned} 𝑑_{𝑥} & = (1 - 𝜎) 𝑟_{𝑥}, 𝑑_{𝑧} = (1 - 𝜎) 𝑟_{𝑧}, 𝑑_{𝜏} = (1 - 𝜎) 𝑟_{𝜏}, \\ 𝑑_{𝜅} & = 𝜅 𝜏 + Δ 𝜅 Δ 𝜏 - 𝜎 𝜇, \\ 𝑑_{𝑠} & = {\begin{cases} 𝑊^{⊤} (𝜆 ⋄ (𝜆 ⚬ 𝜆 - 𝜎 𝜇 𝑒 + 𝜂)) & (symmetric) \\ 𝑠 + 𝜎 𝜇 \nabla 𝑓_{*} (𝑧) + 𝜂 & (nonsymmetric), \end{cases} \end{aligned}

(48)

where $⋄$ denotes the inverse operator of $⚬$ . Computation of a higher-order correction term $𝜂$ is a heuristic technique that is known to accelerate the convergence of IPMs significantly. The choice of this term varies depending on the choice of the scaling matrix $𝐻$ and whether a given constraint is symmetric or not. For symmetric cones we use the Mehrotra correction $𝜂 = (𝑊^{- ⊤} Δ 𝑠) ⚬ (𝑊 Δ 𝑧)$ , while for nonsymmetric cones we compute $𝜂$ using the 3rd-order correction from [4], i.e.

𝜂 = - \frac{1}{2} \nabla^{3} 𝑓_{*} (𝑧) [Δ 𝑧, \nabla^{2} 𝑓_{*} {(𝑧)}^{- 1} Δ 𝑠] .

(49)

We set $𝜎 = {(1 - 𝛼_{aff})}^{3}$ where $𝛼_{aff}$ is the step size of the affine step.

Interactive figure 3: Affine and centering directions

Interactive figure 3 (Affine and centering directions) A primal-only version on the central path example, with the dual variable tied to the iterate by $𝑧 = - 𝜇 \nabla 𝑓 (𝑠)$ , so that the step towards the central path for a target $𝜇$ is the Newton step $Δ 𝑥 = - \nabla^{2} {𝐹 (𝑥)}^{- 1} (𝑞 / 𝜇 + \nabla 𝐹 (𝑥))$ , where $𝐹 (𝑥) = 𝑓 (𝑏 - 𝐴 𝑥)$ . From $𝑥 (𝜇)$ , the affine step (dashed, target $𝜇 = 0$ ) is taken a fraction $𝛼 / 𝛼_{max}$ of the way to the boundary, which leaves the central path and reduces $𝜇$ to $𝜇^{+} = (1 - 𝛼) 𝜇$ . The centering step (dotted) is one damped Newton step for $𝜇^{+}$ .⁸⁸ The next iteration starts where the centering step ends, near but not on the central path. Only the number $𝜇^{+}$ enters this step; $𝑥 (𝜇^{+})$ is unknown to the method and drawn for reference only. This is the classical predictor-corrector pair; Clarabel combines the two steps into one through $𝜎$ .

- {log}_{10} 𝜇

𝛼 / 𝛼_{max}

0.9

5.2. Solver initialization

An essential consideration is how to choose the initial point $𝑣^{0}$ . The initialization strategies depend on whether the objective is linear or quadratic and whether the cones are symmetric or not. However, in all cases the initial scalar values are set as $𝜏^{0} = 𝜅^{0} = 1$ .

5.2.1. Symmetric cones

The initialization strategies below can be thought of as projections of the equality constraints onto cones in the direction of the idempotent $𝑒$ for each symmetric cone.

Quadratic objective. If $𝑃 \neq 0$ the initialization happens by solving the following linear system

[\begin{matrix} 𝑃 & 𝐴^{⊤} \\ 𝐴 & - 𝐼_{𝐸} \end{matrix}] [\begin{matrix} 𝑥 \\ 𝑧 \end{matrix}] = [\begin{matrix} - 𝑞 \\ 𝑏 \end{matrix}]

(50)

with $𝐼_{𝐸} = 𝐼 - 𝐸$ with $𝐸$ being a matrix with ones on the diagonal corresponding to zero cones (i.e. equality constraints) but otherwise zero. The above linear system corresponds to solving the following optimization problem

\begin{aligned} minimize & \frac{1}{2} 𝑥^{⊤} 𝑃 𝑥 + 𝑞^{⊤} 𝑥 + \frac{1}{2} ‖ 𝐴 𝑥 - 𝑏 ‖_{2}^{2} \\ subject to & 𝐸 𝐴 𝑥 = 𝐸 𝑏, \end{aligned}

(51)

From this solution we project $𝑧$ and $𝑠$ back to the interior of the cone

\begin{aligned} 𝑥^{0} & = 𝑥, \\ 𝑠^{0} & = {\begin{cases} - 𝑧, & 𝛼_{p} < - 𝜖 \\ - 𝑧 + (𝜖 + 𝛼_{p}) 𝑒, & otherwise \end{cases}, \\ 𝑧^{0} & = {\begin{cases} 𝑧, & 𝛼_{d} < - 𝜖 \\ 𝑧 + (𝜖 + 𝛼_{d}) 𝑒, & otherwise \end{cases} . \end{aligned}

(52)

where $𝛼_{p} = inf {𝛼 | - 𝑧 + 𝛼 𝑒 \in 𝒦︀}$ and $𝛼_{d} = inf {𝛼 | - 𝑠 + 𝛼 𝑒 \in 𝒦︀}$ .

Linear objective. In the case of $𝑃 = 0$ we instead solve two linear systems

[\begin{matrix} 0 & 𝐴^{⊤} \\ 𝐴 & - 𝐼_{𝐸} \end{matrix}] [\begin{matrix} 𝑥 & 𝑡 \\ 𝑠 & 𝑧 \end{matrix}] = [\begin{matrix} 0 & - 𝑞 \\ 𝑏 & 0 \end{matrix}] .

(53)

From this we set

\begin{aligned} 𝑥^{0} & = 𝑥, \\ 𝑠^{0} & = {\begin{cases} - 𝑠, & 𝛼_{p} < - 𝜖 \\ - 𝑠 + (𝜖 + 𝛼_{p}) 𝑒, & otherwise \end{cases}, \\ 𝑧^{0} & = {\begin{cases} 𝑧, & 𝛼_{d} < - 𝜖 \\ 𝑧 + (𝜖 + 𝛼_{d}) 𝑒, & otherwise \end{cases}, \end{aligned}

(54)

where $𝛼_{p} = inf {𝛼 | - 𝑠 + 𝛼 𝑒 \in 𝒦︀}$ and $𝛼_{d} = inf {𝛼 | 𝑧 + 𝛼 𝑒 \in 𝒦︀}$ .

5.2.2. Nonsymmetric cones

When the composite cone $𝒦︀$ contains any nonsymmetric cone both primal and dual variables are set as $𝑧 = 𝑠 = - \nabla 𝑓_{* (𝑧)}$ (corresponding to $𝜇^{0} = 1$ ). This is equivalent to solving the unconstrained optimization problem

\underset{𝑧}{minimize} \frac{1}{2} ‖ 𝑧 ‖_{2}^{2} + 𝑓_{*} (𝑧),

(55)

which is strictly convex and has a unique solution. It yields $𝑠_{sym}^{0} = 𝑧_{sym}^{0} = 𝑒$ for symmetric cones, and

𝑠_{exp}^{0} = 𝑧_{exp}^{0} \approx (- 1.051383, 0.556409, 1.258967),

(56)

for exponential cones, and

𝑠_{pow}^{0} = 𝑧_{pow}^{0} \approx (\sqrt{1 + 𝛼}, \sqrt{2 - 𝛼}, 0),

(57)

for power cones with parameter $𝛼$ .

5.3. Termination criteria

All termination criteria are based on unscaled problem data and iterates, i.e. after the Ruiz scaling has been reverted.

5.3.1. Feasibility checks

For checks of primal and dual feasibility we introduce the normalized variables $\hat{𝑥} = 𝑥 / 𝜏$ , $\hat{𝑠} = 𝑠 / 𝜏$ , $\hat{𝑧} = 𝑧 / 𝜏$ and define the primal and dual residual as follows

\begin{aligned} 𝑟_{𝑝} & : = 𝐴 \hat{𝑥} + \hat{𝑠} - 𝑏 & (primal residual) \\ 𝑟_{𝑑} & : = 𝑃 \hat{𝑥} + 𝐴^{⊤} \hat{𝑧} + 𝑞 & (dual residual) . \end{aligned}

(58)

Likewise, we define the primal and dual objectives as

\begin{aligned} 𝑔_{𝑝} & : = \frac{1}{2} {\hat{𝑥}}^{⊤} 𝑃 \hat{𝑥} + 𝑞^{⊤} \hat{𝑥} & (primal objective) \\ 𝑔_{𝑑} & : = - \frac{1}{2} {\hat{𝑥}}^{⊤} 𝑃 \hat{𝑥} - 𝑏^{⊤} \hat{𝑧} & (dual objective) . \end{aligned}

(59)

Using the above definitions we can declare convergence if all three of the following hold:

\begin{aligned} ‖ 𝑟_{𝑝} ‖ & < 𝜖_{𝑓} \cdot max {1, ‖ 𝑏 ‖_{\infty} + ‖ \hat{𝑥} ‖ + ‖ \hat{𝑠} ‖} & (primal convergence) \\ ‖ 𝑟_{𝑑} ‖ & < 𝜖_{𝑓} \cdot max {1, ‖ 𝑞 ‖_{\infty} + ‖ \hat{𝑥} ‖ + ‖ \hat{𝑧} ‖} & (dual convergence) \\ | 𝑔_{𝑝} - 𝑔_{𝑑} | & < 𝜖_{𝑓} \cdot max {1, min {| 𝑔_{𝑝} |, | 𝑔_{𝑑} |}} & (gap convergence) \end{aligned}

(60)

We specify a default value of $𝜖_{𝑓} = 10^{- 8}$ as well as a weaker threshold of $𝜖_{𝑓} = 10^{- 5}$ when testing for “near optimality” in cases of early termination (e.g. lack of progress, timeout, iterations limit, etc.).

5.3.2. Infeasibility checks

When testing for infeasibility we do not normalize iterates, but rather work directly with the unscaled variables since infeasibility corresponds to the case where $𝜏 \to 0$ . We declare primal infeasibility if the following holds

\begin{aligned} ‖ 𝐴^{⊤} 𝑧 ‖ & < - 𝜖_{𝑖, 𝑟} \cdot max {1, ‖ 𝑥 ‖ + ‖ 𝑧 ‖} \cdot (𝑏^{⊤} 𝑧) & (primal infeasibility) \\ 𝑏^{⊤} 𝑧 & < - 𝜖_{𝑖, 𝑎} & (primal infeasibility) . \end{aligned}

(61)

Similarly, we define dual infeasibility if the following hold

\begin{aligned} ‖ 𝑃 𝑥 ‖ & < - 𝜖_{𝑖, 𝑟} \cdot max {1, ‖ 𝑥 ‖} \cdot (𝑞^{⊤} 𝑥) & (dual infeasibility) \\ ‖ 𝐴 𝑥 + 𝑠 ‖ & < - 𝜖_{𝑖, 𝑟} \cdot max {1, ‖ 𝑥 ‖ + ‖ 𝑠 ‖} \cdot (𝑞^{⊤} 𝑥) & (dual infeasibility) \\ 𝑞^{⊤} 𝑥 & < - 𝜖_{𝑖, 𝑎} & (dual infeasibility) \end{aligned}

(62)

We set the relative and absolute tolerances as $𝜖_{𝑖, 𝑟} = 𝜖_{𝑖, 𝑎} = 10^{- 8}$ and allow for weaker thresholds to declare “near infeasibility” certificates in cases of early termination.

6. Implementation details

6.1. Solving the KKT system

The main computational effort in the interior-point method is to solve the Karush-Kuhn-Tucker (KKT) system of the form

\underset{𝐾}{\underset{⏟}{[\begin{matrix} 𝑃 & 𝐴^{⊤} \\ 𝐴 & - 𝐻 \end{matrix}]}} \underset{𝑦}{\underset{⏟}{[\begin{matrix} 𝑥 \\ 𝑧 \end{matrix}]}} = \underset{𝑏}{\underset{⏟}{[\begin{matrix} 𝑏_{𝑥} \\ 𝑏_{𝑧} \end{matrix}]}},

(63)

The scaling matrix, however, has the unfortunate property that it becomes ill-conditioned when the iterates are close to the boundary. As such one instead solves the regularized problem

\underset{𝐾 + 𝐷}{\underset{⏟}{[\begin{matrix} 𝑃 + 𝜖 𝐼 & 𝐴^{⊤} \\ 𝐴 & - (𝐻 + 𝜖 𝐼) \end{matrix}]}} \underset{𝑦}{\underset{⏟}{[\begin{matrix} 𝑥 \\ 𝑧 \end{matrix}]}} = \underset{𝑏}{\underset{⏟}{[\begin{matrix} 𝑏_{𝑥} \\ 𝑏_{𝑧} \end{matrix}]}},

(64)

where $𝐷$ denotes the added diagonal part. Note that while solving (64) is more numerically stable, it is not the system that we actually aim to solve i.e. (63). The regularization $𝜖$ is computed as

𝜖 = 𝜖_{constant} + 𝜖_{proportion} ‖ diag (𝐾) ‖_{\infty} .

(65)

6.2. Regularized iterative refinement

As described in [5], the regularization term can be removed by solving a sequence of systems of the form

(𝐾 + 𝐷) Δ 𝑥_{𝑚} = 𝑏 - 𝐾 𝑥_{𝑚}, 𝑥_{𝑚 + 1} = 𝑥_{𝑚} + Δ 𝑥_{𝑚} .

(66)

The sequence ${𝑥_{𝑚}}$ converges to the true solution of $𝐾 𝑥 = 𝑏$ , provided that it exists. To show this first notice that we can rewrite (66) as

(𝐾 + 𝐷) 𝑥_{𝑚 + 1} - 𝐷 𝑥_{𝑚} = 𝑏 .

(67)

Now assume that $𝑥^{0}$ is the solution to the equation $𝐾 𝑥^{0} = 𝑏$ (which is the equation we aim to solve). Now we perform the mathematician’s favorite trick of adding 0 (in this case $𝐷 𝑥^{0} - 𝐷 𝑥^{0}$ ), it follows that

(𝐾 + 𝐷) 𝑥^{0} - 𝐷 𝑥^{0} = 𝑏 .

(68)

We now subtract (68) from (67)

(𝐾 + 𝐷) (𝑥_{𝑚 + 1} - 𝑥^{0}) - 𝐷 (𝑥_{𝑚} - 𝑥^{0}) = 0 .

(69)

Rewriting the above gives us

𝑥_{𝑚 + 1} - 𝑥^{0} = {(𝐾 + 𝐷)}^{- 1} 𝐷 (𝑥_{𝑚} - 𝑥^{0}) .

(70)

The final step is to take norms on both sides resulting in

‖ 𝑥_{𝑚 + 1} - 𝑥^{0} ‖ \leq ‖ {(𝐾 + 𝐷)}^{- 1} 𝐷 ‖ ‖ 𝑥_{𝑚} - 𝑥^{0} ‖ .

(71)

Thus, if $‖ {(𝐾 + 𝐷)}^{- 1} 𝐷 ‖ < 1$ we have that $𝑥_{𝑚} \to 𝑥^{0}$ for $𝑚 \to \infty$ (more precisely, it suffices that the spectral radius of ${(𝐾 + 𝐷)}^{- 1} 𝐷$ is less than one).

6.3. Matrix equilibration

The equilibration in Clarabel is a modified version of Ruiz scaling [6], [7], similar as to what is done in the OSQP package [5]. The equilibration is based on the scaling of the constant part of the KKT system (i.e. the part excluding the primal-dual scaling matrices) as shown below

𝑀 = [\begin{matrix} 𝑃 & 𝐴^{⊤} \\ 𝐴 & 0 \end{matrix}] .

(72)

The aim is to find a diagonal scaling matrix $𝑆$ such that

𝑆 𝑀 𝑆 = [\begin{matrix} 𝐷 & 0 \\ 0 & 𝐸 \end{matrix}] [\begin{matrix} 𝑃 & 𝐴^{⊤} \\ 𝐴 & 0 \end{matrix}] [\begin{matrix} 𝐷 & 0 \\ 0 & 𝐸 \end{matrix}] = [\begin{matrix} 𝐷 𝑃 𝐷 & 𝐷 𝐴^{⊤} 𝐸 \\ 𝐸 𝐴 𝐷 & 0 \end{matrix}],

(73)

is better conditioned. The above can be seen as a transformation of input data to a conic optimization problem

\begin{aligned} minimize & \frac{1}{2} {\bar{𝑥}}^{⊤} \bar{𝑃} \bar{𝑥} + {\bar{𝑞}}^{⊤} \bar{𝑥} \\ subject to & \bar{𝐴} \bar{𝑥} + \bar{𝑠} = \bar{𝑏} \\ \bar{𝑠} \in 𝒦︀, \end{aligned}

(74)

where $\bar{𝑥} = 𝐷^{- 1} 𝑥$ , $\bar{𝑠} = 𝐸 𝑠$ , $\bar{𝑃} = 𝑐 𝐷 𝑃 𝐷$ , $\bar{𝑞} = 𝑐 𝐷 𝑞$ , $\bar{𝐴} = 𝐸 𝐴 𝐷$ , $\bar{𝑏} = 𝐸 𝑏$ , with $𝑐 > 0$ a cost scaling and $𝐸$ chosen such that $𝐸 𝒦︀ = 𝒦︀$ .

The main idea of the equilibration procedure is to scale the rows (and columns) for the matrix $𝑀$ so that they have equal $ℓ_{𝑝}$ -norm. An efficient approximate solution to this problem is the so-called Modified Ruiz scaling described in [5] based on the Ruiz scaling described in [6].

7. Sparse solvers

More often than not the main computational time in an interior-point method is solving (64). Luckily, the linear system is often sparse leading to fast factorization algorithms. As such it is important to use an efficient implementation of a direct sparse solver if one wants the interior-point solver to be computationally efficient. In Clarabel it is possible to utilize the following sparse solvers

(Rust) faer-rs. Multithreaded. Wasm compatible.
(Rust/Julia) Quasi-Definite LDL (QDLDL) based on [8]. Single threaded.
(Rust/Julia) PanuaPardiso (Requires license). Multithreaded.
(Rust/Julia) oneMKL-Pardiso (Only Intel CPUs). Multithreaded.
(Julia) cuDSS (CUDA Direct Sparse Solver). GPU-based.
(Julia) MA57 (Requires license). Multithreaded.
(Julia) CHOLMOD (Standard Julia Solver). Multithreaded.

7.1. Sparsity exploitation

7.1.1. Augmented sparsity

In numerical computations one often meets a matrix of the form sparse-plus-low-rank. The general form of this is the following

𝐻 : = 𝑆 - \sum_{𝑖 = 1}^{𝑛_{𝑉}} 𝑣_{𝑖} 𝑣_{𝑖}^{⊤} + \sum_{𝑖 = 1}^{𝑛_{𝑈}} 𝑢_{𝑖} 𝑢_{𝑖}^{⊤} = 𝑆 - 𝑉 𝑉^{⊤} + 𝑈 𝑈^{⊤},

(75)

where $𝑆$ is sparse and $𝑉 = [\begin{matrix} 𝑣_{1} & \dots & 𝑣_{𝑛_{𝑉}} \end{matrix}]$ and $𝑈 = [\begin{matrix} 𝑢_{1} & \dots & 𝑢_{𝑛_{𝑈}} \end{matrix}]$ with $𝑛_{𝑈}, 𝑛_{𝑉} ≪ 𝑛$ . The problem with the above matrix is that forming the outer products (i.e. assembling the plus-low-rank part) will ruin the sparsity of $𝑆$ and $𝐻$ in turn will be a dense matrix. Given that our aim is to solve a system of the form $𝐻 𝑥 = 𝑏$ lets look at the structure of $𝐻 𝑥$ by performing the product

𝐻 𝑥 = 𝑆 𝑥 - 𝑉 (\underset{- 𝑡_{𝑉}}{\underset{⏟}{𝑉^{⊤} 𝑥}}) + 𝑈 (\underset{𝑡_{𝑈}}{\underset{⏟}{𝑈^{⊤} 𝑥}}) = 𝑆 𝑥 + 𝑉 𝑡_{𝑉} + 𝑈 𝑡_{𝑈} .

(76)

The above together with the constraints $𝑉^{⊤} 𝑥 + 𝑡_{𝑉} = 0$ and $𝑈^{⊤} 𝑥 - 𝑡_{𝑈} = 0$ we can solve $𝐻 𝑥 = 𝑏$ by solving the following equivalent system

\underset{𝐻_{aug}}{\underset{⏟}{[\begin{matrix} 𝑆 & 𝑉 & 𝑈 \\ 𝑉^{⊤} & 𝐼_{𝑛_{𝑉}} & 0 \\ 𝑈^{⊤} & 0 & - 𝐼_{𝑛_{𝑈}} \end{matrix}]}} [\begin{matrix} 𝑥 \\ 𝑡_{𝑉} \\ 𝑡_{𝑈} \end{matrix}] = [\begin{matrix} 𝑏 \\ 0 \\ 0 \end{matrix}] .

(77)

Note that while $𝐻_{aug}$ is sparse it is larger. An important property of the augmented system is that if $𝑆 - 𝑉 𝑉^{⊤} ≻ 0$ (and hence $𝑆 ≻ 0$ ) then $𝐻_{aug}$ is Quasi-Definite and therefore strongly factorizable, meaning that for every permutation $𝑃$ one can find $𝐿$ and $𝐷$ so that $𝑃 𝐻_{aug} 𝑃^{⊤} = 𝐿 𝐷 𝐿^{⊤}$ [9].

7.2. Second-order cones

7.2.1. The Clarabel approach

Clarabel uses the same approach as in ECOS [10]. In short, they utilize that the structure of $𝑊_{𝑘}^{⊤} 𝑊_{𝑘}$ for a second-order cone is augmented-sparse, meaning that it can be split as follows

𝑊_{𝑘}^{⊤} 𝑊_{𝑘} = 𝑆_{𝑘} - 𝑣_{𝑘} 𝑣_{𝑘}^{⊤} + 𝑢_{𝑘} 𝑢_{𝑘}^{⊤},

(78)

where $𝑆_{𝑘} - 𝑣_{𝑘} 𝑣_{𝑘}^{⊤} ≻ 0$ . We can generalize the above for all second-order cones by defining

\begin{aligned} 𝑆 & = blkdiag (𝑆_{1}, 𝑆_{2}, \dots, 𝑆_{𝑛}), \\ 𝑉 & = blkdiag (𝑣_{1}, 𝑣_{2}, \dots, 𝑣_{𝑛}), \\ 𝑈 & = blkdiag (𝑢_{1}, 𝑢_{2}, \dots, 𝑢_{𝑛}), \end{aligned}

(79)

meaning that

𝑊^{⊤} 𝑊 = 𝑆 - 𝑉 𝑉^{⊤} + 𝑈 𝑈^{⊤} .

(80)

As such the KKT of the form

[\begin{matrix} 𝑃 & 𝐴^{⊤} \\ 𝐴 & - 𝑊^{⊤} 𝑊 \end{matrix}] [\begin{matrix} 𝑥 \\ 𝑧 \end{matrix}] = [\begin{matrix} 𝑏_{𝑥} \\ 𝑏_{𝑧} \end{matrix}]

(81)

can equivalently be solved by solving

[\begin{matrix} 𝑃 & 𝐴^{⊤} & 0 & 0 \\ 𝐴 & - 𝑆 & - 𝑉 & - 𝑈 \\ 0 & - 𝑉^{⊤} & - 𝐼 & 0 \\ 0 & - 𝑈^{⊤} & 0 & 𝐼 \end{matrix}] [\begin{matrix} 𝑥 \\ 𝑧 \\ 𝑡_{𝑉} \\ 𝑡_{𝑈} \end{matrix}] = [\begin{matrix} 𝑏_{𝑥} \\ 𝑏_{𝑧} \\ 0 \\ 0 \end{matrix}],

(82)

The approach is simple as it leaves $𝐴$ intact and does not scale $𝑧$ . However, one has to deal with computing $𝑉$ and $𝑈$ , which is not unique. In fact Clarabel and ECOS apply different approaches when calculating $𝑉$ and $𝑈$ . A detail here is that one wants to compute $𝑉$ and $𝑈$ in a way for which $𝑆 - 𝑉 𝑉^{⊤} ≻ 0$ so that the system remains quasi-definite (and therefore strongly factorizable).

7.2.2. The CVXOPT approach

Another approach is that of CVXOPT. Here the diagonal-plus-low-rank structure of the scaling matrix of a second-order cone (and its inverse) is used, i.e. that

𝑊_{𝑘} = 𝐷_{𝑘} + 𝑢_{𝑘} 𝑢_{𝑘}^{⊤}, 𝑊_{𝑘}^{- 1} = 𝐸_{𝑘} + 𝑣_{𝑘} 𝑣_{𝑘}^{⊤},

(83)

where the $𝐷_{𝑘}$ and $𝐸_{𝑘}$ are diagonal. We can again bunch all the scaling matrices together by defining

𝐷 = blkdiag (𝐷_{1}, 𝐷_{2}, \dots, 𝐷_{𝑛}), 𝐸 = blkdiag (𝐸_{1}, 𝐸_{2}, \dots, 𝐸_{𝑛}),

(84)

so that

𝑊 = 𝐷 + 𝑈 𝑈^{⊤}, 𝑊^{- 1} = 𝐸 + 𝑉 𝑉^{⊤} .

(85)

Note that $𝑈$ and $𝑉$ are defined as before, but are different since the $𝑢_{𝑘}$ and $𝑣_{𝑘}$ come from (83) rather than (78). We now show how the rank-update and the known inverse can be used explicitly in order to preserve sparsity. The KKT system of interest using Nesterov-Todd scaling is the following

[\begin{matrix} 𝑃 & 𝐴^{⊤} \\ 𝐴 & - 𝑊^{⊤} 𝑊 \end{matrix}] [\begin{matrix} 𝑥 \\ 𝑧 \end{matrix}] = [\begin{matrix} 𝑏_{𝑥} \\ 𝑏_{𝑧} \end{matrix}]

(86)

By scaling the $𝑧$ -variable and multiplying the last row with $𝑊^{- 1}$ one finds that

[\begin{matrix} 𝑃 & 𝐴^{⊤} (𝐸 + 𝑉 𝑉^{⊤}) \\ (𝐸 + 𝑉 𝑉^{⊤}) 𝐴 & - 𝐼 \end{matrix}] [\begin{matrix} 𝑥 \\ 𝑊 𝑧 \end{matrix}] = [\begin{matrix} 𝑏_{𝑥} \\ 𝑊^{- 1} 𝑏_{𝑧} \end{matrix}] .

(87)

For simplicity, we now define $𝑆 = 𝐴^{⊤} 𝑉$

[\begin{matrix} 𝑃 & 𝐴^{⊤} 𝐸 + 𝑆 𝑉^{⊤} \\ 𝐸 𝐴 + 𝑉 𝑆^{⊤} & - 𝐼 \end{matrix}] [\begin{matrix} 𝑥 \\ 𝑊 𝑧 \end{matrix}] = [\begin{matrix} 𝑏_{𝑥} \\ 𝑊^{- 1} 𝑏_{𝑧} \end{matrix}] .

(88)

We can now introduce new variables $𝛼, 𝛽 \in ℝ^{𝑛_{𝑉}}$ resulting in

[\begin{matrix} 𝑃 & 𝐴^{⊤} 𝐸 & 𝑆 & 0 \\ 𝐸^{⊤} 𝐴 & - 𝐼 & 0 & 𝑉 \\ 𝑆^{⊤} & 0 & 0 & - 𝐼 \\ 0 & 𝑉^{⊤} & - 𝐼 & 0 \end{matrix}] [\begin{matrix} 𝑥 \\ 𝑊 𝑧 \\ 𝛼 \\ 𝛽 \end{matrix}] = [\begin{matrix} 𝑏_{𝑥} \\ 𝑊^{- 1} 𝑏_{𝑧} \\ 0 \\ 0 \end{matrix}] .

(89)

Where we have used that $𝐸 = 𝐸^{⊤}$ to show why the system is also symmetric. Now the above preserves sparsity. For SOCs, $𝐸_{𝑘} = - 𝛾_{𝑘}^{- 1} 𝐽$ for a scalar $𝛾_{𝑘} > 0$ (in CVXOPT, $𝑊_{𝑘} = 𝛾_{𝑘} (2 𝑤_{𝑘} 𝑤_{𝑘}^{⊤} - 𝐽)$ for a vector $𝑤_{𝑘}$ ), so $𝐸 𝐸$ is diagonal with blocks $𝛾_{𝑘}^{- 2} 𝐼$ . If the scaling is normalized so that $𝛾_{𝑘} = 1$ , then $𝐸 𝐸 = 𝐼$ , and we can reduce the above as

[\begin{matrix} 𝑃 & 𝐴^{⊤} & 𝑆 & 0 \\ 𝐴 & - 𝐼 & 0 & 𝐸 𝑉 \\ 𝑆^{⊤} & 0 & 0 & - 𝐼 \\ 0 & 𝑉^{⊤} 𝐸^{⊤} & - 𝐼 & 0 \end{matrix}] [\begin{matrix} 𝑥 \\ 𝐸 𝑊 𝑧 \\ 𝛼 \\ 𝛽 \end{matrix}] = [\begin{matrix} 𝑏_{𝑥} \\ 𝐸 𝑊^{- 1} 𝑏_{𝑧} \\ 0 \\ 0 \end{matrix}] .

(90)

7.3. Rotated second-order cones

By rotating the rotated second-order cone back to the second-order cone form we can exploit the sparsity similarly to the regular second-order cone. This rotation can be done by applying the following transformation

𝑅 = [\begin{matrix} \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} & 0 & \dots & 0 \\ \frac{1}{\sqrt{2}} & - \frac{1}{\sqrt{2}} & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & 0 & \dots & 1 \end{matrix}] .

(91)

7.4. Power cones

The conjugate barrier for certain power cones has a Hessian that, similar to that of second-order cones, can be represented using augmented sparsity [11]. The intuition here is that e.g. the generalized power cone reduces to a second-order cone when $𝑑_{1} = 1$ . Given that the augmented sparsity is only available for the conjugate barrier the following is only relevant in the case of dual scaling.

The generalized power cone is parametrized by $𝛼 \in ℝ_{+ +}^{𝑑_{1}}$ such that

\sum_{𝑖 \in [𝑑_{1}]} 𝛼_{𝑖} = 1

(92)

and is defined as

{𝒦︀}_{gpow} (𝛼, 𝑑_{1}, 𝑑_{2}) = {(𝑢, 𝑤) \in ℝ_{+}^{𝑑_{1}} \times ℝ^{𝑑_{2}} | \prod_{𝑖 = 1, \dots, 𝑑_{1}} 𝑢_{𝑖}^{𝛼_{𝑖}} \geq ‖ 𝑤 ‖}

(93)

with dual cone

{𝒦︀}_{gpow}^{*} (𝛼, 𝑑_{1}, 𝑑_{2}) = {(𝑢, 𝑤) \in ℝ_{+}^{𝑑_{1}} \times ℝ^{𝑑_{2}} | \prod_{𝑖 = 1, \dots, 𝑑_{1}} {(\frac{𝑢_{𝑖}}{𝛼_{𝑖}})}^{𝛼_{𝑖}} \geq ‖ 𝑤 ‖}

(94)

The power mean cone is parametrized by $𝛼 \in ℝ_{+ +}^{𝑑_{1}}$ such that

\sum_{𝑖 \in [𝑑_{1}]} 𝛼_{𝑖} = 1

(95)

and is defined as

{𝒦︀}_{powm} (𝛼, 𝑑_{1}) = {(𝑢, 𝑤) \in ℝ_{+}^{𝑑_{1}} \times ℝ | \prod_{𝑖 = 1, \dots, 𝑑_{1}} 𝑢_{𝑖}^{𝛼_{𝑖}} \geq 𝑤}

(96)

with dual cone

{𝒦︀}_{powm}^{*} (𝛼, 𝑑_{1}) = {(𝑢, 𝑤) \in ℝ_{+}^{𝑑_{1}} \times ℝ_{-} | \prod_{𝑖 = 1, \dots, 𝑑_{1}} {(\frac{𝑢_{𝑖}}{𝛼_{𝑖}})}^{𝛼_{𝑖}} \geq - 𝑤}

(97)

The relative entropy cone is defined as

{𝒦︀}_{rel} (𝑑) = cl {(𝑢, 𝑣, 𝑤) \in ℝ \times ℝ_{+ +}^{𝑑} \times ℝ_{+ +}^{𝑑} | 𝑢 \geq \sum_{𝑖 = 1, \dots, 𝑑} 𝑤_{𝑖} ln (\frac{𝑤_{𝑖}}{𝑣_{𝑖}})}

(98)

with dual cone

\begin{aligned} {𝒦︀}_{rel}^{*} (𝑑) = cl { & (𝑢, 𝑣, 𝑤) \in ℝ_{+ +} \times ℝ_{+ +}^{𝑑} \times ℝ^{𝑑} | \\ 𝑤_{𝑖} \geq 𝑢 (ln (\frac{𝑢}{𝑣_{𝑖}}) - 1) \forall 𝑖 = 1, \dots, 𝑑} \end{aligned}

(99)

7.5. Generalized power cones

As shown in [11] the Hessian of the conjugate barrier for the generalized power cones is augmented sparse with $𝑛_{𝑢} = 1$ and $𝑛_{𝑣} = 2$ , i.e. that

𝐻^{*} (𝑧) = 𝐷 + 𝑝 𝑝^{⊤} - 𝑞 𝑞^{⊤} - 𝑟 𝑟^{⊤} .

(100)

where $𝑧 ≔ (𝑢, 𝑤)$ and $𝐷 - 𝑞 𝑞^{⊤} - 𝑟 𝑟^{⊤} ≻ 0$ with

𝐷 = [\begin{matrix} ⋱ \\ \frac{𝜏_{𝑖} 𝜙}{𝑢_{𝑖} 𝜁} + \frac{1 - 𝛼_{𝑖}}{𝑢_{𝑖}^{2}} \\ ⋱ \\ \frac{2}{𝜁} 𝐼_{𝑑_{2}} \end{matrix}], 𝑝 = [\begin{matrix} 𝑝_{0} \frac{𝜏}{𝜁} \\ 𝑝_{1} \frac{𝑤}{𝜁} \end{matrix}], 𝑞 = [\begin{matrix} 𝑞_{0} \frac{𝜏}{𝜁} \\ 0 \end{matrix}], 𝑟 = [\begin{matrix} 0 \\ 𝑟_{1} \frac{𝑤}{𝜁} \end{matrix}],

(101)

where

\begin{aligned} 𝑝_{0} & = \sqrt{\frac{𝜙 (𝜙 + ‖ 𝑤 ‖^{2})}{2}}, & 𝑝_{1} = - 2 \sqrt{\frac{2 𝜙}{𝜙 + ‖ 𝑤 ‖^{2}}} \\ 𝑞_{0} & = \sqrt{\frac{𝜙 𝜁}{2}}, & 𝑟_{1} = 2 \sqrt{\frac{𝜁}{𝜙 + ‖ 𝑤 ‖^{2}}} \end{aligned}

(102)

and $𝜙 = \prod_{𝑖 \in [𝑑_{1}]} {(\frac{𝑢_{𝑖}}{𝛼_{𝑖}})}^{2 𝛼_{𝑖}}$ , $𝜏_{𝑖} = \frac{2 𝛼_{𝑖}}{𝑢_{𝑖}}$ and $𝜁 = 𝜙 - ‖ 𝑤 ‖^{2}$ .

7.6. Power mean cone

For the power mean cone it is shown in [11] that the Hessian of the conjugate barrier is augmented sparse with $𝑛_{𝑢} = 1$ and $𝑛_{𝑣} = 2$ , i.e. that

𝐻^{*} (𝑧) = 𝐷 + 𝑝 𝑝^{⊤} - 𝑞 𝑞^{⊤} - 𝑟 𝑟^{⊤} .

(103)

where $𝑧 ≔ (𝑢, 𝑤)$ and $𝐷 - 𝑞 𝑞^{⊤} - 𝑟 𝑟^{⊤} ≻ 0$ with

𝐷 = [\begin{matrix} ⋱ \\ \frac{𝜏_{𝑖} 𝜙}{𝑢_{𝑖}} + \frac{1 - 𝛼_{𝑖}}{𝑢_{𝑖}^{2}} \\ ⋱ \\ 𝜃 \end{matrix}], 𝑝 = [\begin{matrix} 𝑝_{0} 𝜏 \\ 𝑝_{1} \frac{1}{𝜁} \end{matrix}], 𝑞 = [\begin{matrix} 𝑞_{0} 𝜏 \\ 0 \end{matrix}], 𝑟 = [\begin{matrix} 0 \\ 𝑟_{1} \frac{1}{𝜁} \end{matrix}],

(104)

where

𝜃 = 1 + 𝑤^{- 2}, 𝑝_{0} = 𝜙, 𝑝_{1} = 1, 𝑞_{0} = \sqrt{𝜁 𝜙}, 𝑟_{1} = 𝜁

(105)

and

$𝜙 = \prod_{𝑖 \in [𝑑_{1}]} {(\frac{𝑢_{𝑖}}{𝛼_{𝑖}})}^{𝛼_{𝑖}}$ , $𝜁 = 𝜙 + 𝑤$ and $𝜏_{𝑖} = \frac{𝛼_{𝑖}}{𝑢_{𝑖} 𝜁}, 𝑖 = 1, \dots, 𝑑$ .

7.7. Relative entropy cone

For the relative entropy cone it is shown in [11] that the Hessian of the conjugate barrier is sparse with the following structure

𝐻^{*} (𝑧) = [\begin{matrix} \nabla_{𝑢, 𝑢}^{2} 𝑓_{*} & \nabla_{𝑢, 𝑣}^{2} 𝑓_{*} & \nabla_{𝑢, 𝑤}^{2} 𝑓_{*} \\ \nabla_{𝑣, 𝑢}^{2} 𝑓_{*} & 𝐷_{𝑣, 𝑣} & 𝐷_{𝑣, 𝑤} \\ \nabla_{𝑤, 𝑢}^{2} 𝑓_{*} & 𝐷_{𝑤, 𝑣} & 𝐷_{𝑤, 𝑤} \end{matrix}]

(106)

where $𝑧 ≔ (𝑢, 𝑣, 𝑤)$ and $𝐻^{* (𝑧)} ≻ 0$

\begin{aligned} 𝐷_{𝑣, 𝑣} & = diag ([\begin{matrix} \frac{𝑢 (𝛾_{1} + 𝑢)}{𝛾_{1}^{2} 𝑣_{1}^{2}} + \frac{1}{𝑣_{1}^{2}} & \dots & \frac{𝑢 (𝛾_{𝑑} + 𝑢)}{𝛾_{𝑑}^{2} 𝑣_{𝑑}^{2}} + \frac{1}{𝑣_{𝑑}^{2}} \end{matrix}]) \\ 𝐷_{𝑤, 𝑤} & = diag ([\begin{matrix} 𝛾_{1}^{- 2} & \dots & 𝛾_{𝑑}^{- 2} \end{matrix}]) \\ 𝐷_{𝑣, 𝑤} = 𝐷_{𝑤, 𝑣} & = diag ([\begin{matrix} {𝑢 (𝛾_{1}^{2} 𝑣_{1})}^{- 1} & \dots & {𝑢 (𝛾_{𝑑}^{2} 𝑣_{𝑑})}^{- 1} \end{matrix}]), \end{aligned}

(107)

with $𝛾_{𝑖} = 𝑤_{𝑖} - 𝑢 (ln (\frac{𝑢}{𝑣_{𝑖}}) - 1) > 0, 𝑖 = 1, \dots, 𝑑$ .

7.8. Semidefinite cones

We start by introducing the vectorization of a symmetric matrix $𝑀$ defined below

svec (𝑀) : = {\begin{aligned} [ & 𝑀_{1, 1}, \sqrt{2} 𝑀_{2, 1}, 𝑀_{2, 2}, \sqrt{2} 𝑀_{3, 1}, \dots, \\ \sqrt{2} 𝑀_{𝑛, 𝑛 - 2}, \sqrt{2} 𝑀_{𝑛, 𝑛 - 1}, 𝑀_{𝑛, 𝑛}] \end{aligned}}^{⊤} .

(108)

Note that the scaling of $\sqrt{2}$ is in order to preserve inner-products i.e. that $tr (𝐿^{⊤} 𝑀) = svec {(𝐿)}^{⊤} svec (𝑀)$ . We equivalently introduce the matrization of a vector (and denote it by $smat$ ) which unsurprisingly has the property $𝑀 = smat (svec (𝑀))$ as

smat (𝑚) : = [\begin{matrix} 𝑚_{1} & 𝑚_{2} / \sqrt{2} & \dots & 𝑚_{𝑛 (𝑛 - 1) / 2 + 1} / \sqrt{2} \\ 𝑚_{2} / \sqrt{2} & 𝑚_{3} & \dots & 𝑚_{𝑛 (𝑛 - 1) / 2 + 2} / \sqrt{2} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 𝑚_{𝑛 (𝑛 - 1) / 2 + 1} / \sqrt{2} & 𝑚_{𝑛 (𝑛 - 1) / 2 + 2} / \sqrt{2} & \dots & 𝑚_{𝑛 (𝑛 + 1) / 2} \end{matrix}] .

(109)

Note that the scaling ensures that we preserve inner-products i.e. that $𝑢^{⊤} 𝑣 = tr (smat (𝑢) smat (𝑣))$ .

7.8.1. Computing dense $𝑊_{𝑘}^{⊤} 𝑊_{𝑘}$

In CVXOPT the Nesterov-Todd scaling of a semidefinite cone is described by the product

𝑊_{𝑘} svec (𝑉) = svec (𝑅_{𝑘}^{⊤} 𝑉 𝑅_{𝑘}),

(110)

However, stating it like this gives no insights into the structure of $𝑊_{𝑘}$ , which is needed when assembling the KKT system (to be precise $𝑊_{𝑘}^{⊤} 𝑊_{𝑘}$ is needed for the KKT assembly). It turns out that the structure comes from the symmetric Kronecker product given by

(𝐺 \otimes_{𝑠} 𝐻) svec (𝑉) = \frac{1}{2} svec (𝐻 𝑉 𝐺^{⊤} + 𝐺 𝑉 𝐻^{⊤}) .

(111)

Using the above we can see that this means that the Nesterov-Todd scaling of the semidefinite cone is

(𝑅^{⊤} \otimes_{𝑠} 𝑅^{⊤}) svec (𝑉) = svec (𝑅^{⊤} 𝑉 𝑅),

(112)

from which it becomes clear that $𝑊_{𝑘} = (𝑅_{𝑘}^{⊤} \otimes_{𝑠} 𝑅_{𝑘}^{⊤})$ . We’re now almost ready to describe $𝑊_{𝑘}^{⊤} 𝑊_{𝑘}$ , but first we need to introduce the following property of the symmetric Kronecker product [12]

(𝐴 \otimes_{𝑠} 𝐵) (𝐶 \otimes_{𝑠} 𝐷) = \frac{1}{2} (𝐴 𝐶 \otimes_{𝑠} 𝐵 𝐷 + 𝐴 𝐷 \otimes_{𝑠} 𝐵 𝐶) .

(113)

Using the above it follows that

\begin{aligned} 𝑊_{𝑘}^{⊤} 𝑊_{𝑘} & = (𝑅_{𝑘} \otimes_{𝑠} 𝑅_{𝑘}) (𝑅_{𝑘}^{⊤} \otimes_{𝑠} 𝑅_{𝑘}^{⊤}) \\ = 𝑅_{𝑘} 𝑅_{𝑘}^{⊤} \otimes_{𝑠} 𝑅_{𝑘} 𝑅_{𝑘}^{⊤} \end{aligned}

(114)

Now how do we write the symmetric Kronecker product as an explicit matrix? First we need to introduce a mapping $𝑄$ that works as follows

𝑄 vec (𝑆) = svec (𝑆), 𝑄^{⊤} svec (𝑆) = vec (𝑆)

(115)

Using $𝑄$ it is now possible to describe the symmetric Kronecker product as

𝐺 \otimes_{𝑠} 𝐻 = \frac{1}{2} 𝑄 (𝐺 \otimes 𝐻 + 𝐻 \otimes 𝐺) 𝑄^{⊤}

(116)

This definition intuitively makes sense. The $𝑄^{⊤}$ on the right is applied first and lifts a $svec$ to $vec$ i.e. from a symmetric vector to a regular vector. We then apply the sum of regular Kronecker products (the sum makes it map symmetric matrices to symmetric matrices) and then pull the product down back to $svec$ using $𝑄$ . Note that in general $𝐺 \otimes_{𝑠} 𝐻$ is not symmetric. Instead, it represents a single half of an otherwise symmetric product. In the case where $𝐻 = 𝐺$ , as is the case of NT-scaling, the product is symmetric and in fact

𝑊_{𝑘}^{⊤} 𝑊_{𝑘} = 𝑅_{𝑘} 𝑅_{𝑘}^{⊤} \otimes_{𝑠} 𝑅_{𝑘} 𝑅_{𝑘}^{⊤} = 𝑄 (𝑅_{𝑘} 𝑅_{𝑘}^{⊤} \otimes 𝑅_{𝑘} 𝑅_{𝑘}^{⊤}) 𝑄^{⊤}

(117)

Cite this post

bibtex


@misc{paltorp2026conicoptimization,
  author = {Paltorp, Mikkel},
  title = {{Conic Optimization}},
  year = {2026},
  month = jan,
  url = {https://mipals.github.io/pdf/2026/01-04-conic-optimization/conic-optimization.pdf},
  note = {Blog post}
}

Bibliography

[1] P. J. Goulart and Y. Chen, “Clarabel: An interior-point solver for conic programs with quadratic objectives.” [Online]. Available: https://arxiv.org/abs/2405.12762
[2] L. Vandenberghe, “The CVXOPT linear and quadratic cone program solvers.” [Online]. Available: https://www.seas.ucla.edu/~vandenbe/publications/coneprog.pdf
[3] L. Tunçel, “Generalization of Primal-Dual Interior-Point Methods to Convex Optimization Problems in Conic Form,” Foundations of Computational Mathematics, vol. 1, no. 3, pp. 229–254, 2001, doi: 10.1007/s002080010009.
[4] J. Dahl and E. D. Andersen, “A primal-dual interior-point algorithm for nonsymmetric exponential-cone optimization,” Mathematical Programming, 2022, doi: 10.1007/s10107-021-01631-4.
[5] B. Stellato, G. Banjac, P. Goulart, A. Bemporad, and S. Boyd, “OSQP: an operator splitting solver for quadratic programs,” Mathematical Programming Computation, vol. 12, no. 4, 2020, doi: 10.1007/s12532-020-00179-2.
[6] D. Ruiz, “A scaling algorithm to equilibrate both rows and columns norms in matrices.” [Online]. Available: https://cerfacs.fr/wp-content/uploads/2017/06/14_DanielRuiz.pdf
[7] R. Takapoui and H. Javadi, “Preconditioning via Diagonal Scaling.” [Online]. Available: https://arxiv.org/abs/1610.03871
[8] T. A. Davis, “Algorithm 849: A concise sparse Cholesky factorization package,” ACM Transactions on Mathematical Software, vol. 31, no. 4, pp. 587–591, 2005, doi: 10.1145/1114268.1114277.
[9] R. J. Vanderbei, “Symmetric Quasidefinite Matrices,” SIAM Journal on Optimization, 1995, doi: 10.1137/0805005.
[10] A. Domahidi, E. Chu, and S. Boyd, “ECOS: An SOCP solver for embedded systems,” in European Control Conference (ECC), EUCA, 2013, pp. 3071–3076. doi: 10.23919/ECC.2013.6669541.
[11] Y. Chen and P. Goulart, “An Efficient IPM Implementation for A Class of Nonsymmetric Cones.” [Online]. Available: https://arxiv.org/abs/2305.12275
[12] K. Schacke, “On the Kronecker Product.” [Online]. Available: https://www.math.uwaterloo.ca/~hwolkowi/henry/reports/kronthesisschaecke04.pdf

Conic Optimization

1. Introduction

2. Primal and dual problems

2.1. The dual problem

2.2. Homogeneous embedding

3. Barrier functions

3.1. The logarithmic barrier function for the three symmetric cones

3.2. Gradients of the barrier function for the three symmetric cones

3.3. Hessians of the barrier function for the three symmetric cones

3.4. The Jordan product for the three symmetric cones

3.5. Logarithmically homogeneous self-concordant barrier

4. The central path

4.1. Scaling matrices

4.1.1. Symmetric cones

4.1.2. Nonsymmetric cones

5. Computing step directions

5.1. The affine and centering directions

5.2. Solver initialization

5.2.1. Symmetric cones

5.2.2. Nonsymmetric cones

5.3. Termination criteria

5.3.1. Feasibility checks

5.3.2. Infeasibility checks

6. Implementation details

6.1. Solving the KKT system

6.2. Regularized iterative refinement

6.3. Matrix equilibration

7. Sparse solvers

7.1. Sparsity exploitation

7.1.1. Augmented sparsity

7.2. Second-order cones

7.2.1. The Clarabel approach

7.2.2. The CVXOPT approach

7.3. Rotated second-order cones

7.4. Power cones

7.5. Generalized power cones

7.6. Power mean cone

7.7. Relative entropy cone

7.8. Semidefinite cones

7.8.1. Computing dense 𝑊𝑘⊤𝑊𝑘

Cite this post

Bibliography

7.8.1. Computing dense $𝑊_{𝑘}^{⊤} 𝑊_{𝑘}$