使用矩阵结构加速Matlab

Question

假设我有一个 N x K 的矩阵 A 和一个 N x P 的矩阵 B。我想进行以下计算，以获得最终的 N x P 矩阵 X。

X(n,p) = B(n,p) - dot(gamma(p,:),A(n,:))

where

gamma(p,k) = dot(A(:,k),B(:,p))/sum( A(:,k).^2 )

在MATLAB中，我的代码如下：

for p = 1:P
    for n = 1:N
        for k = 1:K
            gamma(p,k) = dot(A(:,k),B(:,p))/sum(A(:,k).^2);
        end
        x(n,p) = B(n,p) - dot(gamma(p,:),A(n,:));
    end
end

这段代码使用了三个循环，效率非常低下！有没有好的方法可以加速这段代码？

- Langma

---

我认为你打错了字...s是什么？j又是什么？ - Rody Oldenhuis

另外，不要使用变量名gamma，因为这会覆盖Matlab内置的gamma函数... - bla

Answer 1

在我看来，你可以将gamma计算提升出循环；至少，在gamma计算中我没有看到任何对N的依赖。

因此，可以这样做：

for p = 1:P
    for k = 1:K
        gamma(p,k) = dot(A(:,k),B(:,p))/sum(A(:,k).^2);
    end
end
for p = 1:P
    for n = 1:N
        x(n,p) = B(n,p) - dot(gamma(p,:),A(n,:));
    end
end

我不太熟悉你的代码（或matlab），无法确定你是否可以合并这两个循环，但如果可以：

for p = 1:P
    for k = 1:K
        gamma(p,k) = dot(A(:,k),B(:,p))/sum(A(:,k).^2);
    end
    for n = 1:N
        x(n,p) = B(n,p) - dot(gamma(p,:),A(n,:));
    end
end

- Eric Brown

---

据我所知，这样的合并在任何语言中都应该是可能的 ^_^ - Rody Oldenhuis

@RodyOldenhuis - 我不确定在p的值之间是否存在任何隐含的依赖关系（例如，gamma（p + 1）是否取决于gamma（p））。如果是这样，那么需要使用单独的循环。 - Eric Brown

如果伽马没有出现在右侧，那怎么可能呢？ - Rody Oldenhuis

@RodyOldenhuis - 通过 A 和 B 隐式地传递。（再次强调，我对 Matlab 或其符号一无所知。） - Eric Brown

Answer 2

可以使用 bsxfun 进行除法运算，使用矩阵乘法进行循环操作：

gamma = bsxfun(@rdivide, B.'*A, sum(A.^2));
x = B - A*gamma.';

这里有一个测试脚本

N = 3;
K = 4;
P = 5;

A = rand(N, K);
B = rand(N, P);

for p = 1:P
    for n = 1:N
        for k = 1:K
            gamma(p,k) = dot(A(:,k),B(:,p))/sum(A(:,k).^2);
        end
        x(n,p) = B(n,p) - dot(gamma(p,:),A(n,:));
    end
end

gamma2 = bsxfun(@rdivide, B.'*A, sum(A.^2));
X2 = B - A*gamma2.';

isequal(x, X2)
isequal(gamma, gamma2)

返回

ans =
     1
ans =
     1

- Mohsen Nosratinia

---

+1：绝对更简洁，而且速度也更快。甚至看起来有点像纯向量公式的样子 :) - Rody Oldenhuis

谢谢你的回答。我会尝试学习如何实现bsxfun！ - Langma

Answer 3

编程相关内容：

bxfun很慢... 可以尝试以下代码（可能有转置错误）：

modA = A * (1./sum(A.^2,2)) * ones(1,k);
gamma = B' * modA;
x = B - A * gamma';

- user3672668