Understanding how Numpy does SVD

Question

I have been using different methods to calculate both the rank of a matrix and the solution of a matrix system of equations. I came across the function linalg.svd. Comparing this to my own effort of solving the system with Gaussian Elimination, it appears to be both faster and more precise. I'm trying to understand how this is possible.

As far as I know, the linalg.svd function uses a QR algorithm to calculate the eigenvalues of my matrix. I know how this works mathematically, but I don't know how Numpy manages to do it so quickly and without losing much precision.

So my question: How does the numpy.svd function work, and more specifically, how does it manage to do it fast and accurately (compared to gaussian elimination)?

Answer 1

There are a number of issues in your question.

Do not use Gaussian Elimination (LU factorization) to calculate the numerical rank of a matrix. LU factorization is unreliable for this purpose in floating-point arithmetic. Instead, use a rank-revealing QR decomposition (such as xGEQPX or xGEPQY in LAPACK, where x is C, D, S, or Z, though those routines are difficult to track down; see JedBrown's answer on a related question), or use an SVD (singular value decomposition, such as xGESDD or xGESVD, where x is again C, D, S, or Z). The SVD is a more accurate, reliable algorithm for the determination of numerical rank, but it requires more floating-point operations.

However, for solving a linear system, LU factorization (with partial pivoting, which is the standard implementation in LAPACK) is extremely reliable in practice. There are some pathological cases for which LU factorization with partial pivoting is unstable (see Lecture 22 in Numerical Linear Algebra by Trefethen and Bau for details). QR factorization is a more stable numerical algorithm for solving linear systems, which is probably why it gives you such precise results. However, it requires more floating-point operations than LU factorization by a factor of 2 for square matrices (I believe; JackPoulson may correct me on that). For rectangular systems, QR factorization is a better choice because it will yield least-squares solutions to overdetermined linear systems. SVD can also be used to solve linear systems, but it will be more expensive than QR factorization.

janneb is correct that numpy.linalg.svd is a wrapper around xGESDD in LAPACK. Singular value decompositions proceed in two stages. First, the matrix to be decomposed is reduced to bidiagonal form. The algorithm used to reduce to bidiagonal form in LAPACK is probably the Lawson-Hanson-Chan algorithm, and it does use QR factorization at one point. Lecture 31 in Numerical Linear Algebra by Trefethen and Bau gives an overview of this process. Then, xGESDD uses a divide-and-conquer algorithm to calculate the singular values and left and right singular vectors from the bidiagonal matrix. To get background on this step, you'll need to consult Matrix Computations by Golub and Van Loan, or Applied Numerical Linear Algebra by Jim Demmel.

Finally, you should not confuse singular values with eigenvalues. These two sets of quantities are not the same. The SVD computes the singular values of a matrix. Cleve Moler's Numerical Computing with MATLAB gives a nice overview of the differences between singular values and eigenvalues. In general, there is no obvious relationship between the singular values of a given matrix and its eigenvalues, except in the case of normal matrices, where the singular values are the absolute value of the eigenvalues.

Answer 2

Due to the wording of your question, I am assuming that your matrix is square. LAPACK's SVD routines, such as zgesvd, essentially proceed in three stages for square matrices:

1. Implicitly computing unitary matrices $U_A$ and $V_A$, as products of Householder transforms, such that the general matrix $A$ is reduced to a real, upper bidiagonal matrix $B := U_A^H A V_A$. On exit of this subroutine, the Householder vectors (normalized so that their first entry is implicitly one) for $U_A$ and $V_A$ are respectively stored in the portions of $B$ below and to the right of the main and super-diagonal. This step requires $O(n^3)$ work.
2. A variation on an algorithm for computing the eigenvalue decomposition of a real symmetric tridiagonal matrix is used to compute a bidiagonal SVD. Though the best-known algorithm for real symmetric tridiagonal EVP's (MRRR) is not yet, to my knowledge, stable for the bidiagonal SVD, there is an interesting discussion here. LAPACK currently uses a divide and conquer approach for the bidiagonal SVD. The bidiagonal SVD yields $\{U_B,V_B,\Sigma\}$ such that $B = U_B \Sigma V_B^H$. This step will require $O(n^2)$ work when MRRR is made stable for SVD, but currently requires as much as $O(n^3)$ work.
3. This step is trivial. Since $U_A B V_A^H = A$, $A = (U_A U_B) \Sigma (V_A V_B)^H$, so that the SVD can be explicitly formed after applying the Householder transforms which implicitly define $U_A$ and $V_A$ to $U_B$ and $V_B$, respectively. This step required $O(n^3)$ work.

Answer 3

numpy.linalg.svd is a wrapper around {Z,D}GESDD from LAPACK. LAPACK, in turn, is very carefully written by some of the world's foremost experts in numerical linear algebra. Indeed, it'd be very surprising if someone not intimately familiar with the field would succeed in beating LAPACK (either in speed or accuracy).

As for why QR is better than Gaussian elimination, that is probably more appropriate for https://scicomp.stackexchange.com/