Compute some Eigenvalues first thing in the morning, and nothing worse will happen to you all day.



#20

[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 »


#21

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? :-)


#22

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


#23

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

#24

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 »

#25

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.


#26

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.


#27

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.


#28

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.


#29

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


#30

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.


#31

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

#32

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.


#33

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.


#34

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


#35

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.


#36

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.


#37

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!


#38

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.


Possibly Related Threads...
Thread Author Replies Views Last Post
  An amazing day: Giving a talk at HP about their calculators Geir Isene 9 1,600 12-16-2013, 06:14 PM
Last Post: aurelio
  Using the Prime to solve for eigenvalues Michael de Estrada 28 3,076 10-27-2013, 07:21 AM
Last Post: Tarcisi C
  HP Prime - One thing I can not draw in advanced graphing dg1969 2 526 10-19-2013, 02:00 PM
Last Post: dg1969
  HHC 2013 Day 2 Highlights Eddie W. Shore 6 953 09-23-2013, 04:03 PM
Last Post: Kimberly Thompson
  HHC 2013: Day 1 Highlights Eddie W. Shore 28 2,616 09-23-2013, 03:22 PM
Last Post: Brad Barton
  One thing I really love on HP PRIME... dg1969 0 293 07-31-2013, 06:50 AM
Last Post: dg1969
  HP 50g Takes hours to compute a convergent series Matthew Richards 10 1,086 05-14-2013, 03:30 PM
Last Post: Thomas Klemm
  Happy Mother's Day! Eddie W. Shore 1 396 05-12-2013, 11:35 AM
Last Post: Walter B
  happy fibonacci day 5/8/13 Allen 8 1,104 05-09-2013, 01:48 AM
Last Post: Gerson W. Barbosa
  HP on Saturday Morning Eddie W. Shore 8 816 04-16-2013, 11:16 AM
Last Post: Steve B (Canada)

Forum Jump: