Re: Rodger was right Message #8 Posted by Rodger Rosenbaum on 8 May 2006, 10:10 p.m., in response to message #7 by Palmer O. Hanson, Jr.
Palmer said:
I don't propose to blindly use Valentin's linear equation problem as the sole measure for selection of a solver. I do say that, from a practical standpoint, if one solver can solve his problem successfully, and another cannot, then in the absence of other problems which yield contrary results I would probably decide to use the solver which can solve his problem successfully. Such a solver may be an overkill for many problems; however, such a solver may provide some comfort.
What I'm trying to say is that the results we're getting from the various methods of solution of this highly illconditioned problem are mostly noise and accidents, and we really can't draw conclusions about the performance of any of them as a result of testing them on Valentin's matrix. None of them is trustworthy with such a high condition number.
Consider this: the original A matrix has the number 7.2 as its {1 2} element. If you add 1/90 to this element, the matrix becomes exactly singular. Notice that 1/90 expressed as a decimal is .01111111...; this means that if you simply append the numeral 1 to the 7.2 repeatedly, the modified matrix's determinant approaches zero more and more closely as you add 1's. Here is a table of determinants and condition numbers for the A matrix with a succession of 1's added to the {1 2} element. These values were computed with Mathematica using very high precision; the HP49G+ in approximate mode can only get a few digits right for the condition number.
1,2 element Exact Determinant Condition Number
7.2 1/1E7 3.17397549713E12
7.21 1/1E8 3.16183337880E13
7.211 1/1E9 3.16061908862E14
7.2111 1/1E10 3.16049765882E15
Now, an interesting thing happens. For the case where the {1 2} element of the A matrix is 7.211, and the {1 1} element of the B matrix is 45.311, the LSQ method gets a solution of [1 1 1 1 1 1 1]T, but the B/A and RREF methods get a substantially worse result. Shall we therefore conclude that the LSQ method is much better than the others? No. It is just an accident, numerical happenstance, so to speak, that we get an exactly correct result with the LSQ method for a system with a condition number 100 times larger than the original system.
The same thing happens with the B/A method when applied to the original problem. The fact that it gets [1 1 1 1 1 1 1]T is just an accident. If we transpose the original A and change B to [19.5 58.7 36.3 35.9 47. 33.3 28.9] and solve with B/A, we no longer get [1 1 1 1 1 1 1]T. Why not? It seems like we ought to; it's the same size problem, with the same condition number. It's because the result we get when we solve a system with a condition number of 3E12 (or greater) is contaminated with a large amount of noise, and sometimes a result may be very close to the exact solution just by chance.
If we want a valid determination of the performance of the various methods, we need to solve many systems of similar condition number, so that we aren't fooled by an accidental good result. If we do this, we will find that the HP49G+, like any computing device, is bound by some rules. The rule of thumb is that a matrix computation will lose LOG10(k) digits of accuracy, where k is the condition number of the problem. The INV(A)*B method is the worst, because the inverse is returned to the stack with its elements rounded to 12 digits, and then this reduced accuracy inverse is used to multiply B. The other methods keep intermediate results to 15 digit accuracy, and only round to 12 digits when delivering the final result.
For the case where the {1 2} element of A becomes 7.2111, the HP49G+ thinks the determinant is exactly zero. Make the {1 1} element of B equal to 45.3111. The B/A method fails; the LSQ method returns a result which is mostly error, as does the RREF method. The internal COND command fails when applied to this modified A matrix. The SVD command fails to get the smallest singular value right (make sure flag 54 is set when using the SVD for these sorts of problems). Here, the condition number is about 3E15, and even the internal 15 digit calculations can't do the job.
But, the SVD method can still solve the system, since the smallest singular value is discarded for the solution.
My point is, when the condition number of the system is very high, on the order of 10^N, where N is the number of digits used for the calculations, the ordinary methods of solution just don't work. But the SVD method which I described in other posts still works. I even gave an example of a system whose A matrix was exactly singular, and showed that the SVD method gave a reasonable solution.
In an earlier post, I suggested:
Try reducing the condition number by a factor of 10 by changing the {1 2} element of the A matrix from 7.2 to 7.1 (in the original system), and the {1 1} element of B to 45.2; the change to B is to keep the exact solution as [[ 1 1 1 1 1 1 1 ]]T. Then solve the system by the various methods and compare results.
Reduce the condition some more, by a factor of about 473, by changing the {6 2} element of A from 7.2 to 7.1, and changing the {6 1} element of B from 24.9 to 24.8, then solve this system by the various methods and compare results.
If one does this, the various methods will be seen to give comparable results (except for the INV(A)*B method, as I explained above). Keep in mind that if you solve a system with a condition number of k (with k<10, let's say), you should only expect to get, nominally, 15k accurate digits (on the HP49G+; other calculators may use other than 15 digits internally) in your results. This means that you just aren't going to get a very good solution for a problem with a condition number of 3E12 on the HP49G+; if you do, it's an accident, and it doesn't mean that the method which accidentally got a good result is the best method for solving linear systems in general; or, for solving systems with high condition numbers. For that, one should use the SVD method.
If you want to get the mathematically exact result (which may not have much validity if it's a series of measurements from a physical system), then you will just have to use a system with more digits, such as the HP49G+ in exact mode, or one of the high accuracy programs on the PC, such as Derive, Mathematica, Maple, etc.
