[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Balancing linear systems
From: |
Thomas Shores |
Subject: |
Re: Balancing linear systems |
Date: |
Sun, 26 Oct 2008 20:04:40 -0500 |
User-agent: |
Thunderbird 2.0.0.16 (X11/20080723) |
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
address@hidden
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