0

Suppose a is an 1 x n vector, b is an n x 1 vector, A is an n x n matrix and xx is an m x 1 vector. I want to compute fxx which can be computed as follows:

for i = 1:m
    fxx(i)=a*expm(xx(i)*A)*b;
end

this of course works but it feels like this is a very slow method, is there a better built-in way for dealing with this type of for loop?

Adriaan
  • 17,741
  • 7
  • 42
  • 75
Darkwizie
  • 147
  • 4

1 Answers1

2

The problem is that expm is only implemented for 2D matrices, so you need to call that on 2D slices. The only little bit of possible speed-up is to precalculate xx*A to a 3D matrix:

n=5;
m=10;
a = rand(1,n);
b = rand(n,1);
A = rand(n,n);
xx = rand(m,1);

tmp = bsxfun(@times,permute(xx,[3 2 1]),A); % Pre-calculate the matrix multiplication

fxx = zeros(m,1);
for ii = 1:m
    fxx(ii) = a*expm(tmp(:,:,m))*b;
end

although I doubt it'll save you a lot of time this way. The main speed-up compared to your code, would be to preallocate fxx, as I have done, since that saves MATLAB from allocating and deallocating your array each loop iteration.

In your comments you state that m is generally very large. If that is the case, parfor might help you if you have the parallel computation toolbox. In that case, pre-calculating the 3D matrix might not work, since it'll be very large in RAM, so simply use:

fxx=zeros(m,1);
parfor ii = 1:m
    fxx(ii) = a*expm(xx(ii)*A)*b;
end

For more information on how parfor can help you, and how it can be even faster when slicing your variables properly, see the following question: Saving time and memory using parfor?

Adriaan
  • 17,741
  • 7
  • 42
  • 75
  • I see, the problem is that `m` is generally very large and I have been told to try and prevent using large for loops in matlab, but if there is no better way, then so be it. – Darkwizie Sep 13 '18 at 08:49
  • @Darkwizie that used to be true in the past, but nowadays (since the advent of R2016b) the JIT and general engine of MATLAB have been improved so far that loops are quite fast already. What you might try, if you have the parallel toolbox, is a `parfor` loop, since all calculations are independent. If `m` is very large, it might even be best to only use `parfor`, since the 3D matrix could be too large to fit in memory. – Adriaan Sep 13 '18 at 08:52