Well-conditioned Primal-Dual Interior-point Method for Accurate Low-rank Semidefinite Programming^††thanks: Financial support for this work was provided by NSF CAREER Award ECCS-2047462 and ONR Award N00014-24-1-2671.

By taking advantage of fast matrix-vector products with the data $\mathcal{A}$, an iterative approach can reduce the cost of solving the Newton subproblem (2) to as low as $O(n^{3})$ time and $O(n^{2})$ memory. However, without an effective preconditioner, these improved complexity figures are limited to coarse accuracies of $\mu\approx 10^{-2}$. The key issue is that the scaling point $W$ becomes ill-conditioned as $\operatorname{cond}(W)=\Theta(1/\mu)$. Iterating until the numerical floor yields a numerically exact solution to (2) that is indistinguishable from the true solution under finite precision, but the number of iterations to achieve this grows as $O(1/\mu\log(1/\varepsilon_{\mathrm{mach}}))$. Truncating after a much smaller number of iterations yields an inexact solution to (2) that can slow the convergence of the outer IPM.

Toh and Kojima [45] proposed the first effective preconditioner for SDPs. Their critical insight is the empirical observation that the scaling point $W\in\mathcal{S}_{++}^{n}$ becomes ill-conditioned as $\operatorname{cond}(W)=\Theta(1/\mu)$ only because its eigenvalues split into two well-conditioned clusters:

Based on this insight, they formulated a spectral preconditioner that attempts to counteract the ill-conditioning between the two clusters, and demonstrated promising computational results. However, the empirical effectiveness was not rigorously guaranteed, and its expensive setup cost of $O(mn^{3})$ time made its overall solution time comparable to a direct solution of (2), particularly for problems with $m=\Theta(n^{2})$ constraints.

Inspired by the above, Zhang and Lavaei [53] proposed a spectral preconditioner that is much cheaper to set up, and also rigorously guarantees a high quality of preconditioning under the splitting assumption (A0). Their key idea is to rewrite the scaling point as the low-rank perturbation of a well-conditioned matrix, and then to approximate the well-conditioned part by the identity matrix

Indeed, it is easy to verify that the joint condition number between $W$ and $UU^{T}+\tau I$ is bounded under (A0) even as $\mu\to 0^{+}$. Substituting $W\approx UU^{T}+\tau I$ into (2) reveals a low-rank perturbed problem that admits a closed-form solution via the Sherman–Morrison–Woodbury formula, which they proposed to use as a preconditioner. Assuming (A0) and the availability of fast matrix-vector products that we outline as (A3) below, they proved that their method computes a numerically exact solution in $O(n^{3}r^{3}+n^{3}r^{2}\log(1/\varepsilon_{\mathrm{mach}}))$ time and $O(n^{2}r^{2})$ memory. This same preconditioning idea was later adopted in the solver Loraine [27].

In practice, however, the Zhang–Lavaei preconditioner loses its effectiveness for smaller values of $\mu$, until abruptly failing at around $\mu\approx 10^{-6}$. This is an inherent limitation of preconditioning in finite-precision arithmetic; to precondition a matrix with a diverging condition number $\Theta(1/\mu^{2})$, the preconditioner’s condition number must also diverge as $\Theta(1/\mu^{2})$. As $\mu$ approaches zero, it eventually becomes impossible to solve the preconditioner accurately enough for it to remain effective. This issue is further exacerbated by the use of the Sherman–Morrison–Woodbury formula, which is itself notorious for numerical instability.

Moreover, the effectiveness of both prior preconditioners critically hinges on the splitting assumption (A0), which lacks a rigorous justification. It was informally argued in [45, 44] that the assumpution holds under strict complementarity (1d) and the primal-dual nondegeneracy of Alizadeh, Haeberly, and Overton [2]. But primal-dual nondegeneracy is an excessively strong assumption, as it would require the number of constraints $m$ to be within the range of

In particular, it is not applicable to the diverse range of low-rank matrix recovery problems, like matrix sensing [37, 17], matrix completion [15], phase retrieval [18], that necessarily require $m>nr^{\star}-\frac{1}{2}r^{\star}(r^{\star}-1)$ measurements in order to uniquely recover an underlying rank-$r^{\star}$ ground truth. It is also not applicable to problems with $m=\Theta(n^{2})$ constraints, which as we explained above, are the most difficult SDPs to solve for a fixed value of $n$.

In this paper, we derive a well-conditioned reformulation of the Newton subproblem (2) that costs just $O(n^{3})$ time and $O(n^{2})$ memory to set up. In principle, as the reformulation is already well-conditioned by construction, it allows any iterative method to compute a numerically exact solution in $O(\log(1/\varepsilon_{\mathrm{mach}}))$ iterations for all values of $0<\mu\leq 1$. In practice, the convergence rate can be substantially improved by the use of a well-conditioned preconditioner. Assuming (A0) and the fast matrix-vector products in (A3) below, we provably compute a numerically exact solution in $O(n^{3}r^{3}+n^{3}r^{2}\log(1/\varepsilon_{\mathrm{mach}}))$ time and $O(n^{2}r^{2})$ memory. But unlike Zhang and Lavaei [53], our well-conditioned preconditioner maintains its full effectiveness even at extremely high accuracies like $\mu\approx 10^{-12}$, hence fully addressing a critical weakness of this prior work.

Our main theoretical results are a tight upper bound on the condition number of the reformulation and the preconditioner (3.3), and a tight bound on the number of preconditioned iterations needed to solve (2) to $\epsilon$ accuracy (3.5). These rely on the same splitting assumption (A0) as prior work, which we show is implied by a strengthened version of strict complementarity (1d). In 3.1, we prove that (A0) holds if the primal-dual iterates $(X,S)$ used to compute the scaling point $W$ in (2) satisfy the following centering condition (the spectral norm $\|\cdot\|$ is defined as the maximum singular value):

In particular, the restrictive primal-dual nondegeneracy assumption adopted by prior work [45, 44] is completely unnecessary. By repeating the algebraic derivations in [30], it can be shown that (A1) holds under strict complementarity (1d) if the IPM maintains its primal-dual iterates $(X,S)$ within the “wide” neighborhood $\mathcal{N}_{\infty}$ [48, 35, 6]. In turn, most IPMs that achieve an $\epsilon$ duality gap in $O(\sqrt{n}\log(1/\epsilon))$ iterations—including essentially all theoretical short-step methods [35, 36, 43] but also many practical methods like the popular solver SeDuMi [41]—do so by keeping their iterates within the $\mathcal{N}_{\infty}$ neighborhood [42]. In Section 7.1, we experimentally verify for SDPs satisfying strict complementarity (1d) that SeDuMi [41] keeps its primal-dual iterates $(X,S)$ sufficiently centered to satisfy (A1).

Under (A1), which implies the splitting assumption (A0) adopted in prior work via 3.1, we prove that $O(\log(1/\epsilon))$ preconditioned iterations solve (2) to $\epsilon$ accuracy for all values of $0<\mu\leq 1$. In order for the reformulation and the preconditioner to remain well-conditioned, however, we require a strengthened version of the primal uniqueness assumption (1c). Concretely, we require the linear operator $\mathcal{A}$ to be injective with respect to the tangent space of the manifold of rank-$r^{\star}$ positive semidefinite matrices, evaluated at the unique primal solution $X^{\star}=U^{\star}U^{\star T}$:

This is stronger than primal uniqueness because it is possible for a rank-$r^{\star}$ primal solution to be unique with just $m=\frac{1}{2}r^{\star}(r^{\star}+1)$ constraints [2], but (A2) can hold only when the number of constraints $m$ is no less than the dimension of the tangent space

Fortunately, many rigorous proofs of primal uniqueness, such as for matrix sensing [37, 17], matrix completion [15], phase retrieval [18], actually work by establishing (A2). In problems where (A2) does not hold, it becomes possible for our reformulation and preconditioner to grow ill-conditioned in the limit $\mu\to 0^{+}$. For these cases, the behavior of our method matches that of Zhang and Lavaei [53]; down to modest accuracies of $\mu\approx 10^{-6}$, our preconditioned iterations will still rapidly converge to $\epsilon$ accuracy in $O(\log(1/\epsilon))$ iterations.

Our time complexity figures require the same fast matrix-vector products as previously adopted by Zhang and Lavaei [53] (and implicitly adopted by Toh and Kojima [45]). We assume, for $U,V\in\mathbb{R}^{n\times r}$ and $y\in\mathbb{R}^{m}$, that the two matrix-vector products

can both be performed in $O(n^{2}r)$ time and $O(n^{2})$ storage. This is satisfied in problems with suitable sparsity and/or Kronecker structure. In Section 7.1, we rigorously prove that the robust matrix completion and sensor network localization problems satisfy (A3) due to suitable sparsity in the operator $\mathcal{A}$.

Our well-conditioned reformulation of the Newton subproblem (2) is inspired by Toh [44] and shares superficial similarities. Indeed, the relationship between our reformulation and the Zhang–Lavaei preconditioner [53] has an analogous relationship to Toh’s reformulation with respect to the Toh–Kojima preconditioner [45]. However, the actual methods are very different. First, Toh’s reformulation remains well-conditioned only under primal-dual nondegeneracy, which as we explained above, makes it inapplicable to most low-rank matrix recovery problems, as well as problems with $m=\Theta(n^{2})$ constraints. Second, the effectiveness of Toh’s preconditioner is not rigorously guaranteed; due to its reliance on the PSQMR method, the Faber–Manteuffel Theorem [22] says that it is fundamentally impossible to prove convergence for the preconditioned iterative method. Third, our set-up cost of $O(n^{3}r^{3})$ is much cheaper than Toh’s set-up cost of $O(mn^{3})$; this parallels the reduction in set-up cost achieved by the Zhang–Lavaei preconditioner [53] over the Toh–Kojima preconditioner [45].

We mention another important line of work on the use of partial Cholesky preconditioners [7, 24, 4], which was extended to IPMs for SDPs in [5]. These preconditioners are based on the insight that a small number of ill-conditioned dimensions in the Newton subproblem (2) can be corrected by a low-rank partial Cholesky factorization of its underlying Schur complement. However, to the best of our knowledge, such preconditioners do not come with rigorous guarantees of an improvement.

Outside of SDPs, our reformulation is also inspired by the layered least-squares of Bobrovnikova and Vavasis [9]; indeed, our method can be viewed as the natural SDP generalization. However, an important distinction is that they apply MINRES directly to solve the indefinite reformulation, whereas we combine PCG with an indefinite preconditioner to dramatically improve the convergence rate. This idea of using PCG to solve indefinite equations with an indefinite preconditioner had first appeared in the scientific computing literature [29, 25, 39].

We mention that first-order methods have also been derived to exploit the low-rank structure of the SDP (1) without giving up the convexity [14, 51, 49]. These methods gain their efficiency by keeping all their iterates low-rank, and by performing convex updates on the low-rank iterates $X=UU^{T}$ while maintaining them in low-rank factored form $X_{+}=U_{+}U_{+}^{T}$. While convex first-order methods cannot become stuck at a spurious local minimum, they still require $\mathrm{poly}(1/\epsilon)$ iterations to converge to an $\epsilon$-accurate solution of (1). In contrast, our proposed method achieves the same in just $O(\log(1/\epsilon))$ iterations, for an exponential factor improvement.

Finally, for low-rank SDPs with small treewidth sparsity, chordal conversion methods [23, 28] compute an $\epsilon$ accurate solution in guaranteed $O(n^{1.5}\log(1/\epsilon))$ time and $O(n)$ memory [54, 52]. Where the property holds, chordal conversion methods achieve state-of-the-art speed and reliability. Unfortunately, many real-world problems do not enjoy small treewidth sparsity [32]. For low-rank SDPs that are sparse but non-chordal, our proposed method can solve them in as little as $O(n^{3.5}\log(1/\epsilon))$ time and $O(n^{2})$ memory.

Write $\mathbb{R}^{n\times n}\supseteq\mathcal{S}^{n}\supseteq\mathcal{S}_{+}^{n}\supseteq\mathcal{S}_{++}^{n}$ as the sets of $n\times n$ real matrices, real symmetric matrices, positive semidefinite matrices, and positive definite matrices, respectively. Write $\operatorname{vec}:\mathbb{R}^{n\times n}\to\mathbb{R}^{n^{2}}$ as the column-stacking vectorization, and $A\otimes B$ as the Kronecker product satisfying $(B\otimes A)\operatorname{vec}(X)=\operatorname{vec}(AXB^{T})$. Write $\|\cdot\|$ and $\|\cdot\|_{F}$ as spectral norm and Frobenius norm, respectively. Let $\Psi_{n}$ denote the $n^{2}\times n(n+1)/2$ basis for the set of vectorized real symmetric matrices $\operatorname{vec}(\mathcal{S}^{n})\subseteq\operatorname{vec}(\mathbb{R}^{n\times n})$. We define the symmetric vectorization of $X$ as

We note that $\left\langle X,S\right\rangle=\left\langle\operatorname{vec}_{\mathcal{S}}(X),\operatorname{vec}_{\mathcal{S}}(S)\right\rangle=\left\langle\operatorname{vec}(X),\operatorname{vec}(S)\right\rangle$ for all $X,S\in\mathcal{S}^{n}$. Like the identity matrix, we will frequently suppress the subscript $n$ if the dimensions can be inferred from context. Given $A,B\in\mathbb{R}^{n\times r}$, we define their symmetric Kronecker product as $A\otimes_{\mathcal{S}}B\overset{\mathrm{def}}{=}\frac{1}{2}\Psi_{n}^{T}(A\otimes B+B\otimes A)\Psi_{r}$ so that $(A\otimes_{\mathcal{S}}B)\operatorname{vec}_{\mathcal{S}}(X)=(B\otimes_{\mathcal{S}}A)\operatorname{vec}_{\mathcal{S}}(X)=\operatorname{vec}_{\mathcal{S}}(AXB^{T})$ holds for all $X\in\mathcal{S}^{r}$. (See also [2, 45, 44] for earlier references on the symmetric vectorization and Kronecker product.)

Given a system of equations $Ax=b$ with a symmetric $A$ and right-hand side $b$, we define the corresponding minimum residual (MINRES) iterations with respect to positive definite preconditioner $P$, and initial point $x_{0}$ implicitly in terms of its Krylov optimality condition.

The precise implementation details of MINRES are complicated, and can be found in the standard reference [26, Algorithm 4]. We only note that the method uses $O(n)$ memory over all iterations, and that each iteration costs $O(n)$ time plus the the matrix-vector product $p\mapsto Ap$ and the preconditioner linear solve $r\mapsto P^{-1}r$. In exact arithmetic, MINRES terminates with the exact solution in $n$ iterations. With round-off noise, however, exact termination does not occur. Instead, the behavior of the iterates are better predicted by the following bound.

This paper proposes a well-conditioned reformulation $Ax=b$ for the Newton subproblem (2), whose condition number $\kappa=\operatorname{cond}(A)$ remains bounded over all iterations of the IPM. In principle, 2.2 says that it takes $O(\log(1/\epsilon))$ iterations to solve this reformulation to $\epsilon$ accuracy.

Given a system of equations $Ax=b$ with a symmetric $A$ and right-hand side $b$, we define the corresponding preconditioned conjugate gradient (PCG) iterations with respect to preconditioner $P$ not necessarily positive definite, and initial point $x_{0}$ explicitly in terms of its iterations.

We can verify from 2.3 that the method uses $O(n)$ memory over all iterations, and that each iteration costs $O(n)$ time plus the the matrix-vector product $p\mapsto Ap$ and the preconditioner linear solve $r\mapsto P^{-1}r$. The following is a classical iteration bound for when PCG is used to solve a symmetric positive definite system of equations with a positive definite preconditioner.

In the late 1990s and early 2000s, it was pointed out in several parallel works [29, 25, 39] that PCG can also be used to solve indefinite systems, with the help of an indefinite preconditioner. The proof of 2.5 follows immediately by substituting into 2.3.

Therefore, the indefinite preconditioned PCG (3a) is guaranteed to converge because it is mathematically equivalent to PCG on the underlying positive definite Schur complement system (3b). Nevertheless, (3a) is more preferable when the matrix $C$ is close to singular, because the two indefinite matrices in (3a) can remain well-conditioned even as $\operatorname{cond}(C)\to\infty$. As the two Schur complements become increasingly ill-conditioned, the preconditioning effect in (3b) begins to fail, but the indefinite PCG in (3a) will maintain its rapid convergence as if perfectly conditioned.