Menu
▾
▴
#2596 error in matrix exponentiation
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
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 m (length sp))
(dotimes (i m)
(setq zi (nth i sp))
~~~~~~
Here is another example from 5.33. The
ratmxworkaround does not help here.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.
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.
The latest version of this patch uses
invert_by_luinstead ofinvert, but is still a nice one-liner. This would really help out with SageMath integration if it made its way to a release:Also the
rtest_matrixexptest fails for me here on NixOS.Diff:
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 ensurescalarmatrixpis enabled for the duration of the computation.with
scalarmatrixp: falseit works on 5.46.0, too.I'll locally set
scalarmatrixp: falsefor 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?
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.
Related
Commit: [3683af]
Cause found, patch attached. Two independent defects, both in share/linearalgebra.
1. The residue division is not exact over the eigenvalues.
spectral_repbuilds each projector as a residue of the resolvent at oneeigenvalue, so
rational-residuedivides the characteristic polynomial by(z - zi)withzialready substituted. Onceziis a cubic irrationalitythat division stops being exact: the projectors come out wrong by O(1) and
check-spectral-reprejects them, which is the reported error.For a simple root the projector is
b(mat)/p'(lam)withb(z) = p(z)/(z-lam),and nothing forces
lamto be a radical while that is being built. Held as ageneric root --
tellratputs the arithmetic inQ(...)[lam]/(p(lam))-- everyintermediate stays a polynomial in
lamof degree below n, andP^2 = Pandmat.P = lam*P, proved there, hold at all roots at once. Substituting theeigenvalues afterwards is then exact.
The new branch is only reached when
check-spectral-rephas already failed,every eigenvalue is simple, and the characteristic polynomial carries no
parameter, so it cannot change a result that works today.
2.
matrixexpandmatrixfunran under the defaultdomain: real.There the simplifier rewrites a cube root of an imaginary number onto another
branch:
Both cube to the base, so nothing downstream notices, and the closing
fullratsimpquietly returned a matrix that is notexp(mat): measuredagainst 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_repitself 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:
xequalcannot referee thiscase, because
fullratsimpdoes 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: truereports "No unexpected errors found out of 20,941tests".
The two patches are
git format-patchoutput and apply withgit am. The sametwo commits are also on
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.
file: 0001-spectral_rep-project-with-a-generic-eigenvalue-when-.patch
file: 0002-matrixexp-matrixfun-work-in-DOMAIN-COMPLEX.patch
Thanks for this. Here is an alternative that is, I think, much simpler. Let's replace the spendy test
check-spectral-repwith a check that the putative nilpotent is really a nilpotent. I think the previous check was overkill.Doing so yields:
Thanks for this. Here is an alternative that is, I think, much simpler. Let's replace the spendy test
check-spectral-repwith a check that the putative nilpotent is really a nilpotent. I think the previous check was overkill.Doing so yields:
Fixed by Commit [c5762f] . Closing ticket.
Related
Commit: [c5762f]
Thanks for c5762fdc4. spectral_rep now succeeds on the 3x3 matrix, but the
matrix exponential that comes out of it is wrong, and nothing signals an error.
Checked on master 92979cf7b against a truncated power series:
The new test in rtest_matrixexp.mac checks first(xxx).second(xxx)+third(xxx)
against mat. That is true, but it doesn't look at matrixexp itself.
matrixfun has the same problem:
The cause: spectral_rep binds domain to complex, while matrixexp and matrixfun
simplify their result under the default domain: real. There the simplifier
moves a cube root of an imaginary eigenvalue onto another branch. Binding
domain: complex in matrixexp and matrixfun as well brings the errors down to
8e-14 and 2e-12.
I bind it only when the matrix contains %i. For a real matrix it isn't needed
(the real matrix below comes out right under domain: real), and with it the
parametric tests 74 and 75 in rtest_matrixexp no longer simplify:
matrixfun(lambda([x], x+1/x), mat) gives a 370k-character expression instead
of matrix([x,2,-1],...).
A second case still fails on master, the same matrix without the %i:
The eigenvalues are cubic irrationalities again, and nilpotent-p can't show the
remainder is zero over them. My earlier patch handles this with a fallback for
simple eigenvalues: it builds the projector once at a generic root, in
Q(...)[lam]/(p(lam)) via tellrat, and checks P^2 = P and mat.P = lam*P there
before substituting the roots. It only runs where the current code gives up.
With it the error is 6e-14.
The attached patch applies to current master and contains both changes. It
also adds tests that compare matrixexp and matrixfun(sin) numerically with
the series. I ran these cases against it: the two above, a rotation
matrix([0,1],[-1,0]) with and without the t argument, a Jordan block
matrix([2,1],[0,2]), and matrixfun(lambda([z], z^3), m) against m.m.m. All
agree to 2e-12 or better, and domain is unchanged afterwards. The full
run_testsuite(share_tests=true) on master with the patch (SBCL 2.6.9) reports no
unexpected errors out of 21,687 tests.
Could you reopen this?
Try this experiment: Globally set
domainto complex. With the current linearalgebra, I getShow me what you get.
Last edit: Barton Willis 2 days ago
Thanks for reopening. I get the same number: your session on master 92979cf7b gives 1.3402904179514637e-26. Without
domain : complexthe same session gives 3.1271276229227505, and that is what a user gets by default.Your session also shows something I missed. My patch binds domain: complex inside matrixexp, and then
bfloat(rectform(matrixexp(m)))is right at the default domain: real. But as soon as the result is simplified again under domain: real, it is wrong again (patched, domain: real,xxxas in your session):So the expression itself has different values under the two domains, and binding domain inside matrixexp only moves the problem to the next simplification. I would drop that part of my patch. What does survive is rectform of the result under domain: complex: ratsimp and ev of it under domain: real stay below 1e-23. But the expression grows from 63k to 601k characters, so I doubt it is worth it. How would you like to handle the default domain?
The second case doesn't depend on domain. On master the real matrix fails under both:
With the tellrat fallback from my patch,
ratsimp(matrixexp(r))agrees with the series to 8e-27 under domain: real and 3e-27 under domain: complex. I can send a patch with only that part and its tests.Isn't setting
domaintocomplexand possiblym1pbranchtotat the top-level ofmatrixexpandmatrixfunsufficient? Plus it avoids the extra complexity of tellrat for a bug that is really about the simplification of radicals, not about the matrixexp code.Actually, I've never used tellrat all that much and I'm not sure that I trust it.
get:
Setting
domaintocomplexandm1pbranchtotat the top-level ofmatrixexpandmatrixfunI getWhat do you think of this?