如何在matlab中分析向量外积

Question

在我的 matlab 分析过程中，我注意到一行代码消耗的时间比我想象的要多得多。知道如何让它更快吗？

X = Y(ids_A, ids_A) - (Y(ids_A,k) * Y(k,ids_A))/Y(k,k); 

X 和 Y 是具有相同大小 (dxd) 的对称矩阵，k 是 Y 中单个行/列的索引，ids_A 是所有其他行/列的索引向量（因此 Y(ids_A, k) 是列向量，Y(k,ids_A) 是行向量)

ids_A = setxor(1:d,k); 

谢谢！

Answer 1

您也许可以通过调用bsxfun 来替换外部乘法：

X = Y(ids_A, ids_A) - (bsxfun(@times, Y(ids_A,k), Y(k,ids_A))/Y(k,k));

那么上面的代码是如何工作的呢？我们来看看当一个向量是4个元素，其他3个元素时外积的定义：

^{来源：Wikipedia}

如您所见，外部乘积由元素乘积创建，其中第一个向量 u 水平复制，而第二个向量 v 垂直复制。然后，您可以找到每个元素的元素乘积来产生结果。这是通过bsxfun 雄辩地完成的：

bsxfun(@times, u, v.');

u 是列向量，v.' 是行向量。 bsxfun 自然地复制数据以遵循上述模式，然后我们使用 @times 来执行 element-wise 产品。

Answer 2

我假设你的代码看起来像这样 -

for k = 1:d
    ids_A = setxor(1:d,k);
    X = Y(ids_A, ids_A) - (Y(ids_A,k) * Y(k,ids_A))/Y(k,k);
end

使用给定的代码 sn-p，可以安全地假设您在该循环中以某种方式使用 X。您可以在此类循环开始之前计算所有X 矩阵作为预计算步骤，并且这些计算可以作为矢量化方法执行。

关于代码 sn-p 本身，可以看出您在每次迭代中使用setxor“转义”一个索引。现在，如果您使用矢量化方法，您可以在 one-go 中执行所有这些数学运算，然后删除合并到矢量化方法中但不是预期的元素。这确实是下面列出的基于 bsxfun 的矢量化方法的精髓 -

%// Perform all matrix-multiplications in one go with bsxfun and permute
mults = bsxfun(@times,permute(Y,[1 3 2]),permute(Y,[3 2 1]));

%// Scale those with diagonal elements from Y and get X for every iteration
scaledvals = bsxfun(@rdivide,mults,permute(Y(1:d+1:end),[1 3 2]));
X_vectorized = bsxfun(@minus,Y,scaledvals);

%// Find row and column indices as linear indices to be removed from X_all
row_idx = bsxfun(@plus,[0:d-1]*d+1,[0:d-1]'*(d*d+1));
col_idx = bsxfun(@plus,[1:d]',[0:d-1]*(d*(d+1)));

%// Remove those "setxored" indices and then reshape to expected size
X_vectorized([row_idx col_idx])=[];
X_vectorized = reshape(X_vectorized,d-1,d-1,d);

基准测试

基准代码

d = 50;          %// Datasize
Y = rand(d,d);    %// Create random input
num_iter = 100;   %// Number of iterations to be run for each approach

%// Warm up tic/toc.
for k = 1:100000
    tic(); elapsed = toc();
end

disp('------------------------------ With original loopy approach')
tic
for iter = 1:num_iter
    for k = 1:d
        ids_A = setxor(1:d,k);
        X = Y(ids_A, ids_A) - (Y(ids_A,k) * Y(k,ids_A))/Y(k,k);
    end
end
toc
clear X k ids_A

disp('------------------------------ With proposed vectorized approach')
tic
for iter = 1:num_iter
    mults = bsxfun(@times,permute(Y,[1 3 2]),permute(Y,[3 2 1]));
    scaledvals = bsxfun(@rdivide,mults,permute(Y(1:d+1:end),[1 3 2]));
    X_vectorized = bsxfun(@minus,Y,scaledvals);

    row_idx = bsxfun(@plus,[0:d-1]*d+1,[0:d-1]'*(d*d+1));
    col_idx = bsxfun(@plus,[1:d]',[0:d-1]*(d*(d+1)));

    X_vectorized([row_idx col_idx])=[];
    X_vectorized = reshape(X_vectorized,d-1,d-1,d);
end
toc

结果

案例#1：d = 50

------------------------------ With original loopy approach
Elapsed time is 0.849518 seconds.
------------------------------ With proposed vectorized approach
Elapsed time is 0.154395 seconds.

案例#2：d = 100

------------------------------ With original loopy approach
Elapsed time is 2.079886 seconds.
------------------------------ With proposed vectorized approach
Elapsed time is 2.285884 seconds.

案例#1：d = 200

------------------------------ With original loopy approach
Elapsed time is 7.592865 seconds.
------------------------------ With proposed vectorized approach
Elapsed time is 19.012421 seconds.

结论

人们很容易注意到，在处理大小不超过100 x 100 的矩阵时，所提出的矢量化方法可能是更好的选择 内存饥渴的 bsxfun 让我们慢下来。