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

Question

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?

Answer 1

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}

Answer 2

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.

Answer 3

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.

Answer 4

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.

Answer 5

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]];