HP71 matrix methods (LU decomp?) Message #1 Posted by Karl Schneider on 5 Sept 2010, 8:13 p.m.
A recent post noted a bug in the HP71B (ROM 1BBBB and earlier) in which some leap years were incorrectly identified. JF Garnier responded that it was a known bug, but did not cite his own web page listing them:
http://www.jeffcalc.hp41.eu/emu71/bug71.html
At the bottom of the page are listed some inexact results of matrix arithmetic using the HP71 Math ROM. E.g., the calculated determinant of the matrix
[ 1 2 ]
[ 3 4 ]
is not quite an integer (2.00000000001 instead of 2), and the elements of the calculated inverse similarly are not quite exact. Also, the calculated determinant of the singular matrices
[ 0 1 ] and [ 0 0 ]
[ 0 0 ] [ 1 0 ]
are returned respectively as 0 with a message "WRN: Overflow", and 0 (with no message).
I suspect that the HP71B (1983) and its Math Pac ROM (1984) utilizes linearalgebra methods similar to those of the 1982 HP15C (with its 10digit Coconut processor), adapted for the 12digit Saturn processor. These methods seem to have been refined further for the Saturnbased HP28C (1986), HP42S (1988), HP48S (1990), and HP48G (1993), which do not exhibit all of these minor flaws.
Valentin Albillo has explained that the functionality contained in the HP71 Math ROM was removed from the base HP71 ROM in order to make room for other features (i.e., support for the IEEE floatingpoint standard that was then a proposal and was later adopted). So, time might have been short for incorporating all refinements in the HP15C to the HP71's new Saturn processor.
For reasons of RAM size and execution speed, the HP15C performs an LU decomposition of a square matrix A (larger than 1 x 1) prior to calculating determinants, inverses, and system solutions. The Doolittle algorithm is used, producing a lowertriangular unit matrix L having ones on the main diagonal, while the uppertriangular matrix U has nonzero values. The permutation matrix P performs any necessary row transpositions, and is an identity matrix if none are required:
PA = LU => A = P^{1}LU
det(A) = det(P^{1}) * det(L) * det(U) = (1 or 1) * 1 * det(U)
where det (U) is the product of its maindiagonal entries, and det (P^{1}) is odd for an odd number of row transpositions.
It follows that
A^{1} = (P^{1}LU)^{1} = U^{1}L^{1}P
On the HP15C, L and U are composed "in place", and rowswap information is stored in spare bits (no P matrix is stored). L and U can be inverted in place, and the resulting inverses can also be multiplied together in place.
In the case of a singular or nearsingular input matrix, the HP15C will also modify elements slightly to allow the necessary LU decomposition without dividebyzero or overflow.
The Saturnbased models use a Crout algorithm, which produces an uppertriangular unit matrix U. Otherwise, it appears that all of these models take the same basic approach as that of the HP15C. The HP71B's methods are stated on pp. 8688 of the Math Pac Owner's Manual.
For the matrix
A = [ 1 2 ]
[ 3 4 ]
the following can be extracted from the HP15C calculation of determinant or system solution:
P = [ 0 1 ] L = [ 1 0 ] U = [ 3 4 ]
[ 1 0 ] [ 0.3333333333 1 ] [ 0 0.6666666668 ]
The "LU" command of the HP48G returns the following:
P = [ 0 1 ] L = [ 3 0 ] U = [ 1 1.33333333333 ]
[ 1 0 ] [ 1 0.666666666667 ] [ 0 1 ]
It is easy to see where these methods may produce noninteger or inexact values in determinants and inverses. However, the HP15C, HP42S, and HP48G did produce 'cleaner' results, as listed below:
A^{1} =
exact, HP48G HP71B, HP28C, HP42S, HP48S(?) HP15C
[ 2 1 ] [ 1.99999999998 0.999999999994 ] [ 2 0.9999999999 ]
[ 1.5 0.5 ] [ 1.49999999999 0.499999999997 ] [ 1.5 0.5 ]
As JF Garnier pointed out, and as is documented on page 88 of the HP71B Math Pac Owner's Manual, inverses can be computed more accurately by using the linearsystem command MAT X = SYS(A,I) with an identity matrix as I. In this case, an exact inverse is returned by the HP71B (but not by the HP42S "SIMQ" command).
It seems possible that some algorithmic refinements in the HP15C did not get incorporated into the HP71 Math ROM. For example:
1. Not adjusting elements of the singular matrix A might explain the "overflow" messages that likely occurred during LU decomposition.
2. Some inexactitude in the calculations using elements of the Croutalgorithm L and U matrices might explain the determinant of 2.0000000001 (3. * .66666666667), as well as the slightlyinexact elements of the inverse.
The HP28C, HP42S, and HP48G return integer values for the determinants, and in the process do not warn of overflows or return "0". However, the HP28C and HP42S return the very same inexact matrixinverse values as the HP71B does, while the HP48G returns the exact values.
The HP48G does return "LU error: Overflow" for LU decomposition of the singular matrix
[ 0 1 ]
[ 0 0 ]
apparently becuase no nonzero element is present in the left column. However, a valid Crout LU decomposition can apparently be produced; please see this Forum thread from 2006.
A comment by Rodger Rosenbaum in a thread from 2007:
Quote:
The HP48G and its successors have matrix arithmetic that was reworked by Paul McClellan to use 15 digit arithmetic for internal operations, rounding to 12 digits only occurring after everything is done.
The HP71 and HP48S rounded internal operations to 12 digits at intermediate steps in the carrying out of matrix operations.
Further insights, anyone?
 KS
Edited: 7 Sept 2010, 4:58 a.m. after one or more responses were posted
