Sparse semidefinite programs with guaranteed near-linear time complexity via dualized clique tree conversion

Abstract

Clique tree conversion solves large-scale semidefinite programs by splitting an $n\times n$ matrix variable into up to n smaller matrix variables, each representing a principal submatrix of up to $\omega \times \omega$. Its fundamental weakness is the need to introduce overlap constraints that enforce agreement between different matrix variables, because these can result in dense coupling. In this paper, we show that by dualizing the clique tree conversion, the coupling due to the overlap constraints is guaranteed to be sparse over dense blocks, with a block sparsity pattern that coincides with the adjacency matrix of a tree. We consider two classes of semidefinite programs with favorable sparsity patterns that encompass the MAXCUT and MAX k-CUT relaxations, the Lovasz Theta problem, and the AC optimal power flow relaxation. Assuming that $\omega \ll n$, we prove that the per-iteration cost of an interior-point method is linear O(n) time and memory, so an $\epsilon$-accurate and $\epsilon$-feasible iterate is obtained after $O(\sqrt{n}\log (1/\epsilon ))$ iterations in near-linear $O(n^{1.5}\log (1/\epsilon ))$ time. We confirm our theoretical insights with numerical results on semidefinite programs as large as $n=13{,}659$.

Appendices

A Linear independence and Slater’s conditions

In this section, we prove that (CTC) inherits the assumptions of linear independence and Slater’s conditions from (SDP). We begin with two important technical lemmas.

Lemma 10

The matrix ${\mathbf {N}}$ in (13) has full row rank, that is $\det ({\mathbf {N}}{\mathbf {N}}^{T})\ne 0$.

Proof

We make ${\mathbf {N}}=[{\mathbf {N}}_{i,j}]_{i,j=1}^{\ell }$ upper-triangular by ordering its blocks topologically on T: each nonempty block row ${\mathbf {N}}_{j}$ contains a nonzero block at ${\mathbf {N}}_{j,j}$ and a nonzero block at ${\mathbf {N}}_{j,p(j)}$ where the parent node $p(j)>j$ is ordered after j. Then, the claim follows because each diagonal block ${\mathbf {N}}_{j,j}$ implements a surjection and must therefore have full row-rank. $\square$

Lemma 11

(Orthogonal complement) Let $d=\frac{1}{2}\sum _{j=1}^{\ell }|J_{j}|(|J_{j}|+1)$. Implicitly define the $d\times \frac{1}{2}n(n+1)$ matrix ${\mathbf {P}}$ to satisfy

$$\begin{aligned} {\mathbf {P}}\,\mathrm {svec}\,(X)&=[\mathrm {svec}\,(X[J_{1},J_{1}]);\ldots ;\mathrm {svec}\,(X[J_{\ell },J_{\ell }])]\qquad \forall X\in {\mathbb {S}}^{n}. \end{aligned}$$

Then, (i) ${\mathbf {N}}{\mathbf {P}}=0$; (ii) every $x\in {\mathbb {R}}^{d}$ can be decomposed as $x={\mathbf {N}}^{T}u+{\mathbf {P}}v$.

Proof

For each $x=[\mathrm {svec}\,(X_{j})]_{j=1}^{\ell }\in {\mathbb {R}}^{d}$, Theorem 4 says that there exists a Z satisfying ${\mathbf {P}}\,\mathrm {svec}\,(Z)=x$ if and only if ${\mathbf {N}}x=0$. Equivalently, $x\in {\mathrm {span}}({\mathbf {P}})$ if and only if $x\perp {\mathrm {span}}({\mathbf {N}}^{T})$. The “only if” part implies (i), while the “if” part implies (ii).

$\square$

Define the $m\times \frac{1}{2}n(n+1)$ matrix ${\mathbf {M}}$ as the vectorization of the linear constraints in (SDP), as in

$$\begin{aligned} {\mathbf {M}}\,\mathrm {svec}\,(X)=\begin{bmatrix}\mathrm {svec}\,(A_{1})^{T}\\ \vdots \\ \mathrm {svec}\,(A_{m})^{T} \end{bmatrix}\mathrm {svec}\,(X)=\begin{bmatrix}\mathrm {svec}\,(A_{1})^{T}\mathrm {svec}\,(X)\\ \vdots \\ \mathrm {svec}\,(A_{m})^{T}\mathrm {svec}\,(X) \end{bmatrix}=\begin{bmatrix}A_{1}\bullet X\\ \vdots \\ A_{m}\bullet X \end{bmatrix}. \end{aligned}$$

In reformulating (SDP) into (CTC), the splitting conditions (10) can be rewritten as the following

$$\begin{aligned} c^{T}{\mathbf {P}}&=\mathrm {svec}\,(C)^{T},&{\mathbf {A}}{\mathbf {P}}&={\mathbf {M}}, \end{aligned}$$

(38)

where $c=[\mathrm {svec}\,(C_{j})]_{j=1}^{\ell }$ and ${\mathbf {A}}=[\mathrm {svec}\,(A_{i,j})^{T}]_{i,j=1}^{m,\ell }$ are the data for the vectorized verison of (CTC).

Proof of Lemma 1

We will prove that

$$\begin{aligned} \text {exists }[u;v]\ne 0,\quad {\mathbf {A}}^{T}u+{\mathbf {N}}^{T}v=0\qquad \iff \qquad \text {exists }y\ne 0,\quad \sum _{i=1}^{m}y_{i}A_{i}=0. \end{aligned}$$

($\implies$) We must have $u\ne 0$, because ${\mathbf {N}}$ has full row rank by Lemma 10, and so ${\mathbf {A}}^{T}0+{\mathbf {N}}^{T}v=0$ if and only if $v=0$. Multiplying by ${\mathbf {P}}$ yields ${\mathbf {P}}^{T}({\mathbf {A}}^{T}u+{\mathbf {N}}^{T}v)={\mathbf {M}}^{T}u+0=0$ and so setting $y=u\ne 0$ yields ${\mathbf {M}}^{T}y=0$. () We use Lemma 11 to decompose ${\mathbf {A}}^{T}y={\mathbf {P}}z+{\mathbf {N}}^{T}v$. If $\mathrm {svec}\,(\sum _{i}y_{i}A_{i})={\mathbf {M}}^{T}y={\mathbf {P}}^{T}{\mathbf {A}}^{T}y=0,$ then ${\mathbf {P}}^{T}{\mathbf {P}}z={\mathbf {P}}^{T}{\mathbf {A}}^{T}y-{\mathbf {P}}^{T}{\mathbf {N}}^{T}v=0$ and so ${\mathbf {P}}z=0$. Setting $u=-y\ne 0$ yields ${\mathbf {A}}^{T}u+{\mathbf {N}}^{T}v=0$. $\square$

Proof of Lemma 2

We will prove that

$$\begin{aligned} \text {exists }x\in \mathrm {Int}({\mathcal {K}}),\quad \begin{bmatrix}{\mathbf {A}}\\ {\mathbf {N}}\end{bmatrix}x=\begin{bmatrix}b\\ 0 \end{bmatrix}\qquad \iff \qquad \text {exists }X\succeq 0,\quad {\mathbf {M}}\,\mathrm {svec}\,(X)=b. \end{aligned}$$

Define the chordal completion F as in (8). Observe that ${\mathbf {M}}\,\mathrm {svec}\,(X)={\mathbf {M}}\,\mathrm {svec}\,(Z)$ holds for all pairs of $P_{F}(X)=Z$, because each $A_{i}\in {\mathbb {S}}_{F}^{n}$ satisfies $A_{i}\bullet X=A_{i}\bullet P_{F}(X)$. Additionally, the positive definite version of Theorem 3 is written

$$\begin{aligned} \text {exists }X\succ 0,\quad P_{F}(X)=Z\qquad \iff \qquad {\mathbf {P}}\,\mathrm {svec}\,(Z)\in \mathrm {Int}({\mathcal {K}}). \end{aligned}$$

(39)

This result was first established by Grone et al. [52]; a succinct proof can be found in [15, Theorem 10.1]. ($\implies$) For every x satisfying ${\mathbf {N}}x=0$, there exists Z such that ${\mathbf {P}}\,\mathrm {svec}\,(Z)=x$ due to Lemma 11. If additionally $x\in \mathrm {Int}({\mathcal {K}})$, then there exists $X\succ 0$ satisfying $Z=P_{F}(X)$ due to (39). We can verify that ${\mathbf {M}}\,\mathrm {svec}\,(X)={\mathbf {M}}\,\mathrm {svec}\,(Z)={\mathbf {A}}{\mathbf {P}}\,\mathrm {svec}\,(Z)={\mathbf {A}}x=b$. () For every $X\succ 0,$ there exists Z satisfying $Z=P_{F}(X)$ and ${\mathbf {P}}\,\mathrm {svec}\,(Z)\in \mathrm {Int}({\mathcal {K}})$ due to (39). Set $u={\mathbf {P}}\,\mathrm {svec}\,(Z)$ and observe that $u\in \mathrm {Int}({\mathcal {K}})$ and ${\mathbf {N}}u={\mathbf {N}}{\mathbf {P}}\,\mathrm {svec}\,(Z)=0$. If additionally ${\mathbf {M}}\,\mathrm {svec}\,(X)=b$, then ${\mathbf {A}}u={\mathbf {A}}{\mathbf {P}}\,\mathrm {svec}\,(Z)={\mathbf {M}}\,\mathrm {svec}\,(Z)=b$. $\square$

Proof of Lemma 2

We will prove that

$$\begin{aligned} \text {exists }u,v,\quad c-{\mathbf {A}}^{T}u-{\mathbf {N}}^{T}v\in \mathrm {Int}({\mathcal {K}}_{*})\qquad \iff \qquad \text {exists }y,\quad C-\sum _{i}y_{i}A_{i}\succ 0. \end{aligned}$$

Define the chordal completion F as in (8). Theorem 3 in (39) has a dual theorem

$$\begin{aligned} \text {exists }S\succ 0,\quad S\in {\mathbb {S}}_{F}^{n}\qquad \iff \qquad \text {exists }h\in \mathrm {Int}({\mathcal {K}}_{*}),\quad \mathrm {svec}\,(S)={\mathbf {P}}^{T}h.\qquad \quad \end{aligned}$$

(40)

This result readily follows from the positive semidefinite version proved by Alger et al. [76]; see also [15, Theorem 9.2]. ($\implies$) For each $h=c-{\mathbf {A}}^{T}u-{\mathbf {N}}^{T}v$, define $S=C-\sum _{i}u_{i}A_{i}$ and observe that

$$\begin{aligned} {\mathbf {P}}^{T}h={\mathbf {P}}^{T}(c-{\mathbf {A}}^{T}u-{\mathbf {N}}^{T}v)=\mathrm {svec}\,(C)-{\mathbf {M}}^{T}u-0=\mathrm {svec}\,(S). \end{aligned}$$

If additionally $h\in {\mathcal {K}}_{*}$, then $S\succ 0$ due to (40). () For each $S=C-\sum _{i}y_{i}A_{i}\succ 0$, there exists an $h\in \mathrm {Int}({\mathcal {K}}_{*})$ satisfying $\mathrm {svec}\,(S)={\mathbf {P}}^{T}h$ due to (40). We use Lemma 11 to decompose $h={\mathbf {P}}u+{\mathbf {N}}^{T}v$. Given that $\mathrm {svec}\,(S)={\mathbf {P}}^{T}h={\mathbf {P}}^{T}{\mathbf {P}}u+0$, we must actually have ${\mathbf {P}}u=c-{\mathbf {A}}^{T}y$ since ${\mathbf {P}}^{T}(c-{\mathbf {A}}^{T}y)=\mathrm {svec}\,(C)-{\mathbf {M}}^{T}y=\mathrm {svec}\,(S)$. Hence $h=c-{\mathbf {A}}^{T}y+{\mathbf {N}}^{T}v$ and $h\in \mathrm {Int}({\mathcal {K}}_{*})$. $\square$

B Extension to inequality constraints

Consider the modifying the equality constraint in (11) into an inequality constraint, as in

$$\begin{aligned} \text {minimize }c^{T}x\quad \text {subject to }\quad {\mathbf {A}}x\le b,\quad {\mathbf {N}}x=0,\quad x\in {\mathcal {K}}. \end{aligned}$$

The corresponding dualization reads

$$\begin{aligned} \text {maximize }-c^{T}y\quad \text {subject to }\quad \begin{bmatrix}{\mathbf {A}}\\ {\mathbf {N}}\\ -I \end{bmatrix}y+s=\begin{bmatrix}b\\ 0\\ 0 \end{bmatrix},\quad s\in {\mathbb {R}}_{+}^{m}\times \{0\}^{f}\times {\mathcal {K}}, \end{aligned}$$

where m denotes the number of rows in ${\mathbf {A}}$ and f now denotes the number of rows in ${\mathbf {N}}$. Embedding the equality constraint into a second-order cone, the associated normal equation takes the form

$$\begin{aligned} \left( {\mathbf {D}}_{s}+{\mathbf {A}}^{T}{\mathbf {D}}_{l}{\mathbf {A}}+{\mathbf {N}}^{T}{\mathbf {D}}_{f}{\mathbf {N}}\right) \varDelta y=r, \end{aligned}$$

where ${\mathbf {D}}_{s}$ and ${\mathbf {D}}_{f}$ are comparable as before in (15) and (23), and ${\mathbf {D}}_{l}$ is a diagonal matrix with positive diagonal elements. This matrix has the same sparse-plus-rank-1 structure as (22), and can therefore be solved using the same rank-1 update

$$\begin{aligned} \varDelta y=({\mathbf {H}}+{\mathbf {q}}{\mathbf {q}}^{T})^{-1}r=\left( I-\frac{({\mathbf {H}}^{-1}{\mathbf {q}}){\mathbf {q}}^{T}}{1+{\mathbf {q}}^{T}({\mathbf {H}}^{-1}{\mathbf {q}})}\right) {\mathbf {H}}^{-1}r, \end{aligned}$$

where ${\mathbf {H}}$ and ${\mathbf {q}}$ now read

$$\begin{aligned} {\mathbf {H}}&={\mathbf {D}}_{s}+{\mathbf {A}}^{T}{\mathbf {D}}_{l}{\mathbf {A}}+\sigma {\mathbf {N}}^{T}{\mathbf {N}},&{\mathbf {q}}&={\mathbf {N}}^{T}w. \end{aligned}$$

The matrix ${\mathbf {H}}$ has the same block sparsity graph as the tree graph T, so we can evoke Lemma 5 to show that the cost of computing $\varDelta y$ is again $O(\omega ^{6}n)$ time and $O(\omega ^{4}n)$ memory.

C Interior-point method complexity analysis

We solve the dualized problem (21) by solving its extended homogeneous self-dual embedding

$$\begin{aligned}&\text {min.}\quad (\nu +1)\theta \end{aligned}$$

(41a)

$$\begin{aligned}&\text {s.t.} \quad \begin{bmatrix}0 &{}\quad +{\mathbf {M}}^{T} &{}\quad -{\mathbf {c}}&{}\quad -r_{d}\\ -{\mathbf {M}}&{}\quad 0 &{}\quad +{\mathbf {b}}&{}\quad -r_{p}\\ +{\mathbf {c}}^{T} &{}\quad -{\mathbf {b}}^{T} &{}\quad 0 &{}\quad -r_{c}\\ r_{d}^{T} &{}\quad r_{p}^{T} &{}\quad r_{c} &{}\quad 0 \end{bmatrix}\begin{bmatrix}x\\ y\\ \tau \\ \theta \end{bmatrix}+\begin{bmatrix}s\\ 0\\ \kappa \\ 0 \end{bmatrix} =\begin{bmatrix}0\\ 0\\ 0\\ \nu +1 \end{bmatrix} \end{aligned}$$

(41b)

$$\begin{aligned}&x,s\in {\mathcal {C}}\quad \kappa ,\tau \quad \ge 0, \end{aligned}$$

(41c)

where the data is given in standard form

$$\begin{aligned} {\mathbf {M}}=\begin{bmatrix}0&-\begin{bmatrix}{\mathbf {A}}^{T}&{\mathbf {N}}^{T}\end{bmatrix}&+I\end{bmatrix},\quad {\mathbf {c}}^{T}=\begin{bmatrix}0&\begin{bmatrix}b^{T}&0\end{bmatrix}&0\end{bmatrix},\quad {\mathbf {b}}=-c,\end{aligned}$$

(41d)

$$\begin{aligned} {\mathcal {C}}=\{(x_{0},x_{1}\}\in {\mathbb {R}}^{p+1}:\Vert x_{1}\Vert \le x_{0}\}\times {\mathcal {K}}, \end{aligned}$$

(41e)

and the residual vectors are defined

$$\begin{aligned} r_{d}={\mathbf {1}}_{{\mathcal {C}}}-{\mathbf {c}},\quad r_{p}=-{\mathbf {M}}{\mathbf {1}}_{{\mathcal {C}}}+{\mathbf {b}},\quad r_{c}=1+{\mathbf {c}}^{T}{\mathbf {1}}_{{\mathcal {C}}}. \end{aligned}$$

(41f)

Here, $\nu$ is the order of the cone ${\mathcal {C}}$, and ${\mathbf {1}}_{{\mathcal {C}}}$ is its identity element

$$\begin{aligned} \nu&=1+|J_{1}|+\cdots +|J_{\ell }|,&{\mathbf {1}}_{{\mathcal {C}}}&=[1;0;\ldots ;0;\mathrm {svec}\,(I_{|J_{1}|});\ldots ;\mathrm {svec}\,(I_{|J_{\ell }|})]. \end{aligned}$$

(42)

Problem (41) has optimal value $\theta ^{\star }=0$. Under the primal-dual Slater’s conditions (Assumption 2), an interior-point method is guaranteed to converge to an $\epsilon$-accurate solution with $\tau >0$, and this yields an $\epsilon$-feasible and $\epsilon$-accurate solution to the dualized problem (21) by rescaling $x/\tau$ and $y=y/\tau$ and $s=s/\tau$. The following result is adapted from [38, Lemma 5.7.2] and [77, Theorem 22.7].

Lemma 12

($\epsilon$-accurate and $\epsilon$-feasible) If $(x,y,s,\tau ,\theta ,\kappa )$ satisfies (41b) and (41c) and

$$\begin{aligned} \mu =\frac{x^{T}s+\tau \kappa }{\nu +1}&\le \epsilon ,&\tau \kappa \ge&\gamma \mu \end{aligned}$$

for constants $\epsilon ,\gamma >0$, then the rescaled point $(x/\tau ,y/\tau ,s/\tau )$ satisfies

$$\begin{aligned} \Vert {\mathbf {M}}(x/\tau )-{\mathbf {b}}\Vert&\le K\epsilon ,&\Vert {\mathbf {M}}^{T}(y/\tau )+(s/\tau )-{\mathbf {c}}\Vert&\le K\epsilon ,&\frac{(x/\tau )^{T}(s/\tau )}{\nu +1}\le K\epsilon \end{aligned}$$

where K is a constant.

Proof

Note that (41b) implies $\mu =\theta$ and

$$\begin{aligned} \Vert {\mathbf {M}}(x/\tau )-{\mathbf {b}}\Vert&=\frac{\Vert r_{p}\Vert \mu }{\tau },&\Vert {\mathbf {M}}^{T}(y/\tau )+(s/\tau )-{\mathbf {c}}\Vert&=\frac{\Vert r_{d}\Vert \mu }{\tau },&\\ \frac{(x/\tau )^{T}(s/\tau )}{\nu +1}&=\frac{\mu }{\tau ^{2}}. \end{aligned}$$

Hence, we obtain our desired result by upper-bounding $1/\tau$. Let $(x^{\star },y^{\star },s^{\star },\tau ^{\star },\theta ^{\star },\kappa ^{\star })$ be a solution of (41), and note that for every $(x,y,s,\tau ,\theta ,\kappa )$ satisfying (41b) and (41c), we have the following via the skew-symmetry of (41b)

$$\begin{aligned} (x-x^{\star })^{T}(s-s^{\star })+(\tau -\tau ^{\star })(\kappa -\kappa ^{\star })=0. \end{aligned}$$

Rearranging yields

$$\begin{aligned} (\nu +1)\mu =x^{T}s+\tau \kappa =(x^{\star })^{T}s +x^{T}(s^{\star })+\tau \kappa ^{\star }+\tau ^{\star }\kappa \ge \tau ^{\star }\kappa \end{aligned}$$

and hence

$$\begin{aligned} \kappa \tau ^{\star }\le \mu (\nu +1),\quad \tau \kappa \ge \gamma \mu \qquad \implies \qquad \tau \ge \frac{\gamma }{\nu +1}\tau ^{\star }. \end{aligned}$$

If (SDP) satisfies the primal-dual Slater’s condition, then (CTC) also satisfies the primal-dual Slater’s condition (Lemmas 2& 3). Therefore, the vectorized version (11) of (CTC) attains a solution $({\hat{x}},{\hat{y}},{\hat{s}})$ with ${\hat{x}}^{T}{\hat{s}}=0$, and the following

$$\begin{aligned} x^{\star }=\tau ^{\star }\begin{bmatrix}\Vert {\hat{y}}\Vert \\ {\hat{y}}\\ {\hat{s}} \end{bmatrix},\quad y^{\star }=\tau ^{\star }{\hat{x}},\quad s^{\star }=\tau ^{\star }\begin{bmatrix}0\\ 0\\ {\hat{x}} \end{bmatrix},\quad \tau ^{\star }=\frac{\nu +1}{\Vert {\hat{y}}\Vert +{\mathbf {1}}_{{\mathcal {K}}}^{T}{\hat{s}}+{\mathbf {1}}_{{\mathcal {K}}}^{T}{\hat{x}}+1}, \end{aligned}$$

with $\theta ^{\star }=\kappa ^{\star }=0$ is a solution to (41). This proves the following upper-bound

$$\begin{aligned} \frac{1}{\tau }\le K_{\tau }=\frac{1}{\gamma }\cdot \min _{{\hat{x}},{\hat{y}},{\hat{s}}}\{\Vert {\hat{y}}\Vert +{\mathbf {1}}_{{\mathcal {K}}}^{T}{\hat{s}}+{\mathbf {1}}_{{\mathcal {K}}}^{T}{\hat{x}}+1:({\hat{x}},{\hat{y}},{\hat{s}}) \text { solves }(11)\text { with }{\hat{x}}^{T}{\hat{s}}=0\}. \end{aligned}$$

Setting $K=\max \{\Vert r_{p}\Vert K_{\tau },\Vert r_{d}\Vert K_{\tau },K_{\tau }^{2}\}$ yields our desired result. $\square$

We solve the homogeneous self-dual embedding (41) using the short-step method of Nesterov and Todd [78, Algorithm 6.1] (and also Sturm and S. Zhang [79, Section 5.1]), noting that SeDuMi reduces to it in the worst case; see [61] and [80]. Beginning at the following strictly feasible, perfectly centered point

$$\begin{aligned} \theta ^{(0)}=\tau ^{(0)}=\kappa ^{(0)}=1,\quad y^{(0)}=0,\quad x^{(0)}=s^{(0)}={\mathbf {1}}_{{\mathcal {C}}}, \end{aligned}$$

(43)

with barrier parameter $\mu =1$, we take the following steps

$$\begin{aligned} \mu ^{+}&=\left( 1-\frac{1}{15\sqrt{\nu +1}}\right) \cdot \frac{x^{T}s+\tau \kappa }{\nu +1},\nonumber \\ (x^{+},y^{+},s^{+},\tau ^{+},\theta ^{+},\kappa ^{+})&=(x,y,s,\tau ,\theta ,\kappa )+(\varDelta x,\varDelta y,\varDelta s,\varDelta \tau ,\varDelta \theta ,\varDelta \kappa ). \end{aligned}$$

(44a)

along the search direction defined by the linear system [81, Eqn. 9]

$$\begin{aligned} \begin{bmatrix}0 &{}\quad +{\mathbf {M}}^{T} &{}\quad -{\mathbf {c}}&{}\quad -r_{p}\\ -{\mathbf {M}}&{}\quad 0 &{}\quad +{\mathbf {b}}&{} \quad -r_{d}\\ +{\mathbf {c}}^{T} &{}\quad -{\mathbf {b}}^{T} &{}\quad 0 &{}\quad -r_{c}\\ +r_{p}^{T} &{}\quad +r_{d}^{T} &{}\quad +r_{c} &{}\quad 0 \end{bmatrix}\begin{bmatrix}\varDelta x\\ \varDelta y\\ \varDelta \tau \\ \varDelta \theta \end{bmatrix}+\begin{bmatrix}\varDelta s\\ 0\\ \varDelta \kappa \\ 0 \end{bmatrix}&=0,\end{aligned}$$

(45a)

$$\begin{aligned} s+\varDelta s+\nabla ^{2}F(w)\varDelta x+\mu ^{+}\nabla F(x)&=0,\end{aligned}$$

(45b)

$$\begin{aligned} \kappa +\varDelta \kappa +(\kappa /\tau )\varDelta \tau -\mu ^{+}\tau ^{-1}&=0. \end{aligned}$$

(45c)

Here, F is the usual self-concordant barrier function on ${\mathcal {C}}$

$$\begin{aligned} F([x_{0};x_{1};\mathrm {svec}\,(X_{1});\ldots ;\mathrm {svec}\,(X_{\ell })])=&-\log \left( \frac{1}{2}x_{0}^{2}-\frac{1}{2}x_{1}^{T}x_{1}\right) -\sum _{j=1}^{\ell }\log \det (X_{j}) \end{aligned}$$

(46)

and $w\in \mathrm {Int}({\mathcal {C}})$ is the unique scaling point satisfying $\nabla ^{2}F(w)x=s$, which can be computed from x and s in closed-form. The following iteration bound is an immediate consequence of [78, Theorem 6.4]; see also [79, Theorem 5.1].

Lemma 13

(Short-Step Method) The sequence in (44) arrives at an iterate $(x,y,s,\tau ,\theta ,\kappa )$ satisfying the conditions of Lemma 12 with $\gamma =9/10$ in at most $O(\sqrt{\nu }\log (1/\epsilon ))$ iterations.

The cost of each interior-point iteration is dominated by the cost of computing the search direction in (45). Using elementary but tedious linear algebra, we can show that if

$$\begin{aligned} ({\mathbf {M}}{\mathbf {D}}^{-1}{\mathbf {M}}^{T})\begin{bmatrix}v_{1}&v_{2}&v_{3}\end{bmatrix}=\begin{bmatrix}0&-{\mathbf {b}}&r_{p}\end{bmatrix}-{\mathbf {M}}{\mathbf {D}}^{-1}\begin{bmatrix}d&{\mathbf {c}}&r_{d}\end{bmatrix} \end{aligned}$$

(47a)

where ${\mathbf {D}}=\nabla ^{2}F(w)$ and $d=-s-\mu ^{+}\nabla F(x)$, and

$$\begin{aligned} \begin{bmatrix}u_{1}&u_{2}&u_{3}\end{bmatrix}={\mathbf {D}}^{-1}(\begin{bmatrix}d&c&r_{d}\end{bmatrix}+{\mathbf {M}}^{T}\begin{bmatrix}v_{1}&v_{2}&v_{3}\end{bmatrix}), \end{aligned}$$

(47b)

then

$$\begin{aligned} \left( \begin{bmatrix}-{\mathbf {D}}_{0} &{} -r_{c}\\ r_{c} &{} 0 \end{bmatrix}-\begin{bmatrix}{\mathbf {c}}&{} r_{d}\\ -{\mathbf {b}}&{} r_{p} \end{bmatrix}^{T}\begin{bmatrix}u_{2} &{} u_{3}\\ v_{2} &{} v_{3} \end{bmatrix}\right) \begin{bmatrix}\varDelta \tau \\ \varDelta \theta \end{bmatrix}&=\begin{bmatrix}-d_{0}\\ 0 \end{bmatrix}-\begin{bmatrix}{\mathbf {c}}&{} r_{d}\\ -{\mathbf {b}}&{} r_{p} \end{bmatrix}^{T}\begin{bmatrix}u_{1}\\ v_{1} \end{bmatrix},\end{aligned}$$

(47c)

$$\begin{aligned} \begin{bmatrix}\varDelta x\\ \varDelta y \end{bmatrix}&=\begin{bmatrix}u_{1}\\ v_{1} \end{bmatrix}-\begin{bmatrix}u_{1} &{} u_{2}\\ v_{1} &{} v_{2} \end{bmatrix}\begin{bmatrix}\varDelta \tau \\ \varDelta \theta \end{bmatrix},\end{aligned}$$

(47d)

$$\begin{aligned} \varDelta s&=d-{\mathbf {D}}\varDelta x,\end{aligned}$$

(47e)

$$\begin{aligned} \varDelta \kappa&=d_{0}-{\mathbf {D}}_{0}\varDelta \tau , \end{aligned}$$

(47f)

where ${\mathbf {D}}_{0}=\kappa /\tau$ and $d_{0}=-\kappa +\mu ^{+}\tau ^{-1}$. Hence, the cost of computing the search direction is dominated by the cost of solving the normal equation for three different right-hand sides. Here, the normal matrix is written

$$\begin{aligned} {\mathbf {M}}{\mathbf {D}}^{-1}{\mathbf {M}}^{T}=\mathrm {diag}(W_{1}\otimes _{s}W_{1},\ldots ,W_{\ell } \otimes _{s}W_{\ell })+\begin{bmatrix}{\mathbf {A}}\\ {\mathbf {N}}\end{bmatrix}^{T}(w_{1}w_{1}^{T}+\sigma I)\begin{bmatrix}{\mathbf {A}}\\ {\mathbf {N}}\end{bmatrix}, \end{aligned}$$

where $\sigma =\frac{1}{2}(w_{0}^{2}-w_{1}^{T}w_{1})>0$ and $\otimes _{s}$ denotes the symmetric Kronecker product [82] implicitly defined to satisfy

$$\begin{aligned} (A\otimes _{s}B)\mathrm {svec}\,(X)=\frac{1}{2}\mathrm {svec}\,(AXB^{T}+BXA^{T})\qquad \text {for all }X=X^{T}. \end{aligned}$$

Under the hypothesis on ${\mathbf {A}}$ stated in Theorem 5, the normal matrix satisfies the assumptions of Lemma 5, and can therefore be solved in linear O(n) time and memory.

Proof of Theorem 5

Combining Lemmas 12 and 13 shows that the desired $\epsilon$-accurate, $\epsilon$-feasible iterate is obtained after $O(\sqrt{\nu }\log (1/\epsilon ))$ interior-point iterations. At each iteration we perform the following steps: (1) compute the scaling point w; (2) solve the normal equation (47a) for three right-hand sides; (3) back-substitute (47b)–(47f) for the search direction and take the step in (44). Note from the proof of Lemma 5 that the matrix $[{\mathbf {A}};{\mathbf {N}}]$ has at most $O(\omega ^{2}n)$ rows under Assumption 1, and therefore $\mathrm {nnz}\,({\mathbf {M}})=O(\omega ^{4}n)$ under the hypothesis of Theorem 5. Below, we show that the cost of each step is bounded by $O(\omega ^{6}n)$ time and $O(\omega ^{4}n)$ memory.

Scaling point We partition $x=[x_{0};x_{1};\mathrm {svec}\,(X_{1});\ldots ;\mathrm {svec}\,(X_{\ell })]$ and similarly for s. Then, the scaling point w is given in closed-form [61, Section 5]

$$\begin{aligned}&\begin{bmatrix}u_{0}\\ u_{1} \end{bmatrix}=2^{-3/4}\begin{bmatrix}1 &{} 1\\ -s_{1}/\Vert s_{1}\Vert &{} s_{1}/\Vert s_{1}\Vert \end{bmatrix}\begin{bmatrix}(s_{0}-\Vert s_{1}\Vert )^{1/2}\\ (s_{0}+\Vert s_{1}\Vert )^{1/2} \end{bmatrix},\quad \begin{bmatrix}v_{0}\\ v_{1} \end{bmatrix}=\begin{bmatrix}u_{0} &{} u_{1}^{T}\\ u_{1} &{} \frac{1}{2}(u_{0}^{2}-u_{1}^{T}u_{1})I \end{bmatrix}\begin{bmatrix}x_{0}\\ x_{1} \end{bmatrix},\\&\begin{bmatrix}w_{0}\\ w_{1} \end{bmatrix}=2^{-1/4}\begin{bmatrix}u_{0} &{} u_{1}^{T}\\ u_{1} &{} \frac{1}{2}(u_{0}^{2}-u_{1}^{T}u_{1})I \end{bmatrix}\begin{bmatrix}1 &{} 1\\ -v_{1}/\Vert v_{1}\Vert &{} v_{1}/\Vert v_{1}\Vert \end{bmatrix}\begin{bmatrix}(v_{0}-\Vert v_{1}\Vert )^{-1/2}\\ (v_{0}+\Vert v_{1}\Vert )^{-1/2} \end{bmatrix},\\&W_{j}=S_{j}^{1/2}(S_{j}^{1/2}X_{j}S_{j}^{1/2})^{-1/2}S_{j}^{1/2}\qquad \text {for all }j\in \{1,\ldots ,\ell \}. \end{aligned}$$

Noting that $\mathrm {nnz}\,(w_{1})\le O(\omega ^{2}n)$, $\ell \le n$ and each $W_{j}$ is at most $\omega \times \omega$, the cost of forming $w=[w_{0};w_{1};\mathrm {svec}\,(W_{1});\ldots ;\mathrm {svec}\,(W_{\ell })]$ is at most $O(\omega ^{3}n)$ time and $O(\omega ^{2}n)$ memory. Also, since

$$\begin{aligned} {\mathbf {D}}=\nabla ^{2}F(w)=\mathrm {diag}\left( \begin{bmatrix}w_{0} &{}\quad w_{1}^{T}\\ w_{1} &{}\quad \frac{1}{2}(w_{0}^{2}-w_{1}^{T}w_{1})I \end{bmatrix},W_{1}\otimes _{s}W_{1},\ldots ,W_{\ell }\otimes _{s}W_{\ell }\right) ^{-1}, \end{aligned}$$

the cost of each matrix-vector product with ${\mathbf {D}}$ and ${\mathbf {D}}^{-1}$ is also $O(\omega ^{3}n)$ time and $O(\omega ^{2}n)$ memory.

Normal equation The cost of matrix-vector products with ${\mathbf {M}}$ and ${\mathbf {M}}^{T}$ is $\mathrm {nnz}\,({\mathbf {M}})=O(\omega ^{4}n)$ time and memory. Using Lemma 5, we form the right-hand sides and solve the three normal equations in (47a) in $O(\omega ^{6}n)$ time and $O(\omega ^{4}n)$ memory.

Back-substitution The cost of back substituting (47b)–(47f) and making the step (44) is dominated by matrix-vector products with ${\mathbf {D}}$, ${\mathbf {D}}^{-1}$, ${\mathbf {M}}$, and ${\mathbf {M}}^{T}$ at $O(\omega ^{4}n)$ time and memory.$\square$