matlab中循环的向量化答案

【问题标题】：Vectorization of a loop in matlabmatlab中循环的向量化
【发布时间】：2020-05-22 20:41:44
【问题描述】：

我有一个循环

for i = 2:K
    T(K,i) = ((4^(i-1))*T(K,i-1)-T(K-1,i-1))/(4^(i-1)-1); 
end

其中T 是一个二维矩阵（给定行中的第一个元素，并且上面行中的所有元素都已经存在），K 是一个标量。我试图将这个循环矢量化以使其更快：

i = 2:K;
T(K,i) = ((4.^(i-1)).*T(K,i-1)-T(K-1,i-1))./(4.^(i-1)-1);

它可以编译，但会产生不正确的结果。你能告诉我我哪里出错了吗？

@编辑：我已经写了这个，但结果仍然是错误的

     i = 2:K;
     i2 = 1:(K-1);
     temp1 = T(K,i2)
     temp2 = T(K-1,i2)
     T(K,i) = ((4.^(i2)).*temp1-temp2)./(4.^(i2)-1);

【问题讨论】：

您的循环会在迭代时更改 T 的值。
@bla 所以我应该做一些临时变量来存储原始值？
@bla 我已经编辑了问题

标签： matlab optimization vectorization

【解决方案1】：

首先，让我们重新索引您的循环（使用更少的 i-1 表达式）：

for i=1:K-1
  T(K,i+1) = ( 4^i*T(K,i) - T(K-1,i) ) / (4^i-1); 
end

然后（我现在将省略循环），我们可以分解4^i/(4^i-1)：

  T(K,i+1) = ( T(K,i) - T(K-1,i)/4^i ) * (4^i/(4^i-1));

让我们调用a(i) = (4^i/(4^i-1))、b(i) = - T(K-1,i)/4^i，然后扩展我们得到的第一个术语：

T(K,1) = T(K,1)
T(K,2) = T(K,1)*a(1)           + b(1)*a(1)
T(K,3) = T(K,1)*a(1)*a(2)      + b(1)*a(1)*a(2)      + b(2)*a(2)
T(K,4) = T(K,1)*a(1)*a(2)*a(3) + b(1)*a(1)*a(2)*a(3) + b(2)*a(2)*a(3) + b(3)*a(3)

然后c = [1, a(1), a(1)*a(2), ...] = [1, cumprod(a)]

T(K,i) = (T(K,1) + (b(1)/c(1) + b(2)/c(2) + ... + b(i-1)/c(i-1) ) *  c(i)

所以用d = b ./ c，e = cumsum(d)，总结所有计算：

i=1:K-1;
a = 4.^i./(4.^i-1);
b = -T(K-1,1:end-1) ./ 4.^i;
c = [1, cumprod(a)];
d = b ./ c(1:end-1);
e = cumsum(d);
T(K,2:K) = (T(K,1) + e) .* c(2:end);

为了进一步优化这一点，请注意4^14/(4^14 - 1) 等于 1，当使用双精度计算时，所以实际上 T(K,14:K) 可以大幅优化 - 即，您实际上只需要计算 a、c、 1./c 直到索引 13。（我将把它留作练习）。

【讨论】：