Skip to main content

memoization - 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


t[a,b,c]

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


t[a,b,c]

it stores


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


etcetera.


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!


EDIT


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:


memoize[
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:


ClearAll@memoize
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}]},

{lhs,rhs}/.{Hold[x_],Hold[y_]}:>SetDelayed[x,y]
]]]]

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:


ClearAll@memoize

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] ;
lhs/.x_Pattern:>x[[1]]
]]]
} /. {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..



Answer



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.


Solution


Here is the code:


ClearAll[t];    
t[a_]:= c[a]
t[a_,b_]:= d[a,b]
t[a__]:=

With[{vars=Table[Unique[],{Length[{a}]}]},
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]

]
];

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

Global`t
t[a_]:=c[a]

t[a_,b_]:=d[a,b]

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

t[a__]:=With[{vars=Table[Unique[],{Length[{a}]}]},

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:


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

t[a]
]

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


t[m,n,k]

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


t[a__]:=
With[{vars=Table[Unique[],{Length[{a}]}]},
With[{pts=(Pattern[#1,_]&)/@vars},
pts/.
{x__}:>Set@@Hold[t[x],A * t @@ Most[Most[vars]]+B * t @@ Rest[vars]];
t[a]
]

];

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__].


Comments

Popular posts from this blog

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...

plotting - How to draw lines between specified dots on ListPlot?

I would like to create a plot where I have unconnected dots and some connected. So far, I have figured out how to draw the dots. My code is the following: ListPlot[{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4,13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full] I have thought using ListLinePlot command, but I don't know how to specify to the command to draw only selected lines between the dots. Do have any suggestions/hints on how to do that? Thank you. Answer One possibility would be to use Epilog with Line : ListPlot[ {{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4, 13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full, Epilog -> { Line[ ...