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>