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

2 views (last 30 days)
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.

Answers (3)

Kevin Holly
Kevin Holly on 23 Feb 2022
Can you do the following?
h = rand(2790,3591);
Pn = rand(2790,2790);
size(h'*Pn*h)
ans = 1×2
3591 3591

James Tursa
James Tursa on 23 Feb 2022
If you need the same column on each side of the multiply, then maybe this:
result = sum(h'*Pn.*h',2);
  2 Comments
Matt J
Matt J on 23 Feb 2022
Edited: Matt J on 23 Feb 2022
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.

Sign in to comment.


Jan
Jan on 23 Feb 2022
Edited: Jan on 24 Feb 2022
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.
  1 Comment
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

Sign in to comment.

Products


Release

R2021a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!