Fix coefficients in expm helper function #9705
Conversation
According to [Higham 2005], the leading term of the backward error
function of the [m/m] Padé approximant of the exponential function is
x^{2m+1}. Its corresponding coefficient is C_{2m+1} = (m!)^2 / ( (2m)! *
(2m+1)!), cf. Equations (2.2) and (2.6) in said reference.
So far, the code didn't compute C_{2m+1}, but replaced m <- 2m + 1,
which is probably due to the literature's somewhat surprising notation
(cf. the index i = 2m + 1 in Equation (2.6)).
[Al-Mohy and Higham 2009] actually have a [Matlab reference implementation]
available, where they hard-code the cofficients, which are reproduced
verbatim below. Using m instead of m <- 2m + 1, as is currently done, exactly
reproduces these coefficients
```matlab
% Coefficients of leading terms in the backward error functions h_{2m+1}.
Coeff = [1/100800, 1/10059033600, 1/4487938430976000,...
1/5914384781877411840000, 1/113250775606021113483283660800000000];
```
[Higham 2005]: https://doi.org/10.1137/04061101X
[Al-Mohy and Higham 2009]: https://doi.org/10.1137/09074721X
[Matlab reference implementation]: http://eprints.ma.man.ac.uk/1442/03/expm_new.zip
|
Any way to move this forward? I am cannot assess how bad this incorrect factor is, but given that a probably non-trivial number of people is using matrix exponents in scipy, I would love to get this fixed (or assess the severity of this issue, if somebody has that ability). |
Re the unit tests:As I wrote, I cannot assess the severity of this error. I can only state that a) scipy's Right now, the unit tests available for scipy's expm itself are only testing trivial cases. Furthermore, I do not see how unit tests are useful here. If you want, I can write unit tests with the mathematically correct values, but that seems to be circular. Re licensing:I do not quite understand your point. In how far is their code used to fix scipy's code? Scipy implements the actual equations in the cited papers, with my fix calculating the same numbers that the matlab script hardcodes. There is no sharing of code, no reimplementation, not even a derivation of work. The calculation of the correct terms in the error function of the Padé approximants I did a) by hand because I was using scipy's implementation to learn the algorithm and not being able to get the same factors and b) by implementing them in a script. I have crosschecked the correctness of my factors by comparing them with the matlab code above. The matlab code itself was published alongside their paper, see the preprint server of Manchester University: http://eprints.maths.manchester.ac.uk/1442/ Steps forward:The code is mathematically incorrect. You can follow the derivation of the coefficients in the papers cited. Gautschi 2012 is also useful to read here. Scipy should either not mention that they implement Al-Mohy and Higham 2009, or alternatively specify that the implementation is probably not correct and warn people from using it. |
agreed
this is all fine, thanks for clarifying
yes, this PR should be okay as is. While this is a bugfix and we can therefore change the code, we should put a note in the 1.3.0 release notes when this PR is merged. Other than that, this just needs review. Unfortunately the I'll mark it for the next release; if no one has reviewed in 1-2 weeks I'll get to it at some point (no time now, sorry). |
I don't get how this wouldn't be useful. Maybe it is hard to do, and maybe the PR is fine as is, but certainly not ideal to not have a regression guard for something that is mathematically incorrect. |
Oh, better unit tests would definitely be useful! I just don't know how to approach this, because I would need a much deeper understanding of the methods involved. Maybe you can open an issue separate from this PR and take the lead on this? It would definitely be nice to see a unit test on |
|
Makes sense -- I'll think about this a bit! |
|
@SuperFluffy Thanks! Can you add a release note to: https://github.com/scipy/scipy/pulls?q=is%3Aopen+is%3Apr+milestone%3A1.3.0 |
|
@tylerjereddy Awesome that this got merged, makes me very happy! It looks like I cannot assign labels/milestons to issues. I guess I don't have the rights. :) Can you do that for me? Here is the associated issue: #10074 I hope that that's what you meant me to do. You linked the issue tracker filtered by milestone |
done
No the link was wrong. Tyler's question was to add a note to https://github.com/scipy/scipy/wiki/Release-note-entries-for-SciPy-1.3.0#scipysparse-improvements. I'm not sure that this is relevant for users though, since it's not clear how this is relevant to an end user of |
It would indeed be interesting to ask users who rely heavily on |
|
That's a good idea - more a question for the mailing list than the release notes though. If you would want to write a message to scipy-user, that would be helpful |
According to Higham 2005, the leading term of the backward error function of the [m/m] Padé approximant of the exponential function is
x^{2m+1}. Its corresponding coefficient isC_{2m+1} = (m!)^2 / ( (2m)! * (2m+1)!), cf. Equations (2.2) and (2.6) in said reference.So far, the code didn't compute
C_{2m+1}, but replacedm <- 2m + 1, which is probably due to the literature's somewhat surprising notation (cf. the indexi = 2m + 1in Equation (2.6)).Al-Mohy and Higham 2009 actually have a Matlab reference implementation available, where they hard-code the cofficients, which are reproduced verbatim below. Using m instead of
m <- 2m + 1, as is currently done, exactly reproduces these coefficients: