matlab中循环的矢量化

Question

我有一个循环

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);

Answer 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。（我将把它留作练习）。

matlab中循环的矢量化

Vectorization of a loop in matlab

optimization

matlab

vectorization