我在Matlab中有一个函数,它执行Gram-Schmidt Orthogonalisation并对内部产品应用非常重要的权重(我不认为Matlab的内置函数支持这一点)。
就我所知,这个函数工作得很好,但是在大矩阵上太慢了。
什么是最好的改善方法?
我试图转换成一个MEX文件,但我失去了与我正在使用的编译器的并行化,所以它是那么慢。
我想在GPU上运行它,因为元素乘法是高度并行的(但我更希望实现容易移植)
有人能把这段代码矢量化或使它更快吗我不知道该怎么做优雅。。。
我知道StackOverflow的头脑是惊人的,把这看作一个挑战:)
功能
function [Q, R] = Gram_Schmidt(A, w)
[m, n] = size(A);
Q = complex(zeros(m, n));
R = complex(zeros(n, n));
v = zeros(n, 1);
for j = 1:n
v = A(:,j);
for i = 1:j-1
R(i,j) = sum( v .* conj( Q(:,i) ) .* w ) / ...
sum( Q(:,i) .* conj( Q(:,i) ) .* w );
v = v - R(i,j) * Q(:,i);
end
R(j,j) = norm(v);
Q(:,j) = v / R(j,j);
end
end
其中
A
是复数的m x n
矩阵,w
是实数的m x 1
向量。瓶颈
这是
R(i,j)
的表达式,它是函数最慢的部分(如果符号正确,则不能100%确定):其中
w
是一个非负权重函数。加权内积在Wikipedia的几个页面上,this is one on the weight function和this is one on orthogonal functions都有提到。
繁殖
可以使用以下脚本生成结果:
A = complex( rand(360000,100), rand(360000,100));
w = rand(360000, 1);
[Q, R] = Gram_Schmidt(A, w);
其中
A
和w
是输入。速度与计算
如果使用上述脚本,将得到与以下内容同义的探查器结果:
测试结果
通过使用以下脚本将函数与上面的函数进行比较,可以测试结果:
A = complex( rand( 100, 10), rand( 100, 10));
w = rand( 100, 1);
[Q , R ] = Gram_Schmidt( A, w);
[Q2, R2] = Gram_Schmidt2( A, w);
zeros1 = norm( Q - Q2 );
zeros2 = norm( R - R2 );
其中,
Gram_Schmidt
是前面描述的函数,Gram_Schmidt2
是可选函数。结果zeros1
和zeros2
应该非常接近于零。注:
我试着加快计算
R(i,j)
的速度,但没有成功。。。R(i,j) = ( w' * ( v .* conj( Q(:,i) ) ) ) / ...
( w' * ( Q(:,i) .* conj( Q(:,i) ) ) );
最佳答案
(一)
我第一次尝试矢量化:
function [Q, R] = Gram_Schmidt1(A, w)
[m, n] = size(A);
Q = complex(zeros(m, n));
R = complex(zeros(n, n));
for j = 1:n
v = A(:,j);
QQ = Q(:,1:j-1);
QQ = bsxfun(@rdivide, bsxfun(@times, w, conj(QQ)), w.' * abs(QQ).^2);
for i = 1:j-1
R(i,j) = (v.' * QQ(:,i));
v = v - R(i,j) * Q(:,i);
end
R(j,j) = norm(v);
Q(:,j) = v / R(j,j);
end
end
不幸的是,它比原来的函数慢。
2个)
然后我意识到这个中间矩阵的列是增量构建的,而之前的列是不修改的这是我的第二次尝试:
function [Q, R] = Gram_Schmidt2(A, w)
[m, n] = size(A);
Q = complex(zeros(m, n));
R = complex(zeros(n, n));
QQ = complex(zeros(m, n-1));
for j = 1:n
if j>1
qj = Q(:,j-1);
QQ(:,j-1) = (conj(qj) .* w) ./ (w.' * (qj.*conj(qj)));
end
v = A(:,j);
for i = 1:j-1
R(i,j) = (v.' * QQ(:,i));
v = v - R(i,j) * Q(:,i);
end
R(j,j) = norm(v);
Q(:,j) = v / R(j,j);
end
end
从技术上讲,没有进行主要的矢量化;我只预先计算了中间结果,并将计算移到了内部循环之外。
基于快速基准测试,这个新版本绝对更快:
% some random data
>> M = 10000; N = 100;
>> A = complex(rand(M,N), rand(M,N));
>> w = rand(M,1);
% time
>> timeit(@() Gram_Schmidt(A,w), 2) % original version
ans =
1.2444
>> timeit(@() Gram_Schmidt1(A,w), 2) % first attempt (vectorized)
ans =
2.0990
>> timeit(@() Gram_Schmidt2(A,w), 2) % final version
ans =
0.4698
% check results
>> [Q,R] = Gram_Schmidt(A,w);
>> [Q2,R2] = Gram_Schmidt2(A,w);
>> norm(Q-Q2)
ans =
4.2796e-14
>> norm(R-R2)
ans =
1.7782e-12
编辑:
在注释之后,我们可以重写第二个解决方案以去掉if statmenet,方法是将该部分移动到外部循环的末尾(即在计算新列
QQ
之后,我们计算并存储相应的Q(:,j)
)。这个函数在输出上是相同的,时间也没有那么不同;代码只是短了一点!
function [Q, R] = Gram_Schmidt3(A, w)
[m, n] = size(A);
Q = zeros(m, n, 'like',A);
R = zeros(n, n, 'like',A);
QQ = zeros(m, n, 'like',A);
for j = 1:n
v = A(:,j);
for i = 1:j-1
R(i,j) = (v.' * QQ(:,i));
v = v - R(i,j) * Q(:,i);
end
R(j,j) = norm(v);
Q(:,j) = v / R(j,j);
QQ(:,j) = (conj(Q(:,j)) .* w) ./ (w.' * (Q(:,j).*conj(Q(:,j))));
end
end
注意,我使用了
QQ(:,j)
语法(在最新的MATLAB版本中是新的)这允许我们在GPU上运行未经修改的函数(假设您有并行计算工具箱):% CPU
[Q3,R3] = Gram_Schmidt3(A, w);
与。
% GPU
AA = gpuArray(A);
[Q3,R3] = Gram_Schmidt3(AA, w);
不幸的是,在我的情况下,没有任何更快。实际上,在GPU上运行要比在CPU上运行慢很多倍,但值得一试:)