Vectorization of matrix multiplication with scalar (Scalar is value of other matrix at index i) and expm() operation

Hi,

I can't seem wrap my head on vectorizing a for loop in my MATLAB code. Basically, I have this code:

% Let that:
%Tp = scalar
%N, = scalar (Say 1000)
%Ac, = 4x4 matrix
%pre = 1x4 matrix
%post = 4x2 matrix
%wy1 = N+1x1 matrix (so it would be 1001*1)
%wy2 = N+1x1 matrix (so it would be 1001*1)

% preallocate
delta_ksi=Tp/N;
AcT =Ac';
sum_matrix=zeros(4,1);
Fl=zeros(4,1);

% calculate the sum
for i=1:N
 Fl=delta_ksi*expm(AcT*delta_ksi*i)*post*[wy1(1,i); wy2(1,i);];
 sum_matrix= sum_matrix+Fl;
end

   %value I need
  delta_f_des_ff= pre*sum_matrix;

What I have in mind is to construct a 3D matrix Fl_3D (4 x 1x 1000) and then do array multiplication with i = 1:1000, but I kept getting incompatible dimension error when multiplying with [wy1(1,i); wy2(1,i)] which also use the index i.

Any clue on what is the best approach to do this? Is vectorization still possible? Thanks!

Background:

1. I am trying to troubleshoot a bottleneck code in a Simulink project. Code Profiling shows that the above function hogging most of the runtime.
2. I am hoping vectorizing could solve the performance issue. I also just found out the bottleneck comes from expm() operation.

Any suggestion are welcomed!

Answer 1

[Robert]

I don't think this problem is well suited for vectorization since you are already doing matrix multiplication at each iteration of the loop; however, you can get a very nice speedup by taking the expm call outside of the loop.

Using the property e^(aX) = (e^X)^a you can refactor your code to calculate

E = expm(AcT*delta_ksi);

before the loop and use

E^i

in place of

expm(AcT*delta_ksi*i)

in your loop.

Since you are calculating consecutive exponents, why not do the multiplication yourself for an extra speedup?

E0 = expm(AcT*delta_ksi);
E1 = eye(4);
% calculate the sum
for ii = 1:N
    E1 = E0*E1;
    Fl = E1*post*[wy1(1,ii); wy2(1,ii);];
    sum_matrix = sum_matrix+Fl;
end

Bonus suggestion, ii makes a good loop variable name when you want something short because it doesn't conflict with the complex variable i.

[Arya]

The first part is pretty clear and it gives some improvement to the execution speed. Thank you very much!

but I don't really get the second part: Why are you multiplying eye(4) with expm(AcT*delta_ksi)?

Initially, we have:

expm(AcT*delta_ksi*i)

which is refactored into:

 expm(AcT*delta_ksi)^i

It doesn't add up to me how

expm(AcT*delta_ksi)^i;

becomes equivalent to:

 eye(4)*expm(AcT*delta_ksi)

.. but your code is outputting the desired results. Something is wrong with my understanding.

[Robert]

I should have been more clear about that. I am using eye(4) as the zero-th value of E1. That way, when I multiply E1 (the Identity Matrix) by E0 for the first time, I get the value of E0 back.

In this way, E1 always equals E0^ii.

ii             E1
__            _____
 0*            eye(4)
 1          E0*eye(4)
 2       E0*E0*eye(4)
 3    E0*E0*E0*eye(4)
 4 E0*E0*E0*E0*eye(4)
 ...
 n             E0^ii

* before the loop

Now that I think about it, using 1 would probably be more clear and even easier. We just need some initial value for E1 that won't affect the results inside the loop.

So now it looks like:

E0 = expm(AcT*delta_ksi);
E1 = 1; % a dummy initial value
% calculate the sum
for ii = 1:N
    E1 = E0*E1; % E1 = E0^ii
    Fl = E1*post*[wy1(1,ii); wy2(1,ii);];
    sum_matrix = sum_matrix+Fl;
end