multi array vectorization problem (reformulated)

4 views (last 30 days)
Michal
Michal on 18 May 2022
Edited: Michal on 19 May 2022
How to vectorize (for-loop elimination) the following code?
% init
A = [ -0.8013 -0.4981; -0.2278 -0.9009];
t = 0:0.01:1e2;
% eigen space
[V,D]=eig(A,'vector');
D = reshape(D,1,[]);
% computing
expmAt = zeros([size(A),length(t)]);
for i = 1:length(t)
expmAt(:,:,i) = V.*(exp(D*t(i)))/V;
end
  8 Comments
Michal
Michal on 19 May 2022
@Jan I just completely reformulated my question ...

Sign in to comment.

Accepted Answer

Matt J
Matt J on 19 May 2022
Edited: Matt J on 19 May 2022
runtest(1e2)
Elapsed time is 0.002075 seconds. Elapsed time is 0.000805 seconds.
function runtest(N)
A = [ -0.8013 -0.4981; -0.2278 -0.9009];
t = 0:0.01:N;
b = [2;3];
[V,d]=eig(A,'vector');
tic;
E=exp(d*t);
expmAtb=repmat(eye(size(V)),1,1,numel(t));
expmAtb(logical(expmAtb))=E(:);
expmAtb3=squeeze(pagemtimes( pagemtimes(V,expmAtb), V\b));
toc
tic;
expmAtb4=V*(exp(d*t).*(V\b));
toc
end
  1 Comment
Michal
Michal on 19 May 2022
expmAtb4=V*(exp(d*t).*(V\b));
Simple, nice and fast!

Sign in to comment.

More Answers (2)

Bruno Luong
Bruno Luong on 19 May 2022
A = [ -0.8013 -0.4981; -0.2278 -0.9009];
t = 0:3;
% eigen space
[V,D]=eig(A,'vector');
D = reshape(D,1,[]);
% computing
expmAt = zeros([size(A),length(t)]);
for i = 1:length(t)
expmAt(:,:,i) = V.*(exp(D*t(i)))/V;
end
expmAt
expmAt =
expmAt(:,:,1) = 1 0 0 1 expmAt(:,:,2) = 0.4736 -0.2168 -0.0991 0.4303 expmAt(:,:,3) = 0.2458 -0.1960 -0.0896 0.2066 expmAt(:,:,4) = 0.1358 -0.1376 -0.0629 0.1083
expmAt2 = pagemrdivide(V.*exp(D.*reshape(t,1,1,[])),V)
expmAt2 =
expmAt2(:,:,1) = 1 0 0 1 expmAt2(:,:,2) = 0.4736 -0.2168 -0.0991 0.4303 expmAt2(:,:,3) = 0.2458 -0.1960 -0.0896 0.2066 expmAt2(:,:,4) = 0.1358 -0.1376 -0.0629 0.1083
  5 Comments
Michal
Michal on 19 May 2022
@Bruno Luong I am sorry, you are obviously right. Thanks for help. Your final solution is good!

Sign in to comment.


Catalytic
Catalytic on 19 May 2022
Edited: Catalytic on 19 May 2022
% init
A = [ -0.8013 -0.4981; -0.2278 -0.9009];
t = 0:0.01:1e2;
tic;
% eigen space
[V,D]=eig(A,'vector');
D = reshape(D,1,[]);
% computing
expmAt = zeros([size(A),length(t)]);
for i = 1:length(t)
expmAt(:,:,i) = V.*(exp(D*t(i)))/V;
end
toc
Elapsed time is 0.057550 seconds.
tic
out=theTask(A,t);
toc
Elapsed time is 0.007411 seconds.
[lb,ub]=bounds(expmAt-out,'all')
lb = -1.3878e-16
ub = 1.1102e-16
function expmAt=theTask(A,t)
[V,d]=eig(A,'vector');
E=exp(d*t);
expmAt=repmat(eye(size(A)),1,1,numel(t));
expmAt(logical(expmAt))=E(:);
expmAt=pagemtimes( pagemtimes(V,expmAt), inv(V));
end
  3 Comments
Matt J
Matt J on 19 May 2022
I wouldn't expect the addition of squeeze() to impact timing significantly.
Also, I find it strange that pagemrdivide(A,B) is not optimized for the case where B has only one page.

Sign in to comment.

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by