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

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!

• This is closely related Why does Expand not work within a function? if not simply a duplicate. – Artes Oct 1 '16 at 12:23
• I think this is actually really a better duplicate, the actual question and solution in the other are not really a duplicate IMO... – Albert Retey Oct 2 '16 at 14:49

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

• Thanks @Albert! This is what I was searching for. But now, there is another problem arising. Let me illustrate it. Let us say we want to compute the value $h_{250}(x)$ at the point $x=3$. If we first fix $x=3$ and then run the recurrence until the value $n=250$ we obtain Clear[h2]; x = 3.; h2[0] = 0.; h2[n_] := h2[n] = (h2[n - 1] + Exp[-x^2]/(2^(n)*(n)!)*HermiteH[n, x]*HermiteH[n - 1, x]); h2[250] // AbsoluteTiming {0.011518, 0.809828} – Yasar Sep 30 '16 at 9:33
• While the solution proposed diverges In[911]:= Clear[hh]; hh[0] = Function[0]; hh[n_] := hh[n] = Block[{x}, Function @@ {x, hh[n - 1][x] + Exp[-x^2]/(2^n*n!)*HermiteH[n, x]*HermiteH[n - 1, x]}]; h[n_, x_] := hh[n][x] In[915]:= h[340, zz]; // AbsoluteTiming Out[915]= {3.140634, Null} In[916]:= h[250, zz] /. zz -> 3. Out[916]= -5.7134*10^12 and this is due to the fact that we have now a polynomial expression (times $e^{-x^2}$) with very large coefficients. I am thinking about a way to avoid this divergence... – Yasar Sep 30 '16 at 9:40
• An easy one is h[250,3]//N but of course you will find that the approach will then be less efficient than to do the whole recursion with a numeric value for x inserted from the beginning as in you other example... – Albert Retey Sep 30 '16 at 13:16

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

• Hi @Marius, thanks for your reply! I need something like that rather Table[(h[l, y] +h[l, x]), {l, 0, nmax}, {p, l, nmax}]; with $x$ and $y$ two different points. I am then exporting the table for all couple of points $(x,y)$. Sorry if it was not well stated in the first message. – Yasar Sep 30 '16 at 8:18
• @Yasar I think this is quite a nice answer: the functional dependence on nand xare separated, the rest you can add on your own. Besides, what is h[l,y] and h2[l,x] in your formula? You have to be more precise, not just hope someone can guess your thoughts. – yarchik Sep 30 '16 at 8:23
• @yarchik h[l,x] is the function with index $l$ evaluated in $x$. There was a typo in previous comment. The problem here is that, it seems that the function of the variable $x$ is not put into memory which means that when you call it with the argument $y$ it spent time computing again from the recurrence. Hope it is more clear now. – Yasar Sep 30 '16 at 8:34
• So basically, this $h_{n}(x)$ is a certain polynomial $P_{n}(x)$ times $\mathrm{e}^{-x^2}$ and what I want to put in memory are the polynomial $P_{n}(x)$ for $n<nmax$ as a function of $x$ so that when I call it with any argument (let's say $y$) it doesn't go back to the recurrence which is time consuming... – Yasar Sep 30 '16 at 8:47
• Yasar, what you are describing in your last comment is exactly what the code in my answer does, except that what is stored in each h[n] is not a function, but an expression, involving x. There is no useless recurrence going on :) – Marius Ladegård Meyer Sep 30 '16 at 8:58

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] *)


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

• Great, that is a good example for how you can define arbitrary downvalues on the fly as I mentioned in my answer. I was trying to think about a good example which would be relevant for this question but didn't see this one :-). I think it probably is also worth noting that without the restriction to integer positives in the pattern you wouldn't need the DownValues trick, but of course it makes sense to insist on that restriction. – Albert Retey Oct 1 '16 at 21:06
• Another note: I would strongly argue that the automatic ordering here doesn't conform with the documentation, as h[n_Integer,x_] will be considered as more specific than h[1,x_], but h[n_Integer?Positive,x_] will not. I can hardly find any reason why the additional ?Positive would make the pattern more general... – Albert Retey Oct 1 '16 at 21:08
• @AlbertRetey Thanks. Yes, I was surprised (and fooled) by the automatic ordering at first. I have never felt I completely understood the rules for pattern specificity. – Michael E2 Oct 2 '16 at 3:09
• to completely understand I think one would need to look at the source code, the documentation is certainly not detailed enough for that. As with many of these things we have to live with what we get, so it is good to have those tricks like manually setting DownValues at hand and documented in answers like yours... – Albert Retey Oct 2 '16 at 10:27