我有一个稀疏的对称正定矩阵并且我希望某些行/列中的条目与其他条目中的条目具有非常不同的数量级(高达的因子)。
如果我要使用稀疏 Cholesky 分解来求解涉及的线性系统,用更好缩放的版本是否有任何数值好处?我知道对于像 LU 这样的其他分解,缩放非常重要,但 Cholesky 也受益吗?
我有一个稀疏的对称正定矩阵并且我希望某些行/列中的条目与其他条目中的条目具有非常不同的数量级(高达的因子)。
如果我要使用稀疏 Cholesky 分解来求解涉及的线性系统,用更好缩放的版本是否有任何数值好处?我知道对于像 LU 这样的其他分解,缩放非常重要,但 Cholesky 也受益吗?
缩放可以影响矩阵的条件数,并且对于您可能对矩阵执行的某些操作可能很重要。但是,如果您只是求解方程组,则右侧的缩放/取消缩放最终会抵消 Cholesky 因子的缩放,因此缩放矩阵的 Cholesky 因子几乎是(最多舍入误差)等于原始矩阵的缩放 Cholesky 因子。
扩展答案如下:
一个很好的参考是第 7.3 节
Higham, Nicholas J. 数值算法的准确性和稳定性。暹罗,2002 年。
Higham 指的是 1960 年代 van der Sluis 的一篇论文,其中给出了最佳缩放的结果。的对角线的平方根的倒数进行缩放(因此得到的矩阵具有全为 1 的对角线)几乎是最优的。
我写了一个 MATLAB 脚本(在这个答案的底部)来演示这一点。
的随机对称正定矩阵,然后给出错误的对角线缩放,导致条件数为,然后使用范德将对角线缩放为 1 Sluis 缩放,产生一个条件数为的矩阵。这表明缩放矩阵会影响条件数。
脚本的输出是:
Condition number of original A is 9.895810e+03
Condition number of badly scaled A is 2.307700e+18
Condition number of well scaled A is 9.834918e+03
MATLAB 脚本是:
%
% Reset the RNG's.
%
rand('seed',0);
randn('seed',0);
%
% Basic parameters for the test that can be adjusted.
%
%
% n, the size of the matrix.
%
n=1000;
%
% logcondnum, the log10 of the original condition number of the matrix.
%
logcondnum=4;
%
% range of bad scaling factors. from 10^0 to 10^scalingrange
%
scalingrange=8;
%
% Generate the random matrix.
%
M=randn(n,n);
[Q,R]=qr(M);
%
% For the eigenvalues we'll use random values between 10^0 and 10^logcondnum
%
lambdas=10.^(logcondnum*rand(n,1));
%
% Now, construct A using the eigenvalues and Q.
%
A=Q*diag(lambdas)*Q';
%
% Make A perfectly symmetric.
%
A=(A+A')/2;
%
% Now, construct a bad scaling of A.
%
d1=10.^(scalingrange*rand(n,1));
D1=diag(d1);
AS1=D1*A*D1;
%
% Make it symmetric.
%
AS1=(AS1+AS1')/2;
%
% Now use the van der Sluis diagonal scaling to get AS2 from AS1.
%
d2=1./sqrt(diag(AS1));
D2=diag(d2);
AS2=D2*AS1*D2;
%
% Make it symmetric.
%
AS2=(AS2+AS2')/2;
%
% Output information about the condition numbers.
%
fprintf('Condition number of original A is %e\n',cond(A));
fprintf('Condition number of badly scaled A is %e\n',cond(AS1));
fprintf('Condition number of well scaled A is %e\n',cond(AS2));
fprintf('\n');