Poor SVD reconstruction of singular matrix - Computational Science Stack Exchange most recent 30 from scicomp.stackexchange.com 2019-12-10T22:17:35Z https://scicomp.stackexchange.com/feeds/question/27641 https://creativecommons.org/licenses/by-sa/4.0/rdf https://scicomp.stackexchange.com/q/27641 4 Poor SVD reconstruction of singular matrix myseun https://scicomp.stackexchange.com/users/21311 2017-08-16T14:28:22Z 2017-08-16T16:04:34Z <p>I am trying to calculate the singular value decomposition of <a href="https://pastebin.com/takPtiVD" rel="nofollow noreferrer">this matrix</a> using <code>numpy.linalg.svd</code> .</p> <p>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.</p> <pre><code>u, w, vt = np.linalg.svd(M) np.allclose(M, np.dot(u, np.dot(np.diag(w), vt))) &gt;&gt; False </code></pre> <p>The matrix is singular ( <code>w[-1]</code> is zero), and the largest singular value is much bigger than the other non-zero values (<code>3e+24</code> vs <code>7e-2</code> for the next largest).</p> <p>-</p> <p>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?</p> <p>-</p> <p>The reason I am computing the SVD is a diagnosis for the pseudoinverse of <code>M</code> I calculate using <code>numpy.linalg.pinv</code>, which should return the Moore-Penrose pseudoinverse. As I understand it, this particular pseudoinverse should be symmetric since <code>M</code> is symmetric, but it is not. My assumption is that this is because of problems with the SVD, which <code>numpy.linalg.pinv</code> calculates.</p> https://scicomp.stackexchange.com/questions/27641/-/27642#27642 9 Answer by Federico Poloni for Poor SVD reconstruction of singular matrix Federico Poloni https://scicomp.stackexchange.com/users/4405 2017-08-16T16:04:34Z 2017-08-16T16:04:34Z <p>Algorithms for the SVD, as more or less every classical linear algebra algorithm based on orthogonal transformations, are <em>normwise</em> 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.</p> <p>In particular, this bound in norm does <em>not</em> 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.</p> <p><a href="https://docs.scipy.org/doc/numpy/reference/generated/numpy.allclose.html" rel="noreferrer"><code>numpy.allclose</code></a> compares elementwise, not normwise.</p>