Given a symmetric matrix and a nonnegative number , what is the nearest symmetric matrix whose eigenvalues are all at least ? In other words, how can we project a symmetric matrix onto the set of symmetric positive semidefinite matrices with eigenvalues at least ?
Distance can be measured in any norm, but the most common choice is the Frobenius norm, .
This problem occurs in a very wide range of applications. A typical scenario is that a covariance matrix is approximated in some way that does not guarantee positive semidefinitess, for example by treating blocks of the matrix independently. In machine learning, some methods use indefinite kernels and these can require an indefinite similarity matrix to be replaced by a semidefinite one. When , the problem arises when the matrix is positive definite but ill conditioned and a matrix of smaller condition number is required, or when rounding errors result in small negative eigenvalues and a “safely positive definite” matrix is wanted.
The following theorem gives the solution to the problem for the Frobenius norm.
Theorem (Cheng and Higham, 1998).
Let the symmetric matrix have the spectral decomposition and let . The unique matrix with smallest eigenvalue at least nearest to in the Frobenius norm is given by
The theorem says that there is a unique nearest matrix and that is has the same eigenvectors as with any eigenvalues less than increased to . When , another way to express this result is in terms of the polar decomposition , where is orthogonal and is symmetric positive definite (with , since is symmetric): .
One can pose the same nearness problem for nonsymmetric . For the Frobenius norm there is again a unique nearest symmetric matrix whose eigenvalues are all at least , and it is given by applying the theorem to the symmetric part of (which is the nearest symmetric matrix to ), where
For the -norm the answer is known for nonsymmetric only for , so we take from this point on. We write to denote that the symmetric matrix is positive semidefinite and use a superscript to denote the positive semidefinite square root.
Theorem (Halmos, 1972).
The -norm distance from to the nearest positive semidefinite matrix is
and a nearest matrix is
When is symmetric, so that , the nearest matrix and the distance reduce to
Clearly, is not, in general, unique because when is symmetric, instead of shifting by a multiple of we could just perturb the negative eigenvalue to zero, as with the optimal Frobenius norm perturbation, without changing the -norm of the perturbation. A distinguishing feature of in Halmos’s theorem is that for any other nearest positive semidefinite matrix (Bouldin, 1973, Theorem 4.2); thus has the minimum number of zero eigenvalues over all nearest positive semidefinite matrices.
The Frobenius norm solution is an approximate minimizer of the -norm distance since (Higham, 1988)
Halmos’s theorem simplifies the computation of because it reduces the minimization problem to one dimension. However, the problem is nonlinear and has no closed form solution for general , so iterative methods are needed to compute . Two algorithms are developed in Higham (1988). Both are based on the following properties of the matrix in the theorem: is monotone increasing on , and either , in which case , or has a unique zero . The first algorithm uses the bisection method, determining the sign of by attempting a Cholesky factorization of : the sign is nonnegative if the factorization exists and negative otherwise. This approach is attractive if is required only to low accuracy. For higher accuracy computations it is better to apply a more rapidly converging zero finder to . A hybrid Newton–bisection method is used in Higham (1988).
For a numerical example (in MATLAB), we take the Jordan block
The symmetric part of , , has two negative eigenvalues and a zero eigenvalue:
>> A = gallery('jordbloc',5,0); As = sym(A); eig_B = eig((As+As')/2)' eig_B = [-1/2, 0, 1/2, -3^(1/2)/2, 3^(1/2)/2] >> double(eig_B) ans = -5.0000e-01 0 5.0000e-01 -8.6603e-01 8.6603e-01
The nearest symmetric positive semidefinite matrix to in the Frobenius norm can be computed in MATLAB as
B = (A + A')/2; [Q,D] = eig(B); d = diag(D); X_F = Q*diag(max(d,0))*Q';
We can improve this code by using the implicit expansion feature of MATLAB to avoid forming a diagonal matrix. Since the computed result is not exactly symmetric because of rounding errors, we also need to replace it by the nearest symmetric matrix:
[Q,d] = eig(B,'vector'); X_F = Q*(max(d,0).*Q'); X_F = (X_F + X_F')/2;
X_F = 1.972e-01 2.500e-01 1.443e-01 -4.163e-17 -5.283e-02 2.500e-01 3.415e-01 2.500e-01 9.151e-02 -1.665e-16 1.443e-01 2.500e-01 2.887e-01 2.500e-01 1.443e-01 -4.163e-17 9.151e-02 2.500e-01 3.415e-01 2.500e-01 -5.283e-02 -1.665e-16 1.443e-01 2.500e-01 1.972e-01
which has eigenvalues
-3.517e-17 7.246e-17 9.519e-17 5.000e-01 8.660e-01
Notice that the three eigenvalues that would be zero in exact arithmetic are of order the unit roundoff. A nearest matrix in the -norm, as given by Halmos’s theorem. is
X_2 = 8.336e-01 5.000e-01 1.711e-01 0.000e+00 -1.756e-02 5.000e-01 6.625e-01 5.000e-01 1.887e-01 0.000e+00 1.711e-01 5.000e-01 6.450e-01 5.000e-01 1.711e-01 0.000e+00 1.887e-01 5.000e-01 6.625e-01 5.000e-01 -1.756e-02 0.000e+00 1.711e-01 5.000e-01 8.336e-01
and its eigenvalues are
-7.608e-17 1.281e-01 5.436e-01 1.197e+00 1.769e+00
The distances are
This is a minimal set of references, which contain further useful references within.
- Richard Bouldin, Positive Approximants, Trans. Amer. Math. Soc. 177, 391–403, 1973.
- Sheung Hun Cheng and Nicholas Higham, A Modified Cholesky Algorithm Based on a Symmetric Indefinite Factorization, SIAM J. Matrix Anal. Appl. 19(4), 1097–1110, 1998.
- Nicholas J. Higham, Computing a Nearest Symmetric Positive Semidefinite Matrix, Linear Algebra Appl. 103, 103-118, 1988
- Paul R. Halmos, Positive Approximants of Operators, Indiana Univ. Math. J. 21, 951–960, 1972.
Related Blog Posts
- What Is a Symmetric Positive Definite Matrix? (2020)
- What Is the Nearest Symmetric Matrix? (2020)
- What Is the Polar Decomposition? (2020)