(HP71B) Eigenvalues of Symmetric Real Matrix
|
11-26-2023, 05:16 PM
Post: #1
|
|||
|
|||
(HP71B) Eigenvalues of Symmetric Real Matrix
Algorithm from HP15C Advanced Function Handbook, page 125
This version returned both eigenvectors and eigenvalues of symmetric real matrix 10 DESTROY ALL @ OPTION BASE 1 @ INPUT "N=";N 20 DIM A(N,N),B(N,N),L(N,N),Q(N,N),I0(N,N) @ MAT I0=IDN @ MAT Q=I0 30 MAT INPUT A @ MAT L=TRN(A) @ MAT A=A+L @ MAT A=(.5)*A @ MAT L=A 40 DEF FNF(N,D)=N/(D+(1-2*(D<0))*SQRT(N*N+D*D)) ! = tan(atan(N/D)/2) 100 MAT B=ZER 110 FOR I=2 to N @ FOR J=1 to I-1 120 X=L(I,J) @ IF X THEN B(I,J)=FNF(FNF(2*X,L(I,I)-L(J,J)),1) 130 NEXT J @ NEXT I 140 MAT L=TRN(B) @ MAT B=B-L @ DISP "FNORM(B)=";FNORM(B) 150 MAT L=B+I0 @ MAT L=INV(L) @ MAT L=L+L @ MAT L=L-I0 @ MAT Q=Q*L 160 DISP "Q=" @ MAT DISP Q 170 MAT L=TRN(Q) @ MAT L=L*A @ MAT L=L*Q 180 DISP "L=" @ MAT DISP L >run N=3 A(1,1)? 0,1,2 A(2,1)? 1,2,3 A(3,1)? 2,3,4 Code: FNORM(B)= .605552349352 QT A Q = L → A = Q L QT Werner recently PM me on how this work, but I have no clue. Anyone knows how it work, or reference (besides AFH) to where algorithm comes from? |
|||
11-26-2023, 07:20 PM
Post: #2
|
|||
|
|||
RE: (HP71B) Eigenvalues of Symmetric Real Matrix
(11-26-2023 05:16 PM)Albert Chan Wrote: Anyone knows how it work, or reference (besides AFH) to where algorithm comes from? It looks like Jacobi eigenvalue algorithm to me: Quote:Set this equal to 0, and rearrange: Or maybe a variant of it. |
|||
11-26-2023, 07:41 PM
Post: #3
|
|||
|
|||
RE: (HP71B) Eigenvalues of Symmetric Real Matrix
(11-26-2023 07:20 PM)Thomas Klemm Wrote: It looks like Jacobi eigenvalue algorithm to me: Yes, some Jacobi variant, because Bij entries is tan(θ/2). Question is why tan(θ/2) work. |
|||
11-26-2023, 07:44 PM
(This post was last modified: 11-26-2023 07:45 PM by Werner.)
Post: #4
|
|||
|
|||
RE: (HP71B) Eigenvalues of Symmetric Real Matrix
No, it isn’t Jacobi.
It looks a bit like it, but it transforms the whole matrix in one go, slowly (but ever faster, and ultimately quadratically) converging to a diagonal matrix. I understand the Cayley transform Q=(I-B)*inv(I+B), but not why and how the skew-symmetric matrix B would make Q an estimate for the eigenvectors (and not eigen*values* as written on pg 149 of the AFH, a typo), and why it would converge at all. Jacobi itself is simple, just setting a single off-diagonal element to zero. The resulting formula is similar, but not the same - ![]() Cheers, Werner 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE |
|||
11-26-2023, 09:26 PM
(This post was last modified: 11-29-2023 12:13 PM by Albert Chan.)
Post: #5
|
|||
|
|||
RE: (HP71B) Eigenvalues of Symmetric Real Matrix
HP71B, version 2
1. Suggested by Werner, assume user inputed symmetric square matrix. Remove code: A=(A+AT)/2 2. Bij = tan(θ/2) replace with tan(θ)/2, (|θ|≤pi/4), which use 1 sqrt, instead of 2 This has a slight effect to convergence rate, but likely not noticeable. For OP example, it converged slightly better. Cas> float(eigenvals([[0,1,2],[1,2,3],[2,3,4]])) /* for reference */ [0., 6.87298334621, −0.872983346207] This is 6th iteration, with old Bij = tan(θ/2) Code: FNORM(B)= 7.60886994417E-10 This is 6th iteration, with new Bij = tan(θ)/2 Code: FNORM(B)= 8.7618796429E-11 3. Q correction: I = (I+B)-1 (I+B) 2 (I+B)-1 - I = (I+B)-1 - (I+B)-1 B = (I+B)-1 (I-B) HP71-MathPac-suppl_manual.pdf Wrote:MAT X=INV(A)*B directly computes the solution of the matrix equation A X=B instead of RHS cost the same as matrix inverse. (although in this case, LHS is cheap to compute too) 4. Combine "MAT L=TRN(Q) @ MAT L=L*A" to MAT L=TRN(Q)*A 10 DESTROY ALL @ OPTION BASE 1 @ INPUT "N=";N 20 DIM A(N,N),B(N,N),L(N,N),Q(N,N),I0(N,N) @ MAT I0=IDN @ MAT Q=I0 30 MAT INPUT A @ MAT L=A 40 DEF FNF(N,D)=N/(D+(1-2*(D<0))*SQRT(N*N+D*D)) ! = tan(atan(N/D)/2) 100 MAT B=ZER 110 FOR I=2 to N @ FOR J=1 to I-1 120 X=L(I,J) @ IF X THEN B(I,J)=FNF(2*X,L(I,I)-L(J,J))*.5 130 NEXT J @ NEXT I 140 MAT L=TRN(B) @ MAT B=B-L @ DISP "FNORM(B)=";FNORM(B) 150 MAT L=I0-B @ MAT B=I0+B @ MAT L=INV(B)*L @ MAT Q=Q*L 160 DISP "Q=" @ MAT DISP Q 170 MAT L=TRN(Q)*A @ MAT L=L*Q 180 DISP "L=" @ MAT DISP L Comment: this patch make bij closer to tan(θ/2), with little cost, if desired. Cas> series(tan(x/2),polynom) → 1/2*x+1/24*x^3+1/240*x^5 Cas> series(tan(x)/2,polynom) → 1/2*x+1/6*x^3+1/15*x^5 >120 X=L(I,J) @ IF X THEN X=FNF(2*X,L(I,I)-L(J,J))*.5 @ B(I,J)=X-X*X*X |
|||
11-26-2023, 10:05 PM
Post: #6
|
|||
|
|||
RE: (HP71B) Eigenvalues of Symmetric Real Matrix
(11-26-2023 05:16 PM)Albert Chan Wrote: 40 DEF FNF(N,D)=N/(D+(1-2*(D<0))*SQRT(N*N+D*D)) ! = tan(atan(N/D)/2) I might have jumped a step. The reason FNF = tan(atan(N/D)/2) is based on this identity. (05-31-2021 09:51 PM)Albert Chan Wrote: \(\displaystyle\arctan(x) = 2\arctan\left( {x \over \sqrt{1+x^2}+1} \right)\) x = N/D, but we don't want division by 0, so we scale D away. HP71B has signed zero, I was hoping SGN(0) = 1, but that returned 0. Interestingly, SGN(-0) returned -0 That's why the (1-2*(D<0)) factor, to handle D=0 case --> 1-2*0 = 1 Another way is SGN(D+(D=0)), but probably is just as hard to read. atan(N/D) = 2 atan(N/(D+(1-2*(D<0))*sqrt(N*N+D*D))) Divide by 2, then tan() to both sides, we get FNF |
|||
11-27-2023, 07:01 AM
Post: #7
|
|||
|
|||
RE: (HP71B) Eigenvalues of Symmetric Real Matrix
(11-26-2023 07:41 PM)Albert Chan Wrote: Yes, some Jacobi variant, because Bij entries is tan(θ/2). Question is why tan(θ/2) work. No, the Bij entries are tan(θ/4). The Jacobi formula for determining the Givens rotation is tan(θ/2), and the sign is reversed as well. Cheers, Werner 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE |
|||
11-27-2023, 08:07 PM
Post: #8
|
|||
|
|||
RE: (HP71B) Eigenvalues of Symmetric Real Matrix
Hi, Werner
Sorry for the confusion. I am basing on bij = tan( 1/4 * (2θ = atan(2*aij / (aii - ajj))) ) , i ≠ j, aij ≠ 0 Except edge cases, bij = tan(θ/2) , |θ/2| ≤ pi/8 ≈ 0.3927 From https://en.wikipedia.org/wiki/Cayley_transform#Examples, this may be why AFH code work. LHS matches AFH code of half-angle tan, reversed sign, and zero diagonal. RHS = Cayley transformed LHS, matches rotation matrix. Perhaps same idea can be extended from 2x2 to nxn matrix? |
|||
11-27-2023, 08:51 PM
Post: #9
|
|||
|
|||
RE: (HP71B) Eigenvalues of Symmetric Real Matrix
Aha! Indeed, that is a first step in the right direction!
Werner 41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: 2 Guest(s)