3
$\begingroup$

I have a sparse symmetric positive-definite matrix $M$ and I expect the entries in some rows/columns to have very different orders of magnitude (up to a factor of $10^8$) than the entries in others.

If I am going to solve a linear system involving $M$ using sparse Cholesky decomposition, is there any numerical benefit to replacing $M$ with a better-scaled version $M' = DMD^T$? I know for other decompositions like LU, scaling is very important, but does Cholesky benefit as well?

$\endgroup$
1
$\begingroup$

Scaling can affect the condition number of the matrix, and for some things that you might do with the matrix that could be important. However, if you're only solving a system of equations, the scaling/unscaling of the right hand side ends up canceling out the scaling of the Cholesky factor so that the Cholesky factor of the scaled matrix is nearly (up to round off errors) equal to the scaled Cholesky factor of the original matrix.

Expanded answer below:

A good reference for this is section 7.3 of

Higham, Nicholas J. Accuracy and Stability of Numerical Algorithms. SIAM, 2002.

Higham refers to a paper from the 1960's by van der Sluis in which results on optimal scaling are given. For a symmetric positive definite matrix, scaling by the reciprocal of the square root of the diagonal of $A$ (so the resulting matrix has a diagonal of all ones) is nearly optimal.

I wrote a MATLAB script (at the bottom of this answer) that demonstrates this.

The script generates a random symmetric and positive definite matrix with a condition number of $10^{4}$, then gives a bad diagonal scaling resulting in a condition number of $10^{18}$, and then scales the diagonal to 1 using the van der Sluis scaling, resulting in a matrix with a condition number of $10^{4}$. This demonstrates that scaling the matrix affects the condition number.

The output from the script was:

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

The MATLAB script is:

%
% 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');
$\endgroup$
9
  • 1
    $\begingroup$ But what if you do use pivoting (as one probably should)? $\endgroup$ – Wolfgang Bangerth Jun 4 '20 at 20:19
  • 1
    $\begingroup$ Ah, very interesting observation. Indeed, we don't do pivoting in sparse-land either. $\endgroup$ – Wolfgang Bangerth Jun 4 '20 at 23:23
  • 1
    $\begingroup$ I am not sure I agree with this answer. If you don't pivot dynamically for stability, then scaling is basically irrelevant: instead of getting the factorization $M = LL^*$ you get $DMD^* = (DL)(DL)^*$, but the two matrices called $L$ are identical in exact arithmetic, and the amount of cancellation that you get along the algorithm is the same because you make the exact same subtractions, apart from a multiplicative factor. $\endgroup$ – Federico Poloni Jun 5 '20 at 14:10
  • 1
    $\begingroup$ A bit of computational experimentation shows that in double precision with a badly scaled matrix, the Cholesky factor of the scaled matrix is close to but not exactly equal to the scaled Cholesky factor of the original matrix. As long as the unscaled problem is not so badly conditioned that Cholesky factorization fails, this ill-conditioning of the unscaled matrix might not be a problem for solving systems of equations using the Cholesky factorization. $\endgroup$ – Brian Borchers Jun 8 '20 at 2:58
  • 1
    $\begingroup$ I'll add that I've dealt with situations (in a primal-dual interior point code) where Cholesky factorization of the unscaled matrix fails while Cholesky factorization of the scaled matrix succeeds. This is no surprise, since the system matrix becomes extremely ill-conditioned as the solution approaches optimality. $\endgroup$ – Brian Borchers Jun 8 '20 at 20:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.