Two questions on vectorized matrix vector multiplication

The first question is:

x0 = [3;8]
K = [1 2 3 8 9 10; 2 3 9 -1 2 3];

In general K is of size 2 by n, where n is even.

How can I obtain a matrix `E` such that:

 E = [expm([1 2;2 3]) expm([3 8;9 -1]) expm([9 10; 2 3])]

The second question is.

How can I obtain a matrix J such that:

J = [expm([1 2;2 3]*1)*x0 expm([1 2;2 3]*2)*x0 expm([1 2;2 3]*3)*x0]

For the second question: If you know, J here is the solution of ode without using ode solver, with time as 1:1:100 and initial x0

Answer 1

[Walter Roberson]

[r, n] = size(K);
E = zeros(r, n);
for idx = 1 : 2 : n
  E(:, idx:idx+1) = expm(K(:,idx:idx+1));
end

You could convert this into an arrayfun and so not have an explicit loop, but you cannot do this without some level of looping, because expm() only works on square matrices.

You could also reshape K to 2 x 2 x m and then loop over the pages of m:

 K2 = reshape(K, 2, 2, []);
 m = size(K2,2);
 E = zeros(size(K2));
 for idx = 1 : m
   E(:,:,idx) = expm(K(:,:,m));
 end
 E = reshape(E, 2, []);

That could be advantageous if you were switching to GPU computing, as you would be able to use pagefun(@expm, gpu_E)

[Walter Roberson]

It would be trivial to extend this to multiply by x0 -- you just have to watch out because it gives you a different sized output than what I show here.

[K.E.]

but it seems that i have to multiply by x0 everytime inside the loop.loop? Is there a way to repmat x0 and do vectorization multiplication or kron after the loop?

[Walter Roberson]

No. Remember you are doing algebraic matrix multiplication. You have a 2 x 2 "*" a 2 x 1 getting out a 2 x 1 result, so each pair of columns of your E matrix has to collapse into one column in your J matrix. You cannot just repmat: you have to preserve that behavior. That could, for example, involve transposing and left multiplying and transposing the result, in the general case. Or just doing the equivalent operations directly for the special case of 2 x 2.

In my tests, doing the operations directly is a little slower.

If you have R2016b or later you can put be below into a single script; otherwise you will need to put the three parts into separate files.

x0 = [3;8];
N = 10000;
K = randi([-100 100], 2, N);
timeit( @() mulit(K ,x0), 0)
timeit( @() mulit2(K ,x0), 0)

function E = mulit(K, x0)
    [r, n] = size(K);
    m = n/2;
    E = zeros(r, m);
    for idx = 1 : m
        E(:, idx) = expm(K(:,idx*2-[1 0])) * x0;
    end
end

function E = mulit2(K, x0)
    [r, n] = size(K);
    m = n/2;
    F = zeros(r, n);
    for idx = 1 : 2 : n
        F(:, idx:idx+1) = expm(K(:,idx:idx+1));
    end
    E = F(:,1:2:end) * x0(1) + F(:,2:2:end) * x0(2);
end

[Matt J]

That could, for example, involve transposing and left multiplying and transposing the result, in the general case.

If E is maintained in MxMxN form, then you could also do the multiplication with MTIMESX

   for idx = 1 : m
    E(:,:,idx) = expm(K(:,:,m));
   end
   E = reshape(E, 2, []);

   J=reshape(mtimesx(E,x0),[],m);