multi array vectorization problem (reformulated)
显示 更早的评论
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 个评论
Jan
2022-5-18
Which code do you want to accelerate? What are typical inputs? Das "many" mean 25 or 10e6 ?
Matt J
2022-5-18
Is A positive definite? How do you know V is non-singular?
Matt J
2022-5-18
But not symmetric?
Michal
2022-5-18
Jan
2022-5-18
(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.
采纳的回答
runtest(1e2)
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
更多回答（2 个）
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
expmAt2 = pagemrdivide(V.*exp(D.*reshape(t,1,1,[])),V)
5 个评论
How much it accelerates? Result on TMW online server:
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
tic
expmAt = zeros([size(A),length(t)]);
for i = 1:length(t)
expmAt(:,:,i) = V.*(exp(D*t(i)))/V;
end
toc
tic
expmAt2 = pagemrdivide(V.*exp(D.*reshape(t,1,1,[])),V);
toc
Bruno Luong
2022-5-19
Only one pagemtimes is enough. Also never use squeeze if you want a predictable code.
expmAtb2 = reshape(pagemtimes(V.*exp(D.*reshape(t,1,1,[])) ,V\b),[size(V,1),length(t)]);
Bruno Luong
2022-5-19
编辑：Bruno Luong
2022-5-19
"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.
% 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
tic
out=theTask(A,t);
toc
[lb,ub]=bounds(expmAt-out,'all')
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 个评论
In a case of further generalization ...is the following code the "optimal" solution?
No. It wasn't optimal before the generalization either -
% init
A = [ -0.8013 -0.4981; -0.2278 -0.9009];
t = 0:0.01:1e2;
b = [2;3];
tic;
% eigen space
[V,D]=eig(A,'vector');
D = reshape(D,1,[]);
% computing
expmAtb2 = squeeze(pagemtimes(pagemrdivide(V.*exp(D.*reshape(t,1,1,[])),V),b));
toc
tic
expmAtb3=theTask(A,t,b);
toc
[lb,ub]=bounds(expmAtb2(:)-expmAtb3(:),'all')
function expmAt=theTask(A,t,b)
[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), V\b);
end
Michal
2022-5-19
Matt J
2022-5-19
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.
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Translated by 
