#6441: [with patch, needs review] Charpoly (plus adjoint and det)
----------------------------+-----------------------------------------------
Reporter: spancratz | Owner: somebody
Type: defect | Status: new
Priority: minor | Milestone: sage-4.1
Component: linear algebra | Keywords: charpoly, division-free
Reviewer: | Author: Sebastian Pancratz
Merged: |
----------------------------+-----------------------------------------------
Changes (by davidloeffler):
* component: algebra => linear algebra
Old description:
> [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.
New description:
[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.
--
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/6441#comment:5>
Sage <http://sagemath.org/>
Sage: Creating a Viable Open Source Alternative to Magma, Maple, Mathematica,
and MATLAB
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/sage-trac?hl=en
-~----------~----~----~----~------~----~------~--~---
