2

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

让我们X1成为一个长度为 的向量n1，X2成为一个长度为 的向量n2，F1成为一个N1xn1矩阵，F2成为一个N2xn2矩阵，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 循环。如果没有，有什么想法可以加快代码速度吗？

预先感谢您的帮助。

4

1 回答 1

0

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

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

假设公式等于：

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 文件来尝试挤出一点额外的速度。

于 2012-12-05T10:31:57.017 回答