(DM42) Matrix exponential

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