Slight change in lstsq accuracy for 0.17 - expected?
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? Thanks for any insight, Matthew
On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett <matthew.brett@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#scip.... 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? Ralf
Thanks for any insight,
Matthew _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org https://mail.scipy.org/mailman/listinfo/scipy-dev
Hi, On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers <ralf.gommers@gmail.com> wrote:
On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett <matthew.brett@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#scip.... 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
On Wed, Mar 16, 2016 at 6:05 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers <ralf.gommers@gmail.com> wrote:
On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett <matthew.brett@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#scip... .
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
That difference looks pretty big, looks almost like a cutoff effect. What is the matrix X, in particular, what is its condition number or singular values? Might also be interesting to see residuals for the different functions. <snip> Chuck
On Wed, Mar 16, 2016 at 6:41 PM, Charles R Harris <charlesr.harris@gmail.com
wrote:
On Wed, Mar 16, 2016 at 6:05 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers <ralf.gommers@gmail.com> wrote:
On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett <matthew.brett@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#scip... .
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
That difference looks pretty big, looks almost like a cutoff effect. What is the matrix X, in particular, what is its condition number or singular values? Might also be interesting to see residuals for the different functions.
<snip>
Note that this might be an indication that the result of the calculation has an intrinsically large error bar. Chuck
Hi, On Wed, Mar 16, 2016 at 5:48 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 6:41 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 6:05 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers <ralf.gommers@gmail.com> wrote:
On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett <matthew.brett@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#scip.... 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
That difference looks pretty big, looks almost like a cutoff effect. What is the matrix X, in particular, what is its condition number or singular values? Might also be interesting to see residuals for the different functions.
<snip>
Note that this might be an indication that the result of the calculation has an intrinsically large error bar.
In [6]: np.linalg.svd(X, compute_uv=False) Out[6]: array([ 6.97301629e+07, 5.78702921e+05, 5.39646908e+05, 2.99383626e+05, 3.24329208e+04, 1.29520746e+02, 3.16197688e+01, 2.66011711e+00, 2.18287258e-01, 1.60105276e-01, 1.11807228e-01, 2.12694535e-03, 1.72675581e-03, 9.51576795e-04]) In [7]: np.linalg.cond(X) Out[7]: 73278544935.300507 I wasn't sure what you meant about residuals from different functions - can you explain more? Cheers, Matthew
On Wed, Mar 16, 2016 at 6:16 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Wed, Mar 16, 2016 at 5:48 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 6:41 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 6:05 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers <ralf.gommers@gmail.com> wrote:
On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett <matthew.brett@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#scip.... 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
That difference looks pretty big, looks almost like a cutoff effect. What is the matrix X, in particular, what is its condition number or singular values? Might also be interesting to see residuals for the different functions.
<snip>
Note that this might be an indication that the result of the calculation has an intrinsically large error bar.
In [6]: np.linalg.svd(X, compute_uv=False) Out[6]: array([ 6.97301629e+07, 5.78702921e+05, 5.39646908e+05, 2.99383626e+05, 3.24329208e+04, 1.29520746e+02, 3.16197688e+01, 2.66011711e+00, 2.18287258e-01, 1.60105276e-01, 1.11807228e-01, 2.12694535e-03, 1.72675581e-03, 9.51576795e-04])
In [7]: np.linalg.cond(X) Out[7]: 73278544935.300507
Although, the differences in residual error are at the 12th decimal place comparing the sum of squared error for a full-rank vs not-full-rank matrix: In [2]: import scipy.linalg as spl In [3]: X_redundant = np.c_[X, X[:, 0]] In [5]: np.linalg.cond(X_redundant) Out[5]: 4.6338774550388212e+23 In [8]: print_scalar(sse(X, spl.pinv(X), y)) 72.0232781366634640 In [9]: print_scalar(sse(X_redundant, spl.pinv(X_redundant), y)) 72.0232781366614461 Cheers, Matthew
On Wed, Mar 16, 2016 at 7:16 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Wed, Mar 16, 2016 at 5:48 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 6:41 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 6:05 PM, Matthew Brett <matthew.brett@gmail.com
wrote:
Hi,
On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers <ralf.gommers@gmail.com> wrote:
On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett <matthew.brett@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#scip... .
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
That difference looks pretty big, looks almost like a cutoff effect. What is the matrix X, in particular, what is its condition number or singular values? Might also be interesting to see residuals for the different functions.
<snip>
Note that this might be an indication that the result of the calculation has an intrinsically large error bar.
In [6]: np.linalg.svd(X, compute_uv=False) Out[6]: array([ 6.97301629e+07, 5.78702921e+05, 5.39646908e+05, 2.99383626e+05, 3.24329208e+04, 1.29520746e+02, 3.16197688e+01, 2.66011711e+00, 2.18287258e-01, 1.60105276e-01, 1.11807228e-01, 2.12694535e-03, 1.72675581e-03, 9.51576795e-04])
In [7]: np.linalg.cond(X) Out[7]: 73278544935.300507
Condition number is about 7e10, so you might expect to lose about 10 digits of of accuracy just due to roundoff. The covariance of the fitting error will be `inv(X.T, X) * error`, so you can look at the diagonals for the variance of the fitted parameters just for kicks, taking `error` to be roundoff error or measurement error. This is not to say that the results should not be more reproducible, but should be of interest for the validity of the solution. Note that I suspect that there may be a difference in the default `rcond` used, or maybe the matrix is on the edge and a small change in results changes the effective rank. You can run with the `return_rank=True` argument to check that. Looks like machine precision is used for the default value of rcond, so experimenting with that probable won't be informative, but you might want to try it anyway. Chuck
On Wed, Mar 16, 2016 at 11:01 PM, Charles R Harris < charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 7:16 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Wed, Mar 16, 2016 at 5:48 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 6:41 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 6:05 PM, Matthew Brett <
wrote:
Hi,
On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers <ralf.gommers@gmail.com
wrote:
On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett <matthew.brett@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
> - the travis Ubuntu in this case. > > Is that expected? I haven't investigated deeply, but I believe
matthew.brett@gmail.com> platform 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#scip... .
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
That difference looks pretty big, looks almost like a cutoff effect. What is the matrix X, in particular, what is its condition number or singular values? Might also be interesting to see residuals for the different functions.
<snip>
Note that this might be an indication that the result of the calculation has an intrinsically large error bar.
In [6]: np.linalg.svd(X, compute_uv=False) Out[6]: array([ 6.97301629e+07, 5.78702921e+05, 5.39646908e+05, 2.99383626e+05, 3.24329208e+04, 1.29520746e+02, 3.16197688e+01, 2.66011711e+00, 2.18287258e-01, 1.60105276e-01, 1.11807228e-01, 2.12694535e-03, 1.72675581e-03, 9.51576795e-04])
In [7]: np.linalg.cond(X) Out[7]: 73278544935.300507
Condition number is about 7e10, so you might expect to lose about 10 digits of of accuracy just due to roundoff. The covariance of the fitting error will be `inv(X.T, X) * error`, so you can look at the diagonals for the variance of the fitted parameters just for kicks, taking `error` to be roundoff error or measurement error. This is not to say that the results should not be more reproducible, but should be of interest for the validity of the solution. Note that I suspect that there may be a difference in the default `rcond` used, or maybe the matrix is on the edge and a small change in results changes the effective rank. You can run with the `return_rank=True` argument to check that. Looks like machine precision is used for the default value of rcond, so experimenting with that probable won't be informative, but you might want to try it anyway.
statsmodels is using numpy pinv for the linear model and it works in almost all cases very well *). Using pinv with singular (given the rcond threshold) X makes everything well defined and not subject to numerical noise because of the regularization induced by pinv. Our covariance for the parameter estimates, based on pinv(X), is also not noisy but it's of reduced rank. The main problem I have seen so far is that the default rcond is too small for some cases, and it breaks down with near singular X with rcond a bit larger than 1e-15 in an example with dummy variables. https://github.com/statsmodels/statsmodels/issues/2628 http://stackoverflow.com/questions/32599159/robustness-issue-of-statsmodel-l... If I think a bit longer, then I remember two other cases where pinv with ill conditioned X also doesn't work very well with continuous almost collinear data, one NIST example with badly scaled polynomials where precision goes down to a few decimals and a very small economic dataset where the normal equations (or orthogonality conditions?) are not satisfied at a good precision.. Josef
Chuck
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org https://mail.scipy.org/mailman/listinfo/scipy-dev
On Wed, Mar 16, 2016 at 9:39 PM, <josef.pktd@gmail.com> wrote:
On Wed, Mar 16, 2016 at 11:01 PM, Charles R Harris < charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 7:16 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Wed, Mar 16, 2016 at 5:48 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 6:41 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 6:05 PM, Matthew Brett <
wrote:
Hi,
On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers <
ralf.gommers@gmail.com>
wrote: > > > On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett > <matthew.brett@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
>> - the travis Ubuntu in this case. >> >> Is that expected? I haven't investigated deeply, but I believe
matthew.brett@gmail.com> platform 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#scip... . > 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
That difference looks pretty big, looks almost like a cutoff effect. What is the matrix X, in particular, what is its condition number or singular values? Might also be interesting to see residuals for the different functions.
<snip>
Note that this might be an indication that the result of the calculation has an intrinsically large error bar.
In [6]: np.linalg.svd(X, compute_uv=False) Out[6]: array([ 6.97301629e+07, 5.78702921e+05, 5.39646908e+05, 2.99383626e+05, 3.24329208e+04, 1.29520746e+02, 3.16197688e+01, 2.66011711e+00, 2.18287258e-01, 1.60105276e-01, 1.11807228e-01, 2.12694535e-03, 1.72675581e-03, 9.51576795e-04])
In [7]: np.linalg.cond(X) Out[7]: 73278544935.300507
Condition number is about 7e10, so you might expect to lose about 10 digits of of accuracy just due to roundoff. The covariance of the fitting error will be `inv(X.T, X) * error`, so you can look at the diagonals for the variance of the fitted parameters just for kicks, taking `error` to be roundoff error or measurement error. This is not to say that the results should not be more reproducible, but should be of interest for the validity of the solution. Note that I suspect that there may be a difference in the default `rcond` used, or maybe the matrix is on the edge and a small change in results changes the effective rank. You can run with the `return_rank=True` argument to check that. Looks like machine precision is used for the default value of rcond, so experimenting with that probable won't be informative, but you might want to try it anyway.
statsmodels is using numpy pinv for the linear model and it works in almost all cases very well *).
Using pinv with singular (given the rcond threshold) X makes everything well defined and not subject to numerical noise because of the regularization induced by pinv. Our covariance for the parameter estimates, based on pinv(X), is also not noisy but it's of reduced rank.
What I'm thinking is that changing the algorithm, as here, is like introducing noise relative to the original, so even if the result is perfectly valid, the fitted parameters might change significantly. A change in the 11-th digit wouldn't surprise me at all, depending on the actual X and rhs.
The main problem I have seen so far is that the default rcond is too small for some cases, and it breaks down with near singular X with rcond a bit larger than 1e-15 in an example with dummy variables.
https://github.com/statsmodels/statsmodels/issues/2628
http://stackoverflow.com/questions/32599159/robustness-issue-of-statsmodel-l...
If I think a bit longer, then I remember two other cases where pinv with ill conditioned X also doesn't work very well with continuous almost collinear data, one NIST example with badly scaled polynomials where precision goes down to a few decimals and a very small economic dataset where the normal equations (or orthogonality conditions?) are not satisfied at a good precision..
Yeah, it can make a real mess with polynomial fits ;) Chuck
On Wed, Mar 16, 2016 at 9:49 PM, Charles R Harris <charlesr.harris@gmail.com
wrote:
On Wed, Mar 16, 2016 at 9:39 PM, <josef.pktd@gmail.com> wrote:
On Wed, Mar 16, 2016 at 11:01 PM, Charles R Harris < charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 7:16 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Wed, Mar 16, 2016 at 5:48 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 6:41 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 6:05 PM, Matthew Brett <
wrote: > > Hi, > > On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers < ralf.gommers@gmail.com> > wrote: > > > > > > On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett > > <matthew.brett@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
> >> - the travis Ubuntu in this case. > >> > >> Is that expected? I haven't investigated deeply, but I believe
matthew.brett@gmail.com> platform 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#scip... . > > 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
That difference looks pretty big, looks almost like a cutoff effect. What is the matrix X, in particular, what is its condition number or singular values? Might also be interesting to see residuals for the different functions.
<snip>
Note that this might be an indication that the result of the calculation has an intrinsically large error bar.
In [6]: np.linalg.svd(X, compute_uv=False) Out[6]: array([ 6.97301629e+07, 5.78702921e+05, 5.39646908e+05, 2.99383626e+05, 3.24329208e+04, 1.29520746e+02, 3.16197688e+01, 2.66011711e+00, 2.18287258e-01, 1.60105276e-01, 1.11807228e-01, 2.12694535e-03, 1.72675581e-03, 9.51576795e-04])
In [7]: np.linalg.cond(X) Out[7]: 73278544935.300507
Condition number is about 7e10, so you might expect to lose about 10 digits of of accuracy just due to roundoff. The covariance of the fitting error will be `inv(X.T, X) * error`, so you can look at the diagonals for the variance of the fitted parameters just for kicks, taking `error` to be roundoff error or measurement error. This is not to say that the results should not be more reproducible, but should be of interest for the validity of the solution. Note that I suspect that there may be a difference in the default `rcond` used, or maybe the matrix is on the edge and a small change in results changes the effective rank. You can run with the `return_rank=True` argument to check that. Looks like machine precision is used for the default value of rcond, so experimenting with that probable won't be informative, but you might want to try it anyway.
statsmodels is using numpy pinv for the linear model and it works in almost all cases very well *).
Using pinv with singular (given the rcond threshold) X makes everything well defined and not subject to numerical noise because of the regularization induced by pinv. Our covariance for the parameter estimates, based on pinv(X), is also not noisy but it's of reduced rank.
What I'm thinking is that changing the algorithm, as here, is like introducing noise relative to the original, so even if the result is perfectly valid, the fitted parameters might change significantly. A change in the 11-th digit wouldn't surprise me at all, depending on the actual X and rhs.
The main problem I have seen so far is that the default rcond is too small for some cases, and it breaks down with near singular X with rcond a bit larger than 1e-15 in an example with dummy variables.
https://github.com/statsmodels/statsmodels/issues/2628
http://stackoverflow.com/questions/32599159/robustness-issue-of-statsmodel-l...
If I think a bit longer, then I remember two other cases where pinv with ill conditioned X also doesn't work very well with continuous almost collinear data, one NIST example with badly scaled polynomials where precision goes down to a few decimals and a very small economic dataset where the normal equations (or orthogonality conditions?) are not satisfied at a good precision..
Yeah, it can make a real mess with polynomial fits ;)
The classic rant on these matters is in Numerical Methods that Work. <http://www.amazon.com/Numerical-Methods-that-Work-Spectrum/dp/0883854503> Chuck
On Wed, Mar 16, 2016 at 8:49 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 9:39 PM, <josef.pktd@gmail.com> wrote:
On Wed, Mar 16, 2016 at 11:01 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 7:16 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Wed, Mar 16, 2016 at 5:48 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 6:41 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 6:05 PM, Matthew Brett <matthew.brett@gmail.com> wrote: > > Hi, > > On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers > <ralf.gommers@gmail.com> > wrote: > > > > > > On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett > > <matthew.brett@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#scip.... > > 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
That difference looks pretty big, looks almost like a cutoff effect. What is the matrix X, in particular, what is its condition number or singular values? Might also be interesting to see residuals for the different functions.
<snip>
Note that this might be an indication that the result of the calculation has an intrinsically large error bar.
In [6]: np.linalg.svd(X, compute_uv=False) Out[6]: array([ 6.97301629e+07, 5.78702921e+05, 5.39646908e+05, 2.99383626e+05, 3.24329208e+04, 1.29520746e+02, 3.16197688e+01, 2.66011711e+00, 2.18287258e-01, 1.60105276e-01, 1.11807228e-01, 2.12694535e-03, 1.72675581e-03, 9.51576795e-04])
In [7]: np.linalg.cond(X) Out[7]: 73278544935.300507
Condition number is about 7e10, so you might expect to lose about 10 digits of of accuracy just due to roundoff. The covariance of the fitting error will be `inv(X.T, X) * error`, so you can look at the diagonals for the variance of the fitted parameters just for kicks, taking `error` to be roundoff error or measurement error. This is not to say that the results should not be more reproducible, but should be of interest for the validity of the solution. Note that I suspect that there may be a difference in the default `rcond` used, or maybe the matrix is on the edge and a small change in results changes the effective rank. You can run with the `return_rank=True` argument to check that. Looks like machine precision is used for the default value of rcond, so experimenting with that probable won't be informative, but you might want to try it anyway.
statsmodels is using numpy pinv for the linear model and it works in almost all cases very well *).
Using pinv with singular (given the rcond threshold) X makes everything well defined and not subject to numerical noise because of the regularization induced by pinv. Our covariance for the parameter estimates, based on pinv(X), is also not noisy but it's of reduced rank.
What I'm thinking is that changing the algorithm, as here, is like introducing noise relative to the original, so even if the result is perfectly valid, the fitted parameters might change significantly. A change in the 11-th digit wouldn't surprise me at all, depending on the actual X and rhs.
I believe the errors $e = y - X^+ y$ are in general robust to differences between valid pinv solutions, even in extreme cases: In [2]: print_scalar(sse(X, npl.pinv(X), y)) 72.0232781366616024 In [3]: print_scalar(sse(np.c_[X, X], npl.pinv(np.c_[X, X]), y)) 72.0232781366939463 In [4]: npl.cond(np.c_[X, X]) Out[4]: 6.6247965335098062e+26 Cheers, Matthew
On Fri, Mar 18, 2016 at 3:07 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
On Wed, Mar 16, 2016 at 8:49 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 9:39 PM, <josef.pktd@gmail.com> wrote:
On Wed, Mar 16, 2016 at 11:01 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 7:16 PM, Matthew Brett <
wrote:
Hi,
On Wed, Mar 16, 2016 at 5:48 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 6:41 PM, Charles R Harris <charlesr.harris@gmail.com> wrote: > > > > On Wed, Mar 16, 2016 at 6:05 PM, Matthew Brett > <matthew.brett@gmail.com> > wrote: >> >> Hi, >> >> On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers >> <ralf.gommers@gmail.com> >> wrote: >> > >> > >> > On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett >> > <matthew.brett@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#scip... . >> > 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 > > > That difference looks pretty big, looks almost like a cutoff effect. > What > is the matrix X, in particular, what is its condition number or > singular > values? Might also be interesting to see residuals for the different > functions. > > <snip> >
Note that this might be an indication that the result of the calculation has an intrinsically large error bar.
In [6]: np.linalg.svd(X, compute_uv=False) Out[6]: array([ 6.97301629e+07, 5.78702921e+05, 5.39646908e+05, 2.99383626e+05, 3.24329208e+04, 1.29520746e+02, 3.16197688e+01, 2.66011711e+00, 2.18287258e-01, 1.60105276e-01, 1.11807228e-01, 2.12694535e-03, 1.72675581e-03, 9.51576795e-04])
In [7]: np.linalg.cond(X) Out[7]: 73278544935.300507
Condition number is about 7e10, so you might expect to lose about 10 digits of of accuracy just due to roundoff. The covariance of the fitting error will be `inv(X.T, X) * error`, so you can look at the diagonals for the variance of the fitted parameters just for kicks, taking `error` to be roundoff error or measurement error. This is not to say that the results should not be more reproducible, but should be of interest for the validity of the solution. Note that I suspect that there may be a difference in the default `rcond` used, or maybe the matrix is on the edge and a small change in results changes the effective rank. You can run with the `return_rank=True` argument to check that. Looks like machine
used for the default value of rcond, so experimenting with that
matthew.brett@gmail.com> precision is probable
won't be informative, but you might want to try it anyway.
statsmodels is using numpy pinv for the linear model and it works in almost all cases very well *).
Using pinv with singular (given the rcond threshold) X makes everything well defined and not subject to numerical noise because of the regularization induced by pinv. Our covariance for the parameter estimates, based on pinv(X), is also not noisy but it's of reduced rank.
What I'm thinking is that changing the algorithm, as here, is like introducing noise relative to the original, so even if the result is perfectly valid, the fitted parameters might change significantly. A change in the 11-th digit wouldn't surprise me at all, depending on the actual X and rhs.
I believe the errors $e = y - X^+ y$ are in general robust to differences between valid pinv solutions, even in extreme cases:
In [2]: print_scalar(sse(X, npl.pinv(X), y)) 72.0232781366616024
In [3]: print_scalar(sse(np.c_[X, X], npl.pinv(np.c_[X, X]), y)) 72.0232781366939463
In [4]: npl.cond(np.c_[X, X]) Out[4]: 6.6247965335098062e+26
I think this does not depend on the smallest eigenvalues which should be regularized away in cases like this. The potential problem is with eigenvalues that are just a but too large for rcond to remove them but still have numerical problems. I think it's a similar issue about the threshold as the one that you raised some time ago with determining the matrix rank for near singular matrices. Josef
Cheers,
Matthew _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org https://mail.scipy.org/mailman/listinfo/scipy-dev
On Fri, Mar 18, 2016 at 12:19 PM, <josef.pktd@gmail.com> wrote:
On Fri, Mar 18, 2016 at 3:07 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
On Wed, Mar 16, 2016 at 8:49 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 9:39 PM, <josef.pktd@gmail.com> wrote:
On Wed, Mar 16, 2016 at 11:01 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 7:16 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Wed, Mar 16, 2016 at 5:48 PM, Charles R Harris <charlesr.harris@gmail.com> wrote: > > > On Wed, Mar 16, 2016 at 6:41 PM, Charles R Harris > <charlesr.harris@gmail.com> wrote: >> >> >> >> On Wed, Mar 16, 2016 at 6:05 PM, Matthew Brett >> <matthew.brett@gmail.com> >> wrote: >>> >>> Hi, >>> >>> On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers >>> <ralf.gommers@gmail.com> >>> wrote: >>> > >>> > >>> > On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett >>> > <matthew.brett@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#scip.... >>> > 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 >> >> >> That difference looks pretty big, looks almost like a cutoff >> effect. >> What >> is the matrix X, in particular, what is its condition number or >> singular >> values? Might also be interesting to see residuals for the >> different >> functions. >> >> <snip> >> > > Note that this might be an indication that the result of the > calculation has > an intrinsically large error bar.
In [6]: np.linalg.svd(X, compute_uv=False) Out[6]: array([ 6.97301629e+07, 5.78702921e+05, 5.39646908e+05, 2.99383626e+05, 3.24329208e+04, 1.29520746e+02, 3.16197688e+01, 2.66011711e+00, 2.18287258e-01, 1.60105276e-01, 1.11807228e-01, 2.12694535e-03, 1.72675581e-03, 9.51576795e-04])
In [7]: np.linalg.cond(X) Out[7]: 73278544935.300507
Condition number is about 7e10, so you might expect to lose about 10 digits of of accuracy just due to roundoff. The covariance of the fitting error will be `inv(X.T, X) * error`, so you can look at the diagonals for the variance of the fitted parameters just for kicks, taking `error` to be roundoff error or measurement error. This is not to say that the results should not be more reproducible, but should be of interest for the validity of the solution. Note that I suspect that there may be a difference in the default `rcond` used, or maybe the matrix is on the edge and a small change in results changes the effective rank. You can run with the `return_rank=True` argument to check that. Looks like machine precision is used for the default value of rcond, so experimenting with that probable won't be informative, but you might want to try it anyway.
statsmodels is using numpy pinv for the linear model and it works in almost all cases very well *).
Using pinv with singular (given the rcond threshold) X makes everything well defined and not subject to numerical noise because of the regularization induced by pinv. Our covariance for the parameter estimates, based on pinv(X), is also not noisy but it's of reduced rank.
What I'm thinking is that changing the algorithm, as here, is like introducing noise relative to the original, so even if the result is perfectly valid, the fitted parameters might change significantly. A change in the 11-th digit wouldn't surprise me at all, depending on the actual X and rhs.
I believe the errors $e = y - X^+ y$ are in general robust to differences between valid pinv solutions, even in extreme cases:
In [2]: print_scalar(sse(X, npl.pinv(X), y)) 72.0232781366616024
In [3]: print_scalar(sse(np.c_[X, X], npl.pinv(np.c_[X, X]), y)) 72.0232781366939463
In [4]: npl.cond(np.c_[X, X]) Out[4]: 6.6247965335098062e+26
I think this does not depend on the smallest eigenvalues which should be regularized away in cases like this. The potential problem is with eigenvalues that are just a but too large for rcond to remove them but still have numerical problems.
I think it's a similar issue about the threshold as the one that you raised some time ago with determining the matrix rank for near singular matrices.
My point was (I think you're agreeing) that the accuracy of the pinv does not depend in a simple way on the condition number of the input matrix. I'm not quite sure how to assess the change in accuracy of scipy 0.17 on Linux, but it does seem to me to be significantly different (and worse than) the other options (like numpy, scipy 0.16) Matthew
On Fri, Mar 18, 2016 at 3:56 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
On Fri, Mar 18, 2016 at 12:19 PM, <josef.pktd@gmail.com> wrote:
On Fri, Mar 18, 2016 at 3:07 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
On Wed, Mar 16, 2016 at 8:49 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 9:39 PM, <josef.pktd@gmail.com> wrote:
On Wed, Mar 16, 2016 at 11:01 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 7:16 PM, Matthew Brett <matthew.brett@gmail.com> wrote: > > Hi, > > On Wed, Mar 16, 2016 at 5:48 PM, Charles R Harris > <charlesr.harris@gmail.com> wrote: > > > > > > On Wed, Mar 16, 2016 at 6:41 PM, Charles R Harris > > <charlesr.harris@gmail.com> wrote: > >> > >> > >> > >> On Wed, Mar 16, 2016 at 6:05 PM, Matthew Brett > >> <matthew.brett@gmail.com> > >> wrote: > >>> > >>> Hi, > >>> > >>> On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers > >>> <ralf.gommers@gmail.com> > >>> wrote: > >>> > > >>> > > >>> > On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett > >>> > <matthew.brett@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#scip... . > >>> > 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 > >> > >> > >> That difference looks pretty big, looks almost like a cutoff > >> effect. > >> What > >> is the matrix X, in particular, what is its condition number or > >> singular > >> values? Might also be interesting to see residuals for the > >> different > >> functions. > >> > >> <snip> > >> > > > > Note that this might be an indication that the result of the > > calculation has > > an intrinsically large error bar. > > In [6]: np.linalg.svd(X, compute_uv=False) > Out[6]: > array([ 6.97301629e+07, 5.78702921e+05, 5.39646908e+05, > 2.99383626e+05, 3.24329208e+04, 1.29520746e+02, > 3.16197688e+01, 2.66011711e+00, 2.18287258e-01, > 1.60105276e-01, 1.11807228e-01, 2.12694535e-03, > 1.72675581e-03, 9.51576795e-04]) > > In [7]: np.linalg.cond(X) > Out[7]: 73278544935.300507
Condition number is about 7e10, so you might expect to lose about 10 digits of of accuracy just due to roundoff. The covariance of the fitting error will be `inv(X.T, X) * error`, so you can look at the diagonals for the variance of the fitted parameters just for kicks, taking `error` to be roundoff error or measurement error. This is not to say that the results should not be more reproducible, but should be of interest for the validity of the solution. Note that I suspect that there may be a difference in the default `rcond` used, or maybe the matrix is on the edge and a small change in results changes the effective rank. You can run with the `return_rank=True` argument to check that. Looks like machine precision is used for the default value of rcond, so experimenting with that probable won't be informative, but you might want to try it anyway.
statsmodels is using numpy pinv for the linear model and it works in almost all cases very well *).
Using pinv with singular (given the rcond threshold) X makes everything well defined and not subject to numerical noise because of the regularization induced by pinv. Our covariance for the parameter estimates, based on pinv(X), is also not noisy but it's of reduced rank.
What I'm thinking is that changing the algorithm, as here, is like introducing noise relative to the original, so even if the result is perfectly valid, the fitted parameters might change significantly. A change in the 11-th digit wouldn't surprise me at all, depending on the actual X and rhs.
I believe the errors $e = y - X^+ y$ are in general robust to differences between valid pinv solutions, even in extreme cases:
In [2]: print_scalar(sse(X, npl.pinv(X), y)) 72.0232781366616024
In [3]: print_scalar(sse(np.c_[X, X], npl.pinv(np.c_[X, X]), y)) 72.0232781366939463
In [4]: npl.cond(np.c_[X, X]) Out[4]: 6.6247965335098062e+26
I think this does not depend on the smallest eigenvalues which should be regularized away in cases like this. The potential problem is with eigenvalues that are just a but too large for rcond to remove them but still have numerical problems.
I think it's a similar issue about the threshold as the one that you raised some time ago with determining the matrix rank for near singular matrices.
My point was (I think you're agreeing) that the accuracy of the pinv does not depend in a simple way on the condition number of the input matrix.
That's not the point I wanted to make. The condition number doesn't matter if there is an singular/eigenvalue with rcond 1e-16 or 1e-30 or 1e-100. However, the numerical accuracy might depend on the condition number when there is a singular value with rcond of 2e-15 or 5e-15. I didn't manage to find the default cond or rcond for scipy.linalg.pinv and lstsq. The problem I had was with numpy's pinv. I think there is a grey zone in rcond of eigenvalues that might be bad numerically. (My only experience is with SVD based pinv.) Josef
I'm not quite sure how to assess the change in accuracy of scipy 0.17 on Linux, but it does seem to me to be significantly different (and worse than) the other options (like numpy, scipy 0.16)
Matthew _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org https://mail.scipy.org/mailman/listinfo/scipy-dev
On Fri, Mar 18, 2016 at 4:09 PM, <josef.pktd@gmail.com> wrote:
On Fri, Mar 18, 2016 at 3:56 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
On Fri, Mar 18, 2016 at 12:19 PM, <josef.pktd@gmail.com> wrote:
On Fri, Mar 18, 2016 at 3:07 PM, Matthew Brett <matthew.brett@gmail.com
wrote:
On Wed, Mar 16, 2016 at 8:49 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 9:39 PM, <josef.pktd@gmail.com> wrote:
On Wed, Mar 16, 2016 at 11:01 PM, Charles R Harris <charlesr.harris@gmail.com> wrote: > > > > On Wed, Mar 16, 2016 at 7:16 PM, Matthew Brett > <matthew.brett@gmail.com> > wrote: >> >> Hi, >> >> On Wed, Mar 16, 2016 at 5:48 PM, Charles R Harris >> <charlesr.harris@gmail.com> wrote: >> > >> > >> > On Wed, Mar 16, 2016 at 6:41 PM, Charles R Harris >> > <charlesr.harris@gmail.com> wrote: >> >> >> >> >> >> >> >> On Wed, Mar 16, 2016 at 6:05 PM, Matthew Brett >> >> <matthew.brett@gmail.com> >> >> wrote: >> >>> >> >>> Hi, >> >>> >> >>> On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers >> >>> <ralf.gommers@gmail.com> >> >>> wrote: >> >>> > >> >>> > >> >>> > On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett >> >>> > <matthew.brett@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#scip... . >> >>> > 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 >> >> >> >> >> >> That difference looks pretty big, looks almost like a cutoff >> >> effect. >> >> What >> >> is the matrix X, in particular, what is its condition number or >> >> singular >> >> values? Might also be interesting to see residuals for the >> >> different >> >> functions. >> >> >> >> <snip> >> >> >> > >> > Note that this might be an indication that the result of the >> > calculation has >> > an intrinsically large error bar. >> >> In [6]: np.linalg.svd(X, compute_uv=False) >> Out[6]: >> array([ 6.97301629e+07, 5.78702921e+05, 5.39646908e+05, >> 2.99383626e+05, 3.24329208e+04, 1.29520746e+02, >> 3.16197688e+01, 2.66011711e+00, 2.18287258e-01, >> 1.60105276e-01, 1.11807228e-01, 2.12694535e-03, >> 1.72675581e-03, 9.51576795e-04]) >> >> In [7]: np.linalg.cond(X) >> Out[7]: 73278544935.300507 > > > Condition number is about 7e10, so you might expect to lose about 10 > digits of of accuracy just due to roundoff. The covariance of the > fitting > error will be `inv(X.T, X) * error`, so you can look at the diagonals > for > the variance of the fitted parameters just for kicks, taking `error` > to be > roundoff error or measurement error. This is not to say that the > results > should not be more reproducible, but should be of interest for the > validity > of the solution. Note that I suspect that there may be a difference > in the > default `rcond` used, or maybe the matrix is on the edge and a small > change > in results changes the effective rank. You can run with the > `return_rank=True` argument to check that. Looks like machine > precision is > used for the default value of rcond, so experimenting with that > probable > won't be informative, but you might want to try it anyway.
statsmodels is using numpy pinv for the linear model and it works in almost all cases very well *).
Using pinv with singular (given the rcond threshold) X makes everything well defined and not subject to numerical noise because of the regularization induced by pinv. Our covariance for the parameter estimates, based on pinv(X), is also not noisy but it's of reduced rank.
What I'm thinking is that changing the algorithm, as here, is like introducing noise relative to the original, so even if the result is perfectly valid, the fitted parameters might change significantly. A change in the 11-th digit wouldn't surprise me at all, depending on the actual X and rhs.
I believe the errors $e = y - X^+ y$ are in general robust to differences between valid pinv solutions, even in extreme cases:
In [2]: print_scalar(sse(X, npl.pinv(X), y)) 72.0232781366616024
In [3]: print_scalar(sse(np.c_[X, X], npl.pinv(np.c_[X, X]), y)) 72.0232781366939463
In [4]: npl.cond(np.c_[X, X]) Out[4]: 6.6247965335098062e+26
I think this does not depend on the smallest eigenvalues which should be regularized away in cases like this. The potential problem is with eigenvalues that are just a but too large for rcond to remove them but still have numerical problems.
I think it's a similar issue about the threshold as the one that you raised some time ago with determining the matrix rank for near singular matrices.
My point was (I think you're agreeing) that the accuracy of the pinv does not depend in a simple way on the condition number of the input matrix.
That's not the point I wanted to make. The condition number doesn't matter if there is an singular/eigenvalue with rcond 1e-16 or 1e-30 or 1e-100. However, the numerical accuracy might depend on the condition number when there is a singular value with rcond of 2e-15 or 5e-15. I didn't manage to find the default cond or rcond for scipy.linalg.pinv and lstsq. The problem I had was with numpy's pinv. I think there is a grey zone in rcond of eigenvalues that might be bad numerically.
(My only experience is with SVD based pinv.)
In the singular values array that you showed earlier, you have three small eigenvalues, but their rcond are only around 1e-11. So from my experience with SVD pinv, there shouldn't be any problem with a pinv in this case. Josef
Josef
I'm not quite sure how to assess the change in accuracy of scipy 0.17 on Linux, but it does seem to me to be significantly different (and worse than) the other options (like numpy, scipy 0.16)
Matthew _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org https://mail.scipy.org/mailman/listinfo/scipy-dev
On Fri, Mar 18, 2016 at 1:19 PM, <josef.pktd@gmail.com> wrote:
On Fri, Mar 18, 2016 at 4:09 PM, <josef.pktd@gmail.com> wrote:
On Fri, Mar 18, 2016 at 3:56 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
On Fri, Mar 18, 2016 at 12:19 PM, <josef.pktd@gmail.com> wrote:
On Fri, Mar 18, 2016 at 3:07 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
On Wed, Mar 16, 2016 at 8:49 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 9:39 PM, <josef.pktd@gmail.com> wrote: > > > > On Wed, Mar 16, 2016 at 11:01 PM, Charles R Harris > <charlesr.harris@gmail.com> wrote: >> >> >> >> On Wed, Mar 16, 2016 at 7:16 PM, Matthew Brett >> <matthew.brett@gmail.com> >> wrote: >>> >>> Hi, >>> >>> On Wed, Mar 16, 2016 at 5:48 PM, Charles R Harris >>> <charlesr.harris@gmail.com> wrote: >>> > >>> > >>> > On Wed, Mar 16, 2016 at 6:41 PM, Charles R Harris >>> > <charlesr.harris@gmail.com> wrote: >>> >> >>> >> >>> >> >>> >> On Wed, Mar 16, 2016 at 6:05 PM, Matthew Brett >>> >> <matthew.brett@gmail.com> >>> >> wrote: >>> >>> >>> >>> Hi, >>> >>> >>> >>> On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers >>> >>> <ralf.gommers@gmail.com> >>> >>> wrote: >>> >>> > >>> >>> > >>> >>> > On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett >>> >>> > <matthew.brett@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#scip.... >>> >>> > 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 >>> >> >>> >> >>> >> That difference looks pretty big, looks almost like a cutoff >>> >> effect. >>> >> What >>> >> is the matrix X, in particular, what is its condition number >>> >> or >>> >> singular >>> >> values? Might also be interesting to see residuals for the >>> >> different >>> >> functions. >>> >> >>> >> <snip> >>> >> >>> > >>> > Note that this might be an indication that the result of the >>> > calculation has >>> > an intrinsically large error bar. >>> >>> In [6]: np.linalg.svd(X, compute_uv=False) >>> Out[6]: >>> array([ 6.97301629e+07, 5.78702921e+05, 5.39646908e+05, >>> 2.99383626e+05, 3.24329208e+04, 1.29520746e+02, >>> 3.16197688e+01, 2.66011711e+00, 2.18287258e-01, >>> 1.60105276e-01, 1.11807228e-01, 2.12694535e-03, >>> 1.72675581e-03, 9.51576795e-04]) >>> >>> In [7]: np.linalg.cond(X) >>> Out[7]: 73278544935.300507 >> >> >> Condition number is about 7e10, so you might expect to lose about >> 10 >> digits of of accuracy just due to roundoff. The covariance of the >> fitting >> error will be `inv(X.T, X) * error`, so you can look at the >> diagonals >> for >> the variance of the fitted parameters just for kicks, taking >> `error` >> to be >> roundoff error or measurement error. This is not to say that the >> results >> should not be more reproducible, but should be of interest for the >> validity >> of the solution. Note that I suspect that there may be a >> difference >> in the >> default `rcond` used, or maybe the matrix is on the edge and a >> small >> change >> in results changes the effective rank. You can run with the >> `return_rank=True` argument to check that. Looks like machine >> precision is >> used for the default value of rcond, so experimenting with that >> probable >> won't be informative, but you might want to try it anyway. > > > > statsmodels is using numpy pinv for the linear model and it works > in > almost all cases very well *). > > Using pinv with singular (given the rcond threshold) X makes > everything > well defined and not subject to numerical noise because of the > regularization induced by pinv. Our covariance for the parameter > estimates, > based on pinv(X), is also not noisy but it's of reduced rank.
What I'm thinking is that changing the algorithm, as here, is like introducing noise relative to the original, so even if the result is perfectly valid, the fitted parameters might change significantly. A change in the 11-th digit wouldn't surprise me at all, depending on the actual X and rhs.
I believe the errors $e = y - X^+ y$ are in general robust to differences between valid pinv solutions, even in extreme cases:
In [2]: print_scalar(sse(X, npl.pinv(X), y)) 72.0232781366616024
In [3]: print_scalar(sse(np.c_[X, X], npl.pinv(np.c_[X, X]), y)) 72.0232781366939463
In [4]: npl.cond(np.c_[X, X]) Out[4]: 6.6247965335098062e+26
I think this does not depend on the smallest eigenvalues which should be regularized away in cases like this. The potential problem is with eigenvalues that are just a but too large for rcond to remove them but still have numerical problems.
I think it's a similar issue about the threshold as the one that you raised some time ago with determining the matrix rank for near singular matrices.
My point was (I think you're agreeing) that the accuracy of the pinv does not depend in a simple way on the condition number of the input matrix.
That's not the point I wanted to make. The condition number doesn't matter if there is an singular/eigenvalue with rcond 1e-16 or 1e-30 or 1e-100. However, the numerical accuracy might depend on the condition number when there is a singular value with rcond of 2e-15 or 5e-15. I didn't manage to find the default cond or rcond for scipy.linalg.pinv and lstsq. The problem I had was with numpy's pinv. I think there is a grey zone in rcond of eigenvalues that might be bad numerically.
(My only experience is with SVD based pinv.)
In the singular values array that you showed earlier, you have three small eigenvalues, but their rcond are only around 1e-11. So from my experience with SVD pinv, there shouldn't be any problem with a pinv in this case.
So - is there any downside to recommending numpy.linalg.pinv over scipy.linalg.pinv in general? Cheers, Matthew
On Sat, Mar 19, 2016 at 3:31 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
On Fri, Mar 18, 2016 at 1:19 PM, <josef.pktd@gmail.com> wrote:
On Fri, Mar 18, 2016 at 4:09 PM, <josef.pktd@gmail.com> wrote:
On Fri, Mar 18, 2016 at 3:56 PM, Matthew Brett <matthew.brett@gmail.com
wrote:
On Fri, Mar 18, 2016 at 12:19 PM, <josef.pktd@gmail.com> wrote:
On Fri, Mar 18, 2016 at 3:07 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
On Wed, Mar 16, 2016 at 8:49 PM, Charles R Harris <charlesr.harris@gmail.com> wrote: > > > On Wed, Mar 16, 2016 at 9:39 PM, <josef.pktd@gmail.com> wrote: >> >> >> >> On Wed, Mar 16, 2016 at 11:01 PM, Charles R Harris >> <charlesr.harris@gmail.com> wrote: >>> >>> >>> >>> On Wed, Mar 16, 2016 at 7:16 PM, Matthew Brett >>> <matthew.brett@gmail.com> >>> wrote: >>>> >>>> Hi, >>>> >>>> On Wed, Mar 16, 2016 at 5:48 PM, Charles R Harris >>>> <charlesr.harris@gmail.com> wrote: >>>> > >>>> > >>>> > On Wed, Mar 16, 2016 at 6:41 PM, Charles R Harris >>>> > <charlesr.harris@gmail.com> wrote: >>>> >> >>>> >> >>>> >> >>>> >> On Wed, Mar 16, 2016 at 6:05 PM, Matthew Brett >>>> >> <matthew.brett@gmail.com> >>>> >> wrote: >>>> >>> >>>> >>> Hi, >>>> >>> >>>> >>> On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers >>>> >>> <ralf.gommers@gmail.com> >>>> >>> wrote: >>>> >>> > >>>> >>> > >>>> >>> > On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett >>>> >>> > <matthew.brett@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
>>>> >>> >> `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#scip... . >>>> >>> > 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 >>>> >> >>>> >> >>>> >> That difference looks pretty big, looks almost like a cutoff >>>> >> effect. >>>> >> What >>>> >> is the matrix X, in particular, what is its condition number >>>> >> or >>>> >> singular >>>> >> values? Might also be interesting to see residuals for the >>>> >> different >>>> >> functions. >>>> >> >>>> >> <snip> >>>> >> >>>> > >>>> > Note that this might be an indication that the result of the >>>> > calculation has >>>> > an intrinsically large error bar. >>>> >>>> In [6]: np.linalg.svd(X, compute_uv=False) >>>> Out[6]: >>>> array([ 6.97301629e+07, 5.78702921e+05, 5.39646908e+05, >>>> 2.99383626e+05, 3.24329208e+04, 1.29520746e+02, >>>> 3.16197688e+01, 2.66011711e+00, 2.18287258e-01, >>>> 1.60105276e-01, 1.11807228e-01, 2.12694535e-03, >>>> 1.72675581e-03, 9.51576795e-04]) >>>> >>>> In [7]: np.linalg.cond(X) >>>> Out[7]: 73278544935.300507 >>> >>> >>> Condition number is about 7e10, so you might expect to lose about >>> 10 >>> digits of of accuracy just due to roundoff. The covariance of
>>> fitting >>> error will be `inv(X.T, X) * error`, so you can look at the >>> diagonals >>> for >>> the variance of the fitted parameters just for kicks, taking >>> `error` >>> to be >>> roundoff error or measurement error. This is not to say that the >>> results >>> should not be more reproducible, but should be of interest for
of the the
>>> validity >>> of the solution. Note that I suspect that there may be a >>> difference >>> in the >>> default `rcond` used, or maybe the matrix is on the edge and a >>> small >>> change >>> in results changes the effective rank. You can run with the >>> `return_rank=True` argument to check that. Looks like machine >>> precision is >>> used for the default value of rcond, so experimenting with that >>> probable >>> won't be informative, but you might want to try it anyway. >> >> >> >> statsmodels is using numpy pinv for the linear model and it works >> in >> almost all cases very well *). >> >> Using pinv with singular (given the rcond threshold) X makes >> everything >> well defined and not subject to numerical noise because of the >> regularization induced by pinv. Our covariance for the parameter >> estimates, >> based on pinv(X), is also not noisy but it's of reduced rank. > > > What I'm thinking is that changing the algorithm, as here, is like > introducing noise relative to the original, so even if the result is > perfectly valid, the fitted parameters might change significantly. A > change > in the 11-th digit wouldn't surprise me at all, depending on the > actual > X > and rhs.
I believe the errors $e = y - X^+ y$ are in general robust to differences between valid pinv solutions, even in extreme cases:
In [2]: print_scalar(sse(X, npl.pinv(X), y)) 72.0232781366616024
In [3]: print_scalar(sse(np.c_[X, X], npl.pinv(np.c_[X, X]), y)) 72.0232781366939463
In [4]: npl.cond(np.c_[X, X]) Out[4]: 6.6247965335098062e+26
I think this does not depend on the smallest eigenvalues which should be regularized away in cases like this. The potential problem is with eigenvalues that are just a but too large for rcond to remove them but still have numerical problems.
I think it's a similar issue about the threshold as the one that you raised some time ago with determining the matrix rank for near singular matrices.
My point was (I think you're agreeing) that the accuracy of the pinv does not depend in a simple way on the condition number of the input matrix.
That's not the point I wanted to make. The condition number doesn't matter if there is an singular/eigenvalue with rcond 1e-16 or 1e-30 or 1e-100. However, the numerical accuracy might depend on the condition number when there is a singular value with rcond of 2e-15 or 5e-15. I didn't manage to find the default cond or rcond for scipy.linalg.pinv and lstsq. The problem I had was with numpy's pinv. I think there is a grey zone in rcond of eigenvalues that might be bad numerically.
(My only experience is with SVD based pinv.)
In the singular values array that you showed earlier, you have three small eigenvalues, but their rcond are only around 1e-11. So from my experience with SVD pinv, there shouldn't be any problem with a pinv in this case.
So - is there any downside to recommending numpy.linalg.pinv over scipy.linalg.pinv in general?
The only time that I did a comparison was several years ago with the ill conditioned NIST test case. Some of the scipy pinv where a tiny bit more accurate than the numpy pinv. Smallest singular or eigenvalues were zero +/- noise. My impression was that in edge cases where numerical noise matters, some implementation details or which linalg library is used can have a visible effect. (I still have never read a numerical linear algebra book, and my "incompetence" from a thread many years ago is only tempered by 7 more years of usage.) I'm pretty happy with using np.linalg.pinv, although there are also a few cases of noise results. The main reason I stayed with numpy's pinv for statsmodels was that it was quite a bit faster than scipy pinv, partially because it doesn't do the isfinite check and maybe for some other algorithmic reasons. I haven't compared the two since then, and scipy's pinv has changed since then. statsmodels unit test for these kind of good problems are usually in the rtol=1e-12 to 1e-13 range. Going up to 1e-14 usually failed on some machine during release testing. Problems that are not so nice have lower precision. Josef
Cheers,
Matthew _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org https://mail.scipy.org/mailman/listinfo/scipy-dev
numpy.linalg.pinv uses the same solver as scipy.linalg.pinv2 (SVD) while scipy.linalg.pinv uses the least squares solver. -- Olivier
Hi, On Wed, Mar 16, 2016 at 8:39 PM, <josef.pktd@gmail.com> wrote:
On Wed, Mar 16, 2016 at 11:01 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 7:16 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Wed, Mar 16, 2016 at 5:48 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 6:41 PM, Charles R Harris <charlesr.harris@gmail.com> wrote:
On Wed, Mar 16, 2016 at 6:05 PM, Matthew Brett <matthew.brett@gmail.com> wrote:
Hi,
On Sun, Mar 13, 2016 at 3:43 AM, Ralf Gommers <ralf.gommers@gmail.com> wrote: > > > On Sat, Mar 12, 2016 at 6:16 AM, Matthew Brett > <matthew.brett@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#scip.... > 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
That difference looks pretty big, looks almost like a cutoff effect. What is the matrix X, in particular, what is its condition number or singular values? Might also be interesting to see residuals for the different functions.
<snip>
Note that this might be an indication that the result of the calculation has an intrinsically large error bar.
In [6]: np.linalg.svd(X, compute_uv=False) Out[6]: array([ 6.97301629e+07, 5.78702921e+05, 5.39646908e+05, 2.99383626e+05, 3.24329208e+04, 1.29520746e+02, 3.16197688e+01, 2.66011711e+00, 2.18287258e-01, 1.60105276e-01, 1.11807228e-01, 2.12694535e-03, 1.72675581e-03, 9.51576795e-04])
In [7]: np.linalg.cond(X) Out[7]: 73278544935.300507
Condition number is about 7e10, so you might expect to lose about 10 digits of of accuracy just due to roundoff. The covariance of the fitting error will be `inv(X.T, X) * error`, so you can look at the diagonals for the variance of the fitted parameters just for kicks, taking `error` to be roundoff error or measurement error. This is not to say that the results should not be more reproducible, but should be of interest for the validity of the solution. Note that I suspect that there may be a difference in the default `rcond` used, or maybe the matrix is on the edge and a small change in results changes the effective rank. You can run with the `return_rank=True` argument to check that. Looks like machine precision is used for the default value of rcond, so experimenting with that probable won't be informative, but you might want to try it anyway.
statsmodels is using numpy pinv for the linear model and it works in almost all cases very well *).
Yes, the error for numpy.linalg pinv is a lot less: Sympy high-precision pinv (as before): 72.0232781366602524 Scipy 0.16 pinv (as before, difference at 12th dp): 72.0232781366620856 Scipy 0.17 pinv (as before, difference at 9th dp): 72.0232781390479602 Numpy pinv (difference at 12th dp) In [2]: print_scalar(sse(X, np.linalg.pinv(X), y)) 72.0232781366618440 So numpy pinv is more accurate than scipy pinv for this matrix, for scipy 0.17, and both scipy 0.16 and numpy differ at the 12th dp, as does estimation for a singular matrix (previous email). Incidentally, this is all on scipy / openblas / linux. On OSX: Scipy 0.17.0: 72.0232781366634640 Numpy: 72.0232781366609771 Scipy 0.16.0 72.0232781366634640 So - on OSX - differences also predictably at 12th dp for both scipy versions and numpy. It does seem that scipy 0.17 on Linux is the outlier. Matthew
On Thu, Mar 17, 2016 at 2:30 AM, Olivier Grisel <olivier.grisel@ensta.org> wrote:
Can you check with Atlas on linux then?
Yes, same idea, this on another computer, Debian packaged ATLAS: Sympy high-precision pinv SSE 72.0232781366615313 Local numpy 1.10.4 linalg pinv 72.0232781366614176 Local scipy 0.16.0 linalg pinv 72.0232781366600392 Local scipy 0.17.0 linalg pinv 72.0232781368149233 Difference at 10th dp, but still two orders of magnitude different from the other variations. https://github.com/matthew-brett/scipy-pinv/blob/master/pinv_futz.py Matthew
participants (5)
-
Charles R Harris -
josef.pktd＠gmail.com -
Matthew Brett -
Olivier Grisel -
Ralf Gommers