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.
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
https://mail.scipy.org/mailman/listinfo/numpy-discussion
