dgh@validgh.com (David G. Hough on validgh) (06/18/91)
jbs@ibm.com asserts about solving linear equation with matrices composed of elements from a uniform random distribution: > Problems > of size 1000 appear to typically have condition number between 10**5 > and 10**6. I haven't observed any evidence to support this. If uniform random 1000x1000 matrices had such condition numbers then it would be impossible to solve them accurately with Gaussian elimination. I did report that gaussian elimination without pivoting was somewhat unstable on matrices of that size, but that only reflects on the stability of that method, rather than the underlying condition of the problem, independent of the method. For instance, the normalized residual computed in the linpack program, RESIDn = RESID/( N*NORMA*NORMX*EPS ) grows very slowly; typically 1.5-2 for 100x100 systems, and around 10 for 1000x1000 systems, using IEEE arithmetic and conventional partial pivoting. So the growth of RESID itself is only slightly faster than linear with N. Of course the size of the residual b-Ax doesn't tell you too much about the size of the error in x, in general. Luck is with us in this case, however: for single precision matrices, the measured difference in the computed x's vs. the nominally correct x's whose elements are all exactly one is in the same ratio, e.g. about 1e-4 for 1000x1000 vs. 1e-5 for 100x100. Of course the actual correct answer for this problem isn't a vector of 1's since the right hand side b was computed by rounding the result of computing A*(1 1 1 1 1 ...) and so the computed b differs from the one corresponding to the vector of 1's, and correspondingly the correct x corresponding to the right hand side b actually used is somewhat different - but fortunately that doesn't matter much for this well-conditioned problem. But this is a good illustration of a point I made previously: in general you don't know the error, because that would be tantamount to knowing the answer, and most people don't spend most of their time solving problems with known answers. So the normalized residual is a good indication of whether an algorithm has solved a linear equations problem about as well as can be expected in a particular floating-point precision. If an answer with a small error bound is desired, regardless of the condition of the problem, something else will have to be done. This reminds me of a problem that was posed to me several days in advance of my qualifying examination: A very eminent applied mathematician and one of his students had submitted a paper to a prestigious journal showing experimentally that all the theory about iterative improvement must be wrong. (Iterative improvement is a method for improving the accuracy of the answer x for linear equation problems that are moderately ill-conditioned. When the residual can be computed in higher precision, the accuracy of x can often be improved to nearly full working precision. ) The authors had computed solutions to problems from several collections of test matrices of varying conditions and to a number of random matrices and had seldom seen iterative improvement yield more than a digit of improvement. What was really going on here? I couldn't figure it out before the qualifying exam, but for better or worse I passed anyway. -- David Hough dgh@validgh.com uunet!validgh!dgh na.hough@na-net.ornl.gov