# Best approach for numerically computing the pseudo-inverse of a covariance matrix

What are the reasons to prefer eigenvalue decomposition over singular value decomposition for numerically computing the pseudo-inverse of a symmetric real matrix? In the case when you want to form the pseudo-inverse a covariance matrix and you also have access to the original data matrix, does using SVD make more sense?

Using SciPy as an example, the documentation for scipy.linalg.pinvh states

Compute the (Moore-Penrose) pseudo-inverse of a Hermitian matrix.

Calculate a generalized inverse of a Hermitian or real symmetric matrix using its eigenvalue decomposition and including all eigenvalues with 'large' absolute value.

Why do we prefer eigenvalue decomposition in this case?

Suppose instead of computing the pseudo-inverse of an arbitrary matrix, you start instead with an $m \times n$ matrix $X$ and want to compute the pseudo-inverse of its covariance matrix $C = XX^*$. It seems like there are two obvious and closely related approaches.

Since $C$ is real and symmetric, it can be orthogonally diagonalized: $C = Q \Lambda Q^T$, where $Q$ contains eigenvectors, $\Lambda$ contains eigenvalues, and $QQ^T = I$. We can get the pseudo-inverse of $C$ from the pseudo-inverse of $\Lambda$ as $C^+ = Q \Lambda^+ Q^T$.

Alternatively, start from the singular value decomposition of $X$ as $X = U \Sigma V^*$. Then $C = U \Sigma \Sigma^* U^*$, and we can get the pseudo-inverse again from the pseudo-inverse of $\Sigma \Sigma^*$.

I thought that computing singular values was always a well-conditioned problem, and that this would mean SVD would be a better approach for computing the pseudo-inverse of a covariance matrix. Is that correct? What other concerns might I be overlooking?

• Generally speaking the SVD is faster and more numerically stable. However, $X X^*$ is $m \times m$ while $X$ is $m \times n$. This means that if $m \ll n$ then it is probably easier to do what you can with $C$ rather than to work with $X$. – Ian Feb 7 '15 at 0:45
• Thanks for pointing that out, @Ian. In my case I do expect $n$ to be much larger than $m$, so that’s good to keep in mind. – DGrady Feb 7 '15 at 0:48
• When a matrix is Hermitian, its singular values are the absolute values of its eigenvalues. Using the SVD and the Schur decomposition are effectively equivalent. – Michael Grant Feb 7 '15 at 0:48
• Yes. The symmetric singular value problem is very well behaved. You're gaining nothing by computing the SVD instead, except the extra complexity of having to manage left- and right- singular vectors. – Michael Grant Feb 7 '15 at 0:51
• In fact, LAPACK (for instance) doesn't even supply an SVD for symmetric matrices. You're expected to use the eigenvalue routines instead. – Michael Grant Feb 7 '15 at 0:56

I wouldn't say we "prefer" the eigenvalues to the singular values in the symmetric/Hermitian case. They are effectively the same thing! Or, more accurately, the singular values are nothing more than the absolute values of the eigenvalues.

To see why, let $X$ be any Hermitian matrix. Its Schur (eigenvalue) decomposition is $X=Q\Lambda Q^H$ where $\Lambda$ is diagonal, and $Q^HQ=QQ^H=I$. Its singular value decomposition, on the other hand, is $X=U\Sigma V^H$, where $\Sigma$ is diagonal and nonnegative, and $U^HU=UU^H=I$ and $V^HV=VV^H=I$.

Alternatively, we can write these decompositions in dyadic form: $$X=\sum_{i=1}^n \lambda_i q_i q_i^H = \sum_{i=1}^n \sigma_i u_i v_i^H$$ where $q_i$, $u_i$, and $v_i$ are the $i$th columns of $Q$, $U$, and $V$, respectively, and $\lambda_i$ and $\sigma_i$ are the $i$th diagonal elements of $\Lambda$ and $\Sigma$, respectively.

In this latter form, it is simple to see that, given the Schur decomposition, $$\sigma_i = |\lambda_i|, \quad u_i = q_i, \quad v_i = \begin{cases} v_i & \lambda_i \geq 0 \\ -v_i & \lambda_i < 0 \end{cases} \qquad i=1,2,\dots, n$$ gives us a valid singular value decomposition. Our choice to apply the sign flipping to $v_i$ is arbitrary.

Be assured that in practice, this is precisely what you should be doing to compute the singular value decomposition of a Hermitian matrix. The symmetric eigenvalue problem is well-conditioned, and there is no numerical advantage to applying the a nonsymmetric SVD to a Hermitian matrix. In fact, the computation of the singular values of a non-Hermitian matrix $A$ are sometimes computed by applying a symmetric eigenvalue algorithm to the Hermitian matrix $$\begin{bmatrix} 0 & A \\ A^T & 0 \end{bmatrix}$$ I'll leave it as an exercise for you to prove that this works :-)

EDIT: DGrady's comment alerted me to the fact that I didn't read the second half of the question correctly. In it, he was comparing the merits of computing the SVD of $X$ versus the eigenvalue decomposition of $C=XX^H$. In this case, it is indeed going to be a bit better to compute the SVD of $X$ directly.

• Thanks a lot for your answer! Note that I was comparing the SVD of $X$ (an arbitrary rectangular matrix) to the eigenvalue decomposition of $C=XX^*$ (a Hermitian matrix) — so the singular values of $X$ should be the positive square roots of the eigenvalues of $C$, I think. The goal is to find the pseudo-inverse of $C$, and I was curious if there might be an advantage to approaching this through the SVD of $X$. I appreciate you pointing out that this is not the case, since the eigenvalue problem for Hermitian matrices is well-conditioned. – DGrady Feb 11 '15 at 1:07
• Yes, if you are starting with $X$, then using the singular value decomposition may indeed be better. Squaring up to get $C$ will reduce accuracy, but the subsequent eigenvalue decomposition will be stable. – Michael Grant Feb 18 '15 at 15:50