Picard迭代中矩阵列表的优化计算

Question

目前，我正在使用Mathematica代码来执行Picard迭代。代码本身运行良好，但我正在努力使其更高效。我已经取得了一些成功，但我正在寻求建议。也许不可能再加快速度了，但我的想法已经用完了，我希望那些比我更有编程/数学经验的人能够提出一些建议。我只是张贴迭代本身，但可以提供更多的信息，视需要。

下面的代码根据请求被编辑成完全可执行的

另外，我将它从一段时间改为了Do循环，以使测试更容易，因为不需要收敛。

Clear["Global`*"]

ngrid = 2048;
delr = 4/100;
delk = \[Pi]/delr/ngrid;
rvalues = Table[(i - 1/2) delr, {i, 1, ngrid}];
kvalues = Table[(i - 1/2) delk, {i, 1, ngrid}];
wa[x_] := (19 + .5 x) Exp[-.7 x] + 1
wb[x_] := (19 + .1 x) Exp[-.2 x] + 1
wd = SetPrecision[
   Table[{{wa[(i - 1/2) delk], 0}, {0, wb[(i - 1/2) delk]}}, {i, 1, 
     ngrid}], 26];
sigmaAA = 1;
hcloseAA = {};
i = 1;
While[(i - 1/2)*delr < sigmaAA, hcloseAA = Append[hcloseAA, -1]; i++]
hcloselenAA = Length[hcloseAA];
hcloseAB = hcloseAA;
hcloselenAB = hcloselenAA;
hcloseBB = hcloseAA;
hcloselenBB = hcloselenAA;
ccloseAA = {};
i = ngrid;
While[(i - 1/2)*delr >= sigmaAA, ccloseAA = Append[ccloseAA, 0]; i--]
ccloselenAA = Length[ccloseAA];
ccloselenAA = Length[ccloseAA];
ccloseAB = ccloseAA;
ccloselenAB = ccloselenAA;
ccloseBB = ccloseAA;
ccloselenBB = ccloselenAA;
na = 20;
nb = 20;
pa = 27/(1000 \[Pi]);
pb = 27/(1000 \[Pi]);
p = {{na pa, 0}, {0, nb pb}};
id = {{1, 0}, {0, 1}};
AFD = 1;
AFDList = {};
timelist = {};
gammainitial = Table[{{0, 0}, {0, 0}}, {ngrid}];
gammafirst = gammainitial;
step = 1;
tol = 10^-7;
old = 95/100;
new = 1 - old;

Do[
 t = AbsoluteTime[];
 extractgAA = Table[Extract[gammafirst, {i, 1, 1}], {i, hcloselenAA}];
 extractgBB = Table[Extract[gammafirst, {i, 2, 2}], {i, hcloselenBB}];
 extractgAB = Table[Extract[gammafirst, {i, 1, 2}], {i, hcloselenAB}];
 csolutionAA = (Join[hcloseAA - extractgAA, ccloseAA]) rvalues;
 csolutionBB = (Join[hcloseBB - extractgBB, ccloseBB]) rvalues;
 csolutionAB = (Join[hcloseAB - extractgAB, ccloseAB]) rvalues;
 chatAA = FourierDST[SetPrecision[csolutionAA, 32], 4];
 chatBB = FourierDST[SetPrecision[csolutionBB, 32], 4];
 chatAB = FourierDST[SetPrecision[csolutionAB, 32], 4];
 chatmatrix = 
  2 \[Pi] delr Sqrt[2*ngrid]*
   Transpose[{Transpose[{chatAA, chatAB}], 
      Transpose[{chatAB, chatBB}]}]/kvalues;
 gammahat = 
  Table[(wd[[i]].chatmatrix[[i]].(Inverse[
         id - p.wd[[i]].chatmatrix[[i]]]).wd[[i]] - 
      chatmatrix[[i]]) kvalues[[i]], {i, ngrid}];
 gammaAA = 
  FourierDST[SetPrecision[Table[gammahat[[i, 1, 1]], {i, ngrid}], 32],
    4];
 gammaBB = 
  FourierDST[SetPrecision[Table[gammahat[[i, 2, 2]], {i, ngrid}], 32],
    4];
 gammaAB = 
  FourierDST[SetPrecision[Table[gammahat[[i, 1, 2]], {i, ngrid}], 32],
    4];
 gammasecond = 
  Transpose[{Transpose[{gammaAA, gammaAB}], 
     Transpose[{gammaAB, gammaBB}]}]/(rvalues 2 \[Pi] delr Sqrt[
      2*ngrid]);
 AFD = Sqrt[
    1/ngrid Sum[((gammafirst[[i, 1, 1]] - 
           gammasecond[[i, 1, 1]])/(gammafirst[[i, 1, 1]] + 
           gammasecond[[i, 1, 1]]))^2 + ((gammafirst[[i, 2, 2]] - 
           gammasecond[[i, 2, 2]])/(gammafirst[[i, 2, 2]] + 
           gammasecond[[i, 2, 2]]))^2 + ((gammafirst[[i, 1, 2]] - 
           gammasecond[[i, 1, 2]])/(gammafirst[[i, 1, 2]] + 
           gammasecond[[i, 1, 2]]))^2 + ((gammafirst[[i, 2, 1]] - 
           gammasecond[[i, 2, 1]])/(gammafirst[[i, 2, 1]] + 
           gammasecond[[i, 2, 1]]))^2, {i, 1, ngrid}]];
 gammafirst = old gammafirst + new gammasecond;
 time2 = AbsoluteTime[] - t;
 timelist = Append[timelist, time2], {1}]
Print["Mean time per calculation = ", Mean[timelist]]
Print["STD time per calculation = ", StandardDeviation[timelist]]

只是一些关于事情的笔记

ngrid、delr、delk、rvalue、kvalue只是用来使问题离散的值。典型的是

ngrid = 2048;
delr = 4/100;
delk = \[Pi]/delr/ngrid;
rvalues = Table[(i - 1/2) delr, {i, 1, ngrid}];
kvalues = Table[(i - 1/2) delk, {i, 1, ngrid}];

所有使用的矩阵都是2x2，具有相同的对角线。

恒等矩阵和P矩阵(实际上是密度)是

p = {{na pa, 0}, {0, nb pb}};
id = {{1, 0}, {0, 1}};

我所确定的计算中的主要慢点是FourierDST计算(正向和反变换占计算时间的近40% )，gammahat计算占40%，剩余时间以AFD计算为主。在我的i7处理器上，每个周期的平均计算时间是1.52秒。我的希望是让它在一秒钟之内，但这可能是不可能的。我希望引入一些并行计算，这是用两个ParallelTable命令以及使用ParallelSubmit WaitAll进行的尝试。然而，我发现并行计算中的任何加速都被从主内核到其他内核的通信时间所抵消(至少这是我的假设，因为对新数据的计算所需的时间是重新计算现有数据的两倍。我认为这意味着经济增长放缓是因为发布了新的名单)，我和DistributDefinitions以及SetSharedVariable都玩过，但却无法做到这一点。

我想知道的一件事是，使用Table进行离散计算是否是最好的方法？

我也曾想过，我可能会改写它，以便能够编译它，但我的理解是，只有当你处理机器精度时，我需要更高的精度才能收敛。

谢谢您的任何建议。

Answer 1

我将等待acl建议的代码，但从顶部看，我怀疑这个构造：

Table[Extract[gammafirst, {i, 1, 1}], {i, hcloselenAA}]

可以编写，并将执行得更快，如：

gammafirst[[hcloselenAA, 1, 1]]

但我不得不猜出你数据的形状。

Answer 2

在几行中使用：

FourierDST[SetPrecision[Table[gammahat[[i, 1, 1]], {i, ngrid}], 32], 4];

您可以删除Table

FourierDST[SetPrecision[gammahat[[All, 1, 1]], 32], 4];

如果你真的，真的需要这个SetPrecision，你就不能在计算游戏帽的时候马上做吗？

AFAI可以看出，在计算gammahat时使用的所有数字都是精确的。这可能是故意的，但速度很慢。你可以考虑用近似数字代替。

编辑

随着完整的代码在你的最新编辑，只是添加一个//N到你的第二和第三行削减时间，至少在一半，而不降低数字的准确性很多。如果我比较res={gammafirst，gammasecond，AFD}中的所有数字，那么原始的和添加//N的是res1 - res2 // Flatten // Total ==> 1.88267*10^-13。

删除所有的SetPrecision内容将代码的速度提高了7倍，并且结果似乎具有类似的准确性。