Friday, March 16, 2012

SVD decomposition with numpy

The SVD decomposition is a factorization of a matrix, with many useful applications in signal processing and statistics. In this post we will see
  • how to compute the SVD decomposition of a matrix A using numpy,
  • how to compute the inverse of A using the matrices computed by the decomposition,
  • and how to solve a linear equation system Ax=b using using the SVD.
The SVD decomposition of a matrix A is of the fom


Since U and V are orthogonal (this means that U^T*U=I and V^T*V=I) we can write the inverse of A as (see Solving overdetermined systems with the QR decomposition for the tricks)


So, let's start computing the factorization and the inverse
from numpy import *

A = floor(random.rand(4,4)*20-10) # generating a random
b = floor(random.rand(4,1)*20-10) # system Ax=b

U,s,V = linalg.svd(A) # SVD decomposition of A

# computing the inverse using pinv
pinv = linalg.pinv(A)
# computing the inverse using the SVD decomposition
pinv_svd = dot(dot(V.T,linalg.inv(diag(s))),U.T)

print "Inverse computed by lingal.pinv()\n",pinv
print "Inverse computed using SVD\n",pinv_svd
As we can see, the output shows that pinv and pinv_svd are the equal
Inverse computed by lingal.pinv()
[[ 0.06578301 -0.04663721  0.0436917   0.089838  ]
 [ 0.15243004  0.044919   -0.03681885  0.00294551]
 [ 0.18213058 -0.01718213  0.06872852 -0.07216495]
 [ 0.03976436  0.09867452  0.03387334 -0.04270987]]
Inverse computed using SVD
[[ 0.06578301 -0.04663721  0.0436917   0.089838  ]
 [ 0.15243004  0.044919   -0.03681885  0.00294551]
 [ 0.18213058 -0.01718213  0.06872852 -0.07216495]
 [ 0.03976436  0.09867452  0.03387334 -0.04270987]]
Now, we can solve Ax=b using the inverse:


or solving the system

Multiplying by U^T we obtain


Then, letting c be U^Tb and w be V^Tx, we see


Since sigma is diagonal, we can easily obtain w solving the system above. And finally, we can obtain x solving

Let's compare the results of those methods:

x = linalg.solve(A,b) # solve Ax=b using linalg.solve

xPinv = dot(pinv_svd,b) # solving Ax=b computing x = A^-1*b

# solving Ax=b using the equation above
c = dot(U.T,b) # c = U^t*b
w = linalg.solve(diag(s),c) # w = V^t*c
xSVD = dot(V.T,w) # x = V*w

print "Ax=b solutions compared"
print x.T
print xSVD.T
print xPinv.T
As aspected, we have the same solutions:
Ax=b solutions compared
[[ 0.13549337 -0.37260677  1.62886598 -0.09720177]]
[[ 0.13549337 -0.37260677  1.62886598 -0.09720177]]
[[ 0.13549337 -0.37260677  1.62886598 -0.09720177]]

4 comments:

  1. The whole point of using the SVD to compute the solution to Ax = b is that you may have an ill-conditioned or even rank-deficient matrix - meaning the smallest singular are very small or zero, respectively. These types of matrices arise all the time in practice (for example, you'd run into this in your example if you used rand(100,100) instead of rand(4,4)). In these cases, you can truncate those small singular values (and corresponding columns of U and V) and the SVD lets you compute the pseudo-inverse.

    Solving a system by computing the inverse - which doesn't exist for a rank-deficient matrix, and is very inaccurate for a ill-conditioned matrix - is a very poor numerical method. You should mention this. You say that, as "aspected we have the same solutions," but in fact if you took the norm of the difference of these two solutions versus x, you would find the truncated xSVD solution was much more accurate.

    ReplyDelete
  2. Hi,

    I've tried your method on a set of equations of mine which looks like this :

    [[ 0 2 7 8 16 17]
    [ 2 0 5 6 17 16]
    [ 7 5 0 9 15 13]
    [ 8 6 9 0 18 16]
    [16 17 15 18 0 5]
    [17 16 13 16 5 0]]

    and i have :
    onesarr=numpy.ones(shape=(len(A),1))
    u, s, vh = linalg.svd(A)

    c = dot(u.T,onesarr) # c = U^t*b
    w = linalg.solve(diag(s),c) # w = V^t*c
    xSVD = dot(vh.T,w) # x = V*w

    pinv_svd = dot(dot(vh.T,linalg.inv(diag(s))),u.T)
    xPinv = dot(pinv_svd,onesarr) # solving Ax=b computing x = A^-1*b
    pinv = linalg.pinv(A)
    #print "Inverse computed by lingal.pinv()\n",pinv

    x = linalg.solve(A,onesarr)

    print x.T
    print xSVD.T
    print xPinv.T

    and what i'm getting is something like that :

    [[ 0.02525855 -0.00949907 0.01138364 0.02771869 0.0262207 0.01753129]]
    [[ 0.02525855 -0.00949907 0.01138364 0.02771869 0.0262207 0.01753129]]
    [[ 0.02525855 -0.00949907 0.01138364 0.02771869 0.0262207 0.01753129]]

    but the thing is that it's weird that i have negative number, because as my data suggest that i should not have negative numbers, am i doing something wrong in the calculation or there is something is the basic idea of generating the equations ??

    ReplyDelete
  3. I am reforming Linear Algebra to make the SVD central and accessible early in undergrad courses. Then many more significant modelling and application situations are more immediately available, as well as better theory. The book would be an excellent answer to all these queries. I have completed the first draft: download via
    https://raw.github.com/uoa1184615/LinearAlgebraGit/master/larxxia-a1a.pdf

    ReplyDelete
  4. Is there a way to label the variables for SVD?
    Using Python is great for SVD if you have a plain matrix.
    But how can I employ SVD to determine which values I would keep for a regression? Baseline SVD for Python looks pretty useless, since I just get a bunch of unlabeled numbers for my results.

    ReplyDelete

Note: Only a member of this blog may post a comment.