Using the SVD to solve linear systems Message #1 Posted by Rodger Rosenbaum on 31 Mar 2006, 7:53 p.m.
In this post, I'm going to show how to use the SVD to solve illconditioned systems. I'll be using the HP49G+ to do the computations. I have a special directory for trying out various SVD related schemes. Valentin's original illconditioned system will be solved. I created several variables in this directory, namely AA, BB, A, B, U, V and S. I use AA and BB to hold backups of Valentin's A and B matrices. The A and B variables hold the linear system under investigation and U, V and S hold the SVD of the A matrix. Sometimes A and B will hold AA and BB with small changes, so having AA and BB save one from having to type in Valentin's system again.
Some programs used are as follows (sometimes I leave results on the stack, sometimes not, with no particular organizational scheme in mind!):
>SV (this program computes the Singular Value Decomposition of a matrix on level 1 of the stack, and leaves the results, U V S, on the stack)
<< SVD >>
SV> (this program expects the same objects on the stack that SVD would leave there, namely U V S. It then makes a suitable diagonal matrix from S and multiplies the 3 objects, recreating the matrix whose SVD was in U, V and S.)
<< 3 PICK SIZE OBJ> DROP2 3 PICK SIZE OBJ> DROP2 2 >LIST DIAG> SWAP * * >>
PSU (The product TRANSPOSE(V) * DIAG(1/s) * TRANSPOSE(U) is computed. This program expects an integer N on the stack. It then recalls U, V, and S from the variables and computes a pseudoinverse for the matrix whose SVD is stored in those variables. This pseudoinverse can be rank reduced by means of the integer initially on the stack when PSU is invoked. What it does is to zero N of the reciprocals of the singular values starting with the smallest and working toward the largest. If the matrix SVD stored in the variables is of a square matrix, and if N is 0, then the ordinary inverse is returned.)
<< > N << U TRN V TRN S OBJ> OBJ> DROP >LIST 1E499 ADD INV OBJ> 1 SWAP 2 >LIST
>ARRY DUP SIZE OBJ> DROP N  2 >LIST { 1 1 } SWAP SUB 3 PICK SIZE OBJ> DROP2 3
PICK SIZE OBJ> DROP2 SWAP 2 >LIST DIAG> ROT * * >>
Cond (this program computes the condition number of a matrix on the stack as the ratio of the largest to the smallest singular values. Be sure to spell the name with a capital "C" and lower case "ond" so as not to conflict with the builtin COND command.)
<< SVL DUP 1 GET SWAP DUP SIZE 1 GET GET / >>
MADJ (This program mean adjusts the columns of a matrix on the stack. It computes the mean of each column and subtracts that value from all the elements of the corresponding column.)
<< DUP TRN DUP SIZE OBJ> DROP > c r << r 1 2 >LIST 1 CON * r / c DIAG> c r 2 >LIST
1 CON * TRN  >> >>
CORM (This program computes the columnwise correlation matrix of a matrix on the stack. It uses MADJ. It only works right for a square matrix, but it doesn't check for this.)
<< MADJ DUP TRN SWAP * DUP >DIAG OBJ> OBJ> DROP >LIST SQRT INV OBJ> 1 >LIST >ARRY
DUP SIZE OBJ> DROP DIAG> SWAP OVER * * >>
Three other programs I gave in another post, IREF, SRREF and TST0 should be in the directory, too.
First, type Valentin's system into AA and BB; the 7 x 7 design matrix into AA and the 7 x 1 column matrix into BB; then recall AA and BB and store them into the A and B variables. Recall A and press Cond and see 3.17E12; this is a very illconditioned system, and the usual rule of thumb says that you will lose about LOG10(3.17E12) ~= 12 digits in your solution. Fortunately, the HP49G+ uses 15 digits internally to solve the system with the B/A method, so we actually get a reasonable solution Clear the stack.
Type A >SV and when it's done (takes a while), use the blue left arrow key as a shortcut for store, and store the 3 objects on the stack into the variables S, V, and U in that order. Type U V S SV> and see the reconstituted matrix A. Recall A and type . You can see some tiny residual errors from all the number crunching that has occurred. Clear the stack.
Now type 0 PSU. You will see the inverse of A on the stack. (this assumes that the SVD of A is stored in U, V, and S. If not, then you will see the (pseudo) inverse of whatever matrix had its SVD stored in those variables.) Then recall B to the stack and press *. You have just calculated INVERSE(A) * B, which is a traditional way to solve the system. You won't get the same solution as if you calculate the inverse of A with the 1/x key because the A matrix is so illconditioned that numerical errors are becoming overwhelming, and cause differences in the two methods of computing the inverse. Recall the variable S and look at the last element of that vector (the smallest singular value of the matrix A); it is 1.257E11. This is the distance to the nearest singular matrix; the small size of this number is an indication of severe illconditioning. Clear the stack.
When PSU computes the reciprocals of the singular values to create DIAG(1/s), the last one is very large, namely, 7.955E10. This is many orders of magnitude larger than any of the other singular value reciprocals, and its size is the reason the x = INVERSE(A) * B method has difficulties. If the solution vector x computed by this method is used with a perturbed A matrix to compute A * x = B' and the residual (B'B) is computed, it can be shown that the perturbations can be magnified by a factor about equal to that very large last element in DIAG(1/s).
We can eliminate this difficulty by putting the number 1 on the stack and then executing PSU; then recall B and press *; this replaces that very large reciprocal of the last singular value by zero. We now have a solution for the system that has avoided the illconditioning of the nearness of the matrix A to a singular matrix. This solution has avoided the magnification of small perturbations (errors) in the A matrix at the expense of the size of the residual. If the solution is computed with the B/A method, the residual is 0. As we compute successive solutions with N PSU B * (N starting at 0 and increasing 1 at at time), the error magnification decreases because the largest singular value reciprocal in DIAG(1/s) is smaller, but the residual size increases. The residual size obtained with 1 PSU B * B  ABS is 2E9. We can go further; compute 2 PSU B * B  ABS, dropping the last two singular value reciprocals from the solution, and we get a residual size of 7.6E2. We can go further still, but the size of the residual should not be much larger than the expected error in the elements of the A matrix. Since the A matrix was a series of measurements with only two significant digits of accuracy, the error in each element should be taken to be +.05 (1/2 LSD), so the solution given by 2 PSU B * is probably acceptable.
To see the power of the SVD method of solution, modify the A matrix (in the variable A) so that column 7 is identical to column 3. Now the condition number of the system is infinite and it cannot be solved by the B/A method, or the INVERSE(A)*B method, no matter how many significant digits your calculator or computer has. But type A >SV and then store the U V S objects off the stack into the variables U, V and S. Now type 1 PSU B * and see a solution with a residual size of about 3.7E9, a quite respectable result.
Try modifying the A and B variables as I suggested in earlier posts and experiment with solutions you get with the SVD method vs. the traditional methods.
You can use the CORM command to see which of the columns of a matrix are highly correlated. If you want to see how well the rows are correlated, transpose the matrix first. Recall the AA matrix and transpose it; then type CORM. You will see that the {7 3} element of the correlation matrix is .999035, indicating that rows 3 and 7 are closely correlated. This is a clue as to how Valentin created this matrix. Normally, a random matrix won't have such a high correlation between any 2 rows. Test this by executing {7 7} RANM CORM repeatedly and examining the offdiagonal elements; you won't see .999 often.
