(DM42) Matrix exponential

[Albert Chan]

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

[Werner]

(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"
43▸LBL 01 @ add 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

[Gil]

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?

Thanks for your help.

Regards

[Albert Chan]

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

[Albert Chan]

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

[Albert Chan]

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

Add the 2 lines:

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

[Gjermund Skailand]

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

[Werner]

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
41▸LBL 01 @ add I
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

[Gjermund Skailand]

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

[Albert Chan]

(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

[John Keith]

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

[Gil]

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

[Albert Chan]

(08-28-2023 08:57 AM)Gil Wrote:  Is LOG2 (x)
log [in base 2] of x?

Yes, see Bits of integer