Is there a way to vectorize these loops for speed?

I need to construct symmetric matrix M as shown in the following code:

theta = 0.5;
lambda =  0.1; % usually a number between zero and one
N = 1000   % or larger  number
M = zeros(N,N); % allocate initial value 

for t = 1:N
    firstTerm = sum(theta.^(4.*(t-1:-1:0)));
    secondTerm = 0;
    for k = 1:t-1
        secondTerm = secondTerm + sum(theta.^(2*((2*t)-(2*k)-1:-1:t-k)));
    end
    M(t,t) =  3*lambda^2*firstTerm + 6*lambda^2*secondTerm;
    for s=t+1:N
        thirdTerm = 0;
        for r = t+1:s
            thirdTerm = thirdTerm + sum(theta.^(2*(t+s-r-1:-1:s-r)));
        end
        M(t,s) = 3*lambda^2*theta^(2*(s-t))*firstTerm + 6*lambda^2*theta^(2*(s-t))*secondTerm +  lambda^2*thirdTerm;
        M(s,t) = M(t,s);
    end
end

The code takes very long when N is large; is there any better way to compute M? I was hoping that it might be possible to vectorize the computations?

Observe that the diagonal elements are easy to compute, but the off diagonals required a third term in the sum. However, they can be constructed given the diagonal as shown in the code.

Answer 1

[Matt J]

Summations like

sum(theta.^(2*((2*t)-(2*k)-1:-1:t-k)));

can be replaced with closed-form formulas for geometric series, see for example,

https://en.wikipedia.org/wiki/Geometric_progression#Geometric_series

This will save you the overhead of both constructing the series itself and the sum() command.

[Jan]

@Matt J: Geometric series are not trivial to apply, if rounding must be considered. For terms like 0.5^0 to 0.5^1000 the order of operations matter for a sum. Calculating the sum numerically might cause a huge rounding error compared to the application of the formula for a geometric series. So maybe your suggestion is good, or the only reliable method. +1

[Mohamed Abdalmoaty]

Thanks Matt! I was able to reduce all the sums as well as replacing the inner loop in r with closed form expressions. I think the code is much better now! +1

[Matt J]

Thanks, guys.

Answer 2

[Jan]

You compute expensive powers of theta repeatedly. Better do this once and store the values in a vector. What is the largest power of theta? Let's assume it is 4*N.

theta = 0.5;
N = 1000   % or larger  number
M = zeros(N,N); % allocate initial value 

lambda2 = lambda ^ 2;
tpow = theta .^ (0:4*N);
% Or cheaper:
% tpow = cumprod([1, repmat(theta, 1, N-1)]);
for t = 1:N
  firstTerm = sum(tpow(1 + 4.*(t-1:-1:0)));
  secondTerm = 0;
  for k = 1:t-1
      secondTerm = secondTerm + sum(tpow(1 + 2*((2*t)-(2*k)-1:-1:t-k)));
  end
  M(t,t) =  lambda2*(3*firstTerm + 6*secondTerm);
  for s=t+1:N
      thirdTerm = 0;
      for r = t+1:s
          thirdTerm = thirdTerm + sum(tpow(1 + 2*(t+s-r-1:-1:s-r)));
      end
      M(t,s) = lambda2*(3*tpow(1 + 2*(s-t))*firstTerm + ...
                        6*tpow(1 + 2*(s-t))*secondTerm + thirdTerm);
      M(s,t) = M(t,s);
  end
end

Unfortunately I cannot test this. Please provide the value of lambda.

0.5^4000 is smaller than realmin. Therefore the corresponding terms in the sum are 0 and could be omitted.

I've removed some multiplications by lambda^2 and many square operations of lambda by storing it once before the loop. The idea is to avoid all repeated calculations.

[Mohamed Abdalmoaty]

Thanks for the answer, you can try to use any positive value for lambda. II updated my question with a value. However it seems to me that the computational time, while couple of times shorter than the code in the questions, is still long due to the iterative loops. It is of course necessary to avoid all repeated computations to reduce the time, however this is not the bottle neck and the operations are still in an iterated loops (computing element by element).

Can you let me know if you tested it? When I try, your code doesn't give the same values!

[Jan]

@Mohamed: Please post the timings for the original code and my suggestion. Then show us the size of the difference of the methods. Maybe this is caused by rounding only:

t1 = lambda2*(3*tpow(1 + 2*(s-t))*firstTerm + ...
              6*tpow(1 + 2*(s-t))*secondTerm + thirdTerm)
t2 = lambda2*3*tpow(1 + 2*(s-t))*firstTerm + ...
     lambda2*6*tpow(1 + 2*(s-t))*secondTerm + 
     lambda2*thirdTerm

This can slightly differ, but this is not a bug, but an effect of rounding. But maybe I made a mistake during editing. You can find it by your own using the debugger: Compare the intermediate terms.

[Mohamed Abdalmoaty]

When I compared the results yesterday, there was a big difference which I guess cannot be a rounding error! But, I'll have another look when I can. Thanks