(DM42) Matrix exponential
|
08-11-2023, 09:56 PM
(This post was last modified: 08-12-2023 11:24 AM by Gjermund Skailand.)
Post: #1
|
|||
|
|||
(DM42) Matrix exponential
updated 12/8
A matrix exponential function for square matrix for the DM42, Free42 This is just a simple implementation based argument reduction and series expansion. Exp(M) = I + M + 1/(2!)M^2 + 1/(3!)M^3... 1(n!)M^n I is identity matrix The argument reduction reduces the fnrm to below 1.0 prior to the series expansion. Note that running times will increase for arrays with large inputs. Program runs in 4 stack mode. Flag 24 and 25 shall be cleared. Uses registers R00-R02 R00: counter number of convergence iterations R01: stopping criteria is difference between row norm for two iterations to be less than value set in line 30. Test examples use 1e-32, but a more sensible value is probably 1e-16. R02: scratch, used for argument reduction and squaring Code:
example 1 M = [[ 1 4 ] [ 1 1 ]] MEXP= [[ 10.226708, 19.717657] [ 4.929414, 10.226708]] exact (Matrix exponential - Wikipedia): [[(exp(4)+1)/(2e), (exp(4)-1)/e ] [(exp(4)-1)/(4e), (exp(4)+1)/(2e)]] difference between MEXP(M)- exact = [[ 3.0E-32, 6.0E-32] [ 1.9E-32, 3.0E-32 ]] example 2, calculated in about 1.6 s when running on batteries M = [[1 2 3 ] [4 5 6] [7 8 9]] MEXP(M) = [[ 1.118907e6, 1.374815e6, 1.630724e6] [ 2.533881e6, 3.113415e6, 3.692947e6] [ 3.948856e6, 4.852013e6, 5.755171e6]] ps It might be more sensible to use R1=1E-16, as it cuts down calculation time almost in half, and if in doubt about the result, do a second run with R1=1E-32 to compare the first 16 digits. Best regards Gjermund |
|||
08-11-2023, 11:46 PM
Post: #2
|
|||
|
|||
RE: (DM42) Matrix exponential
I tried for the HP50G a very simple approach :
exp(M) = {(I+M/9000/9000)^9000} ^ 9000 (twice 9000, as powers for n>= 10000 are not allowed in HP50G). However the results with the above matrixes are not so accurate. Is there a simple approach to get better results? Thanks for your help. Regards, Gil |
|||
08-12-2023, 09:27 AM
Post: #3
|
|||
|
|||
RE: (DM42) Matrix exponential
Hi Gil
Thanks for your interest. This method for calculation the matrix exponential with power series expansion is a summation of many terms until convergence, not only 3 first terms. In addition there is a significant risk of overflow in intermediate terms which can be much larger than the converged results. Thus the halving calculations in lines 04 to 17, and repeated squaring (by matrix multiplication) lines 58 to 63 of matrix after convergence is achieved. DM42 and Free42 allows for numbers range up to about +/-1E6145, HP50g to about +/-1E499. Unfortunately I am not able catch overflow error on the DM42 for matrix multiplication (not sure whats happening). I can only get it for element-wise operations, and have updated the program with a final test so overflows are likely catched when flag 24 is cleared. The update is 63 LBL 04 64 FNRM 65 LAST 66 RTN 67 END In addition I will add instructions that for the DM42 or Free42 flag 24 should be cleared to catch out of range errors. On the HP50G: if you are writing a user-rpl program, be sure to set system flag -21 "Overflow --> error" and clear system flag -22 "Infinite --> error". There is also a HPGCC2 implementation for matrix exponential (probably more accurate and 50-100 times faster than any simple user-rpl implementation) at https://www.hpcalc.org/hp49/math/numeric/utilgcc.zip If you are willing to do a rom update to enable HPGCC3 programs, I have an even faster update. Send me a message if that is interesting. For very simple testing, you can use a square matrix with entries only at the diagonal, then (and only then) the matrix exponential is equal to elementwise exponention of diagonal elements. For testing on HP50G, you can e.g. use the EIG function on (symmetric) matrix on the HP50G to diagonalize the matrix. See e.g. Wikipedia on matrix exponential for limitations and how to do that. Best regards Gjermund |
|||
08-12-2023, 10:01 AM
(This post was last modified: 08-12-2023 10:34 AM by Gil.)
Post: #4
|
|||
|
|||
RE: (DM42) Matrix exponential
Thanks for your explanations.
As suggested, I used the EGV command and it works nice. \<< EGV DUP SIZE OBJ\-> DROP \-> ev d \<< DUP INV d IDN 1 d FOR i { i i } ev i GET EXP PUT NEXT SWAP * * \->NUM \>> \>> Or « EGV DUP SIZE OBJ—> DROP —> ev d « DUP INV d IDN 1 d FOR i { i i } ev i GET EXP PUT NEXT SWAP * * —>NUM » » With [[ 1. 2. 3. ] [ 4. 5. 6. ] [ 7. 8. 9. ]], we get: [[ 1118906.69941 1374815.06293 1630724.42645 ] [ 2533881.04188 3113415.03136 3692947.02084 ] [ 3948856.38436 4852012.99979 5755170.61523 ]]. By the way, if I get 3 values 10, 20 and 30, are there commands to build a Matrix [10 0 0 0 20 0 0 0 30], without using the above loop n (n=3)? Regards, Gil |
|||
08-12-2023, 11:38 AM
Post: #5
|
|||
|
|||
RE: (DM42) Matrix exponential
Thanks for confirming/testing by diagonalizing.
Regarding the HP50G, you could do [ 10. 20. 30. ] 3. DIAG--> DIAG--> is found in the Math/Matrix/Make menu On a sidenote, when working with vectors, it may be tempting to use -->V3 (and -->V2) to build vectors, be aware it depends on the Rect/Cyl/Sphere mode. Best regards Gjermund |
|||
08-12-2023, 08:26 PM
Post: #6
|
|||
|
|||
RE: (DM42) Matrix exponential
Thanks, Gjermund.
I changed the program, using "your" DIAG—> command instead of an Identity Matrix. Depending on the size of the Matrix, the execution time now might be much shorter than the previous version. \<< "1 Arg: MAT[n x n]" DROP EGV DUP SIZE OBJ\-> DROP \-> ev d \<< DUP INV 1 d FOR i ev i GET EXP NEXT d \->ARRY d DIAG\-> SWAP * * \->NUM \>> \>> Or « "1 Arg: MAT[n x n]" DROP EGV DUP SIZE OBJ—> DROP —> ev d « DUP INV 1 d FOR i ev i GET EXP NEXT d —>ARRY d DIAG—> SWAP * * —>NUM » » Gil |
|||
08-12-2023, 08:55 PM
(This post was last modified: 08-13-2023 12:43 AM by Gil.)
Post: #7
|
|||
|
|||
RE: (DM42) Matrix exponential
By the way, Gjermund, could you check on your DM42 program the result of EXP (M),
with M [1 2 3 4 5 6 7 8 9 10 11 12 13... 18 19... 24 25... 30 31... 36] ? Element line 6 / col 1 of EXP (M) is with HP50G: 8.87064557977E49 (WolframAlpha rounds it to 8.87064E49), as the result with Maxima is about 8.870639 7724675632217598059059569391342349936428931E49). Thanks in advance for your cooperation. Regards, Gil |
|||
08-13-2023, 10:23 AM
Post: #8
|
|||
|
|||
RE: (DM42) Matrix exponential
Hi Gil
as requested DM42 R01=1e-16 in about 2.5s --> 8.870639772467562572805300897324555e49 DM42 R01=1e-32 in about 3.9s --> 8.870639772467563221759805905956961e49 maxima 8.870639772467563221759805905956939e49 HP50g HPGCC3 8.87063977247E49 (result on real HP50g in HPGCC3 implementation running at 120MHz using standard C double floating point is almost instantly) Best regards Gjermund |
|||
08-13-2023, 10:51 AM
(This post was last modified: 08-13-2023 10:56 AM by Gil.)
Post: #9
|
|||
|
|||
RE: (DM42) Matrix exponential
Amazing how accurate reckons the DM42 with your program — and, of course, the special version for HP50G.
By the way, with my program on emu48 with SamsungA53 phone, the result for exp (M, 6x6) is given in about 0.63 s. Many thanks for your painstaking, Gjermund. Regards, Gil |
|||
08-13-2023, 09:41 PM
Post: #10
|
|||
|
|||
RE: (DM42) Matrix exponential
Hi Gil,
I corrected some bad logic for small values and also done some time-optimizing For convergency test 1e-16, i.e. about 16 correct digits: DM42 on USB calculates the example in 0.67 s DM42 on battery 1.76 s Actually setting convergency test to 1e-12, and R03=8, I can run calculations in 0.6 s on USB power and get 12 correct digits I think it is worth the increased program size. Code:
Best regards Gjermund |
|||
08-13-2023, 09:46 PM
Post: #11
|
|||
|
|||
RE: (DM42) Matrix exponential
I am afraid of being lightning years beyond that level of programming and understanding.
Impressive anyway. Regards, Gil |
|||
08-15-2023, 11:42 PM
(This post was last modified: 08-15-2023 11:50 PM by Gil.)
Post: #12
|
|||
|
|||
RE: (DM42) Matrix exponential
Help, please Gjermund!
Results OK, in exact mode, with my program on emu48 for HP50G for EXP(M7), with M7 = [01. 02. 03. ... 06. 07. 08. ... 14. ... 43. ... 49.], for Matrix with reals 1. 2. ... 49. (with dots at the end). Results not OK, in exact mode, with my program on emu48 for HP50G for EXP(M8), with M8= [01. 02. 03. ... 07. 08. 08. ... 16. ... 57. ... 64.], for Matrix with reals 1. 2. ... 64. (with dots at the end). The result given here is a complex matrix! Two questions: 1) How is this impossible result to be calculated? 2) Do you have this kind of problem on your DM42 and HP50G programs? To get a real answer, I have to delete the dots and enter M8 as [01 02 03 ... 07 08 08 ... 16 ... 57 ... 64] (only integers, without dots). But impossible to calculate EXP (M9), with M9 containing only integers, as not enough memory. Thanks for your help or insight. Regards, Gil . |
|||
08-16-2023, 12:01 PM
(This post was last modified: 08-16-2023 12:07 PM by John Keith.)
Post: #13
|
|||
|
|||
RE: (DM42) Matrix exponential
(08-15-2023 11:42 PM)Gil Wrote: To get a real answer, I have to delete the dots and enter A couple of notes: A 64 x 64 matrix is approaching the limits of memory for the HP 50, especially if the calculation involves other matrices. You don't need to delete the dots to make the matrix exact, just use :: R\->I MAP or \->Q\pi. I am following this thread with great interest. I have experimented with general programs to apply functions to matrices using SVD (singular value decomposition) but the results have been disappointing, especially for EXP and LN. I am hoping that the two of you (and others!) will come up with a good solution. |
|||
08-16-2023, 12:45 PM
Post: #14
|
|||
|
|||
RE: (DM42) Matrix exponential
Thanks for the hint.
What disturbed me was not the conversion from real to integer, but the fact that of getting an impossible complex result for exp (M), with elements of M being real, the decomposition giving already a complex matrix. I am afraid that I can't be of any further help, as my math and programming level is quite a basically one. Calculations with 64×64 Matrixes is impressive anyway, as with the memory content of the my HP50G with EMU48 I hardly can calculate 8×8 or 9×9 Matrixes. Regards, Gil |
|||
08-16-2023, 04:38 PM
Post: #15
|
|||
|
|||
RE: (DM42) Matrix exponential
Dear Gil,
Short answer: there is a nice article "Nineteen dubious ways to calcate the exponential of a matrix, tventyfive years later." Longer answer, mostly incorrect your method for nxn matrix is a bit like solving a n-degree equation, then exponention, and final back multiplication. In the middle step, for unsymmetric matrixes you may end up taking the square root of negative numbers, resulting in imaginary result. Hope this helps best regards Gjermund |
|||
08-16-2023, 09:48 PM
Post: #16
|
|||
|
|||
RE: (DM42) Matrix exponential
(08-12-2023 09:27 AM)Gjermund Skailand Wrote: there is a significant risk of overflow in intermediate terms which can be much larger than the converged results. Perhaps we could do z = expm1(x/2^n), so that overflow is not an issue. Instead, we have underflow issue, but it just flush to a safe zero. expm1(2ε) = (expm1(ε)+1)^2 - 1 = expm1(ε) * (expm1(ε)+2) OP Example: XCas> X := list2mat(range(1,10), 3) XCas> N := max(0, round(logb(maxnorm(X),2)+10)) → 13 XCas> Y := X / 2.^N /* maxnorm 1E-3 or less */ XCas> Z := Y*(1 + Y/2*(1 + Y/3*(1 + Y/4*(1 + Y/5*(1 + Y/6))))) /* Z = expm1(X/2^N) */ XCas> for(k=0; k<N; k++) {Z *= Z+2;} /* Z = expm1(X) */ \(\left(\begin{array}{ccc}1118905.69941&1374815.06293&1630724.42646\\ 2533881.0419&3113414.03138&3692947.02086\\ 3948856.38438&4852012.99982&5755169.61526 \end{array}\right)\) |
|||
08-20-2023, 02:45 PM
(This post was last modified: 08-20-2023 06:57 PM by Albert Chan.)
Post: #17
|
|||
|
|||
RE: (DM42) Matrix exponential
(08-16-2023 09:48 PM)Albert Chan Wrote: Perhaps we could do z = expm1(x/2^n), so that overflow is not an issue. On decimal machine, we like z = expm1(x/10^n), to remove scaling errors. \( (e^x-1)(e^y-1) = e^x e^y - e^x - e^y + 1 = (e^x e^y-1) - (e^x-1) - (e^y-1) \) expm1(x+y) = expm1(x) * expm1(y) + expm1(x) + expm1(y) Using this identity, below code does expm1(log1p(r)*n) = (1+r)^n - 1 This is an important term used in TVM calculations. (Formula arranged to produce Plus42 menus order: N, I, PV, PMT, FV, BEGIN) TVM:0=(1/EXPM1(N*LNP1(I))*(PV+0*PMT+FV)+PV+BEGIN*PMT)*I+PMT Code: function f(r, n) -- expm1(log1p(r)*n), integer n>=0 lua> r = 0.01 lua> f(r,0), f(r,1), f(r,2), f(r,3) 0 0.01 0.0201 0.030301 expm1(log1p(r)) = r, we can build from expm1(r) to expm1(r * n) lua> n = 100 lua> f(expm1(r), n) 1.7182818284590453 lua> expm1(r * n) 1.7182818284590453 |
|||
08-21-2023, 12:11 AM
(This post was last modified: 08-21-2023 12:57 AM by Albert Chan.)
Post: #18
|
|||
|
|||
RE: (DM42) Matrix exponential
EXPM1(Matrix X) for HP71B
10 DESTROY ALL @ OPTION BASE 1 20 INPUT "N=";N @ DIM X(N,N),T(N,N),R(N,N) 30 MAT INPUT X 40 S=10^MAX(0,EXPONENT(MAXAB(X))+4) @ MAT X=(1/S)*X 50 MAT R=(1/4)*X @ MAT T=IDN @ MAT R=R+T 60 MAT R=(1/3)*R @ MAT R=R*X @ MAT R=R+T 70 MAT R=(1/2)*R @ MAT R=R*X @ MAT R=R+T 80 MAT R=R*X @ MAT X=R 90 DISP FNF(S) @ DISP @ MAT DISP R @ END 200 DEF FNF(K) @ IF K<=1 THEN 250 210 DISP FNF(K DIV 2); 220 MAT T=R*R @ MAT R=R+R @ MAT R=R+T 230 IF MOD(K,2)=0 THEN 250 240 MAT T=R*X @ MAT T=T+X @ MAT R=R+T 250 FNF=K @ END DEF >RUN N=3 X(1,1)? 1 X(1,2)? 2 X(1,3)? 3 X(2,1)? 4 X(2,2)? 5 X(2,3)? 6 X(3,1)? 7 X(3,2)? 8 X(3,3)? 9 1 2 4 9 19 39 78 156 312 625 1250 2500 5000 10000 1118905.69939 1374815.06290 1630724.42641 2533881.04184 3113414.03128 3692947.02074 3948856.38427 4852012.99966 5755169.61507 |
|||
08-21-2023, 11:49 AM
Post: #19
|
|||
|
|||
RE: (DM42) Matrix exponential
I might have over-complicated HP71B code with 10^P scalings.
2^P may be good enough. (5^17<1E12, 1/2^17 to 12 decimals is exact) With below code, P>17 only if MAXAB(X) ≥ 2^7.5 ≈ 181 Gain in 10^P scaling down is offset by more computation to then scale back up. With 2^P, we just do P squarings, FNF(K) code not needed. 10 DESTROY ALL @ OPTION BASE 1 20 INPUT "N=";N @ DIM X(N,N),T(N,N),R(N,N) 30 MAT INPUT X 40 P=MAX(0,IROUND(LOG2(MAXAB(X)))+10) @ MAT X=(1/2^P)*X 50 MAT R=(1/4)*X @ MAT T=IDN @ MAT R=R+T 60 MAT R=(1/3)*R @ MAT R=R*X @ MAT R=R+T 70 MAT R=(1/2)*R @ MAT R=R*X @ MAT R=R+T @ MAT R=R*X 80 FOR I=1 TO P @ MAT T=R*R @ MAT R=R+R @ MAT R=R+T @ I; @ NEXT I 90 @@ MAT DISP R @ END >run N=3 X(1,1)? 1 X(1,2)? 2 X(1,3)? 3 X(2,1)? 4 X(2,2)? 5 X(2,3)? 6 X(3,1)? 7 X(3,2)? 8 X(3,3)? 9 1 2 3 4 5 6 7 8 9 10 11 12 13 1118905.69945 1374815.06297 1630724.42651 2533881.04197 3113414.03144 3692947.02095 3948856.38452 4852012.99997 5755169.61548 (08-16-2023 04:38 PM)Gjermund Skailand Wrote: there is a nice article A nasty example from the article, less squaring computations might win over perfect scaling. XCas> approx(EXPM1([[-49,24], [-64, 31]])) /* reference */ -1.73575875814 0.551819099658 -1.47151759909 0.103638240716 HP71B, SCALE BY 10^5 (≈ 2^(5*3.322) = 2^16.6) -1.73575875912 0.551819100440 -1.47151760119 0.103638242380 HP71B, SCALE BY 2^16 or 2^17 -1.73575875814 0.551819099660 -1.47151759909 0.103638240710 |
|||
08-21-2023, 07:03 PM
Post: #20
|
|||
|
|||
RE: (DM42) Matrix exponential
Hi Albert
note the below matrix is row oriented, which c, HP50g, DM42 calculators are using (fortran is column, I don't know about xcas, but I get different result than you) For the 2x2 example , your reference A= [[ -49, 24][-64 ,31 ] ] = [[1, 3] [2, 4]] × [[ -1, 0] [0, -17]] x INV[[ 1, 3][2, 4]] exp(A)=[[1, 3] [2, 4]] × [[ exp(-1) , 0] [0, exp(-17)]] x INV[[ 1, 3][2, 4]] = [[ -0.735759... , 0.551819... ] [-1.471518..., 1.103638...]] I can get full accuracy (31- 32 digits) on DM42. But my implementation on HP50g using double binary floating point fails. br Gjermund |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: