#2596 error in matrix exponentiation

[kcrisman]

This is a regression in 5.30.0 from 5.29.1:

(%i1) m: matrix([%i*%pi]);
(%o1)                             [ %i %pi ]
(%i6) matrixexp(m),keepfloat:true;

Unable to find the spectral representation
 -- an error. To debug this try: debugmode(true);
(%i7) matrixexp(m),keepfloat:false;

Unable to find the spectral representation
 -- an error. To debug this try: debugmode(true);

Point of info - Barton W. has a workaround:

This bug is due to circa January 2013 changes to matrix inversion. A workaround is to set ratmx to true:

Maxima branch_5_30_base_98_g29f9239_dirty http://maxima.sourceforge.net
using Lisp Clozure Common Lisp Version 1.9-r15764  (WindowsX8632)

(%i1) matrixexp(matrix([%i*%pi]));

Unable to find the spectral representation

(%i2) matrixexp(matrix([%i*%pi])), ratmx=true;
(%o2)                                 - 1

[kcrisman]

Potential fix reported downstream on Sage ticket #13973:

~~~~~~~
--- a/share/linearalgebra/matrixexp.lisp 2013-10-07
04:37:12.000000000 +0100
+++ b/share/linearalgebra/matrixexp.lisp 2014-05-16
02:16:09.112011893 +0100
@@ -138,8 +138,8 @@
(print (ratvars = ,$ratvars gcd = '$gcd algebraic = ,$algebraic)) (print(ratfac = ,$ratfac))
(merror "Unable to find the spectrum")))
-

• (setq res ($fullratsimp (ncpower (sub (mult z ($ident n)) mat) -1) z))
  +
• (setq res ($fullratsimp ($invert (sub (mult z ($ident n)) mat) '$crering) z))
  (setq m (length sp))
  (dotimes (i m)
  (setq zi (nth i sp))
  ~~~~~~

[kcrisman]

Here is another example from 5.33. The ratmx workaround does not help here.

(%i2) m:matrix([1,2,3],[3,2,0],[1,2,1]);
                                  [ 1  2  3 ]
                                  [         ]
(%o2)                             [ 3  2  0 ]
                                  [         ]
                                  [ 1  2  1 ]
(%i3) m:%i*m;
                             [  %i   2 %i  3 %i ]
                             [                  ]
(%o3)                        [ 3 %i  2 %i   0   ]
                             [                  ]
                             [  %i   2 %i   %i  ]
(%i4) matrixexp(m),keepfloat:true;

Unable to find the spectral representation
 -- an error. To debug this try: debugmode(true);
(%i5) matrixexp(m),ratmx:true;

Unable to find the spectral representation
 -- an error. To debug this try: debugmode(true);

[kcrisman]

That said, I am not claiming this is a regression, just an example which it doesn't seem ever worked (perhaps for good theoretical reasons, though of course it has such an exp). See http://ask.sagemath.org/question/23784/imaginary-matrix-exponential/ where it came up.

[Michael Orlitzky]

This is about to become important as we begin to allow SageMath users to use the Maxima already installed on their system. Sage has been carrying this custom patch for seven years and has tests for the problem it fixes. Is the patch in kcrisman's comment reasonable? We need a solution that everyone agrees on; otherwise the tests will fail on some machines but not others.

[Michael Orlitzky]

The latest version of this patch uses invert_by_lu instead of invert, but is still a nice one-liner. This would really help out with SageMath integration if it made its way to a release:

diff --git a/share/linearalgebra/matrixexp.lisp b/share/linearalgebra/matrixexp\
.lisp                                                                           
index 218bf35..f2fd468 100644                                                   
--- a/share/linearalgebra/matrixexp.lisp                                        
+++ b/share/linearalgebra/matrixexp.lisp                                        
@@ -138,8 +138,8 @@                                                             
           (print `(ratvars = ,$ratvars gcd = '$gcd algebraic = ,$algebraic))   
           (print `(ratfac = ,$ratfac))                                         
           (merror "Unable to find the spectrum")))                             

-                                                                               
-    (setq res ($fullratsimp (ncpower (sub (mult z ($ident n)) mat) -1) z))     
+                                                                               
+    (setq res ($fullratsimp ($invert_by_lu (sub (mult z ($ident n)) mat) '$cre\
ring) z))                                                                       
     (setq m (length sp))                                                       
     (dotimes (i m)                                                             
       (setq zi (nth i sp))

[Doron Behar]

Also the rtest_matrixexp test fails for me here on NixOS.

[Robert Dodier]

It works as shown with a current version of Maxima (post 5.46), but it needs scalarmatrixp: false, otherwise it bumps into an error about an intermediate result not being a matrix. Presumably to fix that error, somewhere along the way we need to ensure scalarmatrixp is enabled for the duration of the computation.

[Dima Pasechnik]

with scalarmatrixp: false it works on 5.46.0, too.

[Barton Willis]

I'll locally set scalarmatrixp: false for several top-level functions. This will fix the first bug in the report.

The case matrixexp(%i*matrix([1,2,3],[3,2,0],[1,2,1])) fails because Maxima is doesn't simplify some algebraic expressions that almost surely vanish to zero. These simplifications are only needed to verify the answer is correct.

I could turn off this check and Maxima could return a ridiculously lengthy answer for this case. I'm disinclined to turn off the checks.

I'd need to be convinced that we win by using invert_by_lu.

Maybe there is a compromise here, but I don't like the option of allowing the user to optionally turn off the checks. Likely the better way is to return the matrix exponential in terms of the symbolic eigenvalues. That's a big project, I think.

What do you all think about optionally turning off the checks?

[Barton Willis]

The matrixexp(matrix([%i*%pi])) bug is fixed by commit [3683af]. The bug with the 3 x 3 matrix mentioned later on remains unfixed.

I will leave this ticket open.