向量化一个 for 循环，其中输入取决于输出

Question

我有一个关于循环依赖向量化的复杂问题，我想从你那里得到一些帮助。

让我们X1 是一个长度为n1 的向量，X2 是一个长度为n2，F1 的向量是一个N1xn1 矩阵，F2 是一个N2 xn2 矩阵，Q 是N1xN2 矩阵，符号p... 是索引向量。 ntrapz 是梯形数值积分的函数。我想计算矩阵Q如下：

for i1=1:N1
    F1_13tmp=F1(i1, p1_13)';                       % ' 
    F1_13=F1_13tmp(:,ones(n2,1));
    for i2=1:N2
        F2_13 = F2(i2, p2_13);
        Q_13_13 = Q(p1_13, p2_13);
        Q(i1,i2) = Q(i1,i2) + 
              ntrapz(X2(p2_13), ntrapz(X1(p1_13)', Q_13_13.*F1_13).*F2_13);
    end
end

问题是更新Q(i1,i2) 会更改Q_13_13 = Q(p1_13, p2_13) 的值以进行下一次迭代。我想知道我们是否可以矢量化这样的 for 循环。如果没有，有什么加快代码速度的想法吗？

提前感谢您的帮助。

Answer 1

矢量化通常是并行操作的方式。但是，在您的情况下，您似乎有顺序操作，并且这些操作通常不适合矢量化。

我不熟悉您使用的函数，但我能想到的唯一可能加快计算速度的方法是以非依赖形式编写它。请注意，这更像是一个数学练习而不是编程练习，这是一个简单的例子：

假设公式等于：

Y(1)=0.5;
Y(t)=0.3*Y(t-1)+epsilon(t);
% Let as assume epsilon is also known
epsilon = rand(9,1)

计算 Y(10) 的简单方法是执行 10 个计算步骤。

for t=2:10
    Y(t)=0.3*Y(t-1)+epsilon(t);
end

计算 Y(10) 的最快方法是这样做（不确定它是否正确，但应该接近）

Y(10)=Y(1)*0.3^9 + sum(epsilon.^(1:9))

因此得出结论：“矢量化”循环的唯一方法是，您是否可以显着帮助计算机并基本上让它进行不同的计算。如果这是不可能的，您当然可以尝试通过将其制作为 mex 文件来挤出一点额外的速度。