MAINT: Fix the cutoff value inconsistency for pinv2 and pinvh #10067
Conversation
|
Can you explain in the description in what cases this matters for backwards compatibility? |
|
@ilayn If you are comfortable assessing Probably a little too optimistic to hope for any connection though. |
|
@rgommers Done. @tylerjereddy I'll have a look. |
|
@tylerjereddy It looks like the same issue fixed here. For single precision (np.float32, np.complex64 etc.) pinv value is still If OpenBLAS folks touched these and their related LAPACK routines or number representations, they can flip to the left or to the right if they are sitting on the fence. Hence I am not so sure if it is compiler related. Also, unfortunately, |
|
Thanks @ilayn, clear description. Let's leave this open for a while to give people a chance to weigh in. |
|
If there are no more comments on this, I'd like to merge it tomorrow. |
There was a problem hiding this comment.
Looks good to me. Probably worth a note in the backward incompatible changes release notes section, too, though.
Also some tests that show things work better now than on master would help if it's not too much work. Maybe something simple like pinv, pinvh, and pinv2 (good case for parametrize) on a 2x2 double array with 1 and 1e-10 in the diagonal?
|
Ah yes. Good idea. |
|
@larsoner pinvh is basically pinv2 with different solver so I only compared pinv and pinv2. I couldn't replicate a symmetric matrix that exhibits the issue but let me know if you have a better example to test with @tylerjereddy If there is still enough time please consider this to add to 1.3.0 |
|
Sigh... pypy is not having a good day again. @mattip Any suggestions? |
|
This is using NumPy 1.16.3 I assume, which was released yesterday and added a new C-API function. I will issue a PR to go back to using nightly builds rather than official releases, since PyPy fixed when it was merged to NumPy but it is not in a release yet.. |
|
Edit: travis runs seem to use a number of different NumPy versions, including latest. |
|
@tylerjereddy I forgot to add the milestone to this if @larsoner is happy with the changes it is ok to merge. |
There was a problem hiding this comment.
pinvh is basically pinv2 with different solver so I only compared pinv and pinv2. I couldn't replicate a symmetric matrix that exhibits the issue but let me know if you have a better example to test with
What I meant above about "a 2x2 double array with 1 and 1e-10 in the diagonal" is something like the following, which should fail on master for pinvh and pinv2 (and succeed for pinv) but all pass on your PR (at least when I tested it as a script, now generalized to the pytest case, it did I think):
@pytest.mark.parametrize('scale', (1e-20, 1., 1e20))
@pytest.mark.parametrize('pinv_', (pinv, pinvh, pinv2))
def test_auto_rcond(scale, pinv_):
x = np.array([[1, 0], [0, 1e-10]]) * scale
expected = np.diag(1. / np.diag(x))
x_inv = pinv_(x)
assert_allclose(x_inv, expected)
|
I think we should bump the milestone--looks like there's still some active discussion / revision going on here. |
|
Give me an hour and it will be done :) |
|
@larsoner After having sufficient coffee I can give a better response: Now if the cond arg is given it is respected and we only interfere if both are None which is what this PR is fixing. Also thank you for tests indeed it fails on master and passes this PR. Hence please have a look. But I think the terminology is wrong in the inversion context; just as an example https://nl.mathworks.com/help/matlab/ref/rcond.html |
|
OK I thought |
|
Thanks for keeping it sane. I really confused myself with this one 😃 |
|
I restored the milestone to 1.3.0 & did the same for the cognate issue. If this deserves a release note please feel free to add one on the wiki page. |
|
Thanks Tyler, wiki page edited. |
|
Dear member of scipy, Using Thank you for supporting scipy, Gilles |
|
|
I agree with @GillesOrban. @ilayn by reading the lines: https://github.com/scipy/scipy/pull/10067/files#diff-8400a871af5e3ff38643aff36ff192650313e769cdfe80f57b02140def686ee1R1368-R1377 one can see that:
Therefore those two arguments mean the same thing in practice. Another way to frame it: there is no way to pass Finally the docstring does not explain what is the (expected) difference between those two arguments. |
|
Arrrgh, this issue still haunts me 🤦 @ogrisel I think these keywords were incorrect to start with, as you mentioned. I'll come back to that in a sec. However I think there is a discrepancy in that rcond is not a relative condition number as in Back to keywords, ideally we shouldn't have had these wrongly named keywords, shortcut together behind the scenes but I didn't touch it for backwards compatibility. However I do agree that the naming is confusing. And if you want to have the relative cutoff value I would really appreciate a new issue/feature request as here not many people will see it. And more people can chime in. Probably we should have gone for Hopefully this time we can fix it for good and be done with this. |
Fixes #8861
The least square solver bound value is notoriously difficult to select due to the inherent trade-off between selecting the best minimizer for the ||Ax-b|| and the precision of x when b is in the image of the matrix A. Since
pinvis a special least squares problem the usual bound of dtype eps is relaxed to have a factor ofMAX(M, N). Otherwise the problem shown in the linked issue is hit more often than acceptable times.On the other hand SVD based pseudo-inverses enjoy the convenience of having the singular values available for obtaining an upper bound on the error. In this PR we remove the hard-coded values of 1E3 and 1E6 factors for pinv2 and pinvh and tie them to the
max(M,N)*sigma_max*eps. Probably for historical reasons the cutoff values were hardcoded.Again, these values by no means provide a failsafe situation. However, they are kind of the standard values that serve well in most of the cases. Otherwise one can always provide a better
condvalue manually.Here a possible backwards incompatibility can arise if a user relies on the resulting precision too heavily. But in my opinion, it repairs quite more than it might break.