performance - 如何矢量化许多 rank-1 外积的形成？

下面的循环可以矢量化吗？循环中的每次迭代都会形成一个外积，然后将结果对称化并将结果存储为矩阵中的一列。预计 m 较大(例如 1e4)而 s 较小(例如 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

编辑

这是 3 种不同方法的比较:迭代大维度 m 、迭代小维度 s 和基于 bsxfun 的解决方案(公认的最快答案)。

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

下面是时间对比:

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

最佳答案

使用 bsxfun (它是如何完成的，很有趣):

% 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] )';

基准测试结果(我的机器):

原始方法 (迭代大暗淡):

 Elapsed time is 0.217695 seconds.

"new"方法 (迭代较小的暗淡):

 Elapsed time is 0.037538 seconds.

有趣的 bsxfun :

 Elapsed time is 0.021507 seconds.

如您所见，bsxfun 占用了最快循环时间的 ~ 2/3 - 1/2 ...

bsxfun 不是很有趣吗？

关于performance - 如何矢量化许多 rank-1 外积的形成？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/17549434/