(DM42) Matrix exponential
08-21-2023, 07:45 PM
Post: #21
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (DM42) Matrix exponential
(08-21-2023 07:03 PM)Gjermund Skailand Wrote:  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 was doing expm1(A) = exp(A) - 1 (= identity matrix). Adjusted, both agreed.

Quote:But my implementation on HP50g using double binary floating point fails.

This is why example is nasty, if we just simply do e^A = 1 + A + A^2/2! + A^3/3! + ...
The terms blows-up passes true answer. (Matrix A already bigger than final sum!)

A^16/16! ≈ [[6977251.9061, -3488625.95305], [9303002.54147, -4651501.27073]]
A^17/17! ≈ [[-6977251.9061, 3488625.95305], [-9303002.54147, 4651501.27073]]

No. the 2 terms are not exactly negate of the others.
But, catastrophic cancellation made the sum meaningless.
08-23-2023, 07:16 AM (This post was last modified: 08-23-2023 07:38 AM by Werner.)
Post: #22
 Werner Senior Member Posts: 891 Joined: Dec 2013
RE: (DM42) Matrix exponential
(08-21-2023 11:49 AM)Albert Chan Wrote:  1118905.69945      1374815.06297      1630724.42651
2533881.04197      3113414.03144      3692947.02095
3948856.38452      4852012.99997      5755169.61548
(..)
HP71B, SCALE BY 2^16 or 2^17

-1.73575875814      0.551819099660
-1.47151759909      0.103638240710

I can reproduce these exactly on the 42S using the program below.
(well making sure I have P=13 in the first example..)
If you want MEXP, perform an extra M+I

00 { 99-Byte Prgm }
01▸LBL "MXP-1"
02 ENTER
03 FNRM
04 LN
05 2
06 LN
07 ÷
08 10
09 +
10 IP
11 X<0?
12 CLX
13 X<>Y
14 2
15 RCL ST Z
16 +/-
17 Y^X
18 ×
19 RCL ST X
20 4
21 ÷
22 XEQ 01
23 3
24 ÷
25 RCL× ST Y
26 XEQ 01
27 2
28 ÷
29 RCL× ST Y
30 XEQ 01
31 X<>Y
32 ×
33▸LBL 04
34 ENTER
35 +
36 LASTX
37 STOx ST X
38 +
39 DSE ST Y
40 GTO 04
41 RTN
42▸LBL "M+I"
44 EDIT
45 ^
46▸LBL 02
47 ↓
48 ENTER @ no, you can't remove this
49 1
50 +
51 →
52 FC? 77
53 GTO 02
54 EXITALL
55 END

Cheers, Werner

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
08-23-2023, 09:09 AM
Post: #23
 Gil Senior Member Posts: 637 Joined: Oct 2019
RE: (DM42) Matrix exponential
Two questions

Line 40 :
P=MAX(0,IROUND(LOG2(MAXAB(X)))+10)

What is the meaning (perhaps with an example) of the function/command MAXAB(X)?

And LOG2 (x) is log [in base 2] of x?

Regards
08-23-2023, 02:51 PM
Post: #24
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (DM42) Matrix exponential
(08-23-2023 09:09 AM)Gil Wrote:  What is the meaning (perhaps with an example) of the function/command MAXAB(X)?

MAXAB(X) is from JFG "Math Pac 2" - Version 2B Feb. 2021
MAXAB = Array Element Maximum Absolute Value. (Perhaps it should be named MAXABS)
It has a side effect of saving where element is located.

>mat disp x
-5.46969688594
22010.2875425
-22141314.9958

>maxab(x), arow, acol
22141314.9958            3      1
>minab(x), arow, acol
5.46969688594            1      1

Quote:P=MAX(0,IROUND(LOG2(MAXAB(X)))+10)

This is to make sure X/2^P elements size at most ε = 2^(-9.5) ≈ 0.00138

expm1(ε) = ε + ε^2/2! + ε^3/3! + ε^4/4! + ...

If I drop..., expm1 relative error ≈ (ε^5/5!)/ε = ε^4/120 ≤ 3.0E-14, less than machine epsilon.

Note: matrix multiply is a dot product of row and col, it really depends on cell values.
Above is *very* rough estimate of relative errors. (and, only if ε is truly small)
08-23-2023, 07:55 PM
Post: #25
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (DM42) Matrix exponential
Trivia, for A = 2x2 matrix, e^A only need its eigenvalues. (i.e. P, P-1 not needed)

A = P D P-1

D diagonal has eigenvalues of A = c ± g --> D-c has diagonal of ± g

Let X = A-c --> X^2 = g^2      // × identity matrix

sinh(X) = (1 + X^2/3! + X^4/5! + ...) * X = sinh(g)/g * X
cosh(X) = (1 + X^2/2! + X^4/4! + ...)      = cosh(g)      // × identity matrix

e^X = sinh(X) + cosh(X) = sinh(g)/g * X + cosh(g)

XCas> eA := ((sinh(g)/g) * (A-c) + cosh(g)) * e^c

XCas> A := [[a11,a12],[a21,a22]];
XCas> x1, x2 := eigenvalues(A) :;
XCas> c, g := [x1+x2, x1-x2]/2 :;
XCas> simplify(hyp2exp(exp(A) - eA))       // confirmed symbolically

$$\left(\begin{array}{cc}0&0\\0&0\end{array}\right)$$
08-23-2023, 09:39 PM (This post was last modified: 08-24-2023 02:26 AM by Albert Chan.)
Post: #26
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (DM42) Matrix exponential
(08-23-2023 07:55 PM)Albert Chan Wrote:  sinh(X) = (1 + X^2/3! + X^4/5! + ...) * X = sinh(g)/g * X

We can use same idea to get expm1(k*Jn), where Jn = matrix of all ones, dimensions n×n

Identity used: (Jn/n)^(integer_powers) = (Jn/n)

sinh(k*Jn) = (k*Jn) + (k*Jn)^3/3! + (k*Jn)^5/5! + ... = (1 + (k*n)^2/3! + (k*n)^4/5! + ...) * (k*Jn)

sinh(k*Jn) = sinh(k*n) * (Jn/n) = 2*sinh(k*n/2)*cosh(k*n/2) * (Jn/n)
cosh(k*Jn) − 1 = 2*sinh(k/2*Jn)^2 = 2*sinh(k*n/2)*sinh(k*n/2) * (Jn/n)

exp(x*Jn) - 1 = 2*sinh(n/2*x)*exp(n/2*x) * (Jn/n)

exmp1(k*Jn) = expm1(k*n) * (Jn/n)

Any function f with f(0) = 0, will have same pattern. (*)
Taylor series (Jn/n)^(integer_powers) = (Jn/n), which can be factored out.

f(k*Jn) = f(k*n * (Jn/n)) = f(k*n) * (Jn/n)

For examples, these will work: sin, versin, tan, tanh, log1p, ... (and its inverse function)

(*) Exception, if f cannot handle matrix argument.
Example, because Jn is not invertible, this will not work: f(x) := sin(x^2)/x
08-24-2023, 07:59 AM
Post: #27
 Gjermund Skailand Member Posts: 78 Joined: Dec 2013
RE: (DM42) Matrix exponential
Here is another version of Alberts code for the DM42.

There are two tuning parameters, to make use of the 32 digits accuracy of the DM32. Current values should give same result as Werners version.

at line 24, increase from 10 for more reduction of input matrix
at line 34, increase from 4 to e.g. 9, error ~ 0.00138^9/10!= 5.0E-33
If you want to calculate MEXP-1 delete last line (line 60)

Code:
 00 { 103-Byte Prgm } 01▸LBL "MEXP"      @ Extended Albert version 02 FUNC 11            @ has two tuning parameters 03 L4STK               @ uses R00 04 ENTER               @ delete line 60  for MEXP-1 05 DIM? 06 NEWMAT 07 EDIT 08 ← 09▸LBL 00 10 ↓ 11 SIGN 12 → 13 FC? 77 14 GTO 00 15 EXITALL 16 X<>Y 17 ENTER 18 FNRM 19 LN 20 2 21 LN 22 ÷ 23 IP 24 10               @ parameter for halving 25 + 26 X<0? 27 CLX 28 STO 00 29 2 30 X<>Y 31 Y↑X 32 ÷ 33 RCL ST X 34 4                @ parameter for reverse series expansion 35 ÷ 36 LASTX 37 DSE ST X 38 DSE ST X 39 R↓ 40 RCL+ ST Z 41▸LBL 01 42 ISG ST T 43 NOP 44 RCL÷ ST T 45 RCL× ST Y 46 RCL+ ST Z 47 DSE ST T 48 DSE ST T 49 GTO 01 50 X<>Y 51 × 52▸LBL 02 53 ENTER 54 × 55 LASTX 56 STO+ ST X 57 + 58 DSE 00 59 GTO 02 60 + 61 END
Best regards
Gjermund

Attached File(s)
08-24-2023, 01:14 PM
Post: #28
 Werner Senior Member Posts: 891 Joined: Dec 2013
RE: (DM42) Matrix exponential
Hi Gjermund,
adapted mine as well so that the number of expm1 terms could be easily changed.
I switched to 6 i.o. 10 squarings though - for the 42S version (where it matters most), that meant two fewer matrix multiplications (6 terms, 6 squarings give the same accuracy as 4 terms, 10 squarings), and for the DM42 it doesn't matter (6 terms and 14 squarings or 10 terms and 10 squarings).
And, still all in the stack ;-)
(I'd like to make it independent of the calculator's precision, but that proves to be a bit of a challenge, only using the stack ;-)

00 { 98-Byte Prgm }
01▸LBL "E↑M-I"
02 ENTER
03 FNRM
04 LN
05 2
06 LN
07 ÷
08 6
09 +
10 IP
11 X<0?
12 CLX
13 X<>Y
14 2
15 RCL ST Z
16 +/- @ avoids overflow..
17 Y^X
18 ×
19 14 @ 6 for 42S
20 RCL ST Y
21▸LBL 05 @ X := (X/Y+I)*Z, keeping YZT
22 RCL÷ ST Y
23 XEQ 01
24 RCL× ST Z
25 DSE ST Y @ DSE n, skip when n=1
26 DSE ST Y
27 ISG ST Y
28 X=Y? @ aff
29 GTO 05
30 R^
31 X<>Y
32▸LBL 04
33 ENTER
34 +
35 LASTX
36 STO× ST X
37 +
38 DSE ST Y
39 GTO 04
40 RTN
42▸LBL "M+I"
43 EDIT
44 ^
45▸LBL 02
46 ↓
47 SIGN
48 X<0?
49 +/-
50 RCL+ ST L
51 →
52 FC? 77
53 GTO 02
54 EXITALL
55 END

Cheers, Werner

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
08-24-2023, 07:12 PM
Post: #29
 Gjermund Skailand Member Posts: 78 Joined: Dec 2013
RE: (DM42) Matrix exponential
Hi Werner
on DM42 I found that matrix functions are very effective
compared to user operations, so keeping a copy and adding of the Identity matrix really pays off.
Albert's algoritm is much more robust (no surprise there) and about 10% faster than my original implementation.

Best regards
Gjermund
08-26-2023, 10:20 PM
Post: #30
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (DM42) Matrix exponential
(08-23-2023 09:39 PM)Albert Chan Wrote:  exmp1(k*Jn) = expm1(k*n) * (Jn/n)

My old code epsilon has enough safety margin for small matrices.
But, in general, perhaps we need to cover the edge case.

If n×n elements sized about the same, say k, FNORM(X) ≈ k*n
Code:
< 40 P=MAX(0,IROUND(LOG2(MAXAB(X)))+10) @ MAT X=(1/2^P)*X > 40 P=MAX(0,IROUND(LOG2(FNORM(X)))+9) @ MAT X=(1/2^P)*X

Example, expm1(J20) = expm1(20)/20 * J20 = 24258259.7205 * J20

Old Patch: P=10 ---> 24258259.1465 * J20
New Patch: P=13 --> 24258259.7183 * J20
08-27-2023, 04:46 PM
Post: #31
 John Keith Senior Member Posts: 1,039 Joined: Dec 2013
RE: (DM42) Matrix exponential
(08-26-2023 10:20 PM)Albert Chan Wrote:  If n×n elements sized about the same, say k, FNORM(X) ≈ k*n
Code:
< 40 P=MAX(0,IROUND(LOG2(MAXAB(X)))+10) @ MAT X=(1/2^P)*X > 40 P=MAX(0,IROUND(LOG2(FNORM(X)))+9) @ MAT X=(1/2^P)*X

I like this, especially since the HP-28/48 don't have MAXAB but do have ABS (same as FNORM ).
08-28-2023, 08:57 AM
Post: #32
 Gil Senior Member Posts: 637 Joined: Oct 2019
RE: (DM42) Matrix exponential
Is LOG2 (x)
log [in base 2] of x?
08-28-2023, 01:10 PM
Post: #33
 Albert Chan Senior Member Posts: 2,703 Joined: Jul 2018
RE: (DM42) Matrix exponential
(08-28-2023 08:57 AM)Gil Wrote:  Is LOG2 (x)
log [in base 2] of x?

Yes, see Bits of integer
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)