# Nearest positive semidefinite matrix to a symmetric matrix in the spectral norm

So I have a symmetric matrix $$A$$ and I would like to solve the optimization problem, $$\hspace{2.5mm}\text{Minimize}\;\; \|A-S\|_2$$ $$\hspace{-5mm}\text{Subject to}\;\; S\geq0.$$ $$A$$ is given and $$S$$ is the variable.

I can put this straight into my solver (SCSSolver in Convex.jl) but it's too slow for my matrix size (n=500). I'm pretty sure it's the processor and not a lack of RAM, since it taxes my processor to 100% but there is no read/write occurring on my disk while the process is running.

I tried setting it up as a Schur complement problem but that didn't help. Since $$A-S$$ is symmetric I'm really just looking at the largest magnitude eigenvalue, I understand that the power-iteration algorithm is quite efficient, could I somehow incorporate that? What is the standard procedure for speeding this up? Is there a way to make it much faster by restricting to a subset of the positive semidefinite cone? I'm willing to take a hit on the accuracy if necessary.

• Perhaps julia's solver naively precomputes all of the eigenvalues. You could try the Matlab solver (cvx) or the Python solver (cvxpy) in case they happen to be implemented more efficiently. – k20 Feb 5 '15 at 23:00
• What norm is it that you consider? The induced $l_2$ matrix norm, or the Frobenius norm? – Wolfgang Bangerth Feb 6 '15 at 0:32
• I'm asking because if the norm is quadratic, then in essence what you're looking for is just the projection onto the cone of SPD matrices. – Wolfgang Bangerth Feb 6 '15 at 0:33
• @WolfgangBangerth the induced $l_2$ norm (operator/spectral norm) – Thoth Feb 6 '15 at 0:56
• Aw, awkward. No easy way to construct a projection out of that :-( – Wolfgang Bangerth Feb 6 '15 at 2:46

There's no need to use the Schur complement here because $$A$$ and $$S$$ are already symmetric. The conventional formulation of this problem as an SDP is

$$\min t \quad\text{ subject to}\\ A-S+tI \succeq 0 \\ A-S - tI \preceq 0$$

results in an SDP of the form (I'm specific here because there are so many different "standard forms" for SDP)

$$\min c^{T}x\\ F_{0}+x_{1}F_{1}+...+x_{m}F_{m} \succeq 0$$

where

$$x=[t\;\; S_{1,1} \;\; S_{1,2} \;\; S_{1,3} \; \ldots \; S_{n,n}]^{T}$$

has $$1+n(n-1)/2$$ elements, and the matrices $$F_{i}$$ are of sizes $$2n \times 2n$$.

For your problems with $$n=500$$, this means that you have about 125,000 variables $$x_{j}$$, and matrices of size $$1000$$ by $$1000$$. Primal-dual interior-point methods (as implemented in SDPA, SDPT3, SeDuMi, CSDP, etc.) would require the solution of a system of 125,000 equations in 125,000 variables in each iteration, and this system would typically be fully dense.

This is not something you can do on a typical desktop machine but would be within reach on a more powerful server with enough RAM (more than 128 gigabytes.) Solution times would be fairly long (10+ hours)

The "splitting cone solver" that Julia implements appears to me to be an implementation of a first-order ADMM method such as the one described in

The problem with ADMM for SDP is that it doesn't yield very accurate solutions (you'll typically get two or three digits of accuracy at best) and the method requires lots of tuning to achieve good results. It's not nearly as robust (in the sense that it can consistently solve every problem that you throw at it) as the primal-dual interior-point methods.

I can't vouch for the particular implementation of this method in Julia since I've never worked with it. It might be that some other first-order code for SDP could do a better job than SCSSolver.

There is probably some work in the convex optimization literature on this kind of spectral norm minimization that avoids even formulating an SDP, but I'm no expert on that topic.

Are you sure that you really need the best approximation in the matrix 2-norm sense rather than in the Frobenius norm?

Note added much later: It turns out that the approach that computes this projection in the Frobenius norm also works equally well to compute the projection in the spectral norm. First, perform an orthogonal diagonalization of $$A$$ as

$$A=\sum_{i=1}^{n} \lambda_{i} u^{(i)}u^{(i)^{T}}$$

Then

$$S=\sum_{i=1}^{n} \lambda_{i}^{+} u^{(i)}u^{(i)^{T}}$$

When I wrote this answer, I had known that this worked for the Frobenius norm projection, but not that it also works for the 2-norm projection. For non-symmetric $$A$$, or if there are additional constraints, this simple formula isn't applicable and the SDP formulation might still be useful.

• thanks for the response, this is for learning purposes so I don't explicitly need it for anything. I'll probably accept this answer eventually but I'm going to leave the question upon for a bit longer in case someone happens to have seen a way to do this more efficiently. – Thoth Feb 10 '15 at 3:08
• For clarity, the solver is not written in Julia, its actually a C library written by the authors of that paper with a Julia binding. – IainDunning Feb 23 '15 at 0:30
• Now, some years later, SuperSCS and SDPNAL+ both exist. Perhaps one of them would do better. – Mark L. Stone Dec 14 '19 at 20:26

According to the paper

Nicholas J. Higham, MR 943997 Computing a nearest symmetric positive semidefinite matrix, Linear Algebra Appl. 103 (1988), 103--118,

one of the solution is in the following form, $$A_{opt} = S + |\lambda| I$$ where $I$ is the identity matrix, $\lambda$ is the smallest negative eigenvalue of $S$

Since $A$ is symmetric, you can diagonalize it with $A=XDX^t$, and $S=XWX^t$, to get the equivalent formulation $$\text{minimize}\quad\|D-W\|_2, \qquad \text{s.to }W\geq0.$$ So for any positive eigenvalue $d_i$ of $A$, you can use $W_{ii}=d_i$ to subtract it off.

I suspect that the answer might be the absolute value $|\lambda_-|$ of the largest negative eigenvalue of $A$, with $S$ equal to the same matrix as $A$ but with negative eigenvalues zeroed out, and with positive eigenvalues of $A-S$ reduced to be at most as large as $|\lambda_-|$.

• I don't think this works since there is no reason the factorizations of $A$ and $S$ will have the same eigenvector matrices. – Thoth Feb 5 '15 at 19:41
• @BloodPudding I didn't say $W$ is going to be diagonal. Whatever $S$ is, you can construct $W=X^tSX$, then try to solve the problem with a diagonal matrix $D$ and an arbitrary matrix $W$. All I said was that $W$ is probably going to be diagonal. – Kirill Feb 5 '15 at 22:50