加权Gram-Schmidt正交化的MATLAB优化
Jac*_*obD 19 optimization performance matlab linear-algebra vectorization
我在MATLAB中有一个函数,它执行Gram-Schmidt正交化,并对内积进行了非常重要的加权(我不认为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是非负权重函数.加权内积在几个维基百科页面上提到,这是关于权重函数的一个,这是关于正交函数的一个.
再生产
您可以使用以下脚本生成结果:
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) ) ) );
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上慢很多倍,但它值得一试:)
- @John:这正是我答案中的"诡计"; 注意内循环来自`i = 1:j-1`并使用`Q(:,i)`,然后外循环修改当前`j`的`Q(:,j)`.这意味着在内部循环内部,`Q(:,1:j-1)`已经从外部循环的先前迭代计算出来,因此我们可以将该表达式取在内循环之外.我们仍然有内部循环中涉及的向量"v",所以我只是拿到了提名者的右侧,然后用分母预先乘以"QQ",而不是"sum(x.*) y)`我使用内积'x'*y`.说得通？ (2认同)
小智 7
这里有一个很长的讨论,但是,跳到答案.您已通过向量w对R计算的分子和分母进行加权.加权发生在内环上,由三点积,分子中的点Q点w和分母中的Q点Q点w组成.如果你做了一个改变,我认为代码将运行得更快.写入num =(A dot sqrt(w))点(Q dot sqrt(w))并写入den =(Q dot sqrt(w))dot(Q dot sqrt(w)).这将(A点sqrt(w))和(Q dot sqrt(w))乘积计算移出内循环.
我想对Gram Schmidt Orthogonalization的描述进行描述,希望除了提供替代计算解决方案外,还可以进一步了解GSO的优势.
GSO的"目标"是双重的.首先,要启用像Ax = y这样的等式的解,其中A具有比列多得多的行.在测量数据时经常发生这种情况,因为很容易测量比状态数量更多的数据.第一个目标的方法是将A重写为QR,使Q的列正交并归一化,R是三角矩阵.我相信,您提供的算法实现了第一个目标.Q表示A矩阵的基本空间,R表示生成A的每列所需的每个基础空间的幅度.
GSO的第二个目标是按重要性顺序对基础向量进行排序.这是你没有做过的一步.并且,在包括此步骤的同时,可以增加求解时间,结果将根据由A表示的测量中包含的数据来识别x的哪些元素是重要的.
但是,我认为,通过这种实现,解决方案比您提出的方法更快.
Aij = Qij Rij,其中Qj是正交的,Rij是上三角形,Ri,j> i = 0.Qj是A的正交基矢量,Rij是每个Qj参与在A中创建一列.所以,
A_j1 = Q_j1 * R_1,1
A_j2 = Q_j1 * R_1,2 + Q_j2 * R_2,2
A_j3 = Q_j1 * R_1,3 + Q_j2 * R_2,3 + Q_j3 * R_3,3
通过检查,你可以写
A_j1 = ( A_j1 / | A_j1 | ) * | A_j1 | = Q_j1 * R_1,1
然后将Q_j1投射到每个其他列A上以获取R_1,j元素
R_1,2 = Q_j1 dot Aj2
R_1,3 = Q_j1 dot Aj3
...
R_1,j(j>1) = A_j dot Q_j1
然后从A列中减去Q_j1的项目元素(这会将第一列设置为零,因此您可以忽略第一列
for j = 2,n
A_j = A_j - R_1,j * Q_j1
end
现在已经移除了A中的一列,确定了第一个标准正交基矢量Q,j1,并确定了第一个基矢量对每个列R_1,j的贡献,并且第一个基矢量的贡献已经从每列中减去.对剩余的A列重复此过程以获得剩余的Q列和R行.
for i = 1,n
R_ii = |A_i| A_i is the ith column of A, |A_i| is magnitude of A_i
Q_i = A_i / R_ii Q_i is the ith column of Q
for j = i, n
R_ij = | A_j dot Q_i |
A_j = A_j - R_ij * Q_i
end
end
你试图用w来加权A的行.这是一种方法.我会将w归一化,并将效果合并到R.你通过乘以w除去"w"的效果."去除"效果的替代方案是将w的幅度归一化为1.
w = w / | w |
for i = 1,n
R_ii = |A_i inner product w| # A_i inner product w = A_i .* w
Q_i = A_i / R_ii
for j = i, n
R_ij = | (A_i inner product w) dot Q_i | # A dot B = A' * B
A_j = A_j - R_ij * Q_i
end
end
实现w的另一种方法是对w进行归一化,然后将每个A列乘以w进行预乘.这样可以对A行进行干净的加权,并减少乘法次数.
使用以下内容可能有助于加快代码速度
A inner product B = A .* B
A dot w = A' w
(A B)' = B'A'
A' conj(A) = |A|^2
上面的内容可以很容易地在matlab中进行矢量化,就像写的一样.
但是,你缺少A的排名的第二部分,它告诉你哪些状态(A x = y中x的元素)在数据中有显着的表示
排名过程很容易描述,但我会让你弄清楚编程细节.上述过程基本上假设A的列按重要性顺序排列,并且从剩余的所有列中减去第一列,然后从剩余的列中减去第二列等.R的第一行代表A的贡献. Q的第一列到A的每一列.如果你对第一行R贡献的绝对值求和,你得到Q的第一列对矩阵A的贡献的测量.所以,你只需要评估每列的A作为Q的第一列(或下一列),并确定该Q列对A的剩余列的贡献的排名分数.然后选择具有最高等级的A列作为下一个Q列.对此进行编码基本上归结为预先估计A的每个剩余列的下一行R,以便确定哪个排序的R幅度具有最大幅度.具有表示A的原始列顺序的索引向量将是有益的.通过对基础向量进行排序,您最终得到代表A的"主要"基础向量,其通常比A中的列数小得多.
此外,如果对列进行排名,则无需计算R的每一列.当您知道A的哪些列不包含任何有用信息时,保留这些列没有任何实际好处.
在结构动力学中,减少自由度数的一种方法是计算特征值,假设您具有质量和刚度矩阵的代表值.如果您考虑一下,上述方法可用于从测量的响应中"计算"M和K(和C)矩阵,并且还可以识别在数据中显着表示的"测量响应形状".这些比模式形状更难以且可能更重要.因此,您可以通过上述方法解决非常困难的问题,即从测量的响应中估计状态矩阵和所表示的自由度数.如果您阅读了N4SID,他会做类似的事情,除了他使用SVD而不是GSO.我不喜欢N4SID的技术描述,过分关注矢量投影符号,这只是一个点积.
在上述信息中可能有一两个错误,我正在写下这个错误,然后急忙上班.所以,在你实施时检查算法/方程式...祝你好运
回到你的问题,当你用w加权时如何优化算法.这是一个基本的GSO算法,没有排序,与您的函数兼容.
注意,下面的代码是八度,而不是matlab.有一些细微的差别.
function [Q, R] = Gram_Schmidt_2(A, w)
[m, n] = size(A);
Q = complex(zeros(m, n));
R = complex(zeros(n, n));
# Outer loop identifies the basis vectors
for j = 1:n
aCol = A(:,j);
# Subtract off the basis vector
for i = 1:(j-1)
R(i,j) = ctranspose(Q(:,j)) * aCol;
aCol = aCol - R(i,j) * Q(:,j);
end
amp_A_col = norm(aCol);
R(j,j) = amp_A_col;
Q(:,j) = aCol / amp_A_col;
end
end
要获得算法,只需更改一行.但是,你失去了很多速度因为"ctranspose(Q(:,j))*aCol"是一个向量运算但是"sum(aCol.*conj(Q(:,i)).*w)"是一行操作.
function [Q, R] = Gram_Schmidt_2(A, w)
[m, n] = size(A);
Q = complex(zeros(m, n));
R = complex(zeros(n, n));
# Outer loop identifies the basis vectors
for j = 1:n
aCol = A(:,j);
# Subtract off the basis vector
for i = 1:(j-1)
# R(i,j) = ctranspose(Q(:,j)) * aCol;
R(i,j) = sum( aCol .* conj( Q(:,i) ) .* w ) / ...
sum( Q(:,i) .* conj( Q(:,i) ) .* w );
aCol = aCol - R(i,j) * Q(:,j);
end
amp_A_col = norm(aCol);
R(j,j) = amp_A_col;
Q(:,j) = aCol / amp_A_col;
end
end
您可以通过用w的sqrt加权aCol和Q来将其更改回向量运算.
function [Q, R] = Gram_Schmidt_3(A, w)
[m, n] = size(A);
Q = complex(zeros(m, n));
R = complex(zeros(n, n));
Q_sw = complex(zeros(m, n));
sw = w .^ 0.5;
for j = 1:n
aCol = A(:,j);
aCol_sw = aCol .* sw;
# Subtract off the basis vector
for i = 1:(j-1)
# R(i,j) = ctranspose(Q(:,i)) * aCol;
numTerm = ctranspose( Q_sw(:,i) ) * aCol_sw;
denTerm = ctranspose( Q_sw(:,i) ) * Q_sw(:,i);
R(i,j) = numTerm / denTerm;
aCol_sw = aCol_sw - R(i,j) * Q_sw(:,i);
end
aCol = aCol_sw ./ sw;
amp_A_col = norm(aCol);
R(j,j) = amp_A_col;
Q(:,j) = aCol / amp_A_col;
Q_sw(:,j) = Q(:,j) .* sw;
end
end
正如JacobD所指出的,上述功能运行速度不快.可能需要时间来创建其他阵列.三重产品的另一种分组策略是使用conj(Q)对w进行分组.希望这更快......
function [Q, R] = Gram_Schmidt_4(A, w)
[m, n] = size(A);
Q = complex(zeros(m, n));
R = complex(zeros(n, n));
for j = 1:n
aCol = A(:,j);
for i = 1:(j-1)
cqw = conj(Q(:,i)) .* w;
R(i,j) = ( transpose( aCol ) * cqw) ...
/ (transpose( Q(:,i) ) * cqw);
aCol = aCol - R(i,j) * Q(:,i);
end
amp_A_col = norm(aCol);
R(j,j) = amp_A_col;
Q(:,j) = aCol / amp_A_col;
end
end
这是一个驱动程序函数来计时不同的版本.
function Gram_Schmidt_tester_2
nSamples = 360000;
nMeas = 100;
nMeas = 15;
A = complex( rand(nSamples,nMeas), rand(nSamples,nMeas));
w = rand(nSamples, 1);
profile on;
[Q1, R1] = Gram_Schmidt_basic(A);
profile off;
data1 = profile ("info");
tData1=data1.FunctionTable(1).TotalTime;
approx_zero1 = A - Q1 * R1;
max_value1 = max(max(abs(approx_zero1)));
profile on;
[Q2, R2] = Gram_Schmidt_w_Orig(A, w);
profile off;
data2 = profile ("info");
tData2=data2.FunctionTable(1).TotalTime;
approx_zero2 = A - Q2 * R2;
max_value2 = max(max(abs(approx_zero2)));
sw=w.^0.5;
profile on;
[Q3, R3] = Gram_Schmidt_sqrt_w(A, w);
profile off;
data3 = profile ("info");
tData3=data3.FunctionTable(1).TotalTime;
approx_zero3 = A - Q3 * R3;
max_value3 = max(max(abs(approx_zero3)));
profile on;
[Q4, R4] = Gram_Schmidt_4(A, w);
profile off;
data4 = profile ("info");
tData4=data4.FunctionTable(1).TotalTime;
approx_zero4 = A - Q4 * R4;
max_value4 = max(max(abs(approx_zero4)));
profile on;
[Q5, R5] = Gram_Schmidt_5(A, w);
profile off;
data5 = profile ("info");
tData5=data5.FunctionTable(1).TotalTime;
approx_zero5 = A - Q5 * R5;
max_value5 = max(max(abs(approx_zero5)));
profile on;
[Q2a, R2a] = Gram_Schmidt2a(A, w);
profile off;
data2a = profile ("info");
tData2a=data2a.FunctionTable(1).TotalTime;
approx_zero2a = A - Q2a * R2a;
max_value2a = max(max(abs(approx_zero2a)));
profshow (data1, 6);
profshow (data2, 6);
profshow (data3, 6);
profshow (data4, 6);
profshow (data5, 6);
profshow (data2a, 6);
sprintf('Time for %s is %5.3f sec with %d samples and %d meas, max value is %g',
data1.FunctionTable(1).FunctionName,
data1.FunctionTable(1).TotalTime,
nSamples, nMeas, max_value1)
sprintf('Time for %s is %5.3f sec with %d samples and %d meas, max value is %g',
data2.FunctionTable(1).FunctionName,
data2.FunctionTable(1).TotalTime,
nSamples, nMeas, max_value2)
sprintf('Time for %s is %5.3f sec with %d samples and %d meas, max value is %g',
data3.FunctionTable(1).FunctionName,
data3.FunctionTable(1).TotalTime,
nSamples, nMeas, max_value3)
sprintf('Time for %s is %5.3f sec with %d samples and %d meas, max value is %g',
data4.FunctionTable(1).FunctionName,
data4.FunctionTable(1).TotalTime,
nSamples, nMeas, max_value4)
sprintf('Time for %s is %5.3f sec with %d samples and %d meas, max value is %g',
data5.FunctionTable(1).FunctionName,
data5.FunctionTable(1).TotalTime,
nSamples, nMeas, max_value5)
sprintf('Time for %s is %5.3f sec with %d samples and %d meas, max value is %g',
data2a.FunctionTable(1).FunctionName,
data2a.FunctionTable(1).TotalTime,
nSamples, nMeas, max_value2a)
end
在我的旧家用笔记本电脑上,在Octave,结果是
ans = Time for Gram_Schmidt_basic is 0.889 sec with 360000 samples and 15 meas, max value is 1.57009e-16
ans = Time for Gram_Schmidt_w_Orig is 0.952 sec with 360000 samples and 15 meas, max value is 6.36717e-16
ans = Time for Gram_Schmidt_sqrt_w is 0.390 sec with 360000 samples and 15 meas, max value is 6.47366e-16
ans = Time for Gram_Schmidt_4 is 0.452 sec with 360000 samples and 15 meas, max value is 6.47366e-16
ans = Time for Gram_Schmidt_5 is 2.636 sec with 360000 samples and 15 meas, max value is 6.47366e-16
ans = Time for Gram_Schmidt2a is 0.905 sec with 360000 samples and 15 meas, max value is 6.68443e-16
这些结果表明,最快的算法是上面的sqrt_w算法在0.39秒,接着是在0.452秒处将conj(Q)与w(上图)分组,然后在0.905秒处对Amro解的版本2进行分组,然后是问题中的原始算法在0.952,然后是版本5,它交换行/列以查看是否在2.636秒显示行存储(代码未包括).这些结果表明A和Q之间的sqrt(w)分裂是最快的解决方案.但这些结果与JacobD关于sqrt(w)不快的评论不一致.