[Numpy-discussion] About inv/lstsq

[Irvin Probst]

Hi,
I'm not sure if I should send this here or to scipy-user, feel free to 
redirect me there if I'm off topic.

So, there is something I don't understand using inv and lstsq in numpy.

I've built *on purpose* an ill conditioned system to fit a quadric 
a*x**2+b*y**2+c*x*y+d*x+e*y+f, the data points are taken on a narrow 
stripe four times longer than wide. My goal is obviously to find 
(a,b,c,d,e,f) so I built the following matrix:

A[:,0] = data[:,0]**2

A[:,1] = data[:,1]**2

A[:,2] = data[:,1]*data[:,0]

A[:,3] = data[:,0]

A[:,4] = data[:,1]

A[:,5] = 1;


The condition number of A is around 2*1e5 but I can make it much bigger 
if needed by scaling the data along an axis.

I then tried to find the best estimate of X in order to minimize the 
norm of A*X - B with B being my data points and X the vector 
(a,b,c,d,e,f). That's a very basic usage of least squares and it works 
fine with lstsq despite the bad condition number.

However I was expecting to fail to solve it properly using 
inv(A.T.dot(A)).dot(A.T).dot(B) but actually while I scaled up the 
condition number lstsq began to give obviously wrong results (that's 
expected) whereas using inv constantly gave "visually good" results. I 
have no residuals to show but lstsq was just plain wrong (again that is 
expected when cond(A) rises) while inv "worked". I was expecting to see 
inv fail much before lstsq.

Interestingly the same dataset fails in Matlab using inv without any 
scaling of the condition number while it works using \ (mldivide, i.e 
least squares). On octave it works fine using both methods with the 
original dataset, I did not try to scale up the condition number.

So my question is very simple, what's going on here ? It looks like 
Matlab, Numpy and Octave both use the same lapack functions for inv and 
lstsq. As they don't use the same version of lapack I can understand 
that they do not exhibit the same behavior but how can it be possible to 
have lstsq failing before inv(A.T.dot(A)) when I scale up the condition 
number of A ? I feel like I'm missing something obvious but I can't find it.

Thanks.

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20160301/1b8120f8/attachment.html>