Re: Strange behavior of MatrixConditionNumber
- To: mathgroup at smc.vnet.net
- Subject: [mg66557] Re: [mg66541] Strange behavior of MatrixConditionNumber
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Sat, 20 May 2006 04:46:54 -0400 (EDT)
- References: <200605190740.DAA12848@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
The Phantom wrote:
> I get the following in Mathematica 5.2:
>
> <<LinearAlgebra`MatrixManipulation`
>
> ( 5 -8 -9 -7 -4 )
> ( 5 5 -8 3 1 )
> a = ( 3 -2 -6 8 -1 );
> ( 1 -7 9 -5 -2 )
> (-9 9 -4 -7 -5 )
>
> q = SetPrecision[a,99];
>
> N[MatrixConditionNumber[q,2],12]
> Out(1) = 6.50915201983
>
> N[Norm[q,2]*Norm[Inverse[q],2],12]
> Out(2) = 7.54449348901
>
> These two values should be the same, shouldn't they?
Not really. From the comments at the beginning of that package:
(*:Limitations:
The matrix norm is only approximated for the case p=2. The
inverse norm is only approximated for all of the allowable
values of p. Hence, while we use the precision of the input
the process of computing the norms, we only return a machine
number, as at best a few digits are correct. This is adequqate
because typically one is interested only in the order of
magnitude of the norm.
The techniques used, while quite reasonable, may fail for
specialized classes of input. It should also be noted that
any technique used to find the inverse matrix norm (especially
including explicit computation of the inverse matrix) will
fail if the input has inadequate precision. It is easy to
recognize this case, because the log of the condition number
will typically be larger than the precision of the input. In
terms of efficiency and correctness, it is clear from the code
that one can get a good value for either the 1-or-infinity
norm of a matrix. If one is concerned that the matrix may be
pathological for the inverse-norm estimators, perhaps the
safest ploy is to estimate the 2-norm of the inverse; this
uses the power method, and is fairly accurate. But note that
the 1-and-infinity estimators also use an abbreviated
power-method estimate as a reality check on the main method,
so these should not do anything drastically wrong. Note also
that we use a fixed number of iterations in the power-method
computations. One could attain more accuracy (at expense of
run-time) by altering these to use e.g. FixedPoint with a
SameTest option to determine when the desired accuracy has
been met.
*)
> And what's really strange, I get different values for
> N[MatrixConditionNumber[q,2],12] when I execute the program repeatedly.
>
> Does anybody get this strange result?
This is explained by the use of random initial vectors for the iterations.
Specifically, here is a bit of code used in the package (you can see
this yourself if you look at it in a text editor).
EstimatedMatrix2Norm[mat_] := Module[
{j, x, y, len = Length[mat], prec = Precision[mat]},
x = Table[Random[Real, {10, 11}, prec], {len}];
For [j = 1, j <= 5, j++,
x /= Norm[x, Infinity];
y = x;
x = mat.x;
x = Conjugate[Transpose[mat]].x;
x = SetPrecision[x, prec];
];
Sqrt[Norm[x]/Norm[y]]
]
If we change the loop bound from 5 to 500 the estimate for the 2-norm of
q goes up from around 15.3 to about 18, which gets at the first part of
your question. Also the numbers will tend to converge so you'll likely
see results that agree to more places.
As per the comment, one might replace this loop with FixedPoint
iterations. The problem is that can become quite slow. The general need
for condition number estimates is that they be quickly found and in the
right ballpark, not that they be "exact", hence the choice of relatively
few iterations.
Daniel Lichtblau
Wolfram Research
- References:
- Strange behavior of MatrixConditionNumber
- From: The Phantom <phantom@aol.com>
- Strange behavior of MatrixConditionNumber