Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
    
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 Jan 18 '12 at 18:46
    
@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. Jan 18 '12 at 21:47
    
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 Jan 18 '12 at 21:53
    
@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. Jan 18 '12 at 22:00

4 Answers 4

up vote 23 down vote accepted

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}
share|improve this answer
    
Great work! :D I'll wait for other answers, and accept this if there are no better routes. –  J. M. Jan 17 '12 at 23:55
    
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 Bantegui Jan 17 '12 at 23:59
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 Jan 18 '12 at 0:04
    
@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 Bantegui Jan 18 '12 at 0:05
    
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. Jan 18 '12 at 0:14

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.

share|improve this answer
    
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 Feb 6 '12 at 9:59
    
Done :-) Is it okay now? –  faleichik Feb 6 '12 at 10:08
    
+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 Feb 6 '12 at 10:18
    
Yes, I've forgotten to remove this from the previous version. –  faleichik Feb 6 '12 at 10:23
    
Nice code - +1. –  Leonid Shifrin Feb 6 '12 at 11:57

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.

share|improve this answer
    
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 Jan 18 '12 at 16:49
1  
+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 Jan 18 '12 at 16:51
    
@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 Feb 6 '12 at 10:13
    
@Leonid My code was so slow because of the wrong placing of Expand, corrected now. –  Szabolcs Feb 6 '12 at 10:15
    
I see. That new solution is nice indeed. It did not cross my mind to use pure functions here. –  Leonid Shifrin Feb 6 '12 at 11:58

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.

share|improve this answer
    
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. Jan 17 '12 at 23:50
    
I 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 Jan 18 '12 at 0:30

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.