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.
dgesdd
for real valued SVDs. So your real question is probably "how does Lapack dgesdd work?", and that is pretty off topic for stackoverflow. $\endgroup$