|HP-71 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 HP-71B (ROM 1BBBB and earlier) in which some leap years were incorrectly identified. J-F Garnier responded that it was a known bug, but did not cite his own web page listing them:
At the bottom of the page are listed some inexact results of matrix arithmetic using the HP-71 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 HP-71B (1983) and its Math Pac ROM (1984) utilizes linear-algebra methods similar to those of the 1982 HP-15C (with its 10-digit Coconut processor), adapted for the 12-digit Saturn processor. These methods seem to have been refined further for the Saturn-based HP-28C (1986), HP-42S (1988), HP-48S (1990), and HP-48G (1993), which do not exhibit all of these minor flaws.
Valentin Albillo has explained that the functionality contained in the HP-71 Math ROM was removed from the base HP-71 ROM in order to make room for other features (i.e., support for the IEEE floating-point standard that was then a proposal and was later adopted). So, time might have been short for incorporating all refinements in the HP-15C to the HP-71's new Saturn processor.
For reasons of RAM size and execution speed, the HP-15C 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 lower-triangular unit matrix L having ones on the main diagonal, while the upper-triangular matrix U has non-zero values. The permutation matrix P performs any necessary row transpositions, and is an identity matrix if none are required:
PA = LU => A = P-1LU
where det (U) is the product of its main-diagonal entries, and det (P-1) is odd for an odd number of row transpositions.
det(A) = det(P-1) * det(L) * det(U) = (-1 or 1) * 1 * det(U)
It follows that
A-1 = (P-1LU)-1 = U-1L-1P
On the HP-15C, L and U are composed "in place", and row-swap 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 near-singular input matrix, the HP-15C will also modify elements slightly to allow the necessary LU decomposition without divide-by-zero or overflow.
The Saturn-based models use a Crout algorithm, which produces an upper-triangular unit matrix U. Otherwise, it appears that all of these models take the same basic approach as that of the HP-15C. The HP-71B's methods are stated on pp. 86-88 of the Math Pac Owner's Manual.
For the matrix
A = [ 1 2 ]
[ 3 4 ]
the following can be extracted from the HP-15C 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 HP-48G 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 non-integer or inexact values in determinants and inverses. However, the HP-15C, HP-42S, and HP-48G did produce 'cleaner' results, as listed below:
exact, HP-48G HP-71B, HP-28C, HP-42S, HP-48S(?) HP-15C
[ -2 1 ] [ -1.99999999998 0.999999999994 ] [ -2 0.9999999999 ]
[ 1.5 -0.5 ] [ 1.49999999999 -0.499999999997 ] [ 1.5 -0.5 ]
As J-F Garnier pointed out, and as is documented on page 88 of the HP-71B Math Pac Owner's Manual, inverses can be computed more accurately by using the linear-system command MAT X = SYS(A,I) with an identity matrix as I. In this case, an exact inverse is returned by the HP-71B (but not by the HP-42S "SIMQ" command).
It seems possible that some algorithmic refinements in the HP-15C did not get incorporated into the HP-71 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 Crout-algorithm L and U matrices might explain the determinant of 2.0000000001
(3. * .66666666667), as well as the slightly-inexact elements of the inverse.
The HP-28C, HP-42S, and HP-48G return integer values for the determinants, and in the process do not warn of overflows or return "-0". However, the HP-28C and HP-42S return the very same inexact matrix-inverse values as the HP-71B does, while the HP-48G returns the exact values.
The HP-48G does return "LU error: Overflow" for LU decomposition of the singular matrix
[ 0 1 ]
[ 0 0 ]
apparently becuase no non-zero 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:
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?
Edited: 7 Sept 2010, 4:58 a.m. after one or more responses were posted