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

[Ilia Gradina]

Cause found, patch attached. Two independent defects, both in share/linearalgebra.

1. The residue division is not exact over the eigenvalues.
spectral_rep builds each projector as a residue of the resolvent at one
eigenvalue, so rational-residue divides the characteristic polynomial by
(z - zi) with zi already substituted. Once zi is a cubic irrationality
that division stops being exact: the projectors come out wrong by O(1) and
check-spectral-rep rejects them, which is the reported error.

For a simple root the projector is b(mat)/p'(lam) with b(z) = p(z)/(z-lam),
and nothing forces lam to be a radical while that is being built. Held as a
generic root -- tellrat puts the arithmetic in Q(...)[lam]/(p(lam)) -- every
intermediate stays a polynomial in lam of degree below n, and P^2 = P and
mat.P = lam*P, proved there, hold at all roots at once. Substituting the
eigenvalues afterwards is then exact.

The new branch is only reached when check-spectral-rep has already failed,
every eigenvalue is simple, and the characteristic polynomial carries no
parameter, so it cannot change a result that works today.

2. matrixexp and matrixfun ran under the default domain: real.
There the simplifier rewrites a cube root of an imaginary number onto another
branch:

c : (4*sqrt(29)*%i/3 - 244*%i/27)^(1/3)$
float(rectform(c));           =>  1.0644379645853810 - .6145535453890271*%i
float(rectform(ratsimp(c)));  =>  1.2291070907780544*%i

Both cube to the base, so nothing downstream notices, and the closing
fullratsimp quietly returned a matrix that is not exp(mat): measured
against the series to 60 terms the error was 2.7, not 3e-15. Cardano's formula
produces this shape whenever a 3x3 matrix has complex entries. spectral_rep
itself is deliberately left in the default domain -- switching it there makes
spectral_rep(matrix([x,1,0],[1,1,1],[0,1,1])) fail.

Tests. Three statements added to share/linearalgebra/rtest_matrixexp.mac.
They compare against a numerical value on purpose: xequal cannot referee this
case, because fullratsimp does not decide equality over these radicals.

Verification. x86-64 Fedora 44, SBCL 2.6.8, on top of 5.50.0: the
linearalgebra tests go from 440 to 443 with no failures, and the full suite
with share_tests: true reports "No unexpected errors found out of 20,941
tests".

The two patches are git format-patch output and apply with git am. The same
two commits are also on

git://git.code.sf.net/u/ilgrad/maxima  branch matrixexp-2596-3x3

I could not open a merge request for them: the source-branch list on the fork's
"Request Merge" form does not offer any branch pushed after the fork was
created, so the patches go here instead.

[Ilia Gradina]

file: 0001-spectral_rep-project-with-a-generic-eigenvalue-when-.patch

[Ilia Gradina]

file: 0002-matrixexp-matrixfun-work-in-DOMAIN-COMPLEX.patch

[Barton Willis]

Thanks for this. Here is an alternative that is, I think, much simpler. Let's replace the spendy test check-spectral-rep with a check that the putative nilpotent is really a nilpotent. I think the previous check was overkill.

Doing so yields:

(%i85) xxx : spectral_rep(mat)$

(%i86) zzz : first(xxx).second(xxx)+third(xxx) - mat$

(%i87) zzz : rectform(zzz)$

(%i88) zzz : expand(float(zzz));

(%o88) matrix([-(3.885780586188048e-16*%i)-5.551115123125783e-17,
               1.1102230246251565e-16-1.3877787807814457e-16*%i,
               2.7755575615628914e-17],
              [2.220446049250313e-16*%i+5.377642775528102e-17,
               -(3.0531133177191805e-16*%i),-(8.881784197001252e-16*%i)],
              [1.1102230246251565e-16-6.938893903907228e-16*%i,
               6.938893903907228e-17-1.6653345369377348e-16*%i,
               5.0306980803327406e-17*%i-2.220446049250313e-16])

[Barton Willis]

Thanks for this. Here is an alternative that is, I think, much simpler. Let's replace the spendy test check-spectral-rep with a check that the putative nilpotent is really a nilpotent. I think the previous check was overkill.

Doing so yields:

(%i85) xxx : spectral_rep(mat)$

(%i86) zzz : first(xxx).second(xxx)+third(xxx) - mat$

(%i87) zzz : rectform(zzz)$

(%i88) zzz : expand(float(zzz));

(%o88) matrix([-(3.885780586188048e-16*%i)-5.551115123125783e-17,
               1.1102230246251565e-16-1.3877787807814457e-16*%i,
               2.7755575615628914e-17],
              [2.220446049250313e-16*%i+5.377642775528102e-17,
               -(3.0531133177191805e-16*%i),-(8.881784197001252e-16*%i)],
              [1.1102230246251565e-16-6.938893903907228e-16*%i,
               6.938893903907228e-17-1.6653345369377348e-16*%i,
               5.0306980803327406e-17*%i-2.220446049250313e-16])

[Barton Willis]

Fixed by Commit [c5762f] . Closing ticket.