▼
Posts: 2,309
Threads: 116
Joined: Jun 2005
[From my blog]
Actually it is fun, though a bit tedious. My HP49G+ 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 builtin 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 »
▼
Posts: 305
Threads: 17
Joined: Jun 2007
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? :)
▼
Posts: 2,309
Threads: 116
Joined: Jun 2005
I did compute some Eigenvalues, by hand, far too early this morning. If you really want to see them, I can type them up. :)
▼
Posts: 305
Threads: 17
Joined: Jun 2007
It's not so much the eigenvalues that I want to see; it's the program that computes them.
Posts: 2,309
Threads: 116
Joined: Jun 2005
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 »
Posts: 1,755
Threads: 112
Joined: Jan 2005
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 "HP71B  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 5liner
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=N1
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=AU
4 C(I)=DET(A)
@ FOR J=1 TO N
@ D(I,NJ)=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)= x^{5} + 19 x^{4}  79 x^{3}  146 x^{2} + 1153 x  1222
and, as expected in this case, as this is a symmetric matrix,
all five eigenvalues are real (imaginary parts = 0):
x_{1} = 1.49766
x_{2} = 3.36188
x_{3} = 3.55784
x_{4} = 5.67255
x_{5} = 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.
▼
Posts: 305
Threads: 17
Joined: Jun 2007
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.
▼
Posts: 1,792
Threads: 62
Joined: Jan 2005
Rodger 
On my HP71B with Math ROM, the message, "WRN:Overflow" appeared briefly, and then the correct answer of zero appeared. Subsequently, I did a "DESTROY ALL", reentered the matrix, and did "DET" again. No message accompanied the answer of zero. Weird...
It might be that the HP71B Math ROM uses the same basic algorithms as did the HP15C, which introduced the builtin matrix functionality. However, I don't know for certain.
Th HP15C performs an LU decomposition of the matrix prior to calculating the determinant, as described in the HP15C 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 HP15C 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 HP71B 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 HP15C, which naturally returns zero without any warning or error code.
This calculation may be the source of the warning message from the HP71B Math Pac.
 KS
Edited: 17 July 2006, 1:53 a.m.
▼
Posts: 305
Threads: 17
Joined: Jun 2007
Karl said:
The encoded resultant LU matrix calculated by the HP15C 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 HP71B 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 HP15C, 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 resultthe 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.
▼
Posts: 163
Threads: 7
Joined: Jul 2007
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
▼
Posts: 305
Threads: 17
Joined: Jun 2007
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 5by5 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.
▼
Posts: 163
Threads: 7
Joined: Jul 2007
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
Posts: 1,792
Threads: 62
Joined: Jan 2005
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 HP71B, for the DET operation, is doing the same mathematics for the LU decomposition as the two RPLbased models you tried. (I got the same result with my HP49G.)
The HP15C calculates a unit L matrix (which contains ones on the main diagonal). The RPLbased 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 HP15C) 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 computationallyslow HP15C performs the LU decomposition for linearsystem solving: It saves time for subsequent solutions, and is a bit faster yet when L or U is a unit matrix.
The HP15C 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 exactlydetermined 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.
▼
Posts: 305
Threads: 17
Joined: Jun 2007
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 RPLbased 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 10^{450} ]
[ 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.
▼
Posts: 1,792
Threads: 62
Joined: Jan 2005
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 RPLbased 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 HP15C's L matrix is a unit lowertriangular matrix with ones on the main diagonal and also the lowertriangular elements normalized between 1 and +1.
The results of the "LU" operation on the HP48/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 RAMshort HP15C:
 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 RPNbased models is simply a generalpurpose function.
You are indeed correct about the HP71B 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 subthread, 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 HP15C) was the basis of my earlier statement, "...the Doolittle algorithm (used by the HP15C)...".
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 HP15C material I posted was intended only as information for the entire readership. Given the constraints of the HP15C, 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
▼
Posts: 305
Threads: 17
Joined: Jun 2007
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 HP15C's L matrix is a unit lowertriangular matrix with ones on the main diagonal and also the lowertriangular elements normalized between 1 and +1.
The results of the "LU" operation on the HP48/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 RPNbased models is simply a generalpurpose 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 generalpurpose 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 underdetermined systems), because it is numerically favorable. See the books, "Matrix Computations", Golub & Van Loan, and "Matrix Algorithms", G. W. Stewart.
▼
Posts: 1,792
Threads: 62
Joined: Jan 2005
Quote:
I'll ask again, in different words. How does this make the HP15's LU decomposition less than general?
The HP15C's lowerdiagonal L matrix isn't "less than general"  it's a normalized unit lowerdiagonal 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 realvalued general solution of a secondorder 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 HP48 and HP49 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 generalpurpose function on RPN models.
I think you meant "RPL" there. There is no "LU" function built into any RPNbased 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.
▼
Posts: 305
Threads: 17
Joined: Jun 2007
Quote: The HP15C's lowerdiagonal L matrix isn't "less than general"  it's a normalized unit lowerdiagonal 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 nonsingular 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 nonsingular matrices).
Quote: Thus, I would assume that the "LU" function on the HP48 and HP49 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 RPNbased models is simply a generalpurpose function.
I think you meant "RPL" there. There is no "LU" function built into any RPNbased 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!
▼
Posts: 1,792
Threads: 62
Joined: Jan 2005
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 unitmatrix 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 belowdiagonal elements, I also did not see any specific discussion of its purpose in the AFH or the May 1983 HP 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 HP48G*/49G* "LU" function apprarently does not normalize the abovediagonal elements in the unit uppertriangular (U) matrix. Decomposing the 3x3 integer matrix on p. 78 of the HP71 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 nonsingular 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 nonsingular matrices).
I have an HP48G 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 "RPLbased"  specifically in my case, the HP48G and HP49G, which are the only ones I have with the "LU" function.
 KS
Edited: 23 July 2006, 5:59 p.m.
