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)?
Understanding how Numpy does SVD
linear algebra;python;lapack
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.