[SciPy-Dev] Slight change in lstsq accuracy for 0.17 - expected?

[Matthew Brett]

Hi,

On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers <ralf.gommers at gmail.com> wrote:
>
>
> On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett <matthew.brett at gmail.com>
> wrote:
>>
>> Hi,
>>
>> I just noticed that our (neuroimaging) test suite fails for scipy 0.17.0:
>>
>> https://travis-ci.org/matthew-brett/nipy/jobs/115471563
>>
>> This turns out to be due to small changes in the result of
>> `scipy.linalg.lstsq` between 0.16.1 and 0.17.0 - on the same platform
>> - the travis Ubuntu in this case.
>>
>> Is that expected?   I haven't investigated deeply, but I believe the
>> 0.16.1 result (back into the depths of time) is slightly more
>> accurate.  Is that possible?
>
>
> That's very well possible, because we changed the underlying LAPACK routine
> that's used:
> https://github.com/scipy/scipy/blob/master/doc/release/0.17.0-notes.rst#scipy-linalg-improvements.
> This PR also had some lstsq precision issues, so could be useful:
> https://github.com/scipy/scipy/pull/5550
>
> What is the accuracy difference in your case?

Damn - I thought you might ask me that!

The test where it showed up comes from some calculations on a simple
regression problem.

The regression model is $y = X B + e$, with ordinary least squares fit
given by $B = X^+ y$ where $X^+$ is the pseudoinverse of matrix X.

For my particular X, $\Sigma e^2$ (SSE) is smaller using scipy 0.16
pinv, compared to scipy 0.17 pinv:

# SSE using high-precision sympy pinv
72.0232781366615313
# SSE using scipy 0.16 pinv
72.0232781366625261
# SSE using scipy 0.17 pinv
72.0232781390480170

See: https://github.com/matthew-brett/scipy-pinv/blob/master/pinv_futz.py

A difference at the 9th decimal place - but this was enough to give a
difference at the 5th decimal place in a more involved calculation
using the same pinv result:

https://github.com/nipy/nipy/pull/394/files

Does that help?

Cheers,

Matthew