matlab gram-schmidt,性能 – 加权Gram-Schmidt正交化的MATLAB优化

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

请注意，我使用了零(…，’like’，A)语法(最新的MATLAB版本)。这允许我们在GPU上运行未修改的功能(假设你有并行计算工具箱)：

% CPU

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

与

% GPU

AA = gpuArray(A);

[Q3,R3] = Gram_Schmidt3(AA, w);

不幸的是，在我的情况下，这不是更快。实际上，在GPU上运行的速度要比在CPU上慢很多倍，但是值得一提：