CODEWITHSUNDEEP
performancematlabfor-loopmatrixvectorization
Pincopallino
Pincopallino

Reputation: 651

MATLAB way for representing and operating on an array of matrices

I'm writing a program that calculates the time evolution of a two qubit density matrix (that is a 4 x 4 matrix) and I would like to apply certain functions to the density matrix at each time step. For example, a typical function would be trace(rho^2) where rho is the density matrix at a particular time.

In Mathematica I can easily achieve this by creating a list of matrices and using Map to apply the desired function to each matrix. In MATLAB, it seemed obvious to me to implement my time evolution function so that it returns a 4 x 4 x NSTEPS array, but then the only way that came to my mind for applying the function to each 4 x 4 matrix was a for loop, like

for i = 1:size(rho,3)
   p(i) = real(trace(rho(:,:,i)^2));
end

Since I've been taught to avoid for loops as much as possible in MATLAB, what is the MATLAB way (or the most efficient way) of doing such thing?

Upvotes: 1

Views: 179

Answers (2)

Divakar
Divakar

Reputation: 221644

Version 1

p = squeeze(sum(sum(bsxfun(@times,rho,permute(rho,[2 1 3])))))

Version 2

p = sum(reshape(bsxfun(@times,rho,permute(rho,[2 1 3])),[],size(rho,3)))

Benchmarking Code

%%// Original for-loop code
tic
for k = 1:10
    p = zeros(size(rho,3),1);
    for i = 1:size(rho,3)
        p(i) = real(trace(rho(:,:,i)^2));
    end
end
toc
clear p

%%// Version 1
tic
for k = 1:10
    p = squeeze(sum(sum(bsxfun(@times,rho,permute(rho,[2 1 3])))));
end
toc
clear p

%%// Version 2
tic
for k = 1:10
    p = sum(reshape(bsxfun(@times,rho,permute(rho,[2 1 3])),[],size(rho,3)));
end
toc

Benchmarking Results

With rho as 50x50x50

Elapsed time is 0.061841 seconds.
Elapsed time is 0.008360 seconds.
Elapsed time is 0.004700 seconds.

With rho as 200x200x50

Elapsed time is 0.421595 seconds.
Elapsed time is 0.140892 seconds.
Elapsed time is 0.135933 seconds.

With rho as 500x500x50

Elapsed time is 4.973693 seconds.
Elapsed time is 0.899126 seconds.
Elapsed time is 0.894188 seconds.

Upvotes: 1

Dennis Jaheruddin
Dennis Jaheruddin

Reputation: 21563

Assuming you have preallocated your p, I doubt that you can get much of a speedup in any way. Here is how you can make sure the p does not grow every step of the way:

for i = size(rho,3):-1:1
   p(i) = real(trace(rho(:,:,i)^2));
end

Make sure to run a benchmark test, but I would be surprised if this matlab code turned out to be slower than the mathematica equivalent.

Upvotes: 0

Related Questions