Memoization of a function defined by a recurrence relation [duplicate]

Question

I have a function which is defined by the following recurrence relation $$h_{n}(x)=h_{n-1}(x)+\frac{\mathrm{e}^{-x^2}}{2^n n!}H_{n}(x)H_{n-1}(x)$$ with the initial condition $h_{0}(x)=0$ and where the $H_{n}(x)$ are the Hermite polynomial of order $n$.

I am trying to compute numerically a table of these function for different values of $x$. Is there a way to compute these and keep them in memory in order to not have to compute them each time for different values of $x$. In other words, I want to store only results where the recursion variable $n$ changes? (let's say for all $n<N_{max}$).

I tried this

Clear[h];
h[0, x_] = 0.;
h[n_, x_] := (h[n, x] = h[n - 1, x] + 
Exp[-x^2]/(2^n*n!)*HermiteH[n,x]*HermiteH[n - 1, x]);

but it seems that it evaluate it every time you call it with a different argument $x$.

Thanks for your help!

Answer 1

Here is something that might be more elaborated than really necessary for your task (for which I think Marius answer might work well enough) but shows some helpful techniques. It basically does the same thing that Marius solution does: it separates the recursion in n from the functional dependence in x but tries to do this a bit more robust and not change the call syntax of h:

Clear[hh];
hh[0] = Function[0];
hh[n_] := hh[n] = Block[{x},
    Function @@ {x,Simplify[
       hh[n - 1][x] + 
        Exp[-x^2]/(2^n*n!)*HermiteH[n, x]*HermiteH[n - 1, x]
     ]
    }
];
h[n_,x_]:=hh[n][x]

hh is a helper function which returns a Function of x only which it generates on the fly and for which it does the required memoization.

Here are some further explanations:

• the outer Block is there so we can evaluate the function body with the symbol x even when there are definitions for x. That is just a safety measure.
• Applying Function to a list is just a trick to evaluate the function body, there are other techniques to achieve the same thing.
• I have not really tested if it really is achieving what I want, but I think in cases like that it might be more efficient to actually do a Simplify of the function body expression as that will avoid that these expressions become extensively large and the Simplify should be reasonably cheap for each step of the recursion. I find that an interesting aspect as it is an optimization that seems only possible when you can manipulate symbolic expressions. You might want to experiment with that.

Last but not least it seems worth noting that memoization is not something that has explicitly built into Mathematica but rather is a specific application of the possibility to generate additional definitions on the fly. So you can not expect that your original definition could do anything but what you have experienced (there is no specific "memoization" functionality which could detect that the x dependency could be separated). On the other hand there is the possibility to make definitions for whatever you want on the fly, you are not restricted to make definitions for the same function with the given arguments although in the form of standard memoization that is probably the most prominent use case.

As the example in your comment (h[250,3.]) shows when using a numeric argument the above approach will not only be slow but even can return wrong results. To get a correct result what you should do is convert to a numeric value only as late as possible, then you can even get as many correct digits as you need with something like N[h[250,3],50]. But as often in such situations when you are interested in (machine precision) numeric values for numeric arguments only, a "brute force" numeric way can be much faster and still correct. Here is code for a numeric approach which seems to be simple and very fast for that case:

hn[m_, x_] := Total@Table[
  Exp[-x^2]/(2^(n)*(n)!)*HermiteH[n, x]*HermiteH[n - 1, x], {n, m}
]

Answer 2

You mean something like this?

Clear[h, x];
h[0] = 0.;
h[n_] := h[n] = h[n - 1] + 
Exp[-x^2]/(2^n*n!)*HermiteH[n,x]*HermiteH[n - 1, x];

tab = Table[h[en] /. x -> ex, {en, 1, nmax}, {ex, xmin, xmax, dx}]

Answer 3

Original post

This is not an answer to the question, but an extended comment. I was curious if the recursion problem (without memorization) would admit an explicit solution in Mathematica. I can withdraw my contribution if you want.

Here we solve the recursion problem in the "natural way", viz. using RSolve. Also we have included the dependence on the variable x.

The recursion equation including the inital condition is

eq = {hr[0, x] == 0, 
   hr[n, x] == 
    hr[n - 1, x] + 
     Exp[-x^2]/(2^n*n!)*HermiteH[n, x]*HermiteH[n - 1, x]};

The solution

hrs[n_, x_] = hr[n, x] /. RSolve[eq, hr[n, x], n][[1]] // FullSimplify

is an expression using DifferenceRoot[] containing a linear difference equation of 5th order which I don't reproduce here.

It is not ecessary to understand it, since the resulting function can simply be used.

Examples

1) we have operable dependence on the second argument

hrs[1, x]
hrs[2, y]

(* Out[271]= E^-x^2 x

Out[272]= 1/2 E^-y^2 (y + 2 y^3) *)

2) Some more terms

t = Table[hrs[k, z], {k, 0, 4}]

(* Out[270]= 
{0, E^-z^2 z, 1/2 E^-z^2 (z + 2 z^3), 1/3 E^-z^2 z (3 - z^2 + 2 z^4), 
 1/24 E^-z^2 z (15 + 34 z^2 - 20 z^4 + 8 z^6)} 
*)

EDIT

I noticed only now the interesting fact that the functions h are solutions of a linear difference equation of the order 5 in n. And this inlcudes the Hermite polynomials.

At the same time this exercise took off my slight horror of the function DifferenceRoot[].

The main expression under the DifferenceRoot[] is the difference equation and the initial conditions:

$$\left\{(-n-1) (n+2) y(n)\\+(n+2) \left(n+2 x^2+2\right) y(n+1)\\+\left(2 n^2-6 n x^2+11 n-16 x^2+16\right) y(n+2)\\+\left(-2 n^2+6 n x^2-13 n+20 x^2-22\right) y(n+3)\\-(n+4) \left(n+2 x^2+4\right) y(n+4)\\+(n+4) (n+5) y(n+5)=0,\\y(0)=0,y(1)=2 x,y(2)=2 x^3+x,y(3)=\frac{2}{3} x \left(2 x^4-x^2+3\right),\\y(4)=\frac{1}{12} x \left(8 x^6-20 x^4+34 x^2+15\right)\right\}$$

Notice that except for the case n = 1 the initial conditions are not given by Hermite polynomials.

The surprising incorporation of the Hermite polynomials into the recursion relation for h might have to do with the fact that HermiteH[] is itself a solution of a linear difference equation:

ff[n_] = f[n] /. 
  First[RSolve[
     f[n + 1] == 2 x f[n] - 2 n f[n - 1] && f[0] == 1 && f[1] == 2 x, 
     f[n], n] // Simplify]

(* Out[100]= HermiteH[n, x] *)

Answer 4

It's tricky in this case to write a good definition along the lines sought by the OP. Albert Retey's composite method is perhaps an easier way to solve such problems. The difficulty I ran into had to do with the ordering of DownValues. The down value

h[n, y_] = h[n - 1, y] + Exp[-y^2]/(2^n*n!)*HermiteH[n, y]*HermiteH[n - 1, y]

is equivalent to

HoldPattern[h[n, y_]] :> 
  Evaluate[h[n - 1, y] + Exp[-y^2]/(2^n*n!)*HermiteH[n, y]*HermiteH[n - 1, y]]

It needs to be put before the general down value for h[n_Integer?Positive, x_] below, so instead of using Set, we have to prepend it to DownValues[h]. Mathematica won't do this by default, since it considers the general pattern to as specific (or more so) than the one we want to be considered more specific, h[n, y_].

Clear[h];
h[0, x_] = 0;
h[n_Integer?Positive, x_] := (Block[{y},
    PrependTo[DownValues[h],
     HoldPattern[h[n, y_]] :> 
      Evaluate[
       h[n - 1, y] + Exp[-y^2]/(2^n*n!)*HermiteH[n, y]*HermiteH[n - 1, y]]
     ]];
   h[n, x]);

If you run the following, the second Trace will show that a stored expression is used to compute the return value. (The first trace is rather long [image].)

h[2, 2] // Trace
h[2, z] // Trace