CODEWITHSUNDEEP
performancematlabvectorization
Michael
Michael

Reputation: 557

How to vectorize the formation of many rank-1 outer products?

Can the loop below be vectorized? Each iteration in the loop forms an outer product, then symmetrizes and stores the result as a column in a matrix. It's anticipated that m is large (e.g., 1e4) and s is small (e.g., 10).

% U and V are m-by-s matrices
A = zeros(s^2, m); % preallocate
for k = 1:m
    Ak = U(k,:)' * V(k,:);
    Ak = (Ak + Ak')/2;
    A(:, k) = Ak(:);
end

Edit

Here is a comparison of 3 different methods: iterate over the the large dimension m, iterate over the small dimension s, and a bsxfun based solution (the accepted and fastest answer).

s = 5; m = 100000;
U = rand(m, s);
V = rand(m, s);

% Iterate over large dimension
tic
B = zeros(s^2, m);
for k = 1:m
    Ak = U(k,:)' * V(k,:);
    Ak = (Ak + Ak')/2;
    B(:, k) = Ak(:);
end
toc

% Iterate over small dimension
tic
A = zeros(s, s, m);
for i = 1:s
    A(i,i,:) = U(:, i) .* V(:, i);
    for j = i+1:s
        A(i,j,:) = (U(:,i).*V(:,j) + U(:, j).*V(:, i))/2;
        A(j,i,:) = A(i,j,:);
    end
end
A = reshape(A, [s^2, m]);
toc

% bsxfun-based solution
tic
A = bsxfun( @times, permute( U, [1 3 2] ), permute( V, [ 1 2 3 ] ) );
A = .5 * ( A + permute( A, [1 3 2] ) );
B = reshape( A, [m, s^2] )';
toc

Here is a time comparison:

Elapsed time is 0.547053 seconds.
Elapsed time is 0.042639 seconds.
Elapsed time is 0.039296 seconds.

Upvotes: 4

Views: 92

Answers (1)

Shai
Shai

Reputation: 114916

Use bsxfun (it's how its done with lots of fun):

% the outer product
A = bsxfun( @times, permute( U, [1 3 2] ), permute( V, [ 1 2 3 ] ) );
% symmetrization
A = .5 * ( A + permute( A, [1 3 2] ) );
% to vector (per product)
B = reshape( A, [m s^2] )';

Benchmark results (my machine):

  1. Original approach (iterate over large dim): 

     Elapsed time is 0.217695 seconds.

  2. "New" approach (iterate over smaller dim):

     Elapsed time is 0.037538 seconds.

  3. Fun with bsxfun:

     Elapsed time is 0.021507 seconds.

As you can see bsxfun takes ~ 2/3 - 1/2 of the time of the fastest loop...

Isn't it just fun with ?

Upvotes: 1

Related Questions