# HP Forums

Full Version: (DM42) Matrix exponential
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
Pages: 1 2
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
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?

Regards,
Gil
Hi Gil

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

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
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
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
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).

Regards,
Gil
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
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
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 0<? ST T       @ convergence test 56 GTO 02         @ repeat 57 0=? 02         @ skip if no argument reduction 58 GTO 04 59▸LBL 03         @ argument back squaring 60 STO× ST X 61 DSE 02 62 GTO 03 63▸LBL 04 64 FNRM           @ check for /force "out of Range" error on DM42  v. Free 3.0.15 65 LASTX          @ delete line 64,65 for later revisions 66 END

Best regards
Gjermund
I am afraid of being lightning years beyond that level of programming and understanding.

Impressive anyway.

Regards,
Gil

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-15-2023 11:42 PM)Gil Wrote: [ -> ]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.

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.
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
Dear Gil,
there is a nice article
"Nineteen dubious ways to calcate the exponential of a matrix, tventyfive years later."

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-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-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     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
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
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
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
Pages: 1 2
Reference URL's
• HP Forums: https://www.hpmuseum.org/forum/index.php
• :