How can I implement dynamic programming for a function with more than one argument?

54

44

Dynamic programming is a technique for avoiding the repeated computation of the same values in a recursive program. Each value computed is immediately stored. If the value is needed again, it is not computed but simply looked up in the table. (1)

I use orthogonal polynomials a fair bit in my work. Since Mathematica supports only the classical ones, I often have to write my own functions. For instance, the monic Charlier polynomials satisfy the three-term recurrence

$$C_{n+1}^{(a)}(x)=(x-a-n)C_n^{(a)}(x)-an C_{n-1}^{(a)}(x)$$

with $C_0^{(a)}(x)=1$ and $C_1^{(a)}(x)=x-a$.

If I want to be able to use monic Charlier polynomials in Mathematica, I can do this:

CharlierC[0, a_, x_] := 1;
CharlierC[1, a_, x_] := x - a;
CharlierC[n_Integer, a_, x_] := (x - a - n + 1) CharlierC[n - 1, a, x] -
                                 a (n - 1) CharlierC[n - 2, a, x]

The problem with this route, of course, is that the effort expended to generate, say, CharlierC[20, a, x] can't be used for evaluating CharlierC[50, a, x]. For a one-argument recursive function (e.g. Fibonacci), dynamic programming is fine and dandy for saving evaluation effort. For a multiple-argument function, imagine what would happen if one had used the definition CharlierC[n_Integer, a_, x_] := CharlierC[n, a, x] = (* stuff *) and then executed Plot[{CharlierC[5, 1, x], CharlierC[6, x, 2]}, {x, -1, 1}].

Is there a way to reap the benefits of dynamic programming on a multiple-argument function, while storing only results where the recursion variable (n in the Charlier example) changes?

J. M.'s ennui

Posted 2012-01-17T23:27:20.203

Reputation: 115 520

As I am sure you know CharlierC[n, x, a] == (-1)^n HypergeometricU[-n, 1 - n + x, a], so there is no need to use the recursive definition for this one. – Sasha – 2012-01-18T18:46:36.700

Can you give an update of any of the functions available in versions up to 11.1 provide more orthogonal polynomial other that the classic ones? – Jose Enrique Calderon – 2017-04-27T13:16:11.733

@Jose, I don't think anything new came up after version 6. – J. M.'s ennui – 2017-05-10T22:02:52.897

@Sasha: I certainly do, but I also deal with even more complicated polynomials than Charlier polynomials; I just wanted a simple example of an unsupported polynomial system. – J. M.'s ennui – 2012-01-18T21:47:01.197

Out of curiosity, is this related to probability anyhow? Like I just learned that Charlier polynomials are orthogonal with respect to Poisson measure, which makes them well suited for developing distribution-specific quadrature rules

– Sasha – 2012-01-18T21:53:05.587

@Sasha: I'm actually working with Askey-Wilson class polynomials, as well as $q$-analogs of classical orthogonal polynomials, for one of my personal projects. I know about being able to use QHypergeometricPFQ[], but it seems that even a recursive definition can sometimes go a bit faster than using that. But yeah, even the humble Charlier polynomial is fascinating! :)

– J. M.'s ennui – 2012-01-18T22:00:46.620

Answers

41

Yes, there is, although the speed-up is not as dramatic as for 1D memoization:

ClearAll[CharlierC];
CharlierC[0, a_, x_] := 1;
CharlierC[1, a_, x_] := x - a;
CharlierC[n_Integer, a_, x_] :=
  Module[{al, xl},
    Set @@ Hold[CharlierC[n, al_, xl_],          
        Expand[(xl - al - n + 1) CharlierC[n - 1, al, xl] - 
             al (n - 1) CharlierC[n - 2, al, xl]
        ]];
    CharlierC[n, a, x]
];

(Thanks to @Mike Bantegui for pointing out the wastefulness of Simplify, which has now been removed). What you memoize here are function definitions.Expand is used to not accumulate the complexity too fast. The idea is that I first create a new pattern-based definition, using a number of tricks to fool the scoping variable - renaming mechanism but localize pattern variables, and then evaluate this definition.

For example:

In[249]:= CharlierC[20,a,x];//Timing
Out[249]= {0.063,Null}

In[250]:= CharlierC[25,a,x];//Timing
Out[250]= {0.078,Null}

While with clear definitions:

In[260]:= CharlierC[25,a,x];//Timing
Out[260]= {0.094,Null}

Here are a first few generated definitions:

In[262]:= Take[DownValues[CharlierC],4]
Out[262]= 
{HoldPattern[CharlierC[0,a_,x_]]:>1,
 HoldPattern[CharlierC[1,a_,x_]]:>x-a,
 HoldPattern[CharlierC[2,al$4106_,xl$4106_]]:>
      al$4106^2-xl$4106-2 al$4106 xl$4106+xl$4106^2,
 HoldPattern[CharlierC[3,al$4105_,xl$4105_]]:>
   -al$4105^3+2 xl$4105+3 al$4105 xl$4105+3 al$4105^2 xl$4105
       -3 xl$4105^2-3 al$4105 xl$4105^2+xl$4105^3}

Leonid Shifrin

Posted 2012-01-17T23:27:20.203

Reputation: 108 027

@Mike Thanks! I came up with this first about 5 years ago, when I needed it to compute functions with a huge number of terms. In those days, I wouldn't make a mistake with Simplify, but I got rusty since I quit Physics :) – Leonid Shifrin – 2012-01-18T00:38:47.980

@Mike Well, I'd keep that, but as you wish of course. – Leonid Shifrin – 2012-01-18T00:42:43.433

I timed both versions, and yours is significantly slower on my machine. 0.016 seconds (mine) vs 1.185 seconds (yours) for n = 20. Still, very cool idea you have. – Mike Bailey – 2012-01-17T23:59:06.437

1@Mike My code creates definitions for absolutely arbitrary a. Also, some time is spent on Expand and Simplify. It is slow the first run, but instant all the following runs. If your task is to vary a (say, plot when a is in some range), I only need to compute stuff once, and it will work for every value of a. If there is a single numerical a, then your code wins, but then arguably a is a constant rather than a true function argument, so the memoization is effectively 1D. – Leonid Shifrin – 2012-01-18T00:04:46.460

@LeonidShifrin: If you want to vary a and x, you can simply evaluate it and attach it to a variable: poly = CharlierC[100, a, x] and then evaluate it for any arbitrary {a, x}. You just have to be careful to only evaluate CharlierC[n, a, x] for sybmolic a and x for my version, and always attach it to something beforehand. – Mike Bailey – 2012-01-18T00:05:52.947

I also deal with polynomials with more than three arguments, like the Hahn polynomials, and I expect to be varying indices and parameters a lot. So, Leonid's proposal seems to be most expedient for me.

– J. M.'s ennui – 2012-01-18T00:14:28.957

@Mike Good point, I will edit. Simplify was taking the most time, I should have tried Expand alone. As to your first point: since in my case you compute all the previous ones in the process, you need almost no time to compute say CharlierC[101, a, x] fully symbolically then. – Leonid Shifrin – 2012-01-18T00:19:25.687

@LeonidShifrin: Just for the record, the way you did this is very cool and borderline black magic. – Mike Bailey – 2012-01-18T00:25:05.960

25

Nice question. This is my suggested implementation. Evaluate all code at once.

Clear[CharlierC, "CharlierC`*"]

CharlierC (* create symbol in current context *)

Begin["CharlierC`"];

implementation[0] := 1;
implementation[1] := x - a;
implementation[n_Integer] := 
 implementation[n] = Expand[(x - a - n + 1) implementation[n - 1] - 
                            a (n - 1) implementation[n - 2] ]

CharlierC[n_Integer ? NonNegative, ai_, xi_] := 
 Block[{a, x}, implementation[n] /. {a -> ai, x -> xi}]

End[];

The Block isn't strictly necessary as a and x are already created in the CharlierC` context where they should be safe.

This function memoizes the symbolic representation of CharlierC only, and only for each n (not a and x). Then substitutes in whatever variables or numbers you pass in.

Szabolcs

Posted 2012-01-17T23:27:20.203

Reputation: 213 047

Now let me look at the other answers which I didn't check on purpose before posting (so the fun won't be spoiled!) – Szabolcs – 2012-01-18T16:49:03.740

2+1. This is simpler than what I suggested. Andy Ross did the first step just like you, but did not do replacements and deleted his answer, which is a pity. I only later realized that his (and now yours) code is no less general than mine, and perhaps also faster. – Leonid Shifrin – 2012-01-18T16:51:39.887

@Leonid My code is much slower than yours, but it is not completely clear to me why. There was also a new, pure function based solution posted today. Now it seems so obvious that pure functions are the simplest solution here ... – Szabolcs – 2012-02-06T10:13:27.193

2@Leonid My code was so slow because of the wrong placing of Expand, corrected now. – Szabolcs – 2012-02-06T10:15:15.820

I see. That new solution is nice indeed. It did not cross my mind to use pure functions here. – Leonid Shifrin – 2012-02-06T11:58:57.340

22

I would also suggest to use pure functions here:

CharlierC[0] = 1 &;
CharlierC[1] = #2 - #1 &;
CharlierC[n_Integer] := (CharlierC[n] = 
    Evaluate[
      Expand[(#2 - #1 - n + 1) CharlierC[n - 1][#1, #2] - #1 (n - 
           1) CharlierC[n - 2][#1, #2]]] &);
CharlierC[20][a, x] // AbsoluteTiming

(* ==> {0.0312414, 
 a^20 - 121645100408832000 x - 128047474114560000 a x - 
  67580611338240000 a^2 x - 23851980472320000 a^3 x - 
  6335682312960000 a^4 x - 1351612226764800 a^5 x - 
  241359326208000 a^6 x - 37132204032000 a^7 x - 
  5028319296000 a^8 x - 609493248000 a^9 x - 67044257280 a^10 x - 
  6772147200 a^11 x - 634888800 a^12 x - ...
*)

I haven't made comparisons, but this seems to be faster.

faleichik

Posted 2012-01-17T23:27:20.203

Reputation: 12 161

The disadvantage with this solution is that it memoizes for the second argument as well (a). If we evaluate CharlierC[5,a], the cached intermediate results can't be used by CharlierC[5,b]. I think it should be possible to fix this. Can you update the solution? – Szabolcs – 2012-02-06T09:59:44.573

Done :-) Is it okay now? – faleichik – 2012-02-06T10:08:28.803

+1, great solution! But to nitpick a bit more, Evaluate is not needed in CharlierC[1] = Evaluate[#2 - #1] & I placed Expand at the wrong position in my solution which caused a serious performance problem. I corrected it now. – Szabolcs – 2012-02-06T10:18:22.867

Yes, I've forgotten to remove this from the previous version. – faleichik – 2012-02-06T10:23:39.177

Nice code - +1. – Leonid Shifrin – 2012-02-06T11:57:41.833

16

You can utilize Dynamic Programming in the following way:

CharlierC[0, a_, x_] := 1
CharlierC[1, a_, x_] := x - a
CharlierC[n_Integer, a_, x_] := CharlierC[n, a, x] = Expand[Expand[(x - a - n + 1) CharlierC[n - 1, a, x]] - Expand[a (n - 1) CharlieC[n - 2, a, x]]]

It basically creates an in-memory store of previous evaluated function values instead of evaluating it every time you call it.

Mind you, if you actually tried evaluating CharlierC[50, a, x] it would still take a very long time. That's a very long polynomial you're trying to evaluate.


Update

The addition of the Expand[..] forces Mathematica to expand and combine like terms. It's actually possible to evaluate CharlierC[50, a, x].

If you wanted to really use that expression, here's how I would do it:

poly = CharlierC[50, a, x];

a = 0.01;
Plot[poly, {x, ...}];

Here's some timings. In between each run I did a ClearAll and evaluated the definition:

AbsoluteTiming[CharlierC[20, a, x];] (* 0.0250014 *)
AbsoluteTiming[CharlierC[40, a, x];] (* 0.1510087 *)
AbsoluteTiming[CharlierC[60, a, x];] (* 0.5150295 *)

What I would recommend is combining the ideas that Leonid and I have:

ClearAll[CharlierC];
CharlierC[0, a_, x_] := 1;
CharlierC[1, a_, x_] := x - a;
CharlierC[n_Integer, a_, x_] := Module[{al, xl},
   Set @@ Hold[CharlierC[n, al_, xl_],
     Expand[
      Expand[(xl - al - n + 1) CharlierC[n - 1, al, xl]] - 
       Expand[al (n - 1) CharlierC[n - 2, al, xl]]
      ]
     ];
   CharlierC[n, a, x]];

The use of Expand works much faster than Simplify.

Mike Bailey

Posted 2012-01-17T23:27:20.203

Reputation: 1 875

Okay, try executing CharlierC[10, 3, 1] and subsequently ?CharlierC. Messy, no? Using that version of CharlierC[] within something like Plot[] will make things worse. – J. M.'s ennui – 2012-01-17T23:50:35.483

1I did actually use Expand in my initial answer as well - it was Simplify which was slow. And, you only need a single top-level one I think - I also just made an edit. But, again, thanks for your point, and +1. – Leonid Shifrin – 2012-01-18T00:30:27.340

1

Decided to add one more version: what happened is that my question How to memoize with patterns? has been marked as a duplicate of this one, which I agree to. However I believe there are more general situations where it should work. I believe my version covers them.

toPatterns[x_] := 
 x /. Map[# -> Pattern[#, _] &, 
   Union@Select[Level[x, {-1}], Head[#] === Symbol &]]
ClearAll[myCharlierC]
myCharlierC[0, a_, x_] := 1;
myCharlierC[1, a_, x_] := x - a;
myCharlierC[n_Integer, a_, x_] := 
  myCharlierC[n, toPatterns[a], toPatterns[x]] = 
   Expand[(x - a - n + 1) myCharlierC[n - 1, a, 
       x] - (n - 1) a  myCharlierC[n - 2, a, x]];

მამუკა ჯიბლაძე

Posted 2012-01-17T23:27:20.203

Reputation: 1 831