+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (HP71B) Eigenvalues of Symmetric Real Matrix (/thread-20906.html)
(HP71B) Eigenvalues of Symmetric Real Matrix - Albert Chan - 11-26-2023 05:16 PM
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
Q=
.8662568603 .227657549016 .444714619016
-.444714619012 .75699314672 .478738220168
-.227657549018 -.612481359872 .75699314673
L=
-.349003969447 -.425695598709 -1.12609254611
-.425695598709 -.348327850394 -4.47666317212E-2
-1.12609254611 -4.47666317247E-2 6.69733181992
>run 100
FNORM(B)= .595643959403
Q=
.824769255982 -.469983482558 .314437912022
.285388154078 .826016169657 .486056466843
-.488169310695 -.311147575236 .815400460218
L=
-.859525935837 9.56986605065E-2 .200156808749
9.56986605065E-2 1.82880814699E-2 .42420712414
.200156808747 .424207124136 6.84123785436
>run 100
FNORM(B)= 8.95738317424E-2
Q=
.86332901543 -.398599604012 .309485648795
.183632804446 .819374127488 .543051592685
-.470044683262 -.412000479214 .78058542189
L=
-.872825507534 -9.30648141545E-3 -2.14243725402E-2
-9.30648141545E-3 -8.11920389207E-5 1.09063583269E-2
-2.14243725373E-2 .010906358327 6.87290669953
>run 100
FNORM(B)= 7.86910670518E-3
Q=
.859879422537 -.408274692867 .306462320485
.193856724533 .816499388716 .543827471319
-.472257291179 -.408216270356 .78123781753
L=
-.872983344667 3.2546038359E-5 5.1700251917E-5
3.2546038358E-5 1.1264547796E-9 1.26811006234E-4
5.170025078E-5 1.2681100214E-4 6.87298334337
>run 100
FNORM(B)= 2.97899107959E-5
Q=
.859892597624 -.408248289526 .30646052708
.193822654531 .816496581206 .54384383001
-.472247286024 -.408248290825 .781227133338
L=
-.872983346221 -8.99996E-10 -2.063351E-9
-8.99996266439E-10 -7.5924E-19 1.07648489092E-9
-2.05974E-9 1.07239E-9 6.87298334594
>run 100
FNORM(B)= 7.60886994417E-10
Q=
.859892597285 -.40824829046 .306460526788
.193822655517 .816496580921 .543843830086
-.472247286237 -.40824829046 .7812271334
L=
-.872983346221 -1.69E-13 -3.663E-12
-1.69203950392E-13 2.E-24 3.73782958716E-12
1.66E-12 3.73E-12 6.87298334595QT 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?
RE: (HP71B) Eigenvalues of Symmetric Real Matrix - Thomas Klemm - 11-26-2023 07:20 PM
(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:
\(
\tan(2\theta) = \frac{2 S_{ij}}{S_{jj} - S_{ii}}
\)
Or maybe a variant of it.
RE: (HP71B) Eigenvalues of Symmetric Real Matrix - Albert Chan - 11-26-2023 07:41 PM
(11-26-2023 07:20 PM)Thomas Klemm Wrote: It looks like Jacobi eigenvalue algorithm to me:
Quote:Set this equal to 0, and rearrange:
\(
\tan(2\theta) = \frac{2 S_{ij}}{S_{jj} - S_{ii}}
\)
Or maybe a variant of it.
Yes, some Jacobi variant, because Bij entries is tan(θ/2). Question is why tan(θ/2) work.
RE: (HP71B) Eigenvalues of Symmetric Real Matrix - Werner - 11-26-2023 07:44 PM
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
RE: (HP71B) Eigenvalues of Symmetric Real Matrix - Albert Chan - 11-26-2023 09:26 PM
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
Q=
.859892597285 -.40824829046 .306460526788
.193822655517 .816496580921 .543843830086
-.472247286237 -.40824829046 .7812271334
L=
-.872983346221 -1.69E-13 -3.663E-12
-1.69203950392E-13 2.E-24 3.73782958716E-12
1.66E-12 3.73E-12 6.87298334595This is 6th iteration, with new Bij = tan(θ)/2
Code:
FNORM(B)= 8.7618796429E-11
Q=
.408248290464 .859892597282 .306460526793
-.816496580928 .193822655521 .543843830107
.408248290464 -.472247286237 .78122713341
L=
0 0 0
-1.E-15 -.872983346218 1.786E-12
0 3.53E-12 6.872983346243. 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
inverting A then multiplying by B
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
RE: (HP71B) Eigenvalues of Symmetric Real Matrix - Albert Chan - 11-26-2023 10:05 PM
(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
RE: (HP71B) Eigenvalues of Symmetric Real Matrix - Werner - 11-27-2023 07:01 AM
(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
RE: (HP71B) Eigenvalues of Symmetric Real Matrix - Albert Chan - 11-27-2023 08:07 PM
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.
![[Image: 104110a11efbaa15730fbbf044b41db512bf63b9]](https://wikimedia.org/api/rest_v1/media/math/render/svg/104110a11efbaa15730fbbf044b41db512bf63b9)
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?
RE: (HP71B) Eigenvalues of Symmetric Real Matrix - Werner - 11-27-2023 08:51 PM
Aha! Indeed, that is a first step in the right direction!
Werner