Skip to main content

Automatic Asymptotics with Sage

20 Jan 2022 - Tags: sage

Recurrence relations show up all the time in combinatorics and computer science, and even simple objects give recurrences that are difficult or impossible to solve. Thankfully, we can often find a generating function for our objects, and through the power of complex analysis and algebraic geometry, we can use the singularities of the generating function in order to get good asymptotic estimates. Software like Maple can automatically compute asymptotic expansions for a lot of generating functions using the asympt function… Can sage?

The answer is “yes”, which is why I’m writing this post, but the longer answer is “with some work”. Which is the real reason I’m writing this post. I had to slog through a lot of kind of crummy documentation to get this working, and I want to make sure that the next people looking into this have an easier time.

There are two modules (that I can find) which provide features for computing asymptotics:

These both have pros and cons:

The former is slightly easier to work with, and allows generating functions with multiple variables. The downside is that it only works for functions with polynomial denominators1.

The latter has nicer documentation, works for more general functions, and gives more detailed asymptotic expansions. The downside is that it only works for functions of a single variable, and it’s kind of annoying to use because it doesn’t accept symbolic inputs natively. You also have to manually provide a list of singularities2.

Let’s start with the multivariate case, since it’s more complicated. I’ll give a code dump, then explain what’s going on.

The tl;dr is that it works for functions of the form $h / p$ where $h$ is a holomorphic function and $p$ is a polynomial. The number of variables is irrelevant.

This is pretty obviously cribbed from the examples here, but the basic idea is this:

We need the denominator of our generating function to be a factored polynomial, so we build a polynomial ring with the variables in the denominator and factor over that ring.

Then, we put this into FFPD to get an object that the module knows how to deal with.

We choose a multi-index $\alpha$ which controls the “direction” in which we take our asymptotics3. For instance, if our generating function is $f(x,y) = \sum_{m,n} F_{m,n} x^m y^n$, then

If we view the $F_{m,n}$ as being located at the lattice points $(m,n)$ in the plane, then $\alpha$ picks out the direction in which we move.

Next we look at the critical points of the denominator. That is, the points where its gradient is either undefined or vanishing. The former is the singular case, and the latter is the smooth case.

In the one variable case, we know the asymptotics are controlled by the singularity closest to the origin. See chapter $5$ of Generatingfunctionology, for instance.

From skimming these articles, we now have a whole algebraic variety of singularities, and for reasons the asymptotics are now governed by the critical points on this variety. In particular, we’re on the hunt for the critical point which is closest to the origin (which we’ll call the Dominant Singularity).

So now, what does the code do? Well we look for the location of the minimal smooth point and the minimal singular point. Then take the smaller one and run the asymptotic function provided by the module.

It seems to work quite well, and is surprisingly fast too. For a simple example, we can try

which outputs4

\[\frac{1}{5} \, \sqrt{5} \left(\frac{2}{\sqrt{5} - 1}\right)^{r}\]

Of course, \(\frac{z}{1-z-z^2}\) is the generating function for the fibonacci numbers, and we’ve successfully recovered the asymptotic formula (though it requires some massaging to make it look like its usual presentation).

This really shines when it comes to multivariate series, though. For instance, we know that

\[\binom{n}{k} = [x^n y^k] (1 - x(1+y))^{-1}\]

So if we’re interested in the asymptotics of $\binom{3n}{n}$ we can compute

which quickly gives

\[\frac{1}{41472} \, \left(\frac{27}{4}\right)^{r} {\left(\frac{10368 \, \sqrt{6} \sqrt{2}}{\sqrt{\pi} \sqrt{r}} - \frac{1008 \, \sqrt{6} \sqrt{2}}{\sqrt{\pi} r^{\frac{3}{2}}} + \frac{49 \, \sqrt{6} \sqrt{2}}{\sqrt{\pi} r^{\frac{5}{2}}}\right)}\]

Now $\binom{15}{5} = 3003$ and evaluating the above at $r=5$ gives $3002.931$.

But what if we’re interested in generating functions whose numerator isn’t holomorphic, or whose denominator isn’t a polynomial? It’s clear that we should be interested in such things, since the famous catalan numbers already have the generating function

\[\frac{1 - \sqrt{1 - 4z}}{2z}\]

For cases like this (which must be of a single variable), we work with the second asymptotics package here.

This is nice because it returns an asymptotic expansion, which is quite intuitive to work with. The major downside, though, might be surprising:

It only works with (callable) python functions.

This is annoying5, but it really isn’t that much of a hassle to rewrite your function using def. We also need to provide the dominant singularities by hand, which is also a bit of work6.

If and when I get around to automatically finding singularities, I’ll update this code, but for now there’s much less boilerplate here than in the previous case:

The documentation for this function is also quite good, so I’ll stick to one example. For the catalan numbers, we’ll have:

This looks like it has singularities at $0$ (which makes the denominator vanish) and $1/4$ (which makes the $\sqrt{\cdot}$ vanish), but actually it only has the singularity at $1/4$. The singularity at $0$ is removable. So we call

which outputs

\[\frac{1}{\sqrt{\pi}} 4^{n} n^{-\frac{3}{2}} - \frac{9}{8 \, \sqrt{\pi}} 4^{n} n^{-\frac{5}{2}} + \frac{145}{128 \, \sqrt{\pi}} 4^{n} n^{-\frac{7}{2}} - \frac{1155}{1024 \, \sqrt{\pi}} 4^{n} n^{-\frac{9}{2}} + \frac{36939}{32768 \, \sqrt{\pi}} 4^{n} n^{-\frac{11}{2}} + O\!\left(4^{n} n^{-6}\right)\]

A quick massage turns this into the (maybe not so) familiar

\[\frac{4^n}{\sqrt{\pi}} \left ( n^{-3/2} - \frac{9}{8} n^{-5/2} + \frac{145}{128} n^{-7/2} - \frac{1155}{1024} n^{-9/2} + \frac{36939}{32768} n^{-11/2} \pm O \left ( n^{-13/2} \right ) \right )\]

which you can find as formula (20) here.

Lastly, we can combine these different asymptotic approximations into one function. In my init.sage file, which you can find here, I’ve defined the function

This is definitely too brittle to give to other people, but it works fine for my purposes. One day I’d like to make it a bit more robust (this is not how try/except statements are supposed to be used), and ideally make it a bit more automatic (in particular following the discussion of singlevar_asy above). Another option might be to try and get it working with the singularlity_analysis function, which would require a certain amount of parsing… Anyways, as it stands it’s a plenty servicable substitute for maple’s asympt function, and I’m glad I took the time to code it up.

If anyone has any other ideas for making this more usable, or encounters any problems, definitely let me know!

If not, take care and stay warm!

I’ll see you all in the next one ^_^.

  1. You can read about these algorithms and their proofs

    • here (for the smooth case)
    • here (for the singular case)

    There’s also Pemantle and Wilson’s Analytic Combinatorics in Several Variables or Melczer’s An Invitation to Analytic Combinatorics, which both look quite good. Unfortunately I haven’t had time to read either. 

  2. I’m curious if I can automate the search for these singularities… I feel like it shouldn’t be too hard, but I need to think about it. 

  3. As an aside, I keep reading this as “asymptomatic” out of the corner of my eye, which definitely tells you something about my current mental state… 

  4. Sorry I’m not making these evaluate inline like I normally do. You’ll have to trust me (or check yourself!) that the computation really is quite fast.

    The issue is that I don’t want to copy/paste the definition of multivar_asy into each cell. I think that clutters things up.

    I could use linked sage-cells, but linked cells have to be evaluated in the right order, and aren’t allowed to be auto-evaluated. I think that’s a bit user-unfriendly, so I’ve decided to just show the inputs/outputs (which is what really matters in any case). 

  5. And part of me wants to go in and try to fix it 

  6. This is something I really think should be automate-able, though, and I’ll probably fight with it some more when I have the time.