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

[Jan]

Which code do you want to accelerate? What are typical inputs? Das "many" mean 25 or 10e6 ?

[Michal]

This one:

(V.*exp(D*t))/V

Typical size of V is ~ 3 x 3

Typical size of D.*t' is ~ (1e4 x 3)

[Matt J]

Is A positive definite? How do you know V is non-singular?

[Michal]

In my specific application the V is always non-singular.

[Matt J]

But not symmetric?

[Michal]

and non-symmetric

[Jan]

(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]

@Jan I just completely reformulated my question ...

Answer 1

[Matt J]

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

[Michal]

expmAtb4=V*(exp(d*t).*(V\b));

Simple, nice and fast!

Answer 2

[Bruno Luong]

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

[Bruno Luong]

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

Elapsed time is 0.029968 seconds.

tic
expmAt2 = pagemrdivide(V.*exp(D.*reshape(t,1,1,[])),V);
toc

Elapsed time is 0.008685 seconds.

[Michal]

In a case of further generalization:

% init
A = [  -0.8013   -0.4981;  -0.2278   -0.9009];
t = 0:0.01:1e2;
b = [2;3];

% eigen space
[V,D]=eig(A,'vector');
D = reshape(D,1,[]);

% computing
tic
expmAtb = zeros(length(A),length(t));
for i = 1:length(t)
    expmAtb(:,i) = V.*(exp(D*t(i)))/V * b;
end
toc

is the following code the "optimal" solution?

tic
expmAtb2 = squeeze(pagemtimes(pagemrdivide(V.*exp(D.*reshape(t,1,1,[])),V),b));
toc

[Bruno Luong]

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]

"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]

@Bruno Luong I am sorry, you are obviously right. Thanks for help. Your final solution is good!

Answer 3

[Catalytic]

% 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

[Catalytic]

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

Elapsed time is 0.016474 seconds.

tic
expmAtb3=theTask(A,t,b);
toc

Elapsed time is 0.003397 seconds.

[lb,ub]=bounds(expmAtb2(:)-expmAtb3(:),'all')

lb = -8.8818e-16

ub = 4.4409e-16

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]

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]

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.