Dynamic Programming with delayed evaluation



By using dynamical programming, we can save intermediate steps for recursive relations, as in

f[n_]:= f[n] = f[n-1] + f[n-2]

However, this only stores definitions for explicit values of n. But I need to be able to store function definitions, where the number of parameters of the function depends on each step. More specifically, I have something like $$ t^{a_1 \cdots a_n} = A \, t^{a_1 \cdots a_{n-2}} + B \, t^{a_2 \cdots a_n} \\ t^a = C^a \\ t^{ab} = D^{ab} $$

Where all $a_i$ are indices with unspecified values. This recursive definition allows me to express any $t^{a_1 \cdots a_n}$ as a combination of $C$'s and $D$'s. It's easy to implement this in mathematica:

t[a_]:= c[a]
t[a_,b_]:= d[a,b]
t[a__]:= A t@@Most[Most[{a}]] + B t@@Rest[{a}]

but for high $n$, calculation times become very long (and especially because my function is a bit more complicated than this). This is where we normally would use dynamic programming, something like:

t[a__]:= t[a] = A t@@Most[Most[{a}]] + B t@@Rest[{a}]

but this is not what we need, because if I execute


it only makes a definition for a, b and c literally, meaning that t[b,a,c] will have to be recalculated recursively. This is because it is stored as Set, and not a SetDelayed:

t[a,b,c] = A c[a] + B d[b,c]

What I would need is some kind of magic dynamic programming that stores function definitions as in

t[a__]:= Magic.... =A t@@Most[Most[{a}]] + B t@@Rest[{a}]

such that when I execute


it stores

t[a_,b_,c_] := A c[a] + B d[b,c]


I've looked into Leonid's answer for this question, but I can't easily see how to adapt it for a growing number of arguments.

Any ideas? Many thanks in advance!


Thanks to Leonid's answer, I now understand a bit how this can be solved. Unfortunately, the way I'd need it, is a bit more complicated. I have several function definitions, with Condition, PatternTest, Optional and next to the recursive ('growing') variables I also have some that aren't.

At first I thought I'd figure it out myself based on Leonid's answer (hence the accept), but I didn't. The idea is to make a wrapper memoize, which I wrap around my function definition such that it memorises the function definitions. An example in pseudocode:

    t[dot[a__], b_?EvenQ, c_:1]/;b!=c := .... some recursive function of t ...

The first thing I tried is to retrieve all patternames, and replace them with Unique[] ones:

SetAttributes[memoize, HoldAllComplete]
memoize[expr_SetDelayed] :=
    (* first we retrieve the lhs and rhs *)
    With[{funcLHS = Hold[expr] /. Hold[SetDelayed[x_, y_]] :> Hold[x] , 
          funcRHS = Hold[expr] /. Hold[SetDelayed[x_, y_]] :> Hold[y]},
        (* next we retrieve the names of the patterns used *)
        With[{patternNames = First /@ Cases[funcLHS, _Pattern, Infinity]},
            (* we make some locally scoped patterns *)
            With[{varsExt = Table[Unique[], {Length[patternNames]}]},
                (* and we express the lhs and rhs in function of the former locally scoped patterns *)
                With[{lhs = funcLHS /. Rule @@@ Transpose[{patternNames, varsExt}], 
                      rhs = funcRHS /. Rule @@@ Transpose[{patternNames, varsExt}]},


And this works, but it doesn't do anything useful. Now I would need to replace the central part with Leonid's answer somehow, which means I have to inject vars into lhs and rhs, with the additional difficulty that Length[vars] won't equal Length[varsExt]. So this is what I tried:

SetAttributes[memoize, HoldAllComplete]
memoize[expr_SetDelayed] :=
    With[{funcLHS = Hold[expr] /. Hold[SetDelayed[x_, y_]] :> Hold[x] , 
          funcRHS = Hold[expr] /. Hold[SetDelayed[x_, y_]] :> Hold[y]},
        With[{patternNames = First /@ Cases[funcLHS, _Pattern, Infinity]},
            With[{varsExt = Table[Unique[], {Length[patternNames]}]},
                With[{lhs = funcLHS /. Rule @@@ Transpose[{patternNames, varsExt}], 
                      rhs = funcRHS /. Rule @@@ Transpose[{patternNames, varsExt}]},

                    {lhs, Hold[
                        With[{vars = Table[Unique[], {Length[{#}]}] & /@ varsExt},
                            With[{pts = (Pattern[#1, _] &) /@ vars},
                                {lhs/.MapThread[ Pattern[#1, _] -> Sequence @@ (Function[x, Pattern[x, _]] /@ #2) & , {varsExt, vars} ]
                                rhs/.MapThread[ #1 -> Sequence @@ #2 & , {varsExt, vars} ]
                                }/.{Hold[x_],Hold[y_]}:>Set[x,y] ;
                     } /. {Hold[x_], Hold[y_]} :> SetDelayed[x, y] 

But the pattern replacing on lhs doesn't work, and it doesn't take care of Condition, Optionaland PatternTest. So I'm totally stuck, although I have the feeling that it's only a few changes to make it work..


Posted 2013-06-19T17:44:40.220

Reputation: 1 025



I will offer a rather cryptic solution using nested version of the injector pattern, but it should be possible to also rewrite it using more readable methods, if really needed.


Here is the code:

t[a_]:=  c[a]
t[a_,b_]:= d[a,b]
                  {brhs___}:>Set@@Hold[t[x],A t[arhs]+B t[brhs]])

We get:

t[a, b, c, d]

(* A d[a, b] + B (A c[b] + B d[c, d])  *)

and now we have the following definitions generated in the process:



t[$21_,$22_,$23_]=A c[$21]+B d[$22,$23]

t[$17_,$18_,$19_,$20_]=A d[$17,$18]+B (A c[$18]+B d[$19,$20])

   With[{pts=(Pattern[#1,_]&)/@vars},pts/. {x__}:>(Most[Most[vars]]/. {arhs___}:>
      (Rest[vars]/. {brhs___}:>Set@@Hold[t[x],A t[arhs]+B t[brhs]]));t[a]]]

How it works

Basically, dummy symbols were used to programmatically create a pattern for the l.h.s., and then to create the corresponding memoized definition.

In more details, here is what happens: whenever t[x,y,...] is called, the general definition t[a__]:=... is only invoked if there isn't yet a definition for the exact number of input arguments Length[{a}]. In the process of evaluation of t[a], such definition is created, so all subsequent times calls with the same number of arguments will use the new specialized definition, not the general one.

The general idea here can be expressed in the following pseudo-code:

  Module[{helper vars},
    t[var1_,...,varn_] = rhs[var1,...varn];

where we first construct a specialized definition for a fixed number of input arguments, and then use it right away for our specific arguments a.

Now, how it works: let's say we start with the call


so that

Clear[m, n, k];
a = Sequence[m, n, k];

models what a becomes. We first generate a set of dummy variables:

vars = Table[Unique[], {Length[{a}]}]

(* {$3, $4, $5} *)

and the corresponding patterns:

pts = (Pattern[#1, _] &) /@ vars

(* {$3_, $4_, $5_} *)

What we would like now to do is to construct this code:

t[$3_, $4_, $5_] = A * t[$3] + B * t[$4, $5]

In principle, we could just try to use

Evaluate[t @@ pts] = A t @@ Most[Most[vars]] + B t @@ Rest[vars]

to create the definition. However, we must be more careful. Direct evaluation like Evaluate[t @@ pts] will lead to the general definition t[a__] applied now for t[$3_, $4_, $5_] (since patterns are also arguments), and this is something we don't want. In other words, we want to prevent unwanted evaluation during assignment.

This is where injector pattern comes handy. Nested injector pattern is used to collect more than one part of an expression within a single nested rule, and is actually an overkill here. The definition can be constructed with this simpler code:

pts /. {x__} :> Set @@ Hold[t[x], A t @@ Most[Most[vars]] + B t @@ Rest[vars]]; 

The construct pts /. {x__} :> ... is used to inject x (in this case, a sequence of paterns $3_, $4_, $5_ ) into the right-hand side before thr r.h.s. is allowed to evaluate. Here we use it to programmatically construct t[$3_, $4_, $5_] - the l.h.s. of the new definition. The construct Set@@Hold[lhs, rhs] is used to fool the built-in variable renaming mechanism, which would otherwise rename our variables in an attempt to protect them from name collisions (which is exactly what we don't want here).

So, the final code can take a simpler form and look like

               {x__}:>Set@@Hold[t[x],A * t @@ Most[Most[vars]]+B * t @@ Rest[vars]];

where, as explained, the first line of code in the inner With creates a specialized definition, and then this definition is immediately used on the actual arguments a. All subsequent times, when called with the same number of arguments as in a, t will use the specialized defintion automatically, and not the general t[a__].

Leonid Shifrin

Posted 2013-06-19T17:44:40.220

Reputation: 108 027

1wow! Genius, but far beyond my leave as I don't fully understand it (not to say almost not at all).. You create a recursive definition by chaining :>'s, is that right? I am a bit puzzled with the mechanics behind it (why use Setinstead of SetDelayed?), how do you ensure it will not redefine t[a,b,c,d] when you invoke it a second time? Forgive my ignorance, I am trying to understand as to be able to expand it for my problem.. – freddieknets – 2013-06-20T11:53:12.610

Ok sorry, I got it partially now.. Can I chain more replacement rules like this if I want to use more combinations of the variables (f.e. Most[Most[vars]], Rest[vars] but also Most[Rest[vars], vars[[1]]etc)? – freddieknets – 2013-06-20T12:03:23.040

@freddieknets Yes, you can. In fact, the solution could have been simpler. Please see my edit, where I added the explanation section. – Leonid Shifrin – 2013-06-20T12:48:29.217

Thanks a lot for your effort! I understand now how it works. However I am having troubles applying it to my problem, as I would need it more generally. I expanded my question, see my edit, so maybe you can have a look (that is, if you feel like it and have the time for it, I wouldn't want you to waste too much time on my mathematica-programming-insight shortage)? Thanks anyway for your help so far! – freddieknets – 2013-06-21T12:50:58.380

@freddieknets It is not clear from your edit what is your real problem now, or what you are trying to accomplish. It also looks like you are over-generalizing it. Provide a minimal specific self-contained example where you clearly state what you are trying to do and what the difficulty is. Best if you do it as a separate follow-up question. I don't promise that I will have the time to go into this in the next few days, though - but making it a separate question will increase the chances of it to be answered by somebody. – Leonid Shifrin – 2013-06-22T19:48:32.967


You can also try this code. It uses only explicit variables a1,a2,a3... in the definitions. The auxiliary function ComputeDefinition recursively constructs f for given length n of the sequence.


f[{a1_}] = c[a1];
f[{a1_, a2_}] = d[a1, a2];
ComputeDefinition[1] = 1;
ComputeDefinition[2] = 1;
ComputeDefinition[10] = 1;

ComputeDefinition[n_] := (
   f[Table[ToExpression["a" <> ToString[i] <> "_"], {i, 1, n}]] =
   A f[Table[ToExpression["a" <> ToString[i]], {i, 1, n - 2}]] +
   B f[Table[ToExpression["a" <> ToString[i]], {i, 2, n}]];
   ComputeDefinition[n + 1])

(* 1 *)




f[{a1_,a2_,a3_}]=A c[a1]+B d[a2,a3]

f[{a1_,a2_,a3_,a4_}]=A d[a1,a2]+B (A c[a2]+B d[a3,a4])

f[{a1_,a2_,a3_,a4_,a5_}]=A (A c[a1]+B d[a2,a3])+B (A d[a2,a3]+B (A c[a3]+B d[a4,a5]))

f[{t1, t2, t3, t4, t5, t6}]

A (A d[t1, t2] + B (A c[t2] + B d[t3, t4])) + B (A (A c[t2] + B d[t3, t4]) +
B (A d[t3, t4] + B (A c[t4] + B d[t5, t6])))

Vahagn Poghosyan

Posted 2013-06-19T17:44:40.220

Reputation: 532