(DM42) Matrix exponential

[Gjermund Skailand]

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

[Gil]

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

[Gjermund Skailand]

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

[Gil]

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

[Gjermund Skailand]

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

[Gil]

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

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

Thanks in advance for your cooperation.

Regards,
Gil

[Gjermund Skailand]

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

[Gil]

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

[Gjermund Skailand]

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

[Gil]

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

Impressive anyway.

Regards,
Gil

[Gil]

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




.

[John Keith]

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

[Gil]

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

[Gjermund Skailand]

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

[Albert Chan]

(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)\)

[Albert Chan]

(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

[Albert Chan]

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

[Albert Chan]

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

[Gjermund Skailand]

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