# Poor SVD reconstruction of singular matrix

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.

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.
• 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. – Kirill Aug 16 '17 at 21:29
• Deleted previous comment as my check was incorrect. This is not the case, the normwise error $\frac{|| A^+ - (A^+)^* ||}{|| A^+ ||} \approx 0.01$ – myseun Aug 17 '17 at 10:16