CODEWITHSUNDEEP
matlaboptimizationvectorization
MartinYakuza
MartinYakuza

Reputation: 135

Vectorization of a loop in matlab

I have a loop

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

where T is a two dimensional matrix (first element in given row, and all elements in rows above are already there) and K is a scalar. I have tried to vectorize this loop to make it faster like this:

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

it compiles, but it yields improper results. Can you tell me where am I making mistake?

@EDIT: I have written this, but still the result is wrong

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

Upvotes: 0

Views: 47

Answers (1)

chtz
chtz

Reputation: 18807

First, let's re-index your loop (having fewer i-1 expressions):

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

Then (I'll leave out the loop for now), we can factor out 4^i/(4^i-1):

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

Let's call a(i) = (4^i/(4^i-1)), b(i) = - T(K-1,i)/4^i, then expanding the first terms we get:

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)

Then with 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)

So with d = b ./ c, e = cumsum(d), summarizing all calculations: 

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

To further optimize this, note that 4^14/(4^14 - 1) equals 1, when calculated with double-precision, so actually T(K,14:K) could be optimized drastically -- i.e., you actually just need to calculate a, c, 1./c up to index 13. (I'll leave that as an exercise).

Upvotes: 2

Related Questions