Balancing linear systems

[Thomas Shores]

Marco Caliari wrote:
> Dear all,
>
> I would like to discuss the following example. Consider the (badly scaled) 
> matrix A given by
>
> n = 10;
> B = rand(n);
> d = 10.^linspace(-14,4,n)';
> A = diag(d)*B;
>
> and the right hand side b
>
> b = d.*(B*ones(n,1));
>
> It is clear that the exact solution of Ax=b is ones(n,1). If you try
>
> x = A\b
>
> you get the following warning
>
> warning: matrix singular to machine precision, rcond = 8.1733e-20
> warning: attempting to find minimum norm solution
>
> (which is fine to me, the matrix A is really ill-conditioned). On the 
> other hand, if you balance your system
>
> d1 = max(A,[],2);
> A = diag(1./d1)*A;
> b = b./d1;
>
> then x = A\b gives you a good solution. Is there any drawback in 
> balancing, in any case, a linear system? I can't imagine a 
> counter-example.
>
> Best regards,
>
> Marco
> _______________________________________________
> Help-octave mailing list
> Help-octave at cae.wisc.edu
> https://www-old.cae.wisc.edu/mailman/listinfo/help-octave
>   
Scaling a coefficient matrix A by multiplication by a diagonal matrix, 
or more generally, balancing the matrix by multiplying on the left and 
right by a diagonal matrix,  are methods which can improve the condition 
number of A.  You gave an example of row scaling.   Of course one 
objective of scaling is to reduce the condition number of the 
coefficient matrix as much as possible.  Optimal scaling has been known 
for a long time (Bauer, 60s) but the method for obtaining it requires 
knowledge of the inverse of A, so it is impractical. Research on 
practical scaling methods continues to this day. The general advice of 
numerical analysts is that scaling should be applied on an ad hoc basis, 
with particular attention to what the source of the problem indicates 
about the coefficients in A.  You can find balancing routines in LAPACK, 
but I don't think that they are applied as the default in any case.

You can find an extensive discussion in Forsythe and Moler (1967), 
"Computer Solutions to Linear Systems".  One can actually make things 
worse by row scaling.  As far as examples go, it's a bit difficult to 
show in octave because the system solver will use condition number to 
detect near singularity and use least squares/pseudo-inverse routines 
when such is detected.  On paper, an example that illustrates this 
difficulty is

a = [eps/2 1 0; 1 1 4/eps; 0 0 1];
b = [1;2;0];

You can solve this system with paper and pencil symbolically and obtain 
an exact solution very near to [1;1;0].  To force Octave to use Gaussian 
elimination (with pivoting), solve the system by way of the PLU 
factorization and solve the resulting system by back solving:

[l,u,p]=lu(a);
x = u\(l\(p*b)) % since l*u = p*a and a*x = b, l lower, u upper and p a 
permutation matrix
x =

   1.00000
   1.00000
   0.00000

Pretty good.  If you look at the permuation matrix p, you'll see that 
the first and second rows of a were interchanged, as should be with 
partial pivoting.  Next, row scale the problem and compute the solution:

 A = diag(1./max(abs(a),[],2))*a;
 B = b./max(abs(a),[],2);
[L,U,P]=lu(A);
X = U\(L\(P*B))
X =

   2.00000
   1.00000
   0.00000

Pretty bad.  Now look at P and you see that there were no row 
interchanges.  The problem is that scaling caused a bad pivot element to 
be used in Gaussian elimination.  BTW, in spite of that, row scaling did 
actually reduce the condition number of the coefficient matrix:

 cond(a)
ans =  3.2452e+32
 cond(A)
ans =  3.6029e+16

Finally, just for amusement, try calculating the determinant of the 
upper triangular u.  Of course, we *know* what it should be, since the 
matrix is upper triangular (make sure it is via triu):

u
u =

   1.0000e+00   1.0000e+00   1.8014e+16
   0.0000e+00   1.0000e+00  -2.0000e+00
   0.0000e+00   0.0000e+00   1.0000e+00

det(triu(u))
ans = 0
det(diag(diag(u)))
ans =  1.0000


Clearly, octave (LAPACK) is not exploiting triangularity for upper 
triangular calculation.


Thomas Shores