Charpoly (plus adjoint and det)

[5689d088-9d81-41e2-b555-2341cea5bc21]

[Description of the problem]

There are some problems in SAGE 4.0.2 when computing characteristic polynomials, determinants and adjoint matrices over general rings or non-exact fields. A more detailed description follows:

o Firstly, SAGE sometimes fails to compute the characteristic polynomial of a matrix over a ring which is not an integral domain. Here is an example:

```
sage: A = matrix(ZZ.quo(12), 3, [5,8,0,10,2,1,8,7,9])
sage: A
[ 5  8  0]
[10  2  1]
[ 8  7  9]
```
The call "A.charpoly()" results in an error, specifically

"ArithmeticError: The inverse of 10 modulo 12 is not defined."

o Secondly, when computing over non-exact fields like Qp sometimes the computation of the characteristic polynomial results in precision loss.

```
sage: R = Zp(5, prec = 10, type = 'capped-abs', print_mode = 'series')
sage: MS = MatrixSpace(R, 10, 10)
sage: A = MS.random_element()
sage: A.charpoly()
```
Often, this call results in a characteristic polynomial with coefficients of less precision than the base ring.  Sometimes (whenever the Hessenberg algorithm requires too many divisions in Zp with regards to the chosen precision), it also fails with the following error:

ValueError: element valuation cannot be negative.

o Thirdly, in some cases the current implementation of charpoly is ridiculously slow (because of the use of the field of fractions in the Hessenberg algorithm). Consider, for example, the following:

```
sage: R.<t> = QQ[]
sage: A = matrix(R, [[-2*t^2 - t - 2, -t, t^2 - 3/2*t - 1, 1/2*t^2 - 4*t - 1, -t^2 - 9*t + 2], \
                         [-t^2 + 1/2, 7/233*t - 2, -4*t + 1/2, 3*t^2 + 3/2*t + 1, -t^2 - t], \
                         [t - 1, t^2 - 5/7*t - 1, 1/2*t - 1/2, 8*t^2 + t - 2/3, t + 1/2], \
                         [-t^2 + 11, -1/2*t - 1/4, 8*t^2 + t, 1, -1/4*t^2 + 1/2*t + 1/7], \
                         [3/2*t^2 + 5*t, 13/50*t^2 - 11*t + 1, -2*t^2 + t, -1, -1/2*t + 3/2]])
```
(Sorry for the bad formatting.)  Computing the characteristic polynomial in SAGE via A.charpoly takes 91.12s on my laptop (Lenovo T500, Intel Centrino duo core), while the simple implementation of the generic division-free algorithm (attached to the ticket) takes only 0.01s.  During a quick (although not statistically sound) test with matrices in the above matrix space, the new code seemed to be faster by a factor between 1,000 and 10,000.

o For other reasons, namely using the expansion of co-factors, the computation of the determinant of a matrix over a general ring is also noticeably slow whenever the matrix has more than about 7 rows. For example, with

```
sage: A = random_matrix(PolynomialRing(QQ, 't'), 7, 7)
```
a call to ``A.det()`` takes about 0.5s.  For a random 8x8 matrix, the same call takes about 3s.  (And it gets worse very quickly in the expected way.)

o Also, just because it's possible, SAGE ought to be able to compute the characteristic polynomial, determinant and adjoint of square matrices over any ring (commutative with unity). To give an example where this currently fails:

```
sage: R.<a,b> = QQ[]
sage: S.<x,y> = R.quo(b^3)
sage: A = matrix(S, [[2,x^2],[x^3*y,x*y^2]])
sage: A
[    2   x^2]
[x^3*y x*y^2]
```
currently only ``A.det()`` works (and uses expansion by co-factors), but ``A.charpoly()`` throws a ``NotImplementedError`` (because SAGE does not create the univariate polynomial ring over ``S``, because the test whether S is a field passes on an exception as the ideal throws an exception when it is asked about maximality) and ``A.adjoint()`` throws

NotImplementedError: computation of adjoint not implemented in general yet

[Suggested fix]

Implement a generic division-free algorithm for the characteristic polynomial (and then use this to work out the adjoint, and read off the determinant).

[Problematic points]

• Adjusting current code that calls charpoly, det or adjoint:

To quickly mention the latter: in the file matrix2.pyx I simply implemented _adjoint(self), since inheriting classes overwrite this.

The calls to charpoly can choose an algorithm (as before, although it wasn't really used since there wasn't any choice). This problem thus translates to choosing the most sensible defaults (possibly depending on whether something is a number field, or an exact field, etc) in the implementation of charpoly, and to go through all the calls of charpoly in the current code, to check what the desired behaviour is. This question only arises as the Hessenberg algorithm runs in time O(n^3) whereas the division-free algorithm implemented here runs in time O(n^4). Thus the division-free algorithm wouldn't necessarily be the right choice between these two in all cases.

In the case of number fields (PARI) or exact fields (Hessenberg form) or Z/nZ (lift to Z), the choice is pretty clear. The same holds in cases where this wasn't implemented before. I think the only remaining cases for which there may or may not be an obvious (at least to me!) choice are non-exact fields like R or Qp (and their extensions). Then it's basically a choice between speed (and managing precision as the caller by ensuring one has good bounds on the input) or guaranteed precision.

Apart from one instance instance related to modular forms, where I could speak to David Loeffler at SAGE Days 16 and then added the paramater algorithm="hessenberg" to the call, I haven't changed the calls to charpoly anywhere.

The implementation of determinant reflects this somewhat, although it is no real problem since in the patch the logic of determinant hasn't changed. The only difference is that the "last resort" algorithm using expansion by co-factors has been replaced by the generic computation of the characteristic polynomial (from which one can then read off the determinant), i.e., this has no downside (apart from a tiny problem with symbolic expression rings, as we now need to specify a variable name for the characteristic polynomial, see below).

• A strange problem with symbolic rings:

The characteristic polynomial needs to have a variable name specified. This gives an intrinsic problem when computing the determinant (a method which does not require the user to add a variable name as input, and rightfully so shouldn't!) over symbolic rings, as one needs to choose a variable name for the characteristic polynomial different from the already existing symbols or else the computed result will be meaningless. At SAGE Days 16, Martin Albrecht and I briefly tried a quick way around this (by asking the symbolic expression ring for all symbols, and then forming a new one by concatenating them all, thereby always generating a new symbol), however, this idea did not work out. Thus at the moment, the fixed choice of the symbol "A0123456789" appears hardcoded.

PS: I think the inability of SAGE to compute these three quantities over all rings, even for small-ish matrices, could be potentially offputting for beginners, so I was tempted to choose "Priority=major", but then again, it doesn't seem to be a lot of work to fix this. This is my first ticket to I don't really know how to choose the priority and just went with "minor", hope that's OK.

Component: linear algebra

Keywords: charpoly, division-free

Author: Sebastian Pancratz

Reviewer: Rob Beezer, Yann Laigle-Chapuy

Merged: Sage 4.1.2.alpha2

Issue created by migration from https://trac.sagemath.org/ticket/6441

[5689d088-9d81-41e2-b555-2341cea5bc21]

comment:38

Dear Rob,

Yes, you're right. I hadn't even heard about the "-long" tests. Regarding the above problems, I take it I only need to look at the first few problems and not the known failures, right?

I've already checked that they definitely are "just" numerical noise, running the methods in question with the parameter "prec" in {10,20,30} gave answers 0.95, 1.0001, 0.99999994, respectively. However, I think the right way to fix this is not to change the doctests, but to pass the parameter algorithm="hessenberg" at the appropriate place to ensure that the behaviour isn't changed. I'll look into this probably tomorrow.

Sebastian

[1d7ec08f-60ae-4512-91a6-8324c06eab9f]

comment:39

Replying to @sagetrac-spancratz:

I hadn't even heard about the "-long" tests.

Sorry, my fault. First time it got me as well.

I take it I only need to look at the first few problems and not the known failures, right?

Yes, you can expect those other failures, but they are not yours to deal with.

I've already checked that they definitely are "just" numerical noise, running the methods in question with the parameter "prec" in {10,20,30} gave answers 0.95, 1.0001, 0.99999994, respectively. However, I think the right way to fix this is not to change the doctests, but to pass the parameter algorithm="hessenberg" at the appropriate place to ensure that the behaviour isn't changed. I'll look into this probably tomorrow.

Sounds good.

Rob

[5689d088-9d81-41e2-b555-2341cea5bc21]

comment:40

I am now adding two updated patches, rebased against 4.1.2.alpha0.

Sebastian

[5689d088-9d81-41e2-b555-2341cea5bc21]

20090906 - Including "is_field" etc for rings, to facilitate the division-free matrix operations

[5689d088-9d81-41e2-b555-2341cea5bc21]

Attachment: trac_6441_a_rings_412rebase.patch.gz

Attachment: trac_6441_b_df_charpoly_412rebase.patch.gz

20090906 - Including divison-free adjoint, charpoly and det

[1d7ec08f-60ae-4512-91a6-8324c06eab9f]

comment:41

Working with the 20090906 patches on 4.1.2.alpha1.

Builds fine, HTML docs constructed with no warnings, and passes all tests with options: -testall -long

Looks ready to go from here.

Rob

[7c09a680-e216-4024-bb8e-9bfd4aa7f313]

comment:42

Merged patches in this order:

1. trac_6441_a_rings_412rebase.patch
2. trac_6441_b_df_charpoly_412rebase.patch

[7c09a680-e216-4024-bb8e-9bfd4aa7f313]

Merged: Sage 4.1.2.alpha2