4
$\begingroup$

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.

$\endgroup$
9
$\begingroup$

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.

$\endgroup$
  • 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 Aug 16 '17 at 21:29
  • $\begingroup$ Thanks for the information, seems like my matrix is badly behaved so have accepted this answer. $\endgroup$ – myseun Aug 16 '17 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 Aug 17 '17 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$ – Federico Poloni Aug 17 '17 at 20:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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