Modify and speed up a for loop

Hi everyone,

Do you have some idea on how to modify and speed up this for loop?

N usually is among 10 and 100, while K is of the order of 10^4. Moreover, u is an N x K matrix and kernel_vect is a K x 1 vector.

Thanks!

Kernel_appo = zeros(N^2,1);
for k=1:K
   Mat_appo = (u(:,k)*u(:,k)');
   Kernel_appo = Kernel_appo  + kernel_vect(k)*Mat_appo(:);
end

[Adam]

Depending how big k is you can potentially pull out the Mat_appo line similar to the following logic:

u = reshape( 1:20, [4 5] ) % Some small test data - K = 5
Mat_appo = u .* reshape( u', [1 5 4] );
Mat_appo = reshape( permute( Mat_appo, [1 3 2] ), 4^2, [] );

then you should end up with what you have as Mat_appo(:) as the columns of the Mat_appo matrix above.

You can probably simplify those reshsapes and permutes too, that's just how I happened to get to the answer, which usually just involves some trial and error reshaping and permuting rather than working out the absolutely neatest way!

Answer 1

[Cyrus Tirband]

The last line can be brought outside the loop like so

Mat_appo = zeros(N,N,K);
for k=1:K
   Mat_appo(:,:,k) = (u(:,k)*u(:,k)');
end
Kernel_appo=reshape(Mat_appo,N^2,[])*kernel_vect(:);

Since K is much larger than N, I reckon it would be faster to rewrite the for loop as follows:

Mat_appo = zeros(N,N,K);
for i = 1:N
    for j = 1:N
        Mat_appo(i,j,:) = u(i,:)*u(j,:);
    end
end
Kernel_appo=reshape(Mat_appo,N^2,[])*kernel_vect(:);

I haven't tested it, but this code should be equivalent to yours and also faster.

[Stef Fets]

Hi Cyrustd,

Thanks for your answer, but I don’t understand well it…

Maybe I should specify that u is an N x K matrix and kernel_vect is a k x 1 vector. All the quantities are real numbers.

[Cyrus Tirband]

Does the code not run for you? It should do the exact same as your code. It first generates an N x N x K matrix for Mat_appo where the kth N x N matrix is the outerproduct of the kth u vector.

Then, it rearranges Mat_appo to a N^2 x K matrix and does a matrix multiplication with kernel_vect (which is K x 1), the result is a N^2 x 1 matrix called Kernel_appo.