[HP-15C] AM#7 7x7 linear system: results, methods Message #37 Posted by Karl Schneider on 13 Apr 2008, 5:33 a.m., in response to message #11 by Valentin Albillo
Valentin (and all readers) --
I've done some of my own experimenting with this interesting problem, which is to solve an exactly-determined 7th-order linear system based on "Albillo Matrix #7" from Valentin's "Mean Matrices" paper from three years ago.
Quote:
13 x1 + 72 x2 + 57 x3 + 94 x4 + 90 x5 + 92 x6 + 35 x7 = 453
40 x1 + 93 x2 + 90 x3 + 99 x4 + x5 + 95 x6 + 66 x7 = 484
48 x1 + 91 x2 + 71 x3 + 48 x4 + 93 x5 + 32 x6 + 67 x7 = 450
7 x1 + 93 x2 + 29 x3 + 2 x4 + 24 x5 + 24 x6 + 7 x7 = 186
41 x1 + 84 x2 + 44 x3 + 40 x4 + 82 x5 + 27 x6 + 49 x7 = 367
3 x1 + 72 x2 + 6 x3 + 33 x4 + 97 x5 + 34 x6 + 4 x7 = 249
43 x1 + 82 x2 + 66 x3 + 43 x4 + 83 x5 + 29 x6 + 61 x7 = 407
which has the obvious, unique solution:
x1 = x2 = x3 = x4 = x5 = x6 = x7 = 1
The determinant is exactly one, as indicated in Valentin's paper, which also noted that ill-conditioned matrices have "determinants that are far too small for their dimensions and element-size." A determinant of 1.00 might be suitable for solving an exactly-determined 2x2 system, but not necessarily a 7x7 system.
The HP-42S calculates the determinant as 0.0703238892937, and the "solution" as [-2.94919 1.00050 -2.52118 -2.54129 0.99756 4.54124 8.40909]T. I do not know which specific methods it employs.
The HP-15C calculates the determinant as -2.956799886, and the "solution" as [-0.55584 1.00020 -0.38722 -0.39514 0.99904 2.39512 3.91892]T.
As described in the HP-15C Advanced Functions Handbook (AFH): For calculation of determinants, inverses, and linear solutions, the HP-15C performs an LU decomposition of the original matrix A using the Doolittle method. For linear solutions, the LU decomposition will be in-place; for determinants and inverses, it can be output to a different matrix. The form is PA = LU, where P represents the row-transpositions, L is a lower-triangular matrix (normalized so that each of its main-diagonal elements is unity and other elements have magnitudes not exceeding unity), and U is an upper-triangular matrix. L and U are stored in the same space as the original matrix A (omitting the main diagonal of L), with normally-unused bits storing the row-interchange information.
Any "P-1LU" decomposed matrix is indicated with "--" appended to its single-letter identifier in the HP-15C. This altered matrix -- if accurately computed -- will produce accurate determinants and linear solutions, much faster than if the original matrix were used.
It's instructive to know why HP implemented the pioneering matrix functionality of the HP-15C in this manner. (After all, why isn't it considered undesirable to alter the user's painstakingly-entered original matrix such that it can't be usefully edited, or to overwrite a different matrix that was stored in necessary RAM space?)
- First and foremost: The decomposition allows in-place matrix inversion, as the L and U components can be created in-place, inverted separately, and multiplied in-place. This conserves scarce user RAM, which is only 64 registers. Inversion of 6x6, 7x7, and 8x8 matrices are possible within available RAM space using this approach. (Reference: pages 96-98 of the AFH.)
- LU decomposition allows for faster linear solutions, because these can be done using simple back-substitution twice, first using U, then L. If repeated solutions are performed, the extra time spent on LU decomposition will be more than recovered by quicker multiple linear solutions.
- The LU transformation is consistent with practical usage of the matrix functionality. The determinant of a matrix is usually calculated to assess how near-singular it might be. Once the LU decomposition is calculated, the determinant is simply the product of the main-diagonal elements of U (before possible negation by P operations). The LU representation gives the user a fast-performing matrix for subsequent linear solutions, if those are undertaken. However, if the user wants to preserve the entered original matrix, he is responsible for allowing space for it (if possible), and for selecting appropriate result-matrix identifiers as needed.
So, how accurately does the HP-15C compute the LU decomposition of Albillo Matrix #7?
The HP-15C's LU-transformation matrix (not showing the hidden row-interchange data) is as follows:
LU =
48. 91. 71. 48. 93. 32. 67.
0.1458333333 79.72916667 18.64583334 -4.999999998 10.4375 19.33333333 -2.77083331
0.8333333333 0.2153122551 26.81865691 60.07656128 -78.74732166 64.17062974 10.76326104
0.0625 0.8317219754 -0.5199980516 65.39830469 41.55794805 49.28864425 7.7139377741
0.2708333333 0.5939378103 0.9954401520 0.3695366092 121.6443367 -10.24144422 4.935104748
0.8541666667 0.07865168535 -0.6753641545 0.6111292124 -0.631667889 4.893688725 -2.33890765
0.8958333333 0.006009929486 0.08515613588 -0.07776712548 0.07860975206 -0.1245007524 0.00000000074
LU partitioned; U above and L below:
48. 91. 71. 48. 93. 32. 67.
------------- 79.72916667 18.64583334 -4.999999998 10.4375 19.33333333 -2.77083331
0.1458333333 -------------- 26.81865691 60.07656128 -78.74732166 64.17062974 10.76326104
0.8333333333 0.2153122551 --------------- 65.39830469 41.55794805 49.28864425 7.7139377741
0.0625 0.8317219754 -0.5199980516 --------------- 121.6443367 -10.24144422 4.935104748
0.2708333333 0.5939378103 0.9954401520 0.3695366092 ----------------- 4.893688725 -2.33890765
0.8541666667 0.07865168535 -0.6753641545 0.6111292124 -0.631667889 ---------------- 0.00000000074
0.8958333333 0.006009929486 0.08515613588 -0.07776712548 0.07860975206 -0.1245007524 ---------------
At the end of this post are the results from Matlab v6.5, including (same-format) LU decomposition, singular value decomposition (SVD), condition number, inversion, and linear solution.
A careful comparison of the HP-15C's elements of L and U with those of Matlab shows that the HP-15C calculates these results almost as accurately as permitted by its 10 external and 13 internal significant digits. However, for this ill-conditioned matrix, it's not good enough for accurate results. Even Matlab's results using double-precision floating-point mathematics aren't perfect, as a comparison with Valentin's results from Mathematica in the "Mean Matrices" paper will show.
I found only two significant differences between the LU results of HP-15C and Matlab exceeding 7 units in the last place (or, "ulp") of the HP-15C's 10 significant digits:
- Element (7,7) of the HP-15C's U matrix was larger in magnitude by a factor of about 300 (7.4E-10 versus -2.5E-12), and had a different sign. Also note that (7,7) is the only small-magnitude term of this particular U matrix. Only two significant digits were calculated by the HP-15C's for U(7,7); this also was probably a symptom of the near-singularity of this 7x7 "Albillo Matrix".
- Element (7,2) of L was 38 ulp's larger than that of Matlab's result.
-- KS
-------------------------------------------------------------
A =
13 72 57 94 90 92 35
40 93 90 99 1 95 66
48 91 71 48 93 32 67
7 93 29 2 24 24 7
41 84 44 40 82 27 49
3 72 6 33 97 34 4
43 82 66 43 83 29 61
b
b =
453
484
450
186
367
249
407
[L,U,P]=lu(A)
L =
1.0000 0 0 0 0 0 0
0.1458 1.0000 0 0 0 0 0
0.8333 0.2153 1.0000 0 0 0 0
0.0625 0.8317 -0.5200 1.0000 0 0 0
0.2708 0.5939 0.9954 0.3695 1.0000 0 0
0.8542 0.0787 -0.6754 0.6111 -0.6317 1.0000 0
0.8958 0.0060 0.0852 -0.0778 0.0786 -0.1245 1.0000
U =
48.0000 91.0000 71.0000 48.0000 93.0000 32.0000 67.0000
0 79.7292 18.6458 -5.0000 10.4375 19.3333 -2.7708
0 0 26.8187 60.0766 -78.7473 64.1706 10.7633
0 0 0 65.3983 41.5579 49.2886 7.7139
0 0 0 0 121.6443 -10.2414 4.9351
0 0 0 0 0 4.8937 -2.3390
0 0 0 0 0 0 -0.0000
P =
0 0 1 0 0 0 0
0 0 0 1 0 0 0
0 1 0 0 0 0 0
0 0 0 0 0 1 0
1 0 0 0 0 0 0
0 0 0 0 1 0 0
0 0 0 0 0 0 1
format long
[L,U,P]=lu(A)
L =
Columns 1 through 5
1.00000000000000 0 0 0 0
0.14583333333333 1.00000000000000 0 0 0
0.83333333333333 0.21531225503005 1.00000000000000 0 0
0.06250000000000 0.83172197543768 -0.51999805134701 1.00000000000000 0
0.27083333333333 0.59393781029527 0.99544015199493 0.36953660933039 1.00000000000000
0.85416666666667 0.07865168539326 -0.67536415452818 0.61112921260013 -0.63166788876545
0.89583333333333 0.00600992944865 0.08515613582111 -0.07776712546120 0.07860975199350
Columns 6 through 7
0 0
0 0
0 0
0 0
0 0
1.00000000000000 0
-0.12450075225428 1.00000000000000
U =
1.0e+002 *
Columns 1 through 5
0.48000000000000 0.91000000000000 0.71000000000000 0.48000000000000 0.93000000000000
0 0.79729166666667 0.18645833333333 -0.05000000000000 0.10437500000000
0 0 0.26818656911419 0.60076561275150 -0.78747321661876
0 0 0 0.65398304671896 0.41557948068398
0 0 0 0 1.21644336729411
0 0 0 0 0
0 0 0 0 0
Columns 6 through 7
0.32000000000000 0.67000000000000
0.19333333333333 -0.02770833333333
0.64170629736086 0.10763261039979
0.49288644224680 0.07713937740537
-0.10241444215755 0.04935104746538
0.04893688732152 -0.02338980764163
0 -0.00000000000250
P =
0 0 1 0 0 0 0
0 0 0 1 0 0 0
0 1 0 0 0 0 0
0 0 0 0 0 1 0
1 0 0 0 0 0 0
0 0 0 0 1 0 0
0 0 0 0 0 0 1
cond(A)
ans =
3.1740e+012
svd(A)
ans =
399.0161
120.2434
81.6069
65.4718
15.0734
2.0586
0.0000
format long
svd(A)
ans =
1.0e+002 *
3.99016086263324
1.20243367935773
0.81606914636306
0.65471837097341
0.15073408085272
0.02058584440459
0.00000000000126
format short
inv(A)
ans =
1.0e+009 *
0.0001 -0.0005 -2.1289 -0.0365 0.2652 0.0036 2.1298
-0.0000 0.0000 0.0003 0.0000 -0.0000 -0.0000 -0.0003
0.0001 -0.0005 -1.8982 -0.0326 0.2364 0.0032 1.8990
0.0001 -0.0005 -1.9090 -0.0328 0.2378 0.0032 1.9098
0.0000 -0.0000 -0.0013 -0.0000 0.0002 0.0000 0.0013
-0.0001 0.0005 1.9090 0.0328 -0.2378 -0.0032 -1.9098
-0.0001 0.0010 3.9941 0.0686 -0.4975 -0.0067 -3.9957
format long
inv(A)
ans =
1.0e+009 *
Columns 1 through 5
0.00007108329001 -0.00050746920953 -2.12894026190704 -0.03654455920893 0.26516332517146
-0.00000000900016 0.00000006400116 0.00026839087074 0.00000460708361 -0.00003342860666
0.00006337915020 -0.00045246921139 -1.89820387651406 -0.03258382832977 0.23642469462416
0.00006374115677 -0.00045505425830 -1.90904900833100 -0.03276999170826 0.23777547513807
0.00000004400080 -0.00000031400570 -0.00131725090541 -0.00002261141035 0.00016406597746
-0.00006374015676 0.00045504725818 1.90901964679815 0.03276948769911 -0.23777181807170
-0.00013335942020 0.00095206427802 3.99411063093057 0.06856134724808 -0.49747363712801
Columns 6 through 7
0.00355411549992 2.12981303474609
-0.00000044800813 -0.00026850087274
0.00316891750936 1.89898205763645
0.00318702283794 1.90983163553408
0.00000219903991 0.00131779091521
-0.00318697383705 -1.90980226200101
-0.00666788600848 -3.99574804264624
x=inv(A)*b
x =
1.00000000000000
0.99999998509884
1.00000000000000
1.00024414062500
1.00000005960464
0.99987792968750
1.00000000000000
--------------------------------------------------------------
Edited: 16 Apr 2008, 2:26 a.m.
|