A Faster Interior Point Method for Semidefinite Programming

Introduction

Semidefinite programs (SDPs) constitute a class of convex optimization problems that optimize a linear objective over the intersection of the cone of positive semidefinite matrices with an affine space. SDPs generalize linear programs and have a plethora of applications in operations research, control theory, and theoretical computer science [VB96]. Applications in theoretical computer science include improved approximation algorithms for fundamental problems (e.g., Max-Cut [GW95], coloring 3-colorable graphs [KMS94], and sparsest cut [ARV09]), quantum complexity theory [JJUW11], robust learning and estimation [CG18, CDG19, CDGW19], and algorithmic discrepancy and rounding [BDG16, BG17, Ban19]. We formally define SDPs with variable size $n\times n$ and $m$ constraints:

where $\langle A,B\rangle:=\sum_{i,j}A_{i,j}B_{i,j}$ is the trace product.

Two prominent methods for solving SDPs, with runtimes depending logarithmically on the accuracy parameter $\epsilon$, are the cutting plane method and the interior point method.

The cutting plane method maintains a convex set containing the optimal solution. In each iteration, the algorithm queries a separation oracle, which returns a hyperplane that divides the convex set into two subsets. The convex set is then updated to contain the subset with the optimal solution. This process is repeated until the volume of the maintained set becomes small enough and a near-optimal solution can be found. Since Khachiyan proved [Kha80] that the ellipsoid method solves linear programs in polynomial time, cutting plane methods have played a crucial role in both discrete and continuous optimization [GLS81, GV02].

In contrast, interior point methods add a barrier function to the objective and, by adjusting the weight of this barrier function, solve a different optimization problem in each iteration. The solutions to these successive problems form a well-defined central path. Since Karmarkar proved [Kar84] that interior point methods can solve linear programs in polynomial time, these methods have become an active research area. Their number of iterations is usually the square root of the number of dimensions, as opposed to the linear dependence on dimensions in cutting plane methods.

Since cutting plane methods use less structural information than interior point methods, they are slower at solving almost all problems where interior point methods are known to apply. However, SDPs remain one of the most fundamental optimization problems where the state of the art is, in fact, the opposite: the current fastest cutting plane methods[JLSW20] improves upon the runtime of [LSW15] in terms of the dependence on $\log(n/\epsilon)$, while the polynomial factors are the same in both runtimes. of [LSW15, JLSW20] solve a general SDP in time $m(mn^{2}+m^{2}+n^{\omega})$, while the fastest SDP solvers based on interior point methods in the work of [NN92] and [Ans00] achieve runtimes of $\sqrt{n}(m^{2}n^{2}+mn^{\omega}+m^{\omega})$ and $(mn)^{1/4}(m^{4}n^{2}+m^{3}n^{\omega})$, respectively, which are slower in the most common regime of $m\in[n,n^{2}]$ (see Table 1.2). This apparent paradox raises the following natural question:

How fast can SDPs be solved using interior point methods?

1 Our results

We present a faster interior point method for solving SDPs. Our main result is the following theorem, the formal version of which is given in Theorem 4.1.

There is an interior point method that solves a general SDP with variable size $n\times n$ and $m$ constraints in timeWe use $O^{*}$ to hide $n^{o(1)}$ and $\log^{O(1)}(n/\epsilon)$ factors and $\widetilde{O}$ to hide $\log^{O(1)}(n/\epsilon)$ factors, where $\epsilon$ is the accuracy parameter. $O^{*}(\sqrt{n}(mn^{2}+m^{\omega}+n^{\omega}))$.

Our runtime can be roughly interpreted as follows:

$\sqrt{n}$ is the iteration complexity of the interior point method with the log barrier function.

$m^{\omega}$ is the cost of inverting the Hessian of the log barrier.

$n^{\omega}$ is the cost of inverting the slack matrix.

Thus, the terms in the runtime of our algorithm arise as a natural barrier to further speeding up SDP solvers. See Section 1.2.2, 1.2.3, and 1.2.4 for more detail.

Table 1.1 compares our result with previous SDP solvers. The first takeaway of this table and Theorem 1.2 is that our interior point method always runs faster than that in [NN92] and is faster than that in [NN94] and [Ans00] when $m\geq n^{1/13}$. A second consequence is that whenever $m\geq\sqrt{n}$, our interior point method is faster than the current fastest cutting plane method [LSW15, JLSW20]. We note that $n\leq m\leq n^{2}$ is satisfied in most SDP applications known to us, such as classical combinatorial optimization problems over graphs, experiment design problems in statistics and machine learning, and sum-of-squares problems. An explicit comparison to previous algorithms in the cases of $m=n$ and $m=n^{2}$ is shown in Table 1.2.

2 Technique overview

By removing redundant constraints, we can, without loss of generality, assume $m\leq n^{2}$ in the primal formulation of the SDP (1). Thereafter, instead of solving the primal SDP, which has variable size $n\times n$, we solve its dual formulation, which has dimension $m\leq n^{2}$:

Interior point methods solve (2) by minimizing the penalized objective function:

Nesterov and Nemirovski [NN92] use the log barrier function,

where $g(y)$ is the log barrier function defined in (4). They proved that choosing $\phi(y)=\sqrt{n}V(y)$ in (3) makes the interior point method converge in $\widetilde{O}(\sqrt{m}n^{1/4})$ iterations, which is smaller than the $\widetilde{O}(\sqrt{n})$ iteration complexity of [NN92] when $m\leq\sqrt{n}$. They also studied the combined volumetric-logarithmic barrier

and showed that taking $\phi(y)=\sqrt{n/m}\cdot V_{\rho}(y)$ for $\rho=(m-1)/(n-1)$ yields an iteration complexity of $\widetilde{O}((mn)^{1/4})$. when $m\leq n$, this iteration complexity is lower than $\widetilde{O}(\sqrt{n})$ of [NN92]. We refer readers to the much simpler proofs in [Ans00] for these results.

However, the volumetric barrier (and thus the combined volumetric-logarithmic barrier) leads to complicated expressions for the gradient and Hessian that make each iteration costly. For instance, the Hessian of the volumetric barrier is

where $Q(y)$, $R(y)$, and $T(y)$ are $m\times m$ matrices such that for each $(j,k)\in[m]\times[m]$,

where $\otimes$ is the Kronecker product (see Section 2.1for formal definition). Due to the complicated formulas in (1.2.1), efficient computation of Newton step in each iteration of the interior point method is difficult; in fact, each iteration runs slower than the Nesterov-Nemirovski interior point method by a factor of $m^{2}$. Since most applications of SDPs known to us have the number of constraints $m$ be at least linear in $n$, the total runtime of interior point methods based on the volumetric barrier and the combined volumetric-logarithmic barrier is inevitably slow.

2.2 Our techniques

Given the inefficiency of implementing the volumetric and volumetric-logarithmic barriers discussed above, this paper uses the log barrier in (4). We now describe some of our key techniques that improve the runtime of the Nesterov-Nemirovski interior point method [NN92].

As noted in Section 1.2.1, the runtime bottleneck in [NN92] is computing the inverse of the Hessian of the log barrier function, where the Hessian is described in (5). In [NN92], each of these $m^{2}$ entries is computed separately, resulting in a runtime of $O(m^{2}n^{2})$ per iteration.

Instead contrast, we show below how to group these computations using rectangular matrix multiplication. The expression from (5) can be re-written as

Thus far, we have reduced the per iteration cost of $O^{*}(m^{2}n^{2}+mn^{\omega})$ for Hessian computation down to

The fast rectangular matrix multiplication approach noted above, however, is still not very efficient, because the Hessian must be computed from scratch in each iteration of the interior point method. If there are $T$ iterations in total, it then takes time

Specifically, first we prove that whenever $S$ is updated in an iteration, the potential function increases by at most $\widetilde{O}(1)$ (see Lemma 6.2). The proof of this statement crucially uses the structural property of interior point method that slack matrices in consecutive steps are sufficiently close to each other. Formally, for any iteration $t\in[T]$, we show in Theorem 5.1 that the consecutive slack matrices $S_{t}$ and $S_{t+1}$ satisfy

Given the low-rank update on $\widetilde{S}$ described above, we show how to efficiently update the approximate Hessian $\widetilde{H}$, defined as

for each entry $(j,k)\in[m]\times[m]$. The approximate slack matrix $\widetilde{S}$ being a spectral approximation of the true slack matrix $S$ implies that the approximate Hessian $\widetilde{H}$ is also a spectral approximation of the true Hessian $H$ (see Lemma 5.3). This approximate Hessian therefore suffices for our algorithm to approximately follow the central path.

To efficiently update the approximate Hessian $\widetilde{H}$ in (10), we notice that a rank-$r$ update on $\widetilde{S}$ implies a rank-$r$ update on $\widetilde{S}^{-1}$ via the Woodbury matrix identity (see Fact 2.4). The change in $\widetilde{S}^{-1}$ can be expressed as

where $r_{t}$ is the rank of the update on $\widetilde{S}_{t}$. Applying Theorem 1.4 with several properties of fast rectangular matrix multiplication that we prove in Section 3 , we upper bound the runtime in (12) by

which implies Theorem 1.2. In Section 1.2.3 and 1.2.4, we discuss bottlenecks to further improving our runtime.

2.3 Bottlenecks of our interior point method

In most cases, the costliest term in our runtime is the per iteration cost of $mn^{2}$, which corresponds to reading the entire input in each iteration. Our subsequent discussions therefore focus on the steps in our algorithm that require at least $mn^{2}$ time per iteration.

When $y$ is updated in each iteration of our interior point method, we need to compute the true slack matrix $S$ as

Computing $S$ is needed to update the approximate slack matrix $\widetilde{S}$ so that $\widetilde{S}$ remains a spectral approximation to $S$. As $S$ might suffer from full-rank changes, it naturally requires $mn^{2}$ time to compute in each iteration. This is the first appearance of the $mn^{2}$ cost per iteration.

Recall from (3) that our interior point method follows the central path defined via the penalized objective function

for a parameter $\eta>0$ and $\phi(y)=-\log\det S$. In each iteration, to perform the Newton step, the gradient of the penalized objective is computed as

for each coordinate $j\in[m]$. Even if we are given $S^{-1}$, it still requires $mn^{2}$ time to compute (13) for all $j\in[m]$. This is the second appearance of the per iteration cost of $mn^{2}$.

Recall from Section 1.2.2 that updating the approximate slack matrix $S$ by rank $r$ means the time needed to update the approximate Hessian is dominated by computing the term

2.4 LP techniques unlikely to improve SDP runtime

The preceeding discussion of bottlenecks suggests that reading the entire input in each iteration, which takes $mn^{2}$ time per iteration, stands as a natural barrier to further improving the runtime of SDP solvers based on interior point methods.

In the context of linear programming (LP), several recent results [CLS19, BLSS20] yield faster interior point methods that bypass reading the entire input in every iteration. Two techniques crucial to these results are: (1) showing that the Hessian (projection matrix) admits low-rank updates, and (2) speeding computation of the Hessian via sampling.

We now describe these techniques in the context of SDP and argue that they are unlikely to improve our runtime.

We saw in Section 1.2.2 that constructing an approximate slack matrix $\widetilde{S}$ that admits low-rank updates in each iterations leveraged the fact that the true slack matrix $S$ changes “slowly” throughout our interior point method as described in (9). One natural question that follows is whether a similar upper bound can be obtained for the Hessian. If such a result could be proved, then one could maintain an approximate Hessian that admitted low-rank updates, which would speed up the approximate Hessian computation. Indeed, in the context of LP, such a bound for the Hessian can be proved (e.g., [BLSS20, Lemma 47]).

Unfortunately, it is impossible to prove such a statement for the Hessian in the context of SDP. To show this, it is convenient to express the Hessian using the Kronecker product (Section 2.1)as

This large change indicates that we are unlikely to obtain an approximation to the Hessian that admits low-rank updates, which is a key difference between LP and SDP.

Recall from (8) that the Hessian can be computed as

For SDP, however, sampling is unlikely to speed up the Hessian computation. In general, we must sample at least $m$ columns (i.e. $|L|\geq m$) of $\mathcal{B}$ to spectrally approximate $H$ or the computed matrix will not be full rank. However, this requires computing the entries of $S^{-1/2}A_{j}S^{-1/2}$ that correspond to $L\subseteq[n^{2}]$ for all $j\in[m]$, which requires reading all $A_{j}$’s and thus still takes $O(mn^{2})$ time.

3 Related work

Linear Programming is a class of fundamental problems in convex optimization. There is a long list of work focused on fast algorithms for linear programming [Dan47, Kha80, Kar84, Vai87, Vai89b, LS14, LS15, Sid15, Lee16, CLS19, Bra20, BLSS20].

Cutting plane method is a class of optimization methods that iteratively refine a convex set that contains the optimal solution by querying a separation oracle. Since its introduction in the 1950s, there has been a long line of work on obtaining fast cutting plane methods [Sho77, YN76, Kha80, KTE88, NN89, Vai89a, AV95, BV02, LSW15, JLSW20].

As the focus of this paper, cutting plane methods and interior point methods solve SDPs in time that depends logarithmically on $1/\epsilon$, where $\epsilon$ is the accuracy parameter. A third class of algorithms, the first-order methods, solve SDPs at runtimes that depend polynomially on $1/\epsilon$. While having worse dependence on $1/\epsilon$ compared to IPM and CPM, these first-order algorithms usually have better dependence on the dimension. There is a long list of work on first-order methods for general SDP or special classes of SDP (e.g. Max-Cut SDP [AK07, GH16, AZL17, CDST19, LP20, YTF+19], positive SDPs [JY11, PT12, ALO16, JLL+20].)

Preliminaries

2 Useful facts

Given a symmetric matrix $B$, a positive semi-definite matrix $A$ and $\alpha\in$, we have

Matrix Multiplication

The main goal of this section is to derive upper bounds on the time to perform the following two rectangular matrix multiplication tasks (Lemma 3.9, 3.10, and 3.11):

Multiplying a matrix of dimensions $m\times n^{2}$ with one of dimensions $n^{2}\times m$,

Multiplying a matrix of dimensions $n\times mn$ with one of dimensions $mn\times n$.

Besides being crucial to the runtime analysis of our interior point method in Section 7, these results (as well as several intermediate results) might be of independent interest.

We need the following definitions to describe the cost of certain fundamental matrix operations we use.

For any three positive integers $n,m,r$, we have

We refer to Table 3 in [GU18] for the latest upper bounds on $\omega(k)$ for different values of $k$. In particular, we need the following upper bounds in our paper.

2 Technical results for matrix multiplication

We assume that $n^{p}$ and $n^{q}$ are integers for notational simplicity. Consider multiplying an $n\times n^{p}$ matrix with an $n^{p}\times n$ matrix. One can cut the $n\times n^{p}$ matrix into $n^{p-q}$ rectangular blocks of size $n\times n^{q}$ and the $n^{p}\times n$ matrix into $n^{p-q}$ rectangular blocks of size $n^{q}\times n$, and compute the multiplication of the corresponding blocks. This approach takes time $n^{p-q+\omega(q)+o(1)}$, from which the desired inequality immediately follows. ∎

Key to our analysis is the following lemma, which establishes the convexity of $\omega(k)$.

The fast rectangular matrix multiplication time exponent $\omega(k)$ as defined in Definition 3.2 is convex in $k$.

Let $k=\alpha\cdot p+(1-\alpha)\cdot q$ for $\alpha\in(0,1)$. For notational simplicity, we assume that $n^{p}$, $n^{q}$ and $n^{k}$ are all integers. Consider a rectangular matrix of dimensions $n\times n^{k}$. Since $\alpha p\leq k$, we can tile this rectangular matrix with matrices of dimensions $n^{\alpha}\times n^{\alpha p}$. Then, the product of this tiled matrix with another similarly tiled matrix of dimensions $n^{k}\times n$ can be obtained by viewing it as a multiplication of a matrix of dimensions $n/n^{\alpha}\times n^{k}/n^{\alpha p}$ with one of dimensions $n^{k}/n^{\alpha p}\times n^{1/\alpha}$, where each “element” of these two matrices is itself a matrix of dimensions $n^{\alpha}\times n^{\alpha p}$. With this recursion in tow, we obtain the following upper bound.

The final step above follows from denoting $m=n^{\alpha}$ and observing that multiplying matrices of dimensions $n^{\alpha}\times n^{\alpha\cdot p}$ costs, by Definition 3.2, $m^{\omega(p)+o(1)}$, which is exactly $n^{\alpha(\omega(p)+o(1))}$. Applying Definition 3.2 and comparing exponents, this implies that

which proves the convexity of the function $\omega(k)$. ∎

We can upper bound $\omega(1.68568)$ in the following sense

where the first step follows from convexity of $\omega$ (Lemma 3.6), the third step follows from $\omega(1.5)\leq 2.79654$ and $\omega(1.75)\leq 3.02159$ (Lemma 3.4). ∎

Property II is then an immediate consequence of the following inequality, which we prove next:

Define $b=2/a\in(0,\infty)$. Then the desired inequality in (14) can be expressed in terms of $b$ as

Notice that the RHS of (15) is a maximum of two linear functions of $b$ and these intersect at $b^{*}=\omega(1)-1$. By the convexity of $\omega({}\cdot{})$ as proved in Lemma 3.6, it suffices to verify (15) at the endpoints $b\rightarrow 0$, $b\rightarrow\infty$ and $b=b^{*}$. In the case where $b=\delta$ for any $\delta<1$, (15) follows immediately from the observation that $\omega(\delta)<\omega(1)$. We next argue about the case $b\rightarrow\infty$. By Lemma 3.4 we have $\omega(2)\leq 3.252$. Using Lemma 3.5, we have $\omega(b)\leq b-2+\omega(2)$. Combining these two facts implies that for any $b>2$, we have

which again satisfies (15). The final case is $b=b^{*}=\omega(1)-1$, for which (15) is equivalent to

By Lemma 3.4, we have that $\omega(1)-2\in[0,0.372927]$. Then to prove (16), it is sufficient to show that

By the convexity of $\omega({}\cdot{})$ as proved in Lemma 3.6, the upper bound of $\omega(2)\leq 3.251640$ in Lemma 3.4, and recalling that $\omega(1)=t+2$ for $t\in[0,0.372927]$, we have for $k\in$,

In particular, using this inequality for $k=t+1$, we have

which is negative on the entire interval $[0,0.372927]$. This establishes (17) and finishes the proof. ∎

For any two positive integers $n$ and $m$, we have

Case 1: $a\in[1.18647,\infty)$. In this case, we have $\omega(2/a)\leq\omega(2/1.18647)\leq\omega(1.68568)<3$, where the last inequality follows from Claim 3.7. This implies that

Case 2: $a\in(0,1.18647]$. In this case, we have $2/a\in[1.68567,\infty)$. Consider the linear function

An application of Lemma 3.5 then gives, for any $t\geq 2$, the inequality

where the last inequality is by definition of $y(t)$ from (19). Therefore, combining the convexity of $\omega({}\cdot{})$, as proved in Lemma 3.6, with (20), (21), and (22), we conclude that for any $t\in[1.68567,\infty)$, the function $\omega$ is bounded from above by the affine function $y$, expressed as follows.

Combining the results from (18) and (23) finishes the proof of the lemma. ∎

Main Theorem

In this section, we give the formal statement of our main result.

Consider a semidefinite program with variable size $n\times n$ and $m$ constraints (assume there are no redundant constraints):

where $\omega$ is the exponent of matrix multiplication, $X^{*}$ is any optimal solution to the semidefinite program in (24), and $\left\|A_{i}\right\|_{1}$ is the Schatten $1$-norm of matrix $A_{i}$.

The proof of Theorem 4.1 is given in the subsequent sections.

Approximate Central Path via Approximate Hessian

Our main result of this section is the following.

where $X^{*}$ is any optimal solution to the semidefinite program in Definition 1.1, and $\left\|A_{i}\right\|_{1}$ is the Schatten $1$-norm of matrix $A_{i}$. Further, in each iteration of Algorithm 1, the following invariant holds for $\alpha_{H}=1.03$:

At the start of Algorithm 1, Lemma 9.1 is called to modify the semidefinite program to obtain an initial dual solution $y$ for the modified SDP that is close to the dual central path at $\eta=1/(n+2)$. This ensures that the invariant $g_{\eta}(y)^{\top}H(y)^{-1}g_{\eta}(y)\leq\epsilon_{N}^{2}$ holds at the start of the algorithm. Therefore, by Lemma 5.4 and Lemma 5.5, this invariant continues to hold throughout the run of the algorithm. Therefore, after $T=\frac{40}{\epsilon_{N}}\sqrt{n}\log\left(\frac{n}{\delta}\right)$ iterations, the step size $\eta$ in Algorithm 1 grows to $\eta=(1+\frac{\epsilon_{N}}{20\sqrt{n}})^{T}/(n+2)\geq 2n/\delta^{2}$. It then follows from Lemma 5.6 that

Thus when the algorithm stops, the dual solution $y$ has duality gap at most $\delta^{2}$ for the modified SDP. Lemma 9.1 then shows how to obtain an approximate solution to the original SDP that satisfies the guarantees in (25).

where we used the fact that $\Delta_{S}=\sum_{i=1}^{m}(\delta_{y})_{i}A_{i}$. It then follows from Lemma 5.4 and the invariant $g_{\eta}(y)^{\top}H(y)^{-1}g_{\eta}(y)\leq\epsilon_{N}^{2}$ that

where $\alpha_{H}=1.03$. Combining Equation (27) with Inequality (5.1) completes the proof of the theorem. ∎

2 Approximate slack update

3 Closeness of slack implies closeness of Hessian

Then both $\widetilde{H}$ and $H$ are positive semidefinite. For any accuracy parameter $\alpha_{S}\geq 1$, if

As the RHS of (30) and (31) are non-negative, both $\widetilde{H}$ and $H$ are positive semidefinite. Since $\widetilde{S}\preceq\alpha_{S}\cdot S$, we have $S^{-1}\preceq\alpha_{S}\cdot\widetilde{S}^{-1}$ (see Section 2.2), which gives the following inequalities

Combining (32) and (33) with (30) and (31) along with the fact that $v$ can be any arbitrary $n$-dimensional vector finishes the proof of the lemma. ∎

4 Approximate Hessian maintenance

In each iteration of Algorithm 1, for $\alpha_{H}=1.03$, the approximate Hessian $\widetilde{H}(y)$ satisfies that

where $\epsilon_{S}=0.01$ as in Algorithm 2. By definition of operator norm, this implies that in each iteration of Algorithm 1, we have, for $\alpha_{S}=1.011$,

The statement of this lemma then follows from Lemma 5.3. ∎

5 Invariance of Newton step size

The following lemma is standard in the theory of interior point methods (e.g. see [Ren01]).

6 Approximate optimality

The following lemma is also standard in interior point method.

Let $y^{*}$ be an optimal solution to the dual formulation (2). Then we have

Low-rank Update

Crucial to being able to efficiently approximate the Hessian in each iteration is the condition that the rank of the update be not too large. We formalize this idea in the following theorem, essential to the runtime analysis in Section 7.

Let $r_{0}=n$ and $r_{i}$ be the rank of the update to the approximate slack matrix $\widetilde{S}$ when calling Algorithm 2 in iteration $i$ of Algorithm 1. Then, over $T$ iterations of Algorithm 1, the ranks $r_{i}$ satisfy the inequality

and Cauchy-Schwarz inequality. This proves the lemma.

Let $\lambda(M)_{[i]}$ denote the $i$’th (ordered) eigenvalue of a matrix $M$. We then have

where the last inequality is because the first assumption from (36) implies $\nu_{i}\geq 0.98$ for all $i\in[n]$. Plugging (40) into the right hand side of (6), we have

Let $W=U\Sigma V^{\top}$ be the singular value decomposition of $W$, with $U$ and $V$ being $n\times n$ unitary matrices. Because of the invariance of the Frobenius norm under unitary transformation, (40) is then equivalent to

Since $U$ and $V$ are unitary, the matrix $WZW^{\top}=U\Sigma V^{\top}ZV\Sigma U^{\top}$ is similar to $\Sigma V^{\top}ZV\Sigma$, and the matrix $Z^{\prime}=V^{\top}ZV$ is similar to $Z$. Therefore,

where the last inequality is by Fact 2.3. We rewrite the Frobenius norm as

Case 1. There does not exist an $i\leq n/2$ that satisfies the two conditions $y_{[2i]}<\epsilon_{S}$ and $y_{[2i]}<(1-1/10\log n)y_{[i]}$. In this case, we have $r=n/2$. We consider two sub-cases.

Case (b). There exists a minimum index $i\leq n/2$ such that $y_{[2j]}<\epsilon_{S}$ holds for all $j$ in the range $i\leq j\leq n/2$. In this case, for all $j$ in the above range, we have that $y_{[2j]}\geq(1-1/10\log n)y_{[j]}$. In particular, picking $j=i,2i,\cdots$ gives

where $\epsilon_{S}=0.01$ by Table 5.1. Therefore, we can bound, from below, the decrease in potential function as

From Lemma 6.3, we have the following potential decrease:

We note that $\Phi(Z^{(0)})=0$ as we initialized $\widetilde{S}=S$ in the beginning of the algorithm, and that the potential function $\Phi(Z)$ is always non-negative. The theorem then follows by summing up (48) over all $T$ iterations. ∎

Runtime Analysis

Our main result of this section is the following bound on the runtime of Algorithm 1.

The total runtime of Algorithm 1 for solving an SDP with variable size $n\times n$ and $m$ constraints is at most $O^{*}\left(\sqrt{n}\left(mn^{2}+\max(m,n)^{\omega}\right)\right)$, where $\omega$ is the matrix multiplication exponent as defined in Definition 3.2.

To prove Theorem 7.1, we first upper bound the runtime in terms of fast rectangular matrix multiplication times. The iteration complexity of Algorithm 1 is $T=\widetilde{O}(\sqrt{n})$.

The total runtime of Algorithm 1 over $T$ iterations is upper bounded as

The total runtime of Algorithm 1 consists of two parts:

Part 1. The time to compute the approximate Hessian $\widetilde{H}(y)$ (which we abbreviate as $\widetilde{H}$) in Line 11 - 15.

Part 2. The total cost of operations other than computing the approximate Hessian.

We analyze the cost of computing the approximate Hessian $\widetilde{H}$.

We start with computing $\widetilde{H}$ in the first iteration of the algorithm. Each entry of $\widetilde{H}$ involves the computation

It first costs $O^{*}(n^{\omega})$ to invert $\widetilde{S}$. Then the cost of computing the key module of the approximate Hessian, $\widetilde{S}^{-1/2}A_{j}\widetilde{S}^{-1/2}$ for all $j\in[m]$, is obtained by stacking the matrices $A_{j}$ together:

Vectorizing the matrices $\widetilde{S}^{-1/2}A_{j}\widetilde{S}^{-1/2}$ into row vectors of length $n^{2}$, for each $j\in[m]$, and stacking these rows vertically to form a matrix $B$ of dimensions $m\times n^{2}$, one observes that $\widetilde{H}=BB^{\top}$. We therefore have,

Combining (50), (51), and the initial cost of inverting $\widetilde{S}$ gives the following cost for computing $\widetilde{H}$ for the first iteration:

Part 1b. Accumulating low-rank changes over all the iterations

Using this bound over all $T=\widetilde{O}(\sqrt{n})$ iterations, and applying $\sum_{i=0}^{T}\sqrt{r_{i}}\leq\widetilde{O}(\sqrt{n})$ from Theorem 6.1, gives

where we incorporated the bound from (52) into the $i=0$ case.

Observe that there are four operations performed in Algorithm 1 other than computing $\widetilde{H}$:

Part 2a. computing the gradient $g_{\eta}(y)$

Part 2b. inverting the approximate Hessian $\widetilde{H}$

Part 2b. The cost of inverting the approximate Hessian $\widetilde{H}$ is $O(m^{\omega+o(1)})$ per iteration.

The total cost of operations other than computing the Hessian over the $T=\widetilde{O}(\sqrt{n})$ iterations is therefore bounded by

Combining (7) and (7) and using $r_{0}=n$ finishes the proof of the lemma.

We give only the proof of Property I, as the proof of Property II is similar. Let $m=n^{a}$. For each $i\in[T]$, let $r_{i}=n^{b_{i}}$, where $b_{i}\in$. Then

For each number $k\in\{0,1,\cdots,\log n\}$, define the set of iterations

Then our assumption on the sequence $\{r_{1},\cdots,r_{T}\}$ can be expressed as $\sum_{k=0}^{\log n}|I_{k}|\cdot 2^{k/2}\leq O(T\log^{1.5}n)$. This implies that for each $k\{0,1,\cdots,\log n\}$, we have $|I_{k}|\leq O(T\log^{1.5}n/2^{k/2})$. Next, taking the summation of Eq. (59) over all $i\in[T]$, we have

where the fourth step follows from $T=\widetilde{O}(\sqrt{n})$. To bound the exponent on $n$ above, we define the function $g$,

This function is convex in $b_{i}$ due to the convexity of the function $\omega$ (Lemma 3.6). Therefore, over the interval $b_{i}\in$, the maximum of $g$ is attained at one of the end points. We simply evaluate this function at the end points.

Case 1. Consider the case $b_{i}=0$. In this case, we have $g(0)=1/2+a\omega(1/a)$. We consider the following two subcases. Case 1a. If $a\geq 1$, then we have

Case 1b. If $a\in(0,1)$, then we define $k=1/a>1$. It follows from Lemma 3.5 and $\omega>1$, that

Combining both Case 1a and Case 1b, we have that

Case 2 Consider the other case of $b_{i}=1$. In this case, $g(1)=1/2-1/2+a\omega(2/a)=a\omega(2/a)$.

We now finish the proof by combining Case 1 and Case 2 as follows.

In light of Lemma 7.4, the upper bound on runtime given in Lemma 7.2 can be written as

Combining this with 3.10, we have the following upper bound on the total runtime of Algorithm 1:

This finishes the proof of the theorem. ∎

Comparison with Cutting Plane Method

In this section, we prove Theorem 1.3, restated below.

Since $m\geq n$ by assumption, Lemma 3.9 and 3.9 further simplify the runtime to

Initialization

Consider a semidefinite program as in Definition 1.1 of dimension $n\times n$ with $m$ constraints, and assume that it has the following properties.

For any $0<\delta\leq 1$, the following modified semidefinite program

The following are feasible primal and dual solutions:

For any feasible primal and dual solutions $(\overline{X},\overline{y},\overline{S})$ with duality gap at most $\delta^{2}$, the matrix $\widehat{X}=R\cdot\overline{X}_{[n]\times[n]}$, where $\overline{X}_{[n]\times[n]}$ is the top-left $n\times n$ block submatrix of $\overline{X}$, is an approximate solution to the original semidefinite program in the following sense:

where $X^{*}$ is any optimal solution to the original SDP and $\left\|A\right\|_{1}$ denotes the Schatten $1$-norm of a matrix $A$.

Notice that $\overline{X}$ is a feasible primal solution to the modified SDP, and that

where the first step follows because the modified SDP is a maximization problem, and the final step is because $X$ is an optimal solution to the original SDP.

Therefore, we can lower bound the objective value for $\overline{X}_{[n]\times[n]}$ in the original SDP as

where the last inequality follows from (63). By matrix Hölder inequality, we have

where the final step follows from the upper bound of $\theta$ in (64). Summing the above inequality up over all $i\in[m]$ finishes the proof of the lemma. ∎

Acknowledgment

We thank Aaron Sidford for many helpful discussions and Deeksha Adil, Sally Dong, Sandy Kaplan, and Kevin Tian for useful feedback on the writing. We gratefully acknowledge funding from CCF-1749609, CCF-1740551, DMS-1839116, Microsoft Research Faculty Fellowship, and Sloan Research Fellowship. Zhao Song is partially supported by Ma Huateng Foundation, Schmidt Foundation, Simons Foundation, NSF, DARPA/SRC, Google and Amazon.

References

Appendix A Matrix Multiplication: A Tensor Approach

The main goal of this section is to rederive, using tensors, some of the technical results from Section 3. In particular, we use tensors to derive upper bounds on the time to perform the following two rectangular matrix multiplication tasks (Lemma A.12 and A.13):

Multiplying a matrix of dimensions $m\times n^{2}$ with one of dimensions $n^{2}\times m$,

Multiplying a matrix of dimensions $n\times mn$ with one of dimensions $mn\times n$.

Our hope is that these techniques will eventually be useful in further improving the results of this paper.

We recall two definitions to describe the cost of certain fundamental matrix operations, along with their properties.

For any three positive integers $n,m,r$, we have

A.2 Matrix multiplication tensor

The rank of a tensor $T$, denoted as $R(T)$, is the minimum number of simple tensors that sum up to $T$. For any two tensors $S=(S_{i,j,k})_{i,j,k}$ and $T=(T_{a,b,c})_{a,b,c}$, we write $S\leq T$ if there exist three matrices $A,B$ and $C$ (of appropriate sizes) such that $S_{i,j,k}=\sum_{a,b,c}A_{i,a}B_{j,b}C_{k,c}T_{a,b,c}$ for all $i,j,k$. For any $i,j,k$, denote $e_{i,j,k}$ the tensor with $1$ in the $(i,j,k)$-th entry, and elsewhere.

For any three positive integers $a,b,c$, we define

to be the matrix-multiplication tensor corresponding to multiplying a matrix of size $a\times b$ with one of size $b\times c$.

It’s not hard to show that for any $n_{i}$ and $m_{i}$ where $i=1,2,3$, we have

Let $\langle n\rangle=\sum_{i\in[n]}e_{i,i,i}$ be the identity tensor. For any three tensors $S,T_{1}$ and $T_{2}$, if $T_{1}\leq T_{2}$, then we have

Tensor rank is monotone under the relation $\leq$, i.e. if $T_{1}\leq T_{2}$, then we have

For any tensors $T_{1}$ and $T_{2}$, we have

The tensor rank of a matrix multiplication tensor is equal to the cost of multiplying the two correponding sized matrices up to some constant factor, i.e.,

A.3 Implication of matrix multiplication technique

where the last line follows from Lemma A.7. Applying Lemma A.8, we have

Using the definition of $\omega(p)$, we have

Comparing the exponent on both sides completes the proof. ∎

The next lemma establishes the convexity of $\omega(k)$ as a function of $k$.

The fast rectangular matrix multiplication time exponent $\omega(k)$ as defined in Definition A.2 is convex in $k$.

Let $k=\alpha\cdot p+(1-\alpha)\cdot q$ for $\alpha\in(0,1)$. We have

where the last line follows from Lemma A.7. By Lemma A.8, we have

By definition of $\omega(\cdot)$, we have

We only prove the case of $m\geq n$, as the other case where $m<n$ is similar. This is an immediate consequence of Lemma A.11 by taking $a=c=n$, $b=n^{2}$, and $k=\lfloor m/n\rfloor$, where $k$ is a positive integer because $m\geq n$. ∎

Applying the tensor rank on both sides, we have

Let $m=n^{a}$, where $a\in(0,\infty)$. We have

The Property II is then an immediate consequence of the following inequality, which we prove next:

Define $b=2/a\in(0,\infty)$. Then the above desired inequality can be expressed in terms of $b$ as

Notice that the RHS of (15) is a maximum of two linear functions of $b$ and these intersect at $b^{*}=\omega(1)-1$. By the convexity of $\omega({}\cdot{})$ as proved in Lemma A.10, it suffices to verify (15) at the endpoints $b\rightarrow 0$, $b\rightarrow\infty$ and $b=b^{*}$. In the case where $b=\delta$ for any $\delta<1$, (15) follows immediately from the observation that $\omega(\delta)<\omega(1)$. For the case $b\rightarrow\infty$, by Lemma A.3 we have $\omega(2)\leq 3.252$. It then follows from Lemma A.9 that for any $b>2$, we have

The final case is where $b=b^{*}=\omega(1)-1$, for which (15) is equivalent to

By Lemma A.3, we have that $\omega(1)-2\in[0,0.372927]$. Then to prove (66), it is sufficient to show that

By the convexity of $\omega({}\cdot{})$ as proved in Lemma A.10 and the upper bound of $\omega(2)\leq 3.251640$ in Lemma A.3, we have for $k\in$,

In particular, using this inequality for $k=t+1$, we have

which is negative on the entire interval $[0,0.372927]$. This establishes (67) and finishes the proof of the lemma. ∎