如何在 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中单个row/column的索引，ids_A是所有其他[的索引向量=19=]（因此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 执行元素级产品。

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

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

关于代码片段本身，可以看出你是 "escaping" 每次迭代的索引 setxor。现在，如果您要使用矢量化方法，您可以一次性执行所有这些数学运算，然后删除矢量化方法中包含但未包含的元素故意的。这确实是接下来列出的基于 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 拖慢了我们的速度。

如何在 matlab 中分析矢量外积

How to profile a vector outer product in matlab

matlab

performance

vector

profiling

基准测试