Two questions on vectorized matrix vector multiplication

2 次查看（过去 30 天）
K.E.
K.E. 2017-7-9
评论： Matt J 2017-7-9
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

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2017-7-9
[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)
  4 个评论
显示 2更早的评论
Walter Roberson
Walter Roberson 2017-7-9
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
Matt J 2017-7-9
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);

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Arithmetic Operations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by