ENH: Use block structure in pade #14933

DerWeh · 2021-10-29T14:05:16Z

Reference issue

None

What does this implement/fix?

The Pade is solved as the matrix quation:

[A, B] [p, q]^T = 0

The Matrix A, however, has the very simple block structure

A = [I, 0]^T

so the matrix problem can be written in the block structure

|I, B_1| |p| = |0|
|0, B_2| |q| = |0|

The new version of pade solves the second line of this matrix problem

B_2 q = 0

and substitutes the solution.
This is more stable for large degree of the denominator polynomial q.

I personally need this, to use pade for the Fourier-Pade transformation with many Fourier coefficients.

Additional information

Additionally, I did the following refactoring:

use Polynomial instead of old poly1d

This makes it much clearer what happens than using loops.

`Polynomial` is preferred according to Numpy. Furthermore, the order is more intuitive, as the Taylor coefficients are also given in accenting order.

The Padé problem is formulated and solved as [B, A] [q, p]^T = 0 The matrix A, however has a very simple block structure: A = [-I, 0]^T where `I` is the identity matrix. We partition `B` likewise B = [B_1, B_2]^T and get the Block equation [[B_1, -I], [[q], [[0], [B_2, 0]] [p]] = [0]] We can now solve the second line, which is a smaller matrix equation, for `q` and substitute it in the first line for to obtain `p`.

tupui

Thanks @DerWeh! I cannot really comment the maths and will let @ev-br and others have a look 😃 (still tests seem good), still I have a suggestion for the tests and a general remark.

If I read this PR correctly, it will introduce a backward incompatibility. I am not sure how critical this function is, but we usually wait 2 releases before such changes. So for now, we could add a deprecation warning (we have a decorator for that @_deprecated) to warn future changes.

tupui · 2021-10-29T16:27:40Z

scipy/interpolate/tests/test_pade.py

 def test_pade_trivial():
    nump, denomp = pade([1.0], 0)
-    assert_array_equal(nump.c, [1.0])
-    assert_array_equal(denomp.c, [1.0])
+    assert_array_equal(nump.coef, [1.0])


Since you touch the whole file, can you change all assert_array_equal,assert_approx_equal to assert_allclose since it's also recommended to do the change?

Sure, I will change it.

But are you sure about assert_array_equal? This is a more strict test, equal, not just close. Also, I cannot find any mention that it is not recommended anymore.

Oups yes you are right I thought about another one.

For reference the list of non recommended functions is https://numpy.org/doc/stable/reference/routines.testing.html#asserts-not-recommended

DerWeh · 2021-10-29T17:29:41Z

You are right @tupui, somehow I missed that returning Polynomial breaks current behavior in spite of breaking current tests. I only thought about the evaluated rational polynomial...

Should I just revert the change, or would it be smarter to add an optional argument allowing to choose between Polynomial and poly1d, which defaults to poly1d and issues a deprecation warning?

tupui · 2021-10-29T18:21:32Z

No worries 😃

I would not add a parameter. Here we can just have a deprecation cycle. Let's wait for other maintainer's opinion before you do any extra work.

DerWeh · 2021-12-06T20:40:43Z

As we are still waiting for someone to comment on the math, I shipped my own version in the meantime as I need it.
Of course, I am willing to push to scipy, and drop maintaining a personal version.

It is documented here: https://gftools.readthedocs.io/en/0.11.0/generated/gftool.hermpade.pade.html

Further changes:

Use QR decomposition to get Null-Vector instead of fixing one parameter and solving the equation.
This is more stable (so it should typically not matter), and, more importantly, also handles the edge case when q[0]=0.
Add option fast to use less stable but much faster solve_toeplitz. I rarely require this to scan a huge amount of Pade approximants.

I also noted, that the example in the documentation is quite bad. It shows a single point, where it's hard to see the difference between Pade and Taylor expansion (you have to know 4 digits of np.e by heart to know what's better). And the Pade approximant isn't suitable to approximate the exponential function anyway (just give it a try).

ev-br

Looks nice! I'm not a power user of this, so am ready to trust your judgement. The list of planned improvements looks nice, these would be very welcome.

The main sticky point of this PR seems to be np.poly1d vs np.Polynomial. Can't just replace one with the other, so maybe it's time to cook up a replacement function and ditch the current pade.

ev-br · 2022-04-30T16:07:41Z

scipy/interpolate/_pade.py

@@ -29,8 +32,7 @@ def pade(an, m, n=None):
    >>> e_exp = [1.0, 1.0, 1.0/2.0, 1.0/6.0, 1.0/24.0, 1.0/120.0]
    >>> p, q = pade(e_exp, 2)

-    >>> e_exp.reverse()
-    >>> e_poly = np.poly1d(e_exp)
+    >>> e_poly = np.polynomial.Polynomial(e_exp)


That's a backcompat break, as Pamphille already said. Two possibilities: either revert back to poly1d (sigh) or add a separate function to return a Polynomial with your other improvements, and deprecate the older one.
I'd recommend the latter maybe?

I would also favor doing the latter. This further allows us to break with the convention q.coef[0]==1 and use the QR decomposition to get the null-vector instead of solving the system. This is conceptionally cleaner and more stable, albeit slower.
In rare cases, q.coef[0]==1 fails. (Had to exclude these in the property tests of my Pade implementation.)

ev-br · 2022-04-30T16:08:54Z

scipy/interpolate/_pade.py

+        return _Polynomial(an), _Polynomial([1])
+    # first solve the Toeplitz system for q
+    # first row might contain tailing zeros
+    top = r_[an[n+1::-1][:m+1], [0]*(m-n-1)]


Minor, trivial, optional: while we're at, how about spelling np.r_ (and removing the old-fashioned from numpy import r_)

ev-br · 2022-04-30T16:11:45Z

scipy/interpolate/_pade.py

+    if m > 100:  # arbitrary threshold when to use dedicated Toeplitz function
+        p = linalg.matmul_toeplitz((an[:n+1], zeros(m+1)), q)
+    else:
+        p = linalg.toeplitz(an[:n+1], zeros(m+1)) @ q


About the threshold: any guidance of how to choose the threshold?

No guidance, but I can simply test on my machine when matmul_toeplitz becomes faster, assuming memory doesn't matter.

Maybe include the user knob to control the switch. Or better, expose a method= kwarg and discuss them in the docstring notes. Or maybe it's all YAGNI.

ev-br · 2022-04-30T16:12:13Z

scipy/interpolate/_pade.py

+    top = r_[an[n+1::-1][:m+1], [0]*(m-n-1)]
+    an_mat = linalg.toeplitz(an[n+1:], top)
+    # we set q[0] = 1 -> first column is -rhs
+    q = r_[1.0, linalg.solve(an_mat[:, 1:], -an_mat[:, 0])]


Can use solve_toeplitz instead?

EDIT: Oh, I see you mentioned this. So it's the matter of numerical stability. Maybe then add a comment, and hold on on the keyword switch until need arises.

As mentioned, if we rewrite Pade I would replace solve by a QR decomposition.

ev-br · 2022-04-30T16:23:49Z

I also noted, that the example in the documentation is quite bad. It shows a single point, where it's hard to see the difference between Pade and Taylor expansion

A better example would definitely be most welcome!

DerWeh · 2022-05-01T18:05:16Z

Nice timing @ev-br, I just got time to work on the topic again and released an own version of Pade in my package, which I need for the Pade-Fourier algorithm. We can take https://gftools.readthedocs.io/en/0.11.0/generated/gftool.hermpade.pade.html as a reference. I used it quite a lot, but only in the Pade-Fourier context, which considers quite large polynomial degrees. With other use cases, I got little experience.

If we replace interpolate.pade, the question is also if it should be in interpolate. In contrast to all other functions in interpolate (as far as I see), Pade does not take function values as input but Fourier coefficients. There is the N-point Pade algorithm, which would interpolate function values, but we do not implement N-point Pade here! This confused me initially a lot. But I don't really see a good place for Pade either.

As soon as we decide how to proceed, I can readily implement it.

ev-br · 2022-05-01T20:20:56Z

Great that timings overlap :-).

In my own experience, I only used Pade approximants for weak-coupling expansions, where you have a handful of terms for a (hopefully convergent) diagrammatic series and try to resum it. So, not as hard-core as your uses, by quite a margin.
But certainly Fourier series is not the only application?

W.r.t if this should live in interpolate or elsewhere. When originally we moved it out of scipy.misc a few years back, there was a brief discussion if it belongs to interpolate or special. Nobody had a horse at that race, so it ended up where it did. Yes, it's a bit of an odd one out (even if scipy.interpolate is not exclusively function values as an input, see e.g. PPoly), and I suspect it will keep being that anywhere in scipy. So unless you have a strong preference, I'd suggest sticking with the status quo. If you do have one, sure, no problem from my side.

DerWeh · 2022-05-01T22:14:08Z

But certainly Fourier series is not the only application?

Sure, it was just the only case where I needed it so far (and had problems with scipy's current implementation).

I'd suggest sticking with the status quo

I don't have any strong opinion where to put Pade, so we can leave it where it is.

So what should I do now? Add a deprecation warning to scipy.interpolate.pade and create a new function? How should this new function be named?
In this case, it probably is a good idea to close this pull request and I redo it from scratch, as we revert pade to the original version.

ev-br · 2022-05-01T22:25:19Z

Add a deprecation warning to scipy.interpolate.pade and create a new function?

yup, can @np.deprecate (grep the scipy codebase for examples)

How should this new function be named?

Naming is hard :-). PadeApproximant?

ev-br · 2022-05-08T16:28:19Z

it probably is a good idea to close this pull request and I redo it from scratch, as we revert pade to the original version.

Ok, closing then, and looking forward for the replacement PR @DerWeh !

DerWeh added 4 commits October 29, 2021 16:00

STY: fix PEP8 for Padé

ff9d09c

MAINT: use Toeplitz matrix for Padé

3fc0297

This makes it much clearer what happens than using loops.

MAINT: return Polynomial instead of poly1d in Padé

9321a9c

`Polynomial` is preferred according to Numpy. Furthermore, the order is more intuitive, as the Taylor coefficients are also given in accenting order.

DerWeh requested a review from ev-br as a code owner October 29, 2021 14:05

github-actions bot added the scipy.interpolate label Oct 29, 2021

STY: fix lint errors

1228515

tupui reviewed Oct 29, 2021

View reviewed changes

tupui added the enhancement A new feature or improvement label Oct 29, 2021

ev-br requested changes Apr 30, 2022

View reviewed changes

ev-br added the needs-work Items that are pending response from the author label Apr 30, 2022

ev-br added this to In review in scipy.interpolate Apr 30, 2022

ev-br moved this from In review to triaged, waiting-for-contributor in scipy.interpolate May 8, 2022

ev-br closed this May 8, 2022

ev-br moved this from triaged, waiting-for-contributor to Abandoned in scipy.interpolate May 22, 2022

ev-br moved this from Abandoned to triaged, waiting-for-contributor in scipy.interpolate May 22, 2022

ev-br mentioned this pull request Feb 10, 2024

ENH: interpolate: improve calculation of Pade approximants with SVD based method #20064

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

ENH: Use block structure in pade #14933

ENH: Use block structure in pade #14933

DerWeh commented Oct 29, 2021

tupui left a comment

tupui Oct 29, 2021

DerWeh Oct 29, 2021

tupui Oct 29, 2021 •

edited

Loading

tupui Oct 29, 2021

DerWeh commented Oct 29, 2021

tupui commented Oct 29, 2021 •

edited

Loading

DerWeh commented Dec 6, 2021 •

edited

Loading

ev-br left a comment •

edited

Loading

ev-br Apr 30, 2022

DerWeh May 1, 2022

ev-br Apr 30, 2022

ev-br Apr 30, 2022

DerWeh May 1, 2022

ev-br May 1, 2022 •

edited

Loading

ev-br Apr 30, 2022 •

edited

Loading

DerWeh May 1, 2022

ev-br commented Apr 30, 2022

DerWeh commented May 1, 2022

ev-br commented May 1, 2022

DerWeh commented May 1, 2022

ev-br commented May 1, 2022

ev-br commented May 8, 2022

ENH: Use block structure in pade #14933

ENH: Use block structure in pade #14933

Conversation

DerWeh commented Oct 29, 2021

Reference issue

What does this implement/fix?

Additional information

tupui left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tupui Oct 29, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

DerWeh commented Oct 29, 2021

tupui commented Oct 29, 2021 • edited Loading

DerWeh commented Dec 6, 2021 • edited Loading

ev-br left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ev-br May 1, 2022 • edited Loading

Choose a reason for hiding this comment

ev-br Apr 30, 2022 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ev-br commented Apr 30, 2022

DerWeh commented May 1, 2022

ev-br commented May 1, 2022

DerWeh commented May 1, 2022

ev-br commented May 1, 2022

ev-br commented May 8, 2022

tupui Oct 29, 2021 •

edited

Loading

tupui commented Oct 29, 2021 •

edited

Loading

DerWeh commented Dec 6, 2021 •

edited

Loading

ev-br left a comment •

edited

Loading

ev-br May 1, 2022 •

edited

Loading

ev-br Apr 30, 2022 •

edited

Loading