我在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 functionthis is one on orthogonal functions都有提到。
繁殖
可以使用以下脚本生成结果:
A = complex( rand(360000,100), rand(360000,100));
w = rand(360000, 1);
[Q, R] = Gram_Schmidt(A, w);

其中Aw是输入。
速度与计算
如果使用上述脚本,将得到与以下内容同义的探查器结果:
测试结果
通过使用以下脚本将函数与上面的函数进行比较,可以测试结果:
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是可选函数。结果zeros1zeros2应该非常接近于零。
注:
我试着加快计算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上运行慢很多倍,但值得一试:)

10-05 20:16
查看更多