performance - 加权 Gram-Schmidt 正交化的 MATLAB 优化

Question

我在 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 是一个非负权重函数。维基百科的几个页面上提到了加权内积，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) ) ) );

Answer 1

1)

我第一次尝试向量化:

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)

然后我意识到这个中间矩阵QQ的列是增量构建的，之前的没有修改。所以这是我的第二次尝试:

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，方法是将该部分移动到外循环的末尾(即在计算新列 Q(:,j) 之后立即执行) ，我们计算并存储相应的QQ(:,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

请注意，我使用了 zeros(..., 'like',A) 语法(在最近的 MATLAB 版本中是新的)。这允许我们在 GPU 上运行未修改的函数(假设您有并行计算工具箱):

% CPU
[Q3,R3] = Gram_Schmidt3(A, w);

对比

% GPU
AA = gpuArray(A);
[Q3,R3] = Gram_Schmidt3(AA, w);

不幸的是，就我而言，它并没有更快。事实上，在 GPU 上运行比在 CPU 上运行要慢很多倍，但值得一试:)