I am trying to calculate the singular value decomposition of this matrix using numpy.linalg.svd .

However, reconstructing the matrix from the SVD gives a poor reconstruction - in particular, the first row and second column are far smaller in the reconstruction than in the original matrix.

u, w, vt = np.linalg.svd(M)
np.allclose(M, np.dot(u, np.dot(np.diag(w), vt)))

>> False

The matrix is singular ( w[-1] is zero), and the largest singular value is much bigger than the other non-zero values (3e+24 vs 7e-2 for the next largest).


Should I expect that the reconstruction from the SVD be poor for a matrix this badly conditioned? Are their other more stable ways I could calculate the SVD otherwise?


The reason I am computing the SVD is a diagnosis for the pseudoinverse of M I calculate using numpy.linalg.pinv, which should return the Moore-Penrose pseudoinverse. As I understand it, this particular pseudoinverse should be symmetric since M is symmetric, but it is not. My assumption is that this is because of problems with the SVD, which numpy.linalg.pinv calculates.


1 Answer 1


Algorithms for the SVD, as more or less every classical linear algebra algorithm based on orthogonal transformations, are normwise backward stable, i.e., it should be guaranteed that $\frac{\|USV^* - A \|}{\|A\|} = O(u)$, where the norms are Euclidean norms, $u$ is the machine precision, and "$O(u)$" means that the first-order term in $u$ is bounded by a polynomial in the matrix dimensions.

In particular, this bound in norm does not guarantee that all the elements are close: for instance, compare $$ a=\begin{bmatrix} 10^{15}\\1 \end{bmatrix} \quad \text{vs.} \quad b=\begin{bmatrix} 10^{15}\\-1 \end{bmatrix}. $$ These two vectors are normwise very close to each other, so $\frac{\|a-b\|}{\|b\|}$ is very small, but they are not so elementwise.

numpy.allclose compares elementwise, not normwise.

  • 1
    $\begingroup$ For the record, for the matrix in question the value of $\|USV^*-A\|_2/\|A\|_2$ is $4.85\times 10^{-16}$, which is close to machine precision, so numpy computed the SVD accurately. $\endgroup$
    – Kirill
    Commented Aug 16, 2017 at 21:29
  • $\begingroup$ Thanks for the information, seems like my matrix is badly behaved so have accepted this answer. $\endgroup$
    – myseun
    Commented Aug 16, 2017 at 21:35
  • $\begingroup$ Deleted previous comment as my check was incorrect. This is not the case, the normwise error $ \frac{|| A^+ - (A^+)^* ||}{|| A^+ ||} \approx 0.01$ $\endgroup$
    – myseun
    Commented Aug 17, 2017 at 10:16
  • $\begingroup$ @myseun That is strange. Try using an eigendecomposition instead of the SVD. Since your matrix is symmetric, the eigendecomposition is stable and it is (modulo signs) an SVD (assuming your library returns an orthogonal eigenvector matrix, which should be the case since every language is basically just wrapping Lapack), so you can use it in the same way to compute the pseudoinverse. $\endgroup$ Commented Aug 17, 2017 at 20:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.