Delooping variables that are loop-dependent

I have some code that I am looking to considerably speed up.

It is related to a previous question I asked (Why is sum considerably slower that adding each individual element together, when using large loops? - MATLAB Answers - MATLAB Central (mathworks.com).

I'm using A to 'normalise' B which then alters A through some propagator C. D then converts the vector A into a single number that is the signal. Everything outside of the loop can be treated as a constant.

Unlike the previous question I asked, B0 and C are now no longer constants and are dependent on n, so I am unable to pull out the matrix power from the loop.

The expm line is considerably the slowest, followed by the mldivide line.

Any help would be much appreciated.

% Constants w.rt. n
A = rand(16,1);
idx = [1 6 11 16];
D = rand(1,16);

N = 1e5;

signal = zeros(N,1);

tic
for n = 1:N

  % Constants dependent on n
  B0_ = rand(16,1);
  C_ = rand(16,16);

  C = expm(C_*1e-6);
  B0 = -C\B0_;

  Aidx = A(idx);
  Asum = sum(Aidx);
  B = B0*Asum;
  A = B + C*(A-B);
  signal(n) = D*A;

end
toc

Elapsed time is 3.467775 seconds.

[Torsten]

I wonder where you see any potential for code optimization.

Answer 1

[Vatsal]

Hi @Jack,

I understand that you are looking to optimize the code. In the given code variable “A” normalise “B” which then alters “A” through some propagator “C”. Afterwards, "D" converts the vector "A" into a single number representing the signal. To improve the code's performance, I have implemented an alternative approach for matrix exponentiation. Instead of using the "expm" function, I utilized the Taylor series method to calculate the matrix exponentiation. This modification has yielded better results compared to using the "expm" function

I have also included the modified code below, which incorporates the Taylor series method for matrix exponentiation:

A = rand(16,1);
idx = [1 6 11 16];
D = rand(1,16);
N = 1e5;
signal = zeros(N,1);
tic
for n = 1:N
  % Constants dependent on n
  B0_ = rand(16,1);
  C_ = rand(16,16);
  X=C_*1e-6;
  E = zeros(size(X));
  F = eye(size(X));
  k = 1;
  while norm(E+F-E,1) > 0
   E = E + F;
   F = X*F/k;
   k = k+1;
  end
  C = E;
  B0 = -C\B0_;
  Aidx = A(idx);
  Asum = sum(Aidx);
  B = B0*Asum;
  A = B + C*(A-B);
  signal(n) = D*A;
end
toc

For more information and different methods for matrix exponentiation, you can refer to the following link:

https://in.mathworks.com/help/matlab/math/matrix-exponentials.html?searchHighlight=matrix%20exponential&s_tid=srchtitle_support_results_1_matrix%2520exponential

I hope this helps!