Improvements to discriminant computation

[gagern]

Currently (i.e. on sage 6.1.1), the discriminant computation fails with a segmenttation fault when the coefficients are too complicated. For example:

PR.<b,t1,t2,x1,y1,x2,y2> = QQ[]
PRmu.<mu> = PR[]
E1 = diagonal_matrix(PR, [1, b^2, -b^2])
RotScale = matrix(PR, [[1, -t1, 0], [t1, 1, 0], [0, 0, 1]])
E2 = diagonal_matrix(PR, [1, 1, 1+t1^2])*RotScale.transpose()*E1*RotScale
Transl = matrix(PR, [[1, 0, -x1], [0, 1, -y1], [0, 0, 1]])
E3 = Transl.transpose()*E2*Transl
E4 = E3.subs(t1=t2, x1=x2, y1=y2)
det(mu*E3 + E4).discriminant()

Apparently the Python interpreter reaches some limit on the depth of the call stack, since it seems to call symtable_visit_expr from Python-2.7.6/Python/symtable.c:1191 over and over and over again. But this ticket here is not about fixing Python. (If you want a ticket for that, feel free to create one and provide a reference here.)

Instead, one should observe that the above is only a discriminant of a third degree polynomial. As such, it has an easy and well-known formula for the discriminant. Plugging the complicated coefficients into this easy formula leads to the solution rather quickly (considering the immense size of said solution). So perhaps we should employ this approach in the discriminant method.

I've discussed this question on sage-support, and will probably continue the discussion. Right now, I'm unsure about when to choose a new method over the existing one. Should it be based on degree, using hardcoded versions of the well-known discriminants for lower degrees? Or should it rather be based on base ring, thus employing a computation over some simple generators instead of the full coefficients if the base ring is a polynomial ring or some other fairly complex object. Here is a code snippet to illustrate the latter idea:

def generic_discriminant(p):
    """
    Compute discriminant of univariate (sparse or dense) polynomial.
   
    The discriminant is computed in two steps. First all non-zero
    coefficients of the polynomial are replaced with variables, and
    the discriminant is computed using these variables. Then the
    original coefficients are substituted into the result of this
    computation. This should be faster in many cases, particularly
    with big coefficients like complicated polynomials. In those
    cases, the matrix computation on simple generators is a lot
    faster, and the final substitution is cheaper than carrying all
    that information along at all times.
    """
    pr = p.parent()
    if not pr.ngens() == 1:
        raise ValueError("Must be univariate")
    ex = p.exponents()
    pr1 = PolynomialRing(ZZ, "x", len(ex))
    pr2.<y> = pr1[]
    py = sum(v*y^e for v, e in zip(pr1.gens(), ex))
    d = py.discriminant()
    return d(p.coefficients())

Personally I think I prefer a case distinction based on base ring, or a combination. If we do the case distinction based on base ring, then I'd also value some input as to what rings should make use of this new approach, and how to check for them. For example, I just noticed that #15061 discusses issues with discriminants of polynomials over power series. So I guess those might benefit from this approach as well.

CC: @pjbruin

Component: algebra

Keywords: discriminant

Author: Martin von Gagern

Branch: 1892c84

Reviewer: Peter Bruin

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

[gagern]

Branch: u/gagern/ticket/16014

[gagern]

New commits:

28980e9	Added test case to expose the segfault of discriminant.
0758de1	Drive-by fix to TeX notation in docstring.
d5e0307	Use simple coefficients if base is multivariate polynomial.
66a8c19	Compute discriminant of power series polynomial using substitution.

[gagern]

Commit: 66a8c19

[gagern]

comment:3

What I implemented so far is a case distinction based on the base ring, with the new implementation being used for multivariate polynomials and power series.

I implemented the new code in the existing method for now, in part because I wasn't sure whether I should add new methods (e.g. discriminant_by_resultant and discriminant_by_substitution) or rather an algorithm argument for the discriminant method. In the latter case, the code would stay in that method but be reorganized because the substitution version could call the resultant version instead of executing that computation directly.

I'm not sure about this whole ZeroDivisionError stuff from #11782 resp. a224fe3. When doing the substitution approach, there should be no zero divisors, so replacing self by poly in that exception handler is only for the sake of consistency, but not really required. I guess if I were to write code for the substitution approach exclusively, I might consider using // poly[n] instead of * an to avoid the detour to the fraction field.

[pjbruin]

comment:5

This looks good and passes doctests. However, I'm wondering if you couldn't be convinced to adopt John Cremona's proposal on sage-dev to have a separate function generic_univariate_discriminant() (or similar), which caches its output, and hardcodes it for small degrees. This would be useful for two reasons:

• computing the universal discriminant takes a substantial amount of time (several seconds for degree 7, several minutes for degree 8, probably extremely long for higher degrees);
• it would make the discriminant() method cleaner and therefore more attractive to add an algorithm flag and/or a more sophisticated algorithm selection heuristic.
  I would think sage.rings.polynomial.polynomial_element would be a logical place for such a function, although sage.rings.invariant_theory would also be defensible.

Here is a simple function to compute the universal discriminant:

def universal_discriminant(n):
    r"""Return the discriminant of the 'universal' univariate polynomial
    `a_n x^n + \cdots + a_1 x + a_0` in `\ZZ[a_0, \ldots, a_n][x]`."""
    R = PolynomialRing(ZZ, n + 1, 'a')
    S = PolynomialRing(R, 'x')
    return S(list(R.gens())).discriminant()

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

1892c84

Introduce a caching function universal_discriminant.

[sagetrac-git]

Changed commit from 66a8c19 to 1892c84

[gagern]

comment:7

Replying to @pjbruin:

I'm wondering if you couldn't be convinced to adopt John Cremona's proposal on sage-dev to have a separate function generic_univariate_discriminant() (or similar), which caches its output, and hardcodes it for small degrees.

Sure I can be convinced, particularly if two people request something. I've implemented the caching, but not the hardcoding so far. I wonder how small “small” actually should be. Up to 7, things compute fairly quickly so there seems to be little benefit from hardcoding, particularly since we do caching. For higher degrees, computing them for hardcoding takes considerable time, and writing down the hard code will take up quite a lot of space in the source file. Should this really be implemented in Cython code, or should we make some effort to obtain or compute discriminants up to say 10 or so and store them in some kind of persistent data representation instead of source code?

I also wonder how public this function should be. I prefer polynomial_element over invariant_theory simply because it keeps all my changes in one place. I've written a full docstring for the code, and inside that docstring I had to import the function. Is this a reasonable approach? Should I add the function to sage.rings.polynomial.all? Should I prepend the function name with an underscore to mark it non-public?

I also have a question regarding my late import. I noticed the cdef void late_import near the beginning of the module. I could add is_PowerSeriesRing to that. But on the other hand, it seems as if that function only gets called by Polynomial.roots. Which suggests that I might not be able to rely on is_PowertSeriesRing if I were to follow that approach. Do you agree with keeping the import where I have it so far? Should the function late_import be renamed to something like late_import_for_root to make it clearer that this does no general late importing but specific to that method? Or should I add my imports to that function and make sure to call late_import in discriminant as well?

Should I add myself as an author of that module, at the beginning of the file?

[pjbruin]

comment:8

Replying to @gagern:

Replying to @pjbruin:

I'm wondering if you couldn't be convinced to adopt John Cremona's proposal on sage-dev to have a separate function generic_univariate_discriminant() (or similar), which caches its output, and hardcodes it for small degrees.

Sure I can be convinced, particularly if two people request something. I've implemented the caching, but not the hardcoding so far. I wonder how small “small” actually should be. Up to 7, things compute fairly quickly so there seems to be little benefit from hardcoding, particularly since we do caching.

Thanks, this looks very good. You're probably right about not hardcoding the polynomials.

For higher degrees, computing them for hardcoding takes considerable time, and writing down the hard code will take up quite a lot of space in the source file. Should this really be implemented in Cython code, or should we make some effort to obtain or compute discriminants up to say 10 or so and store them in some kind of persistent data representation instead of source code?

Maybe, but I doubt it is very important. If at all, it should be done in a later ticket.

I also wonder how public this function should be. I prefer polynomial_element over invariant_theory simply because it keeps all my changes in one place. I've written a full docstring for the code, and inside that docstring I had to import the function. Is this a reasonable approach? Should I add the function to sage.rings.polynomial.all? Should I prepend the function name with an underscore to mark it non-public?

I actually agree completely with the way you did all these things. I wouldn't add anything to sage.rings.polynomial.all; this would add it to the global namespace, which is already very crowded.

I also have a question regarding my late import. I noticed the cdef void late_import near the beginning of the module. I could add is_PowerSeriesRing to that. But on the other hand, it seems as if that function only gets called by Polynomial.roots. Which suggests that I might not be able to rely on is_PowertSeriesRing if I were to follow that approach. Do you agree with keeping the import where I have it so far? Should the function late_import be renamed to something like late_import_for_root to make it clearer that this does no general late importing but specific to that method? Or should I add my imports to that function and make sure to call late_import in discriminant as well?

Here, too, I think the way you did it is absolutely fine. In fact one can argue that modules should generally be imported locally in the function that needs it. On the other hand, global imports (including the late_import approach) are useful if the imported module is used in many places, or if the function is time-critical enough that a local import would add too much overhead.

Should I add myself as an author of that module, at the beginning of the file?

If you want; as you will have seen, authors for some other additions/changes are listed, but it is certainly not done in a very systematic way througout Sage. In any case, your changes and authorship will be remembered by the revision control system, and using git annotate gives a line-by-line overview of who wrote/changed what.

Could you also fill in your name (i.e. real name, not Trac username) in the "Author" field of this ticket?

I will test this and then set it to positive review.

[pjbruin]

Reviewer: Peter Bruin

[pjbruin]

comment:9

All tests pass and I am happy with the patch as it is. Please add your name and set the ticket to positive review.

[gagern]

Author: Martin von Gagern

[vbraun]

Changed branch from u/gagern/ticket/16014 to 1892c84

[embray]

Changed commit from 1892c84 to none

[embray]

comment:12

FWIW the segfault in Python was fixed in newer Python versions: http://bugs.python.org/issue5765

For some reason (it seems just developer resources) a version of the patch was never backported to 2.7. I might take a stab at that (and open a separate ticket if so). Note the patch doesn't magically allow unlimited stack depth and one can still blow the stack. But with the default settings it will raise a RuntimeError instead, and give one the opportunity to increase the limit with sys.setrecursionlimit.