13
$\begingroup$

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)?

$\endgroup$

migrated from stackoverflow.com Apr 4 '12 at 19:26

This question came from our site for professional and enthusiast programmers.

  • 2
    $\begingroup$ numpy uses the Lapack routine 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$ – talonmies Apr 4 '12 at 17:50
  • $\begingroup$ If your REALLY curious, I would suggest examining the LAPACK source. $\endgroup$ – Joel Cornett Apr 4 '12 at 18:20
  • $\begingroup$ Thank you for your comments, and my apologies I'm offtopic. $\endgroup$ – RobVerheyen Apr 4 '12 at 18:31
  • $\begingroup$ This post is a cross-post from Stack Overflow. Cross-posting is generally discouraged on Stack Exchange sites. The standard protocol for reposting a question on a different site is to close, delete, or migrate the original post before attempting to repost on a different site. (If you migrate the question, it is reposted automatically.) $\endgroup$ – Geoff Oxberry Apr 4 '12 at 19:22
  • $\begingroup$ I'm sorry, I was not aware of the protocol. I hope I can still get an answer. $\endgroup$ – RobVerheyen Apr 4 '12 at 21:00
15
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I think "not related" is pretty strong for the relationship between eigenvalues and singular values. The relationship is pretty obscure unless you know the full Jordan decomposition of your matrix, but you can use one to get estimates of the other if you have information (or are willing to make assumptions) about said Jordan decomposition. $\endgroup$ – Dan Apr 4 '12 at 22:28
  • $\begingroup$ What would you suggest instead? $\endgroup$ – Geoff Oxberry Apr 4 '12 at 23:54
  • $\begingroup$ First of all, thank you for the elaborate answer. I found out I cannot use LU decomposition to determine matrix rank the hard way. Your answer seems to imply that the QR factorisation would in fact be a faster method of solving my problem, correct? Is there a distinct advantage in using SVD? I was well aware of the fact that singular values are not eigenvalues. I was referring to the fact that singular values can be calculated as eigenvalues of the matrix multiplied with it's transpose from the left. I'm sorry that wasn't clear. $\endgroup$ – RobVerheyen Apr 5 '12 at 7:14
  • $\begingroup$ I might add that the matrix I'm solving is actually singular. In fact, the matrix rank is only about half of the size of the matrix. Perhaps this makes some method more preferrable? $\endgroup$ – RobVerheyen Apr 5 '12 at 7:21
  • 1
    $\begingroup$ @RobVerheyen: QR will be slower than LU would be, but considerably more accurate. SVD will be even slower than QR, but SVD is considered the most reliable method for determining numerical rank (for instance, MATLAB uses the SVD in its rank function). There's also a little discretion involved when using either approach; in the SVD approach, the numerical rank is the number of singular values above a specified (usually very small) cutoff. (The QR approach is similar, but replaces singular values with diagonal entries of the R matrix.) $\endgroup$ – Geoff Oxberry Apr 5 '12 at 17:10
8
$\begingroup$

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.
$\endgroup$
7
$\begingroup$

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/

$\endgroup$
  • $\begingroup$ Thank you for the answer and the reference. I'll try it over there. $\endgroup$ – RobVerheyen Apr 4 '12 at 18:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.