(DM42) Matrix exponential - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: General Software Library (/forum-13.html) +--- Thread: (DM42) Matrix exponential (/thread-20290.html) Pages: 1 2 (DM42) Matrix exponential - Gjermund Skailand - 08-11-2023 09:56 PM 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:  @ Matrix exponential by series expansion.  @ scaling by argument reduction   @ uses R0-02, flag 24,25 cleared 00 { 104-Byte Prgm } 01▸LBL "MEXP" 02 FUNC 11 03 L4STK 04 2 05 STO 00 06 RCL ST Y 07 FNRM 08 LN 09 RCL 00 10 LN 11 ÷ 12 IP 13 1 14 STO 01 15 + 16 STO 02 17 Y↑X 18 ÷ 19 ENTER 20 ENTER 21 EDIT 22 ↑ 23▸LBL 01 24 ↓ 25 RCL+ 01 26 → 27 FC? 77 28 GTO 01 29 EXITALL 30 1ᴇ-32   @ 1ᴇ-16 convergence test  31 STO 01 32▸LBL 02 33 R↓ 34 X<>Y 35 RCL÷ 00 36 ISG 00 37 NOP 38 RCL× ST Z 39 X<>Y 40 RCL+ ST Y 41 STO- ST L 42 LASTX 43 RNRM 44 RCL- 01 45 X>0? 46 GTO 02 47 R↓ 48 0=? 02 49 GTO 04 50▸LBL 03 51 ENTER 52 × 53 DSE 02 54 GTO 03 55▸LBL 04 56 FNRM 57 LASTX 58 END 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 RE: (DM42) Matrix exponential - Gil - 08-11-2023 11:46 PM 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 RE: (DM42) Matrix exponential - Gjermund Skailand - 08-12-2023 09:27 AM 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 RE: (DM42) Matrix exponential - Gil - 08-12-2023 10:01 AM 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 RE: (DM42) Matrix exponential - Gjermund Skailand - 08-12-2023 11:38 AM 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 RE: (DM42) Matrix exponential - Gil - 08-12-2023 08:26 PM 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 RE: (DM42) Matrix exponential - Gil - 08-12-2023 08:55 PM 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 RE: (DM42) Matrix exponential - Gjermund Skailand - 08-13-2023 10:23 AM 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 RE: (DM42) Matrix exponential - Gil - 08-13-2023 10:51 AM 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 RE: (DM42) Matrix exponential - Gjermund Skailand - 08-13-2023 09:41 PM 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:  @ Matrix exponential, square matrix @ calculation with argument reduction and series expansion @ convergence check uses RNRM difference between iterations, @ uses R00-R03 @ R00: counter, number of iterations @ R01: convergence limit,  set at line 37,  typical R01=1e-16 to R01=1e-32 @ R02: number of halvings and squarings used for argument reduction @ R03: number of initial iterations before convergence check (saves time),  @   suggest R03=10 when R01=1e-16  @   suggest R03=20 when R01=1e-32  @ Optional CF 24, CF 25 for overflow check 00 { 121-Byte Prgm } 01▸LBL "MEXP"  02 FUNC 11     @ 1 argument required, square matrix 03 L4STK 04 10 05 STO 03 06 CLX 07 2 08 STO 00 09 RCL ST Y    @ calculate and determine any argument reduction 10 FNRM 11 LN 12 RCL 00 13 LN 14 ÷ 15 IP 16 3           @ reduce the argument a little more 17 + 18 1 19 STO 01 20 + 21 X<0? 22 CLX         @ if negative, don't use 23 STO 02 24 Y↑X 25 ÷           @ reduce argument 26 ENTER 27 ENTER 28 EDIT        @ add the identity matrix directly 29 ↑ 30▸LBL 01 31 ↓ 32 RCL+ 01 33 → 34 FC? 77 35 GTO 01 36 EXITALL 37 1ᴇ-16        @ convergence crtiteria 38 STO 01 39 R↓ 40▸LBL 02       @ main iteration  41 X<>Y 42 RCL÷ 00 43 ISG 00 44 NOP 45 RCL× ST Z 46 X<>Y 47 RCL+ ST Y 48 DSE 03 49 GTO 02 50 STO- ST L 51 LASTX 52 RNRM 53 RCL- 01 54 R↓ 55 0I 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. RE: (DM42) Matrix exponential - Gil - 08-16-2023 12:45 PM 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 RE: (DM42) Matrix exponential - Gjermund Skailand - 08-16-2023 04:38 PM 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 RE: (DM42) Matrix exponential - Albert Chan - 08-16-2023 09:48 PM (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=0     if n<=1 then return r*n end     local R = f(r, floor(n/2))     R = R * (R+2)               -- expm1(log1p(r)*(n-n%2))     return n%2==0 and R or r*R + r + R end 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 RE: (DM42) Matrix exponential - Albert Chan - 08-21-2023 12:11 AM 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 RE: (DM42) Matrix exponential - Albert Chan - 08-21-2023 11:49 AM 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 "Nineteen dubious ways to calcate the exponential of a matrix, tventyfive years later." 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 RE: (DM42) Matrix exponential - Gjermund Skailand - 08-21-2023 07:03 PM 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