The Museum of HP Calculators

HP Forum Archive 16

[ Return to Index | Top of Index ]

Compute some Eigenvalues first thing in the morning, and nothing worse will happen to you all day.
Message #1 Posted by Eric Smith on 13 July 2006, 1:58 p.m.

[From my blog]

Actually it is fun, though a bit tedious. My HP-49G+ calculator, or Mathematica or Maple on my laptop, can of course determine Eigenvalues and Eigenvectors automatically, but the point of today’s Linear Algebra class (and Monday’s quiz) is to learn how to do it ourselves.

Computing an matrix of cofactors (and its transpose, the adjunct) is pretty tedious too. It doesn’t appear that the HP 49G+ has a built-in function to do that, so I wrote these during class.

COFACTORS:
« DUP SIZE LIST-> DROP DROP
  -> a n
  « 1 n FOR i
      1 n FOR j
        a i ROW- DROP j COL- DROP DET
        -1 i j + ^ *
      NEXT
    NEXT
  n n 2 ->LIST ->ARRY
  »
»

ADJOINT: « COFACTORS TRAN »

      
Re: Compute some Eigenvalues first thing in the morning, and nothing worse will happen to you all day.
Message #2 Posted by Rodger Rosenbaum on 13 July 2006, 4:37 p.m.,
in response to message #1 by Eric Smith

It's not "adjunct", it's "adjoint".

Your subject line suggests you're going to compute some eigenvalues, but instead you computed the adjoint.
When do we get to see the computation of eigenvalues? :-)

            
Re: Compute some Eigenvalues first thing in the morning, and nothing worse will happen to you all da
Message #3 Posted by Eric Smith on 13 July 2006, 5:00 p.m.,
in response to message #2 by Rodger Rosenbaum

I did compute some Eigenvalues, by hand, far too early this morning. If you really want to see them, I can type them up. :-)

                  
Re: Compute some Eigenvalues first thing in the morning, and nothing worse will happen to you all da
Message #4 Posted by Rodger Rosenbaum on 13 July 2006, 11:08 p.m.,
in response to message #3 by Eric Smith

It's not so much the eigenvalues that I want to see; it's the program that computes them.

      
Minors, cofactors, and adjoint
Message #5 Posted by Eric Smith on 13 July 2006, 5:44 p.m.,
in response to message #1 by Eric Smith

Here's a slightly improved version, using a separate subroutine to compute the minor of a matrix element, and avoiding unnecessary conversion of matrix elements to floating point:

MINOR:
« UNROT ROW- DROP SWAP COL- DROP DET »

COFACTORS: « DUP SIZE DUP 1 GET R->I -> a s n « 1 n FOR i 1 n FOR j a i j MINOR -1 i j + ^ * NEXT NEXT s ->ARRY » »

ADJOINT: « COFACTORS TRAN »

      
Re: Compute some Eigenvalues first thing in the morning, and nothing worse will happen to you all day.
Message #6 Posted by Valentin Albillo on 14 July 2006, 8:29 a.m.,
in response to message #1 by Eric Smith

Hi, Eric:

Eric posted:

"Actually it is fun, though a bit tedious."

    Yes, I concur, I do it every morning as well and sometimes also before going to sleep. This is how I do it, extracted from my Datafile article "HP-71B - Math ROM Baker's Dozen - Vol 2", freely downloadable from here:

       "The general case of finding all eigenvalues of any square matrix
        is hardly any more difficult, as demonstrated by this 5-liner 
        that will accept any square matrix and will promptly compute 
        both the coefficients of its Characteristic Polynomial and all
        its eigenvalues, real and/or complex:

    1 DESTROY ALL @ OPTION BASE 0 @ FIX 5 @ INPUT "N=";N @ M=N-1

    2 DIM A(M,M),U(M,M),C(N),D(N,N) @ COMPLEX E(N) @ MAT INPUT A

    3 MAT D=CON @ MAT U=IDN @ FOR I=0 TO N @ IF I THEN MAT A=A-U

    4 C(I)=DET(A) @ FOR J=1 TO N @ D(I,N-J)=I^J @ NEXT J @ NEXT I

    5 MAT C=SYS(D,C) @ MAT DISP C @ MAT E=PROOT(C) @ MAT DISP E

    Once the matrix has been entered, all it does is use DET (Determinant) to find N+1 values of the Characteristic Polynomial, then MAT ... SYS is used to explicitly find (and display) the coefficients of this polynomial [...] then all its roots (the eigenvalues) real and/or complex are computed at once by PROOT and displayed as well.

    Let’s test it with the same 5x5 matrix as before:

    >RUN [ENTER] N=5 [ENTER] A(0,0)? 5,1,2,0,4,1,4,2,1,3,2,2,5,4,0,0,1,4,1,3,4,3,0,3,4 [ENTER]

    -1.00000 19.00000 -79.00000 -146.00000 1153.0000 -1222.00000

    (1.49766, 0.00000) (3.36188, 0.00000) (-3.55784,0.00000) (5.67255,-0.00000) (12.02575,0.00000)

    so the Characteristic Polynomial is :

    P(x)= -x5 + 19 x4 - 79 x3 - 146 x2 + 1153 x - 1222

    and, as expected in this case, as this is a symmetric matrix, all five eigenvalues are real (imaginary parts = 0):

    x1 = 1.49766 x2 = 3.36188 x3 = -3.55784 x4 = 5.67255 x5 = 12.02575

    You can verify them by checking that their product equals the matrix determinant (DET(A) = -1222). To verify the Characteristic Polynomial, you can apply it to the matrix itself: the resulting value should be a zero matrix, i.e.: P(A) = ZER"

Best regards from V.

            
Re: Compute some Eigenvalues first thing in the morning, and nothing worse will happen to you all day.
Message #7 Posted by Rodger Rosenbaum on 15 July 2006, 4:21 a.m.,
in response to message #6 by Valentin Albillo

I was very surprised when this program failed when given the classical small example of a defective matrix, namely:

[[ 0 1 ]
 [ 0 0 ]]

This appears to be due to a bug in the HP71 DET command.

The HP71 correctly computes the determinant of:

  [[ 1 0 ]   [[ 0 0 ]      [[ 0 0 ]         [[ 0 1 ]
   [ 0 0 ]],  [ 1 0 ]] and  [ 0 1 ], but not [ 0 0 ]]

It gets an overflow error for the last one.

It gets the same overflow error for larger matrices in the same pattern.

                  
HP-71B DET calculation for imperfect matrices
Message #8 Posted by Karl Schneider on 16 July 2006, 6:14 p.m.,
in response to message #7 by Rodger Rosenbaum

Rodger --

On my HP-71B with Math ROM, the message, "WRN:Overflow" appeared briefly, and then the correct answer of zero appeared. Subsequently, I did a "DESTROY ALL", re-entered the matrix, and did "DET" again. No message accompanied the answer of zero. Weird...

It might be that the HP-71B Math ROM uses the same basic algorithms as did the HP-15C, which introduced the built-in matrix functionality. However, I don't know for certain.

Th HP-15C performs an LU decomposition of the matrix prior to calculating the determinant, as described in the HP-15C Advanced Functions Handbook available on the MoHPC CD/DVD sets.

The exact LU decomposition of

          A               L     *    U
      [[ 0  1 ]   =   [[ 1  0 ]  [[ 0  1 ]    
       [ 0  0 ]]       [ 0  1 ]]  [ 0  0 ]] 

= I * A

  • All elements above the diagonal of L are zero; all elements below are between -1 and +1 inclusive, and all elements on the diagonal are 1.

  • All elements below the diagonal of U are zero.

The encoded resultant LU matrix calculated by the HP-15C shows 10-99 and 10-10 as the diagonal elements of U. These are very close to, but not quite equal to zero. I'd assume that those values on the HP-71B have negative exponents of greater magnitude.

The determinant of an untransposed LU matrix is simply the product of the diagonal elements of the U matrix. This would be 10-109, an "underflow" on the HP-15C, which naturally returns zero without any warning or error code.

This calculation may be the source of the warning message from the HP-71B Math Pac.

-- KS

Edited: 17 July 2006, 1:53 a.m.

                        
Re: HP-71B DET calculation for imperfect matrices
Message #9 Posted by Rodger Rosenbaum on 19 July 2006, 4:16 a.m.,
in response to message #8 by Karl Schneider

Karl said:

The encoded resultant LU matrix calculated by the HP-15C shows 10-99 and 10-10 as the diagonal elements of U. These are very close to, but not quite equal to zero. I'd assume that those values on the HP-71B have negative exponents of greater magnitude.

The determinant of an untransposed LU matrix is simply the product of the diagonal elements of the U matrix. This would be 10-109, an "underflow" on the HP-15C, which naturally returns zero without any warning or error code.

But the result you got on the HP15C would give an underflow; how does that explain the overflow message on the HP71?

  I tried to calculate the LU decomposition of [[ 0 1 ]
                                                [ 0 0 ]]

on the HP48G and got no result--the calculator simply returned an error of "overflow". The HP49G+, on the other hand, did return a result (deleting the permutation matrix) of:

      L                   U
  [[ 0 0 ]     [[ 1  9.99999999999E499 ]
   [ 0 0 ]]     [ 0          1         ]]

I imagine that something like this may have happened on the HP71; the 9.99999999999E499 is obviously an overflow, but it hasn't affected the diagonal elements.

                              
Re: HP-71B DET calculation for imperfect matrices
Message #10 Posted by Werner on 19 July 2006, 8:26 a.m.,
in response to message #9 by Rodger Rosenbaum

Rodger,
the 48G and 49 work the same way for matrices containing real numbers; no doubt flag -21 (Overflow->error) is set on your 48 and clear on the 49+.
Cheers, Werner

                                    
Re: HP-71B DET calculation for imperfect matrices
Message #11 Posted by Rodger Rosenbaum on 19 July 2006, 10:34 a.m.,
in response to message #10 by Werner

Quote:
Rodger,
the 48G and 49 work the same way for matrices containing real numbers; no doubt flag -21 (Overflow->error) is set on your 48 and clear on the 49+.
Cheers, Werner


Hi, Werner,

You haven't posted for two months. Did you take a long vacation?

Do you mean to say: "(the LU decomposition on) the 48G and 49 work the same way for matrices containing real numbers;"

Parenthetical added by me. Because the unqualified statement is apparently not true; you said a couple of months ago:

"The 48G contained a bug in the eigenvalue routines, for a 5-by-5 Hilbert matrix the result of EGV/EGVL/SCHUR/SRAD was 'Undefined Result'. I corrected it, the 49G gives the correct result."

which would be an example of a case where the 48G and 49G *don't* work the same way for matrices containing real numbers.

As for flag -21, you are quite right; it was set in my HP48G.

Now that you're back, can you tell me why the 49G+ does this:

If you execute (in exact mode) [ 1 2 3 4 5 ] VANDERMONDE, you get a matrix full of exact numbers (type 28).

But, if you execute (in approximate mode) { 1. 2. 3. 4. 5. ] VANDERMONDE, you get a matrix of approximate numbers, except for the first column, which is all exact numbers.

                                          
Re: HP-71B DET calculation for imperfect matrices
Message #12 Posted by Werner on 20 July 2006, 5:55 a.m.,
in response to message #11 by Rodger Rosenbaum

Hi Rodger!
2 months? I have no idea, but it seems I do have to be a little more careful (or precise) in what I post ;-) Perhaps it should have been: the numerical (approximate) matrix routines are the same, save for some bugs that have been corrected. Or even: the 49G working in approx mode on type 3 or 4 matrices produces the same results as the 48G, except in a few cases where bugs have been corrected. Phew.
And about the VANDERMONDE command: that would be a bug. I don't own a 49G+, but the 49G produces the same result.
Werner

                              
DET and LU calculations; HP-15C methods
Message #13 Posted by Karl Schneider on 20 July 2006, 1:01 a.m.,
in response to message #9 by Rodger Rosenbaum

Rodger --

[Aagh! When quoting, it's better to use the "quote" or "italic" formatting feature instead of boldface. Thank you.]

You're probably right that the HP-71B, for the DET operation, is doing the same mathematics for the LU decomposition as the two RPL-based models you tried. (I got the same result with my HP-49G.)

The HP-15C calculates a unit L matrix (which contains ones on the main diagonal). The RPL-based models tried to compute general LU matrices, which contained overflow values in the U matrix.

The basics about LU decomposition and the Doolittle algorithm (used by the HP-15C) can be found here:

http://en.wikipedia.org/wiki/LU_decomposition

and here:

http://mathworld.wolfram.com/LUDecomposition.html

The Wikipedia article's discussion of "Solving linear equations" offers a tidy explanation why the computationally-slow HP-15C performs the LU decomposition for linear-system solving: It saves time for subsequent solutions, and is a bit faster yet when L or U is a unit matrix.

The HP-15C also performs the LU decomposition for determinants, and places the decomposition in the RESULT matrix. Why? Normally, the determinant is calculated in order to determine whether a matrix is singular. Once that answer is obtained, the encoded LU matrix may be used to obtain faster solutions for exactly-determined systems. If free memory space is limited, the RESULT can be specified with the same identifier as the original input matrix, which can be recalculated by inverting the encoded LU matrix twice.

Regards,

-- KS

Edited: 20 July 2006, 2:48 a.m.

                                    
Re: DET and LU calculations; HP-15C methods
Message #14 Posted by Rodger Rosenbaum on 20 July 2006, 3:13 p.m.,
in response to message #13 by Karl Schneider

Why is it "better" to use "quote" or "italic" rather than "bold"? It seems to me that it is a matter of personal preference. I find that the "bold" method looks ok. How about if we have a deal; I won't ask you to change your posting style if you don't ask me to change mine. :-)

I don't understand what you're getting at when you say, "The RPL-based models tried to compute general LU matrices, which contained overflow values in the U matrix."

Isn't the HP15's Doolittle LU decomposition "general"?

The HP15 Advanced Functions Handbook explains those particulars that you found in Wikipedia which apply to the HP15. It also mentions that the HP15 uses the Doolittle method. And, it explains the rationale for keeping the LU decomposition in the RESULT matrix, etc., etc. I've known all this since I bought my HP15 over 20 years ago.

The HP71 Math Pac book says that the HP71 uses the Crout method; the HP48G AUR says that the HP48G also uses the Crout method.

I would expect a Doolittle algorithm to return a result of:

                      A               L     *    U
                  [[ 0  1 ]   =   [[ 1  0 ]  [[ 0  1 ]    
                   [ 0  0 ]]       [ 0  1 ]]  [ 0  0 ]]

but another decomposition which one might expect a Crout algorithm to return would be:

                A               P     *     L    *    U
            [[ 0  1 ]   =   [[ 0 1 ]   [[ 0  0 ]  [[ 1  0 ]    
             [ 0  0 ]]       [ 1 0 ]]   [ 0  1 ]]  [ 0  1 ]]

Either one would be correct; unfortunately, the HP49G+ and HP48G don't do so well.

However, if the HP49G+ is used to compute the LU decomposition of the following matrix:

                             L              U
A = [[ 10-450 1 ]   =   [[ 10-450 0 ] *  [[ 1  10450 ]
     [  0    0 ]]       [  0    0 ]]    [ 0   1   ]]

we can see what's happening. If we multiply L*U, the {1 1} element of L times the {1 2} element of U gives us the 1 we need in the original A matrix, but as the {1 1} element of A gets smaller, eventually we get overflow in U.

The Wikipedia article discusses the conditions for an invertible matrix to have a unique LU decompostion, but says nothing about the uniqueness of a singular matrix's LU decomposition, if one exists.

The matrix B:

     [[ n 1 ]
      [ 0 0 ]]

has two decompositions (where U is a unit matrix) that I can see:

         P     *     L    *     U
     [[ 0 1 ]   [[ 0  0 ]  [[ 1  0 ]    
      [ 1 0 ]    [ n  1 ]]  [ 0  1 ]]
and

         P     *     L    *     U
     [[ 1 0 ]   [[ n  o ]  [[ 1  1/n ]    
      [ 0 1 ]    [ 0  0 ]]  [ 0   1  ]]

The second one has a problem when n -> 0. It's too bad the HP49G+ doesn't recognize the case where the matrix is singular and return a correct decomposition where possible.

                                          
Re: DET and LU calculations; HP-15C methods
Message #15 Posted by Karl Schneider on 21 July 2006, 2:56 a.m.,
in response to message #14 by Rodger Rosenbaum

Quote:
It seems to me that it is a matter of personal preference. I find that the "bold" method looks ok.

I find that it shouts at my eyeballs, particularly if the quoted text is lengthy. I prefer to use bold for extreme emphasis of a single word, or to emphasize a letter or character for notational purposes. So do most participants here, I'd say.

Quote:
I don't understand what you're getting at when you say, "The RPL-based models tried to compute general LU matrices, which contained overflow values in the U matrix."

Isn't the HP15's Doolittle LU decomposition "general"?


As I had stated, The HP-15C's L matrix is a unit lower-triangular matrix with ones on the main diagonal and also the lower-triangular elements normalized between -1 and +1.

The results of the "LU" operation on the HP-48/49 calculators appear not to condition the elements outside the zeroed portions, which can take any value, including zero.

Of course, the intent was different in the two applications, as you well know. The particular implementation of LU decomposition made the following possible on the RAM-short HP-15C:

  • Encoded L and U matrices that take no more space than the original matrix.
  • Inversion "in place", to the same original matrix.
  • Faster system solutions.

"LU" on the RPN-based models is simply a general-purpose function.

You are indeed correct about the HP-71B Math ROM. Page 86 of the manual (which I hadn't read thoroughly, and it doesn't have an index) states that DET, INV, and SYS perform the LU decomposition as an intermediary step, using the Crout factorization. I'm not sure how to access the LU results, or why the step is beneficial for DET, or why it is needed for INV if memory space is not an issue. Clearly, a straightforward DET calculation would have yielded zero without an "overflow warning" for the imperfect 2x2 matrix. [That was the original subject of this sub-thread, remember? :-)]

Quote:
The HP15 Advanced Functions Handbook explains those particulars that you found in Wikipedia which apply to the HP15. It also mentions that the HP15 uses the Doolittle method. And, it explains the rationale for keeping the LU decomposition in the RESULT matrix, etc., etc.

Yes, indeed. Verbiage on p. 97 of the AFH (which I acquired only three years ago, twenty years after I had bought my HP-15C) was the basis of my earlier statement, "...the Doolittle algorithm (used by the HP-15C)...".

I finally saw a description of the Doolittle method, well, yesterday in the Wikipedia article...

Quote:
I've known all this since I bought my HP15 over 20 years ago.

Everything you know and have, I wouldn't presume to know. The links and HP-15C material I posted was intended only as information for the entire readership. Given the constraints of the HP-15C, HP really did a remarkable job implementing practical matrix functionality.

I'd say that we've done a pretty good job of identifying the cause of the overflow messages and examining the various approaches taken the different models for this problem.

Best,

-- KS

                                                
Re: DET and LU calculations; HP-15C methods
Message #16 Posted by Rodger Rosenbaum on 21 July 2006, 7:21 a.m.,
in response to message #15 by Karl Schneider

Quote:
I prefer to use bold for extreme emphasis of a single word, or to emphasize a letter or character for notational purposes. So do most participants here, I'd say.

I started doing it because I saw other people doing it, so I know some other people like it, but I have no way of knowing what most prefer. And I wouldn't suggest that you abandon your preferred style because of an effect it had on my eyeballs.

Quote:
As I had stated, The HP-15C's L matrix is a unit lower-triangular matrix with ones on the main diagonal and also the lower-triangular elements normalized between -1 and +1.

The results of the "LU" operation on the HP-48/49 calculators appear not to condition the elements outside the zeroed portions, which can take any value, including zero.


I'll ask again, in different words. How does this make the HP15's LU decomposition less than general?

Quote:
Of course, the intent was different in the two applications, as you well know.

No, I don't well know it, because it isn't the whole truth. The three reasons you give for the use of the LU decomposition in the HP15 are not the only reasons; see my next comments.

Quote:
"LU" on the RPN-based models is simply a general-purpose function.

Quote:
I'm not sure how to access the LU results (from the HP71 Math Rom), or why the step is beneficial for DET, or why it is needed for INV if memory space is not an issue.

The LU decomposition is not just a general-purpose function on RPN models. The DET and INV functions on the HP48G use the LU decomposition. When solving the system Ax = B on the HP48G with B and A on the stack and then pressing the / key, the LU decomposition is used. The designers gave access to the decomposition, probably because they had a lot more ROM than the HP71 Math Pac did, even though you would almost never want to see it, except for educational purposes. Even had they not provided an LU keyword, the LU decompostion would still have been an important feature of the linear algebra in the HP48G, as it was in the HP71's Math Pac.

The LU decomposition is the preferred method used by modern numerical software, AFAIK, for computing the determinant and inverse, and for solving linear systems (other than over- and under-determined systems), because it is numerically favorable. See the books, "Matrix Computations", Golub & Van Loan, and "Matrix Algorithms", G. W. Stewart.

                                                      
Re: DET and LU calculations; HP-15C methods
Message #17 Posted by Karl Schneider on 22 July 2006, 12:31 a.m.,
in response to message #16 by Rodger Rosenbaum

Quote:
I'll ask again, in different words. How does this make the HP15's LU decomposition less than general?

The HP-15C's lower-diagonal L matrix isn't "less than general" -- it's a normalized unit lower-diagonal matrix, with ones on the main diagonal and all "lower" elements between -1 and +1, as stated. The extra computation that puts the resultant L into this form is both necessary and beneficial for subsequent computations.

I intended "general" to describe a form of result that not only meets the fundamental requirement(s), but also admits all valid possibilities. (An example would be the real-valued general solution of a second-order differential equation, including terms of exponential, sine, and cosine.) Perhaps the related word "generic" is a better fit.

Thus, I would assume that the "LU" function on the HP-48 and HP-49 attempt to produce a "generic" LU decomposition -- one that meets only the basic requirement defining triangularity. Perhaps those solutions -- made accessible by the "LU" function -- are the exact same ones produced by the same microcode for DET, INV, and system solution. Or, maybe not. I haven't delved into it.

Quote:
The LU decomposition is not just a general-purpose function on RPN models.

I think you meant "RPL" there. There is no "LU" function built into any RPN-based model.

As for the rest of your response, I'll take you at your word. It seems reasonable, although you'll have to admit that the LU decomposition made a simple problem:

det ([0 1 ; 0 0]) in Matlab syntax

rather cumbersome indeed!

-- KS

Edited: 22 July 2006, 2:13 a.m.

                                                            
Re: DET and LU calculations; HP-15C methods
Message #18 Posted by Rodger Rosenbaum on 22 July 2006, 6:10 a.m.,
in response to message #17 by Karl Schneider

Quote:
The HP-15C's lower-diagonal L matrix isn't "less than general" -- it's a normalized unit lower-diagonal matrix, with ones on the main diagonal and all "lower" elements between -1 and +1, as stated. The extra computation that puts the resultant L into this form is both necessary and beneficial for subsequent computations.

Why is this form necessary and beneficial? (I'm referring to the scaling of the elements below the diagonal in the L matrix to be between -1 and +1.) I know of no such benefit; the Advanced Functions Handbook doesn't seem to mention it. Am I overlooking some such discussion? (If you're just referring to the making of the L matrix into a unit matrix, then yes, that is a good thing.)

Quote:
I intended "general" to describe a form of result that not only meets the fundamental requirement(s), but also admits all valid possibilities.

The LU decomposition on the HP48 has a unit U matrix (see the AUR), rather than the unit L matrix on the HP15, a result of its use of the Crout algorithm.

As the Wikipedia article says, if we require the diagonal of the L (or U) matrix to be all ones, then the decompostion (of a non-singular matrix) is unique, and the phrase "all valid possibilities" isn't really germane; there's only one possibility. Since the U matrix of the HP48 decomposition is a unit matrix, there is only one decomposition (for non-singular matrices).

Quote:
Thus, I would assume that the "LU" function on the HP-48 and HP-49 attempt to produce a "generic" LU decomposition -- one that meets only the basic requirement defining triangularity.

This is not true, and I guess it's your assumption that it was true that led you to speak of the HP48 providing a "general" LU decomposition; I thought you knew.

Quote:
Quote:
"LU" on the RPN-based models is simply a general-purpose function.

I think you meant "RPL" there. There is no "LU" function built into any RPN-based model.

That's a direct quote from a previous posting of yours. I assumed it was a typo, but since I was cutting and pasting your quotes, I just left it alone. Your were referring to the HP48G/HP49/HP49G+ series calculators, right? I assumed so, since only they have an LU keyword.

It should be noted that the HP48S doesn't have an LU keyword, only the HP48G/G+/GX and HP49G/G+, so sometimes when Karl or I was careless and only said HP48, you know what we meant!

                                                                  
Re: DET and LU calculations; HP-15C & HP-48G methods
Message #19 Posted by Karl Schneider on 22 July 2006, 10:44 a.m.,
in response to message #18 by Rodger Rosenbaum

Quote:
Why is this form necessary and beneficial? (I'm referring to the scaling of the elements below the diagonal in the L matrix to be between -1 and +1.) ... If you're just referring to the making of the L matrix into a unit matrix, then yes, that is a good thing.)

Yes, I was primarily referring to the unit-matrix aspect of the format of L. This simplifies subsequent calculations, and also enables the L and U matrices to be packed into a matrix of the same size as the original, which is absoutely essential from a practical standpoint.

As for the normalization of the below-diagonal elements, I also did not see any specific discussion of its purpose in the AFH or the May 1983 H-P Journal article. I assume that it doesn't just "fall out" from the decomposition and unitization, and that there indeed was a specific purpose for doing it. With all elements of L between -1 and +1, it might reduce possible rounding errors when computing Ly = b -- just speculation on my part.

(The HP-48G*/49G* "LU" function apprarently does not normalize the above-diagonal elements in the unit upper-triangular (U) matrix. Decomposing the 3x3 integer matrix on p. 78 of the HP-71 Math ROM manual yielded a U matrix with one value above +1. That matrix is [5 3 2 ; 7 1 3 ; 6 4 9] by rows.)

Quote:
The LU decomposition on the HP48 has a unit U matrix (see the AUR), rather than the unit L matrix on the HP15, a result of its use of the Crout algorithm.

As the Wikipedia article says, if we require the diagonal of the L (or U) matrix to be all ones, then the decompostion (of a non-singular matrix) is unique, and the phrase "all valid possibilities" isn't really germane; there's only one possibility. Since the U matrix of the HP48 decomposition is a unit matrix, there is only one decomposition (for non-singular matrices).


I have an HP-48G with user's manual but no AUR (Advanced User's Reference?), so I'll accept your points at your word.

Your remaining comments are also correct -- I meant "RPL-based" -- specifically in my case, the HP-48G and HP-49G, which are the only ones I have with the "LU" function.

-- KS

Edited: 23 July 2006, 5:59 p.m.


[ Return to Index | Top of Index ]

Go back to the main exhibit hall