[HP15C] 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 exactlydetermined 7thorder 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 illconditioned matrices have "determinants that are far too small for their dimensions and elementsize." A determinant of 1.00 might be suitable for solving an exactlydetermined 2x2 system, but not necessarily a 7x7 system.
The HP42S 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 HP15C 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 HP15C Advanced Functions Handbook (AFH): For calculation of determinants, inverses, and linear solutions, the HP15C performs an LU decomposition of the original matrix A using the Doolittle method. For linear solutions, the LU decomposition will be inplace; for determinants and inverses, it can be output to a different matrix. The form is PA = LU, where P represents the rowtranspositions, L is a lowertriangular matrix (normalized so that each of its maindiagonal elements is unity and other elements have magnitudes not exceeding unity), and U is an uppertriangular matrix. L and U are stored in the same space as the original matrix A (omitting the main diagonal of L), with normallyunused bits storing the rowinterchange information.
Any "P^{1}LU" decomposed matrix is indicated with "" appended to its singleletter identifier in the HP15C. 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 HP15C in this manner. (After all, why isn't it considered undesirable to alter the user's painstakinglyentered 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 inplace matrix inversion, as the L and U components can be created inplace, inverted separately, and multiplied inplace. 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 9698 of the AFH.)
 LU decomposition allows for faster linear solutions, because these can be done using simple backsubstitution 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 nearsingular it might be. Once the LU decomposition is calculated, the determinant is simply the product of the maindiagonal elements of U (before possible negation by P operations). The LU representation gives the user a fastperforming 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 resultmatrix identifiers as needed.
So, how accurately does the HP15C compute the LU decomposition of Albillo Matrix #7?
The HP15C's LUtransformation matrix (not showing the hidden rowinterchange 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 (sameformat) LU decomposition, singular value decomposition (SVD), condition number, inversion, and linear solution.
A careful comparison of the HP15C's elements of L and U with those of Matlab shows that the HP15C calculates these results almost as accurately as permitted by its 10 external and 13 internal significant digits. However, for this illconditioned matrix, it's not good enough for accurate results. Even Matlab's results using doubleprecision floatingpoint 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 HP15C and Matlab exceeding 7 units in the last place (or, "ulp") of the HP15C's 10 significant digits:
 Element (7,7) of the HP15C's U matrix was larger in magnitude by a factor of about 300 (7.4E10 versus 2.5E12), and had a different sign. Also note that (7,7) is the only smallmagnitude term of this particular U matrix. Only two significant digits were calculated by the HP15C's for U(7,7); this also was probably a symptom of the nearsingularity 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.
