multi array vectorization problem (reformulated)

1 visualización (últimos 30 días)
Michal
Michal el 18 de Mayo de 2022
Editada: Michal el 19 de Mayo de 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 comentarios
Jan
Jan el 18 de Mayo de 2022
(V.*exp(D*t))/V, Typical size of V is ~ 3 x 3, Typical size of D.*t' is ~ (1e4 x 3)
I'm confused. Does "D.*t' is [1e4 x 3]" mean, that D*t is [3 x 1e4] ? Why is it .* in one case and * in the other?
The best idea is to provide Matlab code, which creates the typical input data. This avoids ambiguities.
Michal
Michal el 19 de Mayo de 2022
Editada: Michal el 19 de Mayo de 2022
@Jan I just completely reformulated my question ...

Iniciar sesión para comentar.

Respuesta aceptada

Matt J
Matt J el 19 de Mayo de 2022
Editada: Matt J el 19 de Mayo de 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 comentario
Michal
Michal el 19 de Mayo de 2022
Editada: Michal el 19 de Mayo de 2022
expmAtb4=V*(exp(d*t).*(V\b));
Simple, nice and fast!

Iniciar sesión para comentar.

Más respuestas (2)

Bruno Luong
Bruno Luong el 19 de Mayo de 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 comentarios
Bruno Luong
Bruno Luong el 19 de Mayo de 2022
Editada: Bruno Luong el 19 de Mayo de 2022
"In a case of further generalization"
That's why one should not ask too simplified question at the first place.
The best solution always depends on the detail/size of the problem, and the MATLAB version one is using.
Michal
Michal el 19 de Mayo de 2022
Editada: Michal el 19 de Mayo de 2022
@Bruno Luong I am sorry, you are obviously right. Thanks for help. Your final solution is good!

Iniciar sesión para comentar.


Catalytic
Catalytic el 19 de Mayo de 2022
Editada: Catalytic el 19 de Mayo de 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 comentarios
Michal
Michal el 19 de Mayo de 2022
Your result array expmAtb3 has the size (2,1,length(t)) so you need to add squeeze command:
expmAt=squeeze(pagemtimes( pagemtimes(V,expmAtb__), V\b));
Matt J
Matt J el 19 de Mayo de 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.

Iniciar sesión para comentar.

Categorías

Más información sobre Automotive en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by