Dynamic Programming in Mathematica: how to automat

Question

In Mathematica 8.0, suppose I have some constants:

a:=7
b:=9
c:=13
d:=.002
e:=2
f:=1

and I want to use them to evaluate some interlinked functions

g[0,k_]:=0
g[t_,0]:=e
g[t_,k_]:=g[t-1,k]*a+h[t-1,k-1]*b

h[0,k_]:=0
h[t_,0]:=f
h[t_,k_]:=h[t-1,k]*c+g[t-1,k-1]*d

But this is really slow and in need of dynamic programming, or else you get an exponential slowdown:

g[0, k_] := 0
g[t_, 0] := e
g[t_, k_] := g[t, k] = g[t - 1, k]*a + h[t - 1, k - 1]*b

h[0, k_] := 0
h[t_, 0] := f
h[t_, k_] := h[t, k] = h[t - 1, k]*c + g[t - 1, k - 1]*d

Now it's really fast, but if we ever want to change the constants(say, to use this in a Manipulate function) we have to Clear g and h each time. If we had complex interdependencies, it might be really annoying to clear them all out every single time we wanted a value from g and h.

Is there an easy way to run g and h in a Module or Block or similar, so that I can get a fresh result back each time it's evaluated without the exponential slowdown? Or even a fast way to build up a table of results for both g and h in a nice way? As said, I want to be able to compute g and h in a Manipulate function.

Answer 1

Here is one way, using Block as you requested:

ClearAll[defWrap];
SetAttributes[defWrap, HoldFirst];
defWrap[fcall_] :=
  Block[{g, h},
     (* Same defintions with memoization as you had, but within Block*)

     g[0, k_] := 0;
     g[t_, 0] := e;
     g[t_, k_] := g[t, k] = g[t - 1, k]*a + h[t - 1, k - 1]*b;   
     h[0, k_] := 0;
     h[t_, 0] := f;
     h[t_, k_] := h[t, k] = h[t - 1, k]*c + g[t - 1, k - 1]*d;

     (* Our function call, but within a dynamic scope of Block *)
     fcall];

We will use this to give definitions to f and h as

ClearAll[g, h];
g[tt_, kk_] := defWrap[g[tt, kk]];
h[tt_, kk_] := defWrap[h[tt, kk]];

We call now:

In[1246]:= g[20,10]//Timing
Out[1246]= {0.,253809.}

In[1247]:= h[20,10]//Timing
Out[1247]= {6.50868*10^-15,126904.}

There are no global memoized definitions left after each call - Block takes care to destroy them just before the execution exits Block. In particular, here I will change the parameters and call them again:

In[1271]:= 
a:=1
b:=2
c:=3
d:=.01
e:=4
f:=5

In[1279]:= g[20,10]//Timing
Out[1279]= {0.015,0.808192}

In[1280]:= h[20,10]//Timing
Out[1280]= {0.,1.01024}

An alternative to this scheme would be to explicitly pass all parameters like a,b,c,d,e,f to functions, extending their formal parameter lists (signatures), but this has a disadvantage that the older memoized definitions corresponding to different past parameter values would not be automatically cleared. Another problem with that approach is that the resulting code will be more fragile, w.r.t changes in the number of parameters, etc.

EDIT

However, if you want to build a table of results, this could be somewhat faster since you do it once and for all, and in this case you do want to keep all memoized definitions. So, here is the code:

ClearAll[g, h];
g[0, k_, _] := 0;
g[t_, 0, {a_, b_, c_, d_, e_, f_}] := e;
g[t_, k_, {a_, b_, c_, d_, e_, f_}] := 
     g[t, k, {a, b, c, d, e, f}] = 
        g[t - 1, k, {a, b, c, d, e, f}]*a +  h[t - 1, k - 1, {a, b, c, d, e, f}]*b;

h[0, k_, _] := 0;
h[t_, 0, {a_, b_, c_, d_, e_, f_}] := f;
h[t_, k_, {a_, b_, c_, d_, e_, f_}] := 
     h[t, k, {a, b, c, d, e, f}] = 
        h[t - 1, k, {a, b, c, d, e, f}]*c +  g[t - 1, k - 1, {a, b, c, d, e, f}]*d;

You call it, passing the parameters explicitly:

In[1317]:= g[20,10,{a,b,c,d,e,f}]//Timing
Out[1317]= {0.,253809.}

(I was using the original parameters). You can check that the memoized definitions remain in the global rule base, in this method. Next time you call a function with exact same parameters, it will fetch the memoized definition rather than recompute. Apart from the problems with this approach that I outlined above, you should also watch for the memory usage, since nothing gets cleared.

Answer 2

Memoization Using Auxiliary Symbol

The memoization technique exhibited in the question can be modified so that the definitions of g and h do not need to be re-established whenever the cache needs to be cleared. The idea is to store the memoized values on an auxiliary symbol instead of directly on g and h:

g[0,k_] = 0;
g[t_,0] = e;
g[t_,k_] := memo[g, t, k] /. _memo :> (memo[g, t, k] = g[t-1,k]*a+h[t-1,k-1]*b)

h[0,k_] = 0;
h[t_,0] = f;
h[t_,k_] := memo[h, t, k] /. _memo :> (memo[h, t, k] = h[t-1,k]*c+g[t-1,k-1]*d)

The definitions are essentially the same as the original memoized versions of g and h except that a new symbol, memo, has been introduced. With these definitions in place, the cache can be cleared simply using Clear@memo -- there is no need to redefine g and h anew. Better still, the cache can be localized by placing memo in Block, thus:

Block[{memo, a = 7, b = 9, c = 13, d = 0.002, e = 2, f = 1}
, Table[g[t, k], {t, 0, 100}, {k, 0, 100}]
]

The cache is discarded when the block is exited.

Memoization Using Advice

The original and revised memoization techniques required invasive changes within the function g and h. Sometimes, it is convenient to introduce memoization after the fact. One way to do this would be to use the technique of advising -- a kind of functional programming analog to subclassing in OO programming. A particular implementation of function advice appears regularly in the pages of StackOverflow. However, that technique is also invasive. Let's consider an alternate technique of adding advice to g and h without altering their global definitions.

The trick will be to temporarily redefine g and h within a Block. The redefinitions will first check the cache for the result and, failing that, call the original definitions from outside the block. Let's go back to the original definitions of g and h that are blissfully unaware of memoization:

g[0,k_]:=0
g[t_,0]:=e
g[t_,k_]:=g[t-1,k]*a+h[t-1,k-1]*b

h[0,k_]:=0
h[t_,0]:=f
h[t_,k_]:=h[t-1,k]*c+g[t-1,k-1]*d

The basic schema for this technique looks like this:

Module[{gg, hh}
, copyDownValues[g, gg]
; copyDownValues[h, hh]
; Block[{g, h}
  , m:g[a___] := m = gg[a]
  ; m:h[a___] := m = hh[a]
  ; (* ... do something with g and h ... *)
  ]
]

The temporary symbols gg and hh are introduced to hold the original definitions of g and h. Then g and h are locally rebound to new caching definitions that delegate to the original definitions as necessary. Here is definition of copyDownValues:

ClearAll@copyDownValues
copyDownValues[from_Symbol, to_Symbol] :=
  DownValues[to] =
    Replace[
      DownValues[from]
    , (Verbatim[HoldPattern][from[a___]] :> d_) :> (HoldPattern[to[a]] :> d)
    , {1}
    ]

In an effort to keep this post short(er), this "copy" function is concerned only with down-values. A general advice facility also needs to account for up-values, subvalues, symbol attributes and so on.

This general pattern is easy, if tedious, to automate. The following macro function memoize does this, presented with little comment:

ClearAll@memoize
SetAttributes[memoize, HoldRest]
memoize[symbols:{_Symbol..}, body_] :=
  Module[{pairs, copy, define, cdv, sd, s, m, a}
  , pairs = Rule[#, Unique[#, Temporary]]& /@ symbols
  ; copy = pairs /. (f_ -> t_) :> cdv[f, t]
  ; define = pairs /. (f_ -> t_) :> (m: f[a___]) ~sd~ (m ~s~ t[a])
  ; With[{ temps = pairs[[All, 2]]
         , setup1 = Sequence @@ copy
         , setup2 = Sequence @@ define }
    , Hold[Module[temps, setup1; Block[symbols, setup2; body]]] /.
        { cdv -> copyDownValues, s -> Set, sd -> SetDelayed }
    ] // ReleaseHold
  ]

After much ado, we are now in a position to externally impose memoization upon the non-caching versions of g and h:

memoize[{g, h}
, Block[{a = 7, b = 9, c = 13, d = .002, e = 2, f = 1}
  , Table[g[t, k], {t, 0, 100}, {k, 0, 100}]
  ]
]

Putting it all together, we can now create a responsive Manipulate block:

Manipulate[
  memoize[{g, h}
  , Table[g[t, k], {t, 0, tMax}, {k, 0, kMax}] //
      ListPlot3D[#, InterpolationOrder -> 0, PlotRange -> All, Mesh -> None] &
  ]
, {{tMax, 10}, 5, 25}
, {{kMax, 10}, 5, 25}
, {{a, 7}, 0, 20}
, {{b, 9}, 0, 20}
, {{c, 13}, 0, 20}
, {{d, 0.002}, 0, 20}
, {{e, 2}, 0, 20}
, {{f, 1}, 0, 20}
, LocalizeVariables -> False
, TrackedSymbols -> All
]

The LocalizeVariables and TrackedSymbols options are artifacts of the dependencies that g and h have on the global symbols a through f.