From: Ian on
Hi,

Suppose K is an array with size(K) = [n n L].

I'd really appreciate help trying to vectorize the following loop to solve for A:

A = zeros(n^2, n^2);
for j = 1:L
A = A + kron(conj(K(:,:,j)), K(:,:,j));
end

(context: this converts a set of kraus operators (K) into the Liouville superoperator form (A) in quantum mechanics)

Thanks for any help,
Ian
From: Bruno Luong on
Ian <ian.hincks(a)gmail.com> wrote in message <1148922034.56639.1280946950547.JavaMail.root(a)gallium.mathforum.org>...
> Hi,
>
> Suppose K is an array with size(K) = [n n L].
>
> I'd really appreciate help trying to vectorize the following loop to solve for A:
>
> A = zeros(n^2, n^2);
> for j = 1:L
> A = A + kron(conj(K(:,:,j)), K(:,:,j));
> end

K1 = reshape(conj(K),[1 n 1 n L]);
K2 = reshape(K,[n 1 n 1 L]);
A = bsxfun(@times,K1,K2);
A = sum(reshape(A,[n^2 n^2 L]),3)

% Bruno
From: Ian on
Here's a solution:

KK = K(:,:);
n = size(KK,1);
w = size(KK,2);

q = [1:n];
p = [1:w];
r = [1:w/n];

a1 = reshape( q(ones(1,n),:), [1, n^2]);
a2 = reshape( p(ones(1,n),:), [1, w*n]);
b1 = repmat(q, [1 n]);
b2 = reshape( p, [n, w/n] );
b2 = reshape( b2(:, reshape( r(ones(1,n),:), [1, w] )), [1 w*n]);

A = sum( reshape(conj(KK(a1,a2)) .* KK(b1,b2), n^2, n^2, w/n), 3);
From: Ian on
Thanks Bruno,

That's definitely more elegant than the solution I contrived and posted just now. All of these little matlab tricks can be overwhelming.

Ian
From: Bruno Luong on
Ian <ian.hincks(a)gmail.com> wrote in message <1056167390.62207.1281022752965.JavaMail.root(a)gallium.mathforum.org>...
> All of these little matlab tricks can be overwhelming.

Hi Ian,

Overwhelming and rather fun. That's why many of us are hanging around here... ;-)

Bruno