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.

Answer 1

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]