Matrix Exponential (expm) might give deadly wrong results

[rdbyk]

Describe the bug

This seems to be related to a rather old M***** bug cf. https://blogs.mathworks.com/cleve/2012/07/23/a-balancing-act-for-the-matrix-exponential/#7401fe8a-5a7d-40df-92d7-1ae34f45adf2

To Reproduce

Expected behavior

cf. https://blogs.mathworks.com/cleve/2012/07/23/a-balancing-act-for-the-matrix-exponential/#7401fe8a-5a7d-40df-92d7-1ae34f45adf2

Screenshots cf. above

Desktop (please complete the following information):

• OS: Ubuntu 20.04

Additional context might be interesting, to have a look at https://nhigham.com/2016/03/08/the-improved-matlab-functions-expm-and-logm/ https://nhigham.com/2020/03/27/update-of-catalogue-of-software-for-matrix-functions/

[Nelson-numerical-software]

https://blogs.mathworks.com/cleve/2012/07/23/a-balancing-act-for-the-matrix-exponential/#7401fe8a-5a7d-40df-92d7-1ae34f45adf2

[Nelson-numerical-software]

contribution is welcome :)

[Nelson-numerical-software]

format long g
a = 2e10;
b = 4e8/6;
c = 200/3;
d = 3;
e = 1e-8;
A = [0, e, 0; -(a+b), -d, a; c, 0, -c]
E = expm(A)

expected result:

[ 0.446849, 1.5404410^-9, 0.462811, -5.7430710^6, -0.015283, -4.5265410^6 0.447723, 1.542710^-9, 0.463481]

[Nelson-numerical-software]

current algo call eigen library https://github.com/nelson-lang/nelson/blob/0c5a5a0c3d6ee6d02c8daca24b73c947fe6b68e4/modules/linear_algebra/src/cpp/ExpMatrix.cpp#L53

need to report bug or change code https://github.com/cran/expm/blob/master/src/expm.c (?)

[Nelson-numerical-software]

Temp workaround:

expm(complex(A))

[Nelson-numerical-software]

bug "dedicated" to Nick Higham (RIP)

[rdbyk]

This matrix exponential stuff is a pain ... the Eigen matrix exponential function seems to be from the "unsupported" section ...

But why don't use just the same approach for real matrices as you did for the complex case, i.e.

[V,D]=eig(A) E = Vdiag(exp(diag(D)))inv(V)

and maybe just implement it as a Nelson function for the time being and not a builtin to ease future modifications ?

At least this seems to be faster, than the workaround "expm(complex(A))" ...

Best wishes!

[Nelson-numerical-software]

need to be refactored for 1.6

[Nelson-numerical-software]

format long g
a = 2e10;
b = 4e8/6;
c = 200/3;
d = 3;
e = 1e-8;
A_IN = [0, e, 0; -(a+b), -d, a; c, 0, -c];
N = 3;
NDIAG = 9;
DELTA = 1.0;
BALANC = 'S';
[A_OUT, MDIG, IDIG, IWARN, INFO] = slicot_mb05od(BALANC, NDIAG, DELTA, A_IN);
A_OUT

result is good but slower

With A = [1 -1 -1;0 1 -1; 0 0 1];

 [V,D]=eig(A)
E = V*diag(exp(diag(D)))*inv(V)

result is wrong in this case

Need more investigation