<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=utf-8">
</head>
<body text="#000000" bgcolor="#FFFFFF">
Hi,<br>
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.<br>
<br>
So, there is something I don't understand using inv and lstsq in
numpy.<br>
<br>
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:<br>
<br>
<meta name="qrichtext" content="1">
<p style=" margin-top:0px; margin-bottom:0px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px; -qt-user-state:0;"><!--StartFragment-->A[:,0] = data[:,0]**2</p>
<p style=" margin-top:0px; margin-bottom:0px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px; -qt-user-state:0;">A[:,1] = data[:,1]**2</p>
<p style=" margin-top:0px; margin-bottom:0px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px; -qt-user-state:0;">A[:,2] = data[:,1]*data[:,0]</p>
<p style=" margin-top:0px; margin-bottom:0px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px; -qt-user-state:0;">A[:,3] = data[:,0]</p>
<p style=" margin-top:0px; margin-bottom:0px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px; -qt-user-state:0;">A[:,4] = data[:,1]</p>
<p style=" margin-top:0px; margin-bottom:0px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px; -qt-user-state:0;">A[:,5] = 1;<!--EndFragment--></p>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<style type="text/css">
p, li { white-space: pre-wrap; }
</style><br>
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.<br>
<br>
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.<br>
<br>
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.<br>
<br>
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.<br>
<br>
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.<br>
<br>
Thanks.<br>
<br>
</body>
</html>