Matlab Vectorization - Time evolution of a matrix

Question

I have tried searching for an answer on this, but can't find one that specifically addresses my issue. Although vectorizing in MATLAB must draw many questions in, the problem I am having is less general than typical examples I have found on the web. My background is more in C++ than MATLAB, so this is an odd concept for me to get my head around.

I am trying to evolve a Hamiltonian matrix from it's initial state (being a column vector where all elements but the last is a 0, and the last is a 1) to a final state as time increases. This is achieved by sequentially applying a time evolution operator U to the state. I also want to use the new state at each time interval to calculate an observable property.

I have achieved this, as can be seen in the code below. However, I need to make this code as efficient as possible, and so I was hoping to vectorize, rather than rely on for loops. However, I am unsure of how to vectorize this code. The problem I have is that on each iteration of the for loop, the column vector psi should change its values. Each new psi is then used to calculate my observable M for each interval of time. I am unsure of how to track the evolution of psi such that I can end up with a row vector for M, giving the outcome of the application of each new psi.

time = tmin:dt:tmax;
H = magic(2^N)
X = [0,1;1,0]

%%% INITIALISE COLUMN VECTOR
init = sparse(2^N,1);
init(2^N) = 1;

%%% UNITARY TIME EVOLUTION OPERATOR
U = expm(-1i*H*dt);

%%% TIME EVOLVUTION
for num = 1:length(time)
    psi = U*init;
    init = psi;

    %%% CALCULATE OBSERVABLE
    M(num) = psi' * kron(X,speye(2^(N-1))) * psi

end

Any help would be greatly appreciated.

What is X in M(num) = psi' * kron(X,speye(2^(N-1))) * psi? — rahnema1
Your example code is incomplete and I'm not sure it's correct. The loop doesn't seems to depend on previous iterations and X is undefined. You may be prematurely optimizing. Get your code working and verified with a loop first. Also, for loops are not necessarily slow and integration problems like this are often best solved with a loop (though there are ways to remove the loop for simple cases like this). — horchler
@horchler - I've kept the code incomplete primarily for readability, although of course it is not really testable in its current format. This main thing missing is H, which for the purposes of this can just be thought of as a matrix with dimensions (2^N, 2^N). The code does work, and I have verified it. The loop does depend on previous iteration - the operator U is applied to init to give psi. The value of psi is then assigned to init, and thus on each iteration of the loop, psi changes, and so too then does M. — Epsilon0

pkic pkic · Accepted Answer · 2017-10-19T11:40:04

I have quickly come up with the following partially vectorized code:

time = tmin:dt:tmax;
H = magic(2^N);
X = [0,1;1,0];

%%% INITIALISE COLUMN VECTOR
init = sparse(2^N,1);
init(2^N) = 1;

%%% UNITARY TIME EVOLUTION OPERATOR
U = expm(-1i*H*dt);

%%% TIME EVOLVUTION
% preallocate psi
psi = complex(zeros(2^N, length(time)));

% compute psi for all timesteps
psi(:,1) = U*init;
for num = 2:length(time)
    psi(:,num) = U*psi(:, num-1);    
end

% precompute kronecker product (if X is constant through time)
F = kron(X,speye(2^(N-1)));

%%% CALCULATE OBSERVABLE
M = sum((psi' * F) .* psi.', 2);

However, it seems that the most computationally intensive part of your problem is computation of the psi. For that I can't see any obvious way to vectorize as it depends on the value computed in the previous step.

This line: M = sum((psi' * F) .* psi.', 2);

is a little Matlab trick to compute psi(:,i)'*F*psi(:,i) in a vectorized way.

Matlab Vectorization - Time evolution of a matrix

1 Answers