Re: Balancing linear systems

[Jaroslav Hajek]

On Mon, Oct 27, 2008 at 2:04 AM, Thomas Shores <address@hidden> wrote:
> 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.
>

Hello,

Just FYI, I just checked in a patch
(http://hg.savannah.gnu.org/hgweb/octave/rev/9813c07ca946) that makes
det exploit matrix_type, i.e. it will multiply diagonal elements for
triangular matrices / use Cholesky decomposition for SPD/HPD matrices.

enjoy,

-- 
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz