Vectorising Multiplying each column from a matrix with a square matrix and the tranpose of this column

Hello,

I have an operation that I'm doing with a loop but would like to vectorise if possible.

I have a matrix h (2790x3591) and a matrix Pn(2790x2790).

I need to implement the following operation: h_column' * Pn * h_column, for all columns of h.

Is there a way to vectorise this operation? It takes about 45s to do with the for loop.

NB: The matrices are complex with very small imaginary parts, so that adds a bit of delay because the multiplications are more involved.

Answer 1

[Kevin Holly]

Can you do the following?

h = rand(2790,3591);
Pn = rand(2790,2790);

size(h'*Pn*h)

ans = 1×2

3591 3591

[Kleanthis-Marios Papadopoulos]

I am not looking for direct matrix multiplication.

I believe I found the answer.

I need to take diag(h'*Pn*h).

Answer 2

[James Tursa]

If you need the same column on each side of the multiply, then maybe this:

result = sum(h'*Pn.*h',2);

[Kleanthis-Marios Papadopoulos]

I don't need the same column on each side because the operation is column of h' (turn it into row)*Pn* h.

This isn't what I'm looking for, although it is close.

[Matt J]

I need to take diag(h'*Pn*h).

No, that will be super slow and inefficient. This should be what you want,

result = sum(h.*(Pn*h),1);

James' version was almost equivalent except for a missing transpose,

result = sum(h'*Pn'.*h',2);

but no need to incur the expense of all those transpositions, as far as I can see.

Answer 3

[Jan]

Then please post you code. It is much faster on my computer:

h = rand(2790, 3591);
Pn = rand(2790, 2790);

tic
R = zeros(1, 3591);
for k = 1:3591
   R(k) = h(:, k).' * Pn * h(:, k);
end
toc

Elapsed time is 8.494673 seconds.

tic  % Matt J and James:
R2 = sum(h'*Pn'.*h',2);
toc

Elapsed time is 0.508067 seconds.

[Kleanthis-Marios Papadopoulos]

This operation is part of a function which is part of another function, so the code here on its own won't make much sense.

This is the snippet that I'm concerned about

for j=1:Nsc

hj= findhVectorised(theta,l_ik,uk,J,Fjvec(j),Nsc,array,userpn);

holderonj=[holderonj;hj];

%dotp=[dotp; hj.*conj(hj)];

%numer=[numer; sum(dotp)]

numer(j,:)=sum(hj.*conj(hj),1);

% for iter=1:length(hj)

% den(j,iter)= hj(:,iter)'*Pn*hj(:,iter);

% end

% den2= hj'*Pn*hj;

temp=diag(hj'*Pn*hj);

temp=transpose(temp);

den(j,:)=temp;

end