🚀 go-pugleaf

In <1993Jan11.083253.9521@dutrun2.tudelft.nl> rcpshdb@dutrun2.tudelft.nl (Han de Bruijn) writes:

>The condition number of a matrix  A  is defined as the absolute value of
>its largest eigenvalue divided by its smallest eigenvalue. Right?

>I'm at loss how this is related to the loss of decimals when solving the
>equations system:
>                      A.x = b

>Can someone out there provide a simple argument?


A perturbation "db" in the right hand side (due to rounding errors for
example) results in a change "dx" of the solution vector.  
They are related by:

                |dx| / |x|  <=  cond(A)  |db| / |b|

The proof I am looking at right now assumes that A is symmetric...
-----------------------------------------------------------------------------
Fred Walsteijn                                | Internet: walsteyn@fys.ruu.nl
Institute for Marine and Atmospheric Research | FAX:      31-30-543163
Utrecht University, The Netherlands           | Phone:    31-30-533169

In article <1993Jan11.083253.9521@dutrun2.tudelft.nl> rcpshdb@dutrun2.tudelft.nl (Han de Bruijn) writes:

>The condition number of a matrix  A  is defined as the absolute value of
>its largest eigenvalue divided by its smallest eigenvalue. Right?
>
>I'm at loss how this is related to the loss of decimals when solving the
>equations system:
>                      A.x = b
>
>Can someone out there provide a simple argument?

You did say simple right?  Consider that a singular matrix has smallest
magnitude eigenvalue of zero.  The closer to singular in some sense.
Also, dividing by a small number gives a larger condition number.

andy
pitonyak@cis.ohio-state.edu

In article <1993Jan11.083253.9521@dutrun2.tudelft.nl> rcpshdb@dutrun2.tudelft.nl (Han de Bruijn) writes:
>The condition number of a matrix  A  is defined as the absolute value of
>its largest eigenvalue divided by its smallest eigenvalue. Right?
>
>I'm at loss how this is related to the loss of decimals when solving the
>equations system:
>                      A.x = b
>
>Can someone out there provide a simple argument?

Backward error analysis of Gaussian Elimination shows that the error
introduced during the factorization because of finite precision floating
point arithmetic can be modelled as by making small normwise relative
perturbations to the matrix and the right-hand side, where "small" means
a few times the computer's machine epsilon.

Suppose you have a non-singular linear system A.x = b and you perturb the
matrix and right-hand side.  How much is the solution perturbed; i.e.,
how big is dx if (A+dA).(x+dx) = b+db ?

	(A + dA).(x + dx) = b + db

	(A + dA).dx =  b + db - A.x - dA.x

		    = b + db - b - dA.x

		    = db - dA.x

	dx = (A + dA)^-1.(db - dA.x)

	   =A^-1.(I+A^-1.dA)^-1.(db - dA.x)

	||dx|| <= ||A^-1.(I+A^-1.dA)^-1||.||(db - dA.x)||


Assume ||A^-1||.||dA|| < 1.  Then (I+A^-1.dA)^-1 exists and

	||(I+A^-1.dA)^-1|| <= (1-||A^-1||.||dA||)^-1.

Hence
	||dx|| <= ||A^-1||.||(I+A^-1.dA)^-1||.||(db - dA.x)||

	       	      ||A^-1||
	       <= ----------------- . (||dA||.||x|| + ||db||)
		  1-||A^-1||.||dA||

	||dx||        ||A^-1||        ||dA||.||x|| + ||db||
	------  = ----------------- . ---------------------
	||x||	  1-||A^-1||.||dA||          ||x||

	||dx||        ||A^-1||                  ||db||
	------  = ----------------- . (||dA|| + ------)
	||x||	  1-||A^-1||.||dA||             ||x||

	||dx||      ||A||.||A^-1||     ||dA||      ||db||
	------  = ----------------- . (------ + -----------)
	||x||	  1-||A^-1||.||dA||    ||A||	||A||.||x||

	||dx||      ||A||.||A^-1||     ||dA||   ||db||
	------ <= ----------------- . (------ + -------)
	||x||	  1-||A^-1||.||dA||    ||A||	||A.x||

	||dx||      ||A||.||A^-1||     ||dA||   ||db||
	------ <= ----------------- . (------ + ------)
	||x||	  1-||A^-1||.||dA||    ||A||	||b||

Define kappa(A) = ||A||.||A^-1||.  This says that if dA is not too big, the
normwise relative perturbation in the solution is bounded by the condition
number times the sum of normwise relative perturbations of the matrix and
right-hand side.

Finally, ||A||_2 is given by the largest magnitude of the eigenvalues of A
and ||A^-1||_2 is given by the smallest magnitude of the eigenvalues of A,
if A is a normal matrix, but not in general.

----------------------------------------------------------------------

Dave Dodson		                             dodson@convex.COM
Convex Computer Corporation      Richardson, Texas      (214) 497-4234

In article <1993Jan11.083253.9521@dutrun2.tudelft.nl> rcpshdb@dutrun2.tudelft.nl (Han de Bruijn) writes:
>The condition number of a matrix  A  is defined as the absolute value of
>its largest eigenvalue divided by its smallest eigenvalue. Right?
>
>I'm at loss how this is related to the loss of decimals when solving the
>equations system:
>                      A.x = b
>
>Can someone out there provide a simple argument?

David Dodson has provided a nice derivation of the relevant
perturbation theory.  However, some additional comments are needed to
get to loss of decimals.  (I'll use LaTeX notation here, with
apologies to those who don't read it.)

If we replace $A$ by $A+E$, we find that the change $h$ in the
solution satisfies (up to second order terms in $E$)
\[
\|h\|/\|X\| \leq \kappa(A)\|E\|/\|A\|,
\]
where
\[
\kappa(A) = \|A\| \|A^{-1}\|
\]
is the condition number of $A$.  Now if $\|E\|/\|A\| = 10^{-t}$, then
the elements of $A$ will be perturbed in their $t$th decimal digits
(some qualifications about balance in the sizes of the elements of $E$
and $A$ are needed here).  Hence if $\kappa(A) = 10^{k}$, then
$\|h\|/\|X\| \leq 10^{k-t}$, so that the components of $x$ are
perturbed in their $(t-k)$th digits.  In ordinary English, if $A$ is
known to $t$ digits and $\kappa(A) = 10^{k}$, then $x$ will be known
to $t-k$ digits.

In particular, if a stable algorithm like Gaussian elimination is used
to solve the system, then $t$ is the number of digits used in the
computation.  Thus, if you try to solve a system with condition
number $10^{k}$ using Gaussian elimination, you can expect a loss
of about $k$ digits in the solution.

All this is a useful oversimplification.  For the complete story, see
an introductory text on matrix computations.

Pete Stewart

It is actually even more complex than just looking at the eignvalues.  ONe must
also examine whether the eigenvectors of the ill conditioning are being excited. 
It is a trivial example, but the largest positive definite ill conditioning on a
VAX in single precision is approximately 10E76.  The following matirx is the
worst case conditioning, but is highly stable during Gauss elimination:
 
                                 1.e38         0
              A         =                      
                                   0          1.e-38

This effect is particularly noticable in thin shell problems using dinite element
analysis, since the out of plane eigenvectors are not excited by in plane
bending, evne though the ill conditioning is very high (frequently on the order
of 10e9).
-- 

Regards.
 
Don Rolph  a722756@pan.mc.ti.com WD3 MS10-13 (508)-699-1263

In article <1993Jan11.214251.12368@news.eng.convex.com>
dodson@convex.COM (Dave Dodson) writes:
>In article <1993Jan11.083253.9521@dutrun2.tudelft.nl>
>rcpshdb@dutrun2.tudelft.nl (Han de Bruijn) writes:
>>The condition number of a matrix  A  is defined as the absolute value of
>>its largest eigenvalue divided by its smallest eigenvalue. Right?
>>
>>I'm at loss how this is related to the loss of decimals when solving the
>>equations system:
>>                      A.x = b
>>
>>Can someone out there provide a simple argument?

  [Dave Dodson's good explaination deleted]

One thing that no one has remarked on: The original poster's reference
to *eigenvalues*.

    As Dave Dodson points out, the condition number of A is
norm(A) * norm(A-inverse).  No mention of eigenvalues.

    In fact, the condition number of a matrix has (almost) nothing to
do with the eigenvalues, as the example

        ( 1  t )                        -1   ( 1  -t )
    A = (      )   for t >> 1 shows:   A   = (       )
        ( 0  1 )                             ( 0   1 )

so the eigenvalues are both 1, but norm(A) = norm(A-inverse)
which are approximately  t, so the condition number is approx.
t**2.

    What the condition number *is* related to are the *singular values*.
For more info, look in a good book on numerical linear algebra
(e.g., Golub & van Loan's _Matrix_Computations_) on the "singular
value decomposition" (or SVD).  The SVD, once understood, makes
a lot of this a whole lot simpler to understand.


-- 
Alan McKenney        E-mail:  mckenney@cims.nyu.edu         (INTERNET)
Courant Institute,NYU,USA     ...!cmcl2!cims.nyu.edu!mckenney   (UUCP)

mckenney@cims.nyu.edu (Alan McKenney) writes:


>In article <1993Jan11.214251.12368@news.eng.convex.com>
>dodson@convex.COM (Dave Dodson) writes:
>>In article <1993Jan11.083253.9521@dutrun2.tudelft.nl>
>>rcpshdb@dutrun2.tudelft.nl (Han de Bruijn) writes:
>>>The condition number of a matrix  A  is defined as the absolute value of
>>>its largest eigenvalue divided by its smallest eigenvalue. Right?
>>>
>>>I'm at loss how this is related to the loss of decimals when solving the
>>>equations system:
>>>                      A.x = b
>>>
>>>Can someone out there provide a simple argument?

>  [Dave Dodson's good explaination deleted]

>One thing that no one has remarked on: The original poster's reference
>to *eigenvalues*.

>    As Dave Dodson points out, the condition number of A is
>norm(A) * norm(A-inverse).  No mention of eigenvalues.

>    In fact, the condition number of a matrix has (almost) nothing to
>do with the eigenvalues, as the example

>        ( 1  t )                        -1   ( 1  -t )
>    A = (      )   for t >> 1 shows:   A   = (       )
>        ( 0  1 )                             ( 0   1 )

>so the eigenvalues are both 1, but norm(A) = norm(A-inverse)
>which are approximately  t, so the condition number is approx.
>t**2.

>    What the condition number *is* related to are the *singular values*.

The singular values of A are the square roots of the eigenvalues of
the symmetric matrix AA'.  These are of course the eigenvalues of A
itself whenever A is symmetric and positive semidefinite.  So saying
that the condition number has (almost) nothing to do with the
eigenvalues is not accurate either.  The original poster may have had
the symmetric positive semidefinite situation in mind.

The eigenvalues of AA' also have an exact relationship to the norm of
A if we define the latter appropriately.  For example, if we define
the norm of A as:
                   |Ax|
	|A| = max  ---- 
               x    |x|

Then we get 

|A| = square root of the largest eigenvalue of AA' and

|Inverse(A)| = 1 / (square root of the smallest eigenvalue of AA')

and the condition number is the product of these two norms.
--
T. Scott Thompson              email:  thompson@atlas.socsci.umn.edu
Department of Economics        phone:  (612) 625-0119
University of Minnesota        fax:    (612) 624-0209

thompson@atlas.socsci.umn.edu (T. Scott Thompson) writes:

>The singular values of A are the square roots of the eigenvalues of
>the symmetric matrix AA'.  These are of course the eigenvalues of A
                      ^^^
This should have been A'A.  Many apologies for any confusion.

>|A| = square root of the largest eigenvalue of AA' and
                                                ^^^ 
>|Inverse(A)| = 1 / (square root of the smallest eigenvalue of AA')
                                                               ^^^
Same problem and correction in these two statements.
--
T. Scott Thompson              email:  thompson@atlas.socsci.umn.edu
Department of Economics        phone:  (612) 625-0119
University of Minnesota        fax:    (612) 624-0209

Dave Dodson <dodson@convex.COM> writes:

[original question deleted]
[Dodson's reply deleted except for the following excerpts.]

>Suppose you have a non-singular linear system A.x = b and you perturb the
>matrix and right-hand side.  How much is the solution perturbed; i.e.,
>how big is dx if (A+dA).(x+dx) = b+db ?

>	||dx||      ||A||.||A^-1||     ||dA||   ||db||
>	------ <= ----------------- . (------ + ------)
>	||x||	  1-||A^-1||.||dA||    ||A||	||b||

Dodson's proof assumes that A^(-1) exists, which is not necessary.
Here is a similar but somewhat shorter derivation that only assumes
that there is a unique solution to Ax = b.  It derives the bound using
a first order differential approximation to the deviations, which
simplifies matters considerably.

If Ax = b has a unique solution then A has full column rank, in which
case an explicit formula for the solution is

	x = (A'A)^(-1)A'b = Cb

where C = (A'A)^(-1)A'.  Note that C'C = (A'A)^(-1).

Let |y| denote the usual Euclidean norm for a vector y, and define the
norm of a matrix Z by:

	          |Zy|
	|Z| = max ---- = Sqrt( maximum eigenvalue of Z'Z )
	       y   |y|

This norm satisfies: |ZW| <= |Z| |W|.  Applying the definition we get

	|C| = Sqrt( maximum eigenvalue of (A'A)^(-1) )
	    = 1 / Sqrt( minimum eigenvalue of A'A )

First order perturbations to the linear equation Ax = b must satisfy:

	A(dx) + (dA)x = db      (i.e. we drop the term (dA).(dx) )
or
	A(dx) = db - (dA)x	
or
	dx = C ( db - (dA)x )

Thus	|dx| / |x| <= |C| ( |db|/|x| + |dA| )
	           <= |C| |A| [ |db|/(|A| |x|) + |dA|/|A| ]

However since b = Ax we have |b| <= |A| |x|.  Substitution then gives

	|dx| / |x| <= |C| |A| ( |db|/|b| + |dA|/|A| )

The condition number is |C| |A| = |(A'A)^(-1)A| |A| which is the
square root of the ratio of the largest and smallest eigenvalues of
A'A.

This argument is not quite as useful for a formal analysis of
precision loss as as Dodson's, since it applies to first order
differentials rather than finite differences.  I find it somewhat
"simpler", to understand, however.
--
T. Scott Thompson              email:  thompson@atlas.socsci.umn.edu
Department of Economics        phone:  (612) 625-0119
University of Minnesota        fax:    (612) 624-0209