Speed up iteration of pair of fully coupled recursion equations

Question

I am trying to do several hundred iterations of these fully coupled recursion equations (but will use twenty iterations for the example).

ClearAll["Global`*"]

x1[t_] := x1[t] = (1 - m) x1[t - 1] + m x2[t - 1]
x2[t_] := x2[t] = (1 - m) x2[t - 1] + m x1[t - 1]

x1[0] := 1
x2[0] := 0

ListPlot[Table[{{t, x1[t]}, {t, x2[t]}} /. m -> 0.01, {t, 0, 
20}]] // AbsoluteTiming

Memoization sped it up from 27.933541 seconds to 7.681030 seconds but I think it takes so long because it is still recursively calculating all the values from the other equation. I have read several posts about partially coupled recursion equations (e.g. Solve pair of recurrence relations, How do I use RSolve to solve a system of recurrence relations?) but ultimately, I won't be able to use RSolve because my equations will be too complicated and I cannot find a solution to speed up iterations of fully coupled recursion equations.

Is there a way to link x1[t] and x2[t] together for memoization? Or is a For loop the way to go? Thanks!

score 8 · Accepted Answer · 2016-09-16T07:25:43.867

I think that your code already memoized properly.

Adding a Print on the most right-hand side of your definitions:

ClearAll[x1, x2]

x1[t_] := x1[t] = (
    Print["Computation of x1 for t = ", t]; 
    (1 - m) x1[t - 1] + m x2[t - 1]
);

x2[t_] := x2[t] = (
    Print["Computation of x2 for t = ", t]; 
    (1 - m) x2[t - 1] + m x1[t - 1]
);

x1[0] := 1
x2[0] := 0

and evaluating the table with an iterator that goes up to 3, we get

Table[{{t, x1[t]}, {t, x2[t]}} /. m -> 0.01, {t, 0, 3}]

During evaluation of In[..]:= Computation of x1 for t = 1
During evaluation of In[..]:= Computation of x2 for t = 1
During evaluation of In[..]:= Computation of x1 for t = 2
During evaluation of In[..]:= Computation of x2 for t = 2
During evaluation of In[..]:= Computation of x1 for t = 3
During evaluation of In[..]:= Computation of x2 for t = 3

(* {
    {{0, 1}, {0, 0}}, {{1, 0.99}, {1, 0.01}}, 
    {{2, 0.9802}, {2, 0.0198}}, {{3, 0.970596}, {3, 0.029404}}
   } 
*)

so the computation is done once for each of the x's at each step, which shows that values at previous steps are indeed stored and used.

A possibility to speed up the computation would be to avoid using ReplaceAll and define your functions with m as argument:

ClearAll[x1, x2]

x1[t_, m_] := x1[t, m] = (1 - m) x1[t - 1, m] + m x2[t - 1, m]
x2[t_, m_] := x2[t, m] = (1 - m) x2[t - 1, m] + m x1[t - 1, m]

x1[0, m_] := 1
x2[0, m_] := 0

ListPlot[Table[{{t, x1[t, 0.01]}, {t, x2[t, 0.01]}}, {t, 0, 20}]] // AbsoluteTiming

On my computer, the timing goes from about 7 seconds, to less than 0.05 second.

Thanks! The following also speeds it up on a similar order of magnitude. Subbing in a series of parameters in the recursion definition, rather than taking a parameter value as a function argument, will be helpful to me (and possibly other users) down the road as I get into more complicated equations. parameters = {m->0.01} but can be parameters = {m->0.01, s->0.1, alpha ->0.5}
`x1[t_] := x1[t] = (1 - m) x1[t - 1] + m x2[t - 1]/.parameters x2[t_] := x2[t] = (1 - m) x2[t - 1] + m x1[t - 1]/.parameters x1[0] := 1 x2[0] := 0` — biologyUser, Sep 15 '16 at 22:00
I agree, this code avoids one to rewrite all the definitions. It is better. In case you have many parameters, and more complicated expressions (with more leaves), it is possible that the argument-approach performs better, but I doubt there will be an important difference in efficiency with your new code. (This could make another series of tests! :-)) — , Sep 15 '16 at 22:09
Faster, shorter, and easier to read: DiscretePlot[{x1[t, 0.01], x2[t, 0.01]}, {t, 0, 20}, Filling -> None] — Karsten7, Sep 16 '16 at 08:03

score 5 · Answer 2 · answered Sep 16 '16 at 14:04

RSolve will provide an exact closed-form solution that is faster

Clear[soln]

soln[m_, t_] = RSolve[{
    x1[t] == (1 - m) x1[t - 1] + m x2[t - 1],
    x2[t] == (1 - m) x2[t - 1] + m x1[t - 1],
    x1[0] == 1, x2[0] == 0}, {x1[t], x2[t]}, t][[1]]

Plot[
 Evaluate[
  {x1[t], x2[t]} /. soln[1/100, t]],
 {t, 0, 10},
 PlotLegends -> "Expressions"]

Chris K · Answer 3 · 2016-09-16T13:11:56.873

4

Could you not use RecurrenceTable?

m = 0.01;
res = RecurrenceTable[{
  x1[t] == (1 - m) x1[t - 1] + m x2[t - 1],
  x2[t] == (1 - m) x2[t - 1] + m x1[t - 1],
  x1[0] == 1, x2[0] == 0}, {x1, x2}, {t, 20}]

is superfast. Here's an ugly way to recreate your plot:

ListPlot[{
  Transpose[{Table[t, {t, 0, 20}], Transpose[res][[1]]}], 
  Transpose[{Table[t, {t, 0, 20}], Transpose[res][[2]]}]
}]

Mathematica graphics

edited Sep 16 '16 at 13:11

answered Sep 16 '16 at 02:56

Chris K

20,207
3
39
74

Speed up iteration of pair of fully coupled recursion equations

3 Answers3