[SciPy-user] Python - Matlab least squares difference

[Robert Kern]

A. Rachel Grinberg wrote:
> Hi,
> 
> I noticed a difference between the linear least square solutions in Python 
> and Matlab. I understand that if the system is underdetermined the solution 
> is not going to be unique, nevertheless I would like figure out the 
> algorithm Matlab is using, since it seems "better" to me. For example, let's 
> say I have a Matrix
>       2 1 0
> A = 1 1 1 and b = (0,1)'
> 
> While Matlab yields (0,0,1)' as the solution to A\b, scipy's result for 
> linalg.linear_least_squares(A,b) is
> 
> array([[-0.16666667],
>        [ 0.33333333],
>        [ 0.83333333]])
> 
> Any ideas, how Matlab's A\b is implemented?

Not sure, but it's wrong (more or less). In underdetermined linear least squares
problems, it is conventional to choose the solution that has minimum L2-norm.

>>> A = array([[2., 1, 0], [1, 1, 1]])
>>> A
array([[ 2.,  1.,  0.],
       [ 1.,  1.,  1.]])
>>> b = array([0., 1.])
>>> linalg.lstsq(A, b)
(array([-0.16666667,  0.33333333,  0.83333333]), array([], dtype=float64), 2,
array([ 2.6762432 ,  0.91527173]))
>>> x = _[0]
>>> dot(x, x)
0.83333333333333426
>>> dot([0., 0, 1], [0., 0, 1])
1.0

Implementations that do something else have some 'splainin' to do. Of course,
A\b may not quite be "do linear least squares" in Matlab but something else. I
don't know. FWIW, Octave gives me the same answer as numpy/scipy.

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
 that is made terrible by our own mad attempt to interpret it as though it had
 an underlying truth."
  -- Umberto Eco