在处理递归函数时如何提高(mathematica)性能？

Question

背景.我想打印一张31 ^(1/2)的会聚表.我对表进行了以下递归定义.(交换31 ^(1/2)黄金比率,下表将包含斐波那契系列).

 cf := ContinuedFraction
 tf := TableForm
 p[-1] = 0; p[0] = 1; q[-1] = 1; q[0] = 0;
 a[k_] := cf[Sqrt[31], k][[k]]
 p[k_] := a[k]*p[k - 1] + p[k - 2]
 q[k_] := a[k]*q[k - 1] + q[k - 2]
 s[n_] := Timing[Table[{k, a[k], p[k], q[k]}, {k, 8, 8 n, 8}]] // tf

时间以指数方式快速增长.我不得不alt +.(中止)s [4].

问题:如何在处理递归函数时提高(mathematica)性能？

Answer 1

从快速(不彻底,承认它)看你的代码,它看起来像两个,p并q以两个先前的值递归地定义.这意味着要计算其n值p,~2^n需要进行评估(每一步都会使数字加倍).所以是的,复杂性是指数级的,无论使用Mathematica还是其他任何语言.

如果你坚持使用问题的递归公式(例如为了简单),那么减少性能损失的最简单方法是使用memoization,即做类似的事情

p[k_] := p[k] = a[k]*p[k - 1] + p[k - 2]

Clear[p]在重新定义之前不要忘记.

简而言之,memoization意味着让函数记住每个输入的计算结果,因此后续的评估更快.从两个先前的值(p_(n)和p_(n-1))计算两个值(p_(n + 1)和p_(n))可能更快但更复杂,那么复杂性将是线性的而不是指数.

_{我希望这有帮助.我现在没有Mathematica来测试.}

Answer 2

这是一个小小的进一步改进.由于这是二次无理,您还可以更直接地计算a [k]系数.

In[499]:= Clear[a, p, q, cf]
cf = ContinuedFraction[Sqrt[31]];
cf2len = Length[cf[[2]]];
a[1] = cf[[1]];
a[k_] := cf[[2, Mod[k - 1, cf2len, 1]]]
p[-1] = 0; p[0] = 1; q[-1] = 1; q[0] = 0;
p[k_] := p[k] = a[k]*p[k - 1] + p[k - 2]
q[k_] := q[k] = a[k]*q[k - 1] + q[k - 2]
s[n_] := Timing[Table[{k, a[k], p[k], q[k]}, {k, 8, 8 n, 8}];]

In[508]:= s[1000]

Out[508]= {0.12, Null}

In[509]:= Clear[a, p, q, cf]
cf := ContinuedFraction
p[-1] = 0; p[0] = 1; q[-1] = 1; q[0] = 0;
a[k_] := a[k] = cf[Sqrt[31], k][[k]]
p[k_] := p[k] = a[k]*p[k - 1] + p[k - 2]
q[k_] := q[k] = a[k]*q[k - 1] + q[k - 2]
s[n_] := Timing[Table[{k, a[k], p[k], q[k]}, {k, 8, 8 n, 8}];]

In[516]:= s[1000]

Out[516]= {6.08, Null}

你也可以以封闭的形式获得[k],虽然它不是非常漂亮.

In[586]:= Clear[a];
asoln[k_] = 
 FullSimplify[
  a[k] /. First[
    RSolve[Join[
      Table[a[k] == cf[[2, Mod[k - 1, cf2len, 1]]], {k, 
        cf2len}], {a[k] == a[k - 8]}], a[k], k]], Assumptions -> k > 0]

Out[587]= (1/(8*Sqrt[2]))*(4*(Cos[(k*Pi)/4] + Sin[(k*Pi)/4])*
        (-2*Sqrt[2] + (5 + 2*Sqrt[2])*Sin[(k*Pi)/2]) + 
      Sqrt[2]*(25 - 9*Cos[k*Pi] + 26*Sin[(k*Pi)/2] - 9*I*Sin[k*Pi]))

我不知道这是否可以用来直接解决p [k]和q [k].RSolve似乎无法做到这一点.

---编辑---

正如其他人所提到的,从头到尾构建列表可能更清晰.这是p [k]的处理,使用上面的memoization和NestList.

Clear[a, p, q, cf]
cf = ContinuedFraction[Sqrt[31]];
cf2len = Length[cf[[2]]];
a[1] = cf[[1]];
a[k_] := cf[[2, Mod[k - 1, cf2len, 1]]]

p[-1] = 0; p[0] = 1;
p[k_] := p[k] = a[k]*p[k - 1] + p[k - 2]
s[n_] := Timing[Table[p[k], {k, n}];]

In[10]:= s[100000]

Out[10]= {1.64, Null}

In[153]:= s2[n_] := Timing[ll = Module[{k = 0},
      NestList[(k++; {#[[2]], a[k]*#[[2]] + #[[1]]}) &, {0, 1}, 
       n]][[All, 2]];]

In[154]:= s2[100000]

Out[154]= {0.78, Null}

除了速度更快之外,第二种方法并没有保留大量的定义.并且你并不真正需要它们来生成更多元素,因为这个迭代可以使用最后一个元素中的一对来恢复(确保它们从0和1开始以模8开始).

我会提到可以获得p [k]的封闭形式.我发现将解决方案分解为8(即cf2len)片段并通过重复将它们链接起来很方便.幕后推理来自基本的生成函数操作.我对一个方程和一个初始条件进行了一些略微特殊的处理,以便巧妙地确定a [1]不是重复序列的一部分.

In[194]:= func = Array[f, cf2len];
args = Through[func[n]];
firsteqns = {f[2][n] == a[2]*f[1][n] + f[cf2len][n - 1], 
   f[1][n] == a[9]*f[cf2len][n - 1] + f[cf2len - 1][n - 1]};
resteqns = 
  Table[f[j][n] == a[j]*f[j - 1][n] + f[j - 2][n], {j, 3, cf2len}];
inits = {f[8][0] == 1, f[1][1] == 5};
eqns = Join[firsteqns, resteqns, inits];

In[200]:= 
soln = FullSimplify[args /. First[RSolve[eqns, args, n]], 
   Assumptions -> n > 0];

In[201]:= FullSimplify[Table[soln, {n, 1, 3}]]

Out[201]= {{5, 6, 11, 39, 206, 657, 863, 1520}, {16063, 17583, 33646, 
  118521, 626251, 1997274, 2623525, 4620799}, {48831515, 53452314, 
  102283829, 360303801, 1903802834, 6071712303, 7975515137, 
  14047227440}}

快速检查:

In[167]:= s2[16]; ll

Out[167]= {1, 5, 6, 11, 39, 206, 657, 863, 1520, 16063, 17583, 33646, \
118521, 626251, 1997274, 2623525, 4620799}

我们现在可以从中定义一个函数.

In[165]:= 
p2[k_Integer] := soln[[Mod[k, cf2len, 1]]] /. n -> Ceiling[k/cf2len]

In[166]:= Simplify[p2[4]]

Out[166]= 39

我并不认为这是特别有用的,只是想看看我是否能真正得到一些工作.

---结束编辑---

Daniel Lichtblau