Solve Singular Value Decomposition (SVD) in Python

Question

I amtrying to translate an IDL program to Python. I have to solve the outcome from SVD which I achieve in the following way

from scipy.linalg import svd
A = [[1,2,3],[4,5,6]]
b = [4,4,5]

u,w,v = svd(A)

And this works fine and is translated nicely from IDL. The next step is IN IDL(!)

x = svsol(u,w,v,b)

The u in python and IDL are almost the same (and for the other matrix's as well). The only difference is the dimensions, where IDL's matrix's is larger, but has a lot of zeros. It looks like Python's matrix's are more compressed in that sence.

Does anyone know something similarly for Python.

If anyone needs it, here is the manual for svsol.

Can't you just use: [`scipy.linalg.solve(A,b)`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.solve.html)? — Andy Hayden, Sep 25 '12 at 10:01
No, I am looking for something similar to `svsol` in IDL. The method I am following has to use this way. — Daniel Thaagaard Andreasen, Sep 25 '12 at 10:24

score 4 · Accepted Answer · answered Oct 04 '12 at 21:21

With SVDC and SVSOL in IDL you solve a linear least squares problem by SVD decomposition. This is done in numpy by the numpy.linalg.lstsq function. (No need to compute first the SVD decomposition and then back solve.)

>>> import numpy as np
>>> A = np.array([[1,2,3],[4,5,6]])
>>> b = np.array([4,4])
>>> x, _, _, _ = np.linalg.lstsq(A,b)
>>> x
array([-2.,  0.,  2.])
>>> np.dot(A,x)
array([ 4.,  4.])

Please note that the length of b has to be the same as the number of rows of A, so your example is wrong. Just to make shure that I'm correctly interpreting the IDL semantics, here is the example in the svsol reference manual:

>>> A = np.array(
... [[1.0, 2.0, -1.0, 2.5],
...  [1.5, 3.3, -0.5, 2.0],
...  [3.1, 0.7,  2.2, 0.0],
...  [0.0, 0.3, -2.0, 5.3],
...  [2.1, 1.0,  4.3, 2.2],
...  [0.0, 5.5,  3.8, 0.2]])
>>> B = np.array([0.0, 1.0, 5.3, -2.0, 6.3, 3.8])
>>> x, _, _, _ = np.linalg.lstsq(A,B)
>>> print x
[ 1.00095058  0.00881193  0.98417587 -0.01009547]

That seems to work really great. I think I'll just use the least square method then. Just of curiosity, can you tell my the advantages of SVD (if there is any) over `lstsq`. Thank you really much for your help; it is very preciated! — Daniel Thaagaard Andreasen, Oct 05 '12 at 09:27
@DanielThaagaardAndreasen internally `np.linalg.lstsq` calls LAPACK [`dgelsd`](http://www.netlib.org/lapack/double/dgelsd.f) which in turns is SVD based, so conceptually there are no differences. This said usually `dgelsd` turns out to be faster than full SVD and subsequent back solve. For more info have a a look at the LAPACK [user guide](http://www.netlib.org/lapack/lug/node27.html) — Stefano M, Oct 05 '12 at 11:13

Solve Singular Value Decomposition (SVD) in Python

1 Answers1