What does MPolynomial_libsingular.reduce() do?

[Luca De Feo]

Hi everyone,

Here's a Sage session:

sage: A.<x,y> = QQ[]
sage: (x+y).reduce([(x-y), (x+y)])
0
sage: (x-y).reduce([(x-y), (x+y)])
-2*y

The docstring says reduce computes "the normal form of self w.r.t. I,
i.e. [...] the remainder of this polynomial with respect to the
polynomials in I".

Does anyone have any idea how this normal form is defined? It doesn't
seem to depend on the order of the polynomials in I.

From the source code, I can only tell it calls Singular's kNF, but I
can't find any doc for it. Maybe this function should be underscored?

Cheers,
Luca

[Daniel Krenn]

On 2017-10-16 18:41, Luca De Feo wrote:
> Here's a Sage session:
>
> sage: A.<x,y> = QQ[]
> sage: (x+y).reduce([(x-y), (x+y)])
> 0
> sage: (x-y).reduce([(x-y), (x+y)])
> -2*y
>
> The docstring says reduce computes "the normal form of self w.r.t. I,
> i.e. [...] the remainder of this polynomial with respect to the
> polynomials in I".
>
> Does anyone have any idea how this normal form is defined? It doesn't
> seem to depend on the order of the polynomials in I.

It computes the polynomial "modulo" the given ideal (i.e. compute a
Groebner basis of the ideal and reduce the given polynomial by this basis).

My guess: If only a list of polynomials is given, then it is assumed
that these form a Groebner basis, which seems not to be the case.

>>From the source code, I can only tell it calls Singular's kNF, but I
> can't find any doc for it. Maybe this function should be underscored?

Once we know what it does with lists, the documentation should be made
precise.

I am against underscoring, as for ideals as parameter, this is a
standard operation.

[Martin R. Albrecht]

Hi there,

this is already documented:

“ Return the normal form of self w.r.t. "I", i.e. return the

remainder of this polynomial with respect to the polynomials in

"I". If the polynomial set/list "I" is not a (strong) Groebner
basis the result is not canonical.
”

Cheers,
Martin

--

_pgp: https://keybase.io/martinralbrecht
_www: https://martinralbrecht.wordpress.com
_jab: martinr...@jabber.ccc.de
_otr: 47F43D1A 5D68C36F 468BAEBA 640E8856 D7951CCF

[Luca De Feo]

On Mon, Oct 16, 2017 at 7:35 PM, 'Martin R. Albrecht' via sage-devel
<sage-...@googlegroups.com> wrote:
> Hi there,
>
> this is already documented:
>
> “ Return the normal form of self w.r.t. "I", i.e. return the
> remainder of this polynomial with respect to the polynomials in
> "I". If the polynomial set/list "I" is not a (strong) Groebner
> basis the result is not canonical.
> ”

Can you tell from this documentation what the function will compute
prior to running it? I can't.

I agree with Daniel: this function does something useful and sensible
when I is an ideal, so it shouldn't be underscored.

But I have no idea of what it does when I is a list, except give an
undefined result congruent to self modulo the ideal generated by I. So
I again agree with Daniel: if we can figure out what this function
does, we should document it better. And I would go as far as adding
that if we can't figure it out, we should forbid list input.

Luca

[Travis Scrimshaw]

Can you tell from this documentation what the function will compute
prior to running it? I can't.

It takes I as the generators of the ideal and uses that as the reduction set.

I agree with Daniel: this function does something useful and sensible
when I is an ideal, so it shouldn't be underscored.

But I have no idea of what it does when I is a list, except give an
undefined result congruent to self modulo the ideal generated by I. So
I again agree with Daniel: if we can figure out what this function
does, we should document it better. And I would go as far as adding
that if we can't figure it out, we should forbid list input.

What it does is probably do the reduction using the list in reverse order for this case. As previously mentioned, because it is not a Gröbner basis, there is no guarantee of a canonical result. So IMO it does what the documentation says it does.

Best,

Travis

 

[Luca De Feo]

> It takes I as the generators of the ideal and uses that as the reduction
> set.

That's not a definition. I'm in front of a class asking what this
function does, and I'm unable to give a mathematical definition of
what Sage means by "reduction" modulo something that's not a Groebner
basis.

> What it does is probably do the reduction using the list in reverse order
> for this case.

"Probably" is not a mathematical definition. Besides, I think it's
more complicated than that.

Am I the only one who's regularly embarassed explaining Sage's quirks
to an audience of beginners (or not beginners)?

Luca

[Daniel Krenn]

+1 for doing something.

What about the following fix: When the input is a list/tuple, we check
if it is a Groebner basis or not. If it is, do the computation, if not,
print a warning or raise an error.

Testing if something is a Groebner basis could be done by converting the
list to an object of
<class
'sage.rings.polynomial.multi_polynomial_sequence.PolynomialSequence_generic'>
and use its method .is_groebner()

Daniel

[Luca De Feo]

> +1 for doing something.
>
> What about the following fix: When the input is a list/tuple, we check
> if it is a Groebner basis or not. If it is, do the computation, if not,
> print a warning or raise an error.

Sounds reasonable.

Other options would be:

- Just refuse list input.
- Always compute a Groebner basis.

Since it seems that David Roe wrote the original code, his opinion
would be very valuable.

Luca

[Martin R. Albrecht]

Hi there,

AFAIK if you do that you prevent high-level implementation of
Gröbner basis algorithms in Sage which call reduce,
i.e. polynomial division with remainders, on S-polynomials wrt to
the current basis.

Cheers,
Martin

Daniel Krenn <kr...@aon.at> writes:

[Martin R. Albrecht]

Hi there,

I’m fairly certain that I wrote that code or was at least involved
since I wrote most of the first version of the libsingular stuff,
David’s commit output by `git blame` is a merge commit.

The algorithm in question is described in “Ideals, Varieties, and
Algorithms” by David Cox, John Little and Donal O’Shea in Section
"§3 A Division Algorithm in k[x_1 , … , x_n]".

I didn’t check, though, if under the hood Singular does some
additional tricks which might modify the output when I is not a
GB.

Cheers,
Martin

PS: Checking if something is a GB calls reduce to check. See
`basis_is_groebner`.

[Travis Scrimshaw]

On Tuesday, October 17, 2017 at 4:50:20 AM UTC-5, Luca De Feo wrote:

> It takes I as the generators of the ideal and uses that as the reduction
> set.

That's not a definition. I'm in front of a class asking what this
function does, and I'm unable to give a mathematical definition of
what Sage means by "reduction" modulo something that's not a Groebner
basis.

My understanding is it does the (naïve?) reduction algorithm:

reducing = True
while reducing:
    reducing = False
    for expr in I:
        if expr.leading_monomial().divides(cur.leading_monomial()):
            cur -= cur.leading_term() / expr.leading_term() * expr
            reducing = True

So this is well-defined, but it doesn't guarantee a unique result: it depends on the ordering of I.

I agree with Martin, that this could have unintended consequences with doing things in quotient rings because IIRC this is the key function that is called for every operation. So having something that checks that the input is a GB could slow things down. IMO, it also makes sense to expose this to users.

Best,

Travis

[Simon King]

Hi Luca,

On 2017-10-16, Luca De Feo <de...@lix.polytechnique.fr> wrote:
> On Mon, Oct 16, 2017 at 7:35 PM, 'Martin R. Albrecht' via sage-devel
><sage-...@googlegroups.com> wrote:
>> Hi there,
>>
>> this is already documented:
>>
>> “ Return the normal form of self w.r.t. "I", i.e. return the
>> remainder of this polynomial with respect to the polynomials in
>> "I". If the polynomial set/list "I" is not a (strong) Groebner
>> basis the result is not canonical.
>> ”
>
> Can you tell from this documentation what the function will compute
> prior to running it? I can't.

When reading `normal form` and `Groebner basis` in the same sentence,
the meaning should be clear to anybody who took a course in commutative
algebra. So, the question is: Whom should documentation be addressed to?

I do *not* think that documentation should always be addressed to
non-experts.

Here, we talk about methods of polynomials respectively polynomial
ideals. Thus, a user having a closer look into it is, to some extent,
assumed to know about commutative algebra.

In addition, if the documentation mentions a technical notion such
as `Groebner basis`, and it doesn't ring a bell to the user, then
of course the user should either conclude right away that the
method might not be what (s)he was looking for, or should look up
the technical notion before drawing that conclusion.

Just for comparison: In the method galois_conjugates of number field
elements, it says right in the beginning
Return all Gal(Qbar/Q)-conjugates of this number field element
in the field K.
No explanation what that actually means. Again "of course", if the
user doesn't know what Gal(Qbar/Q) means and just knows from
calculus that there is a conjugation in CC (which is why (s)he
stumbled upon the method), (s)he should certainly not expect to
understand every detail of the documentation.

And likewise for polynomial reduction.

Best regards,
Simon

[Simon King]

On 2017-10-17, Daniel Krenn <kr...@aon.at> wrote:
> What about the following fix: When the input is a list/tuple, we check
> if it is a Groebner basis or not.

Too expensive.

Cheers,
Simon

[Simon King]

On 2017-10-17, Luca De Feo <de...@lix.polytechnique.fr> wrote:
>> It takes I as the generators of the ideal and uses that as the reduction
>> set.
>
> That's not a definition. I'm in front of a class asking what this
> function does, and I'm unable to give a mathematical definition of
> what Sage means by "reduction" modulo something that's not a Groebner
> basis.

Why do they expect you to be able to? Unless it is a class on
commutative algebra.

>> What it does is probably do the reduction using the list in reverse order
>> for this case.
>
> "Probably" is not a mathematical definition. Besides, I think it's
> more complicated than that.

I agree on that part of your statement. When implementing polynomial
reduction, the result (if you reduce by a list of polynomials that
do not form a Gröbner basis) will depend on the order in which
the reductions are done - hence, the outcome is not uniquely
determined by the underlying mathematics but by implementation
details. Not uncommon, if the result isn't unique.

> Am I the only one who's regularly embarassed explaining Sage's quirks
> to an audience of beginners (or not beginners)?

Are your students the only ones that are embarassed if the teacher
tells them that they are looking at a method which can only
be understood with a background in some specific subfield of maths
(such as number fields, commutative algebra, group theory, topology)?
My students certainly wouldn't be embarassed.

Cheers,
Simon

[john_perry_usm]

Perhaps it would be clearer if the wording were changed to, "return a remainder," because there can be more than one, and strike the words between and including "the normal form" and "i.e.", because in many places "normal form" implies uniqueness. For instance, speaking of Cox, Little, and O'Shea, at least their second edition mentions that normal forms are unique (p. 80 and p. 85, Exercises 1 and 4). But not all texts take this route (see below).

Otherwise, I agree that the function should work with lists of polynomials that are not Gröbner bases; there are valuable reasons for it. For instance, when teaching commutative algebra one can illustrate easily that the remainder is not unique by using a basis is not Gröbner.

Singular itself does this (see pp. 51-52 of "A Singular Introduction to Commutative Algebra") and the text even speaks of a non-unique "NF [normal form] w.r.t. a non-standard basis", illustrating how the results are different when the basis is not standard (aka Gröbner) and how there is a unique canonical result when the basis is standard. So in fact on p. 44, "After weakening the requirements, there is no more uniqueness statement for the normal form" and on p. 46 the definition of normal form is applied to a list of polynomials, not to an ideal, basis of an ideal, or even standard (Gröbner) basis.

[Luca De Feo]

> The algorithm in question is described in “Ideals, Varieties, and
> Algorithms” by David Cox, John Little and Donal O’Shea in Section "§3 A
> Division Algorithm in k[x_1 , … , x_n]".

It is not:

sage: A.<x,y> = PolynomialRing(QQ, order="lex")
sage: (x*y^2 + 1).reduce([x*y + 1, y+1])
x + 1
sage: (x*y^2 + 1).reduce([y+1, x*y + 1])
x + 1

whereas Example 1 in §3 gives remainder 2;

sage: A.<x,y> = PolynomialRing(QQ, order="lex")
sage: (x^2*y + x*y^2 + y^2).reduce([x*y - 1, y^2 - 1])
2*x + 1
sage: (x^2*y + x*y^2 + y^2).reduce([y^2 - 1, x*y - 1])
2*x + 1

whereas Example 2 gives remainder x + y + 1.

> My understanding is it does the (naïve?) reduction algorithm:
>
> reducing = True
> while reducing:
> reducing = False
> for expr in I:
> if expr.leading_monomial().divides(cur.leading_monomial()):
> cur -= cur.leading_term() / expr.leading_term() * expr
> reducing = True

It does not:

sage: def travis_red(f, I):
....: cur = f
....: reducing = True
....: while reducing:
....: reducing = False
....: for expr in I:
....: if expr.lm().divides(cur.lm()):
....: cur -= cur.lt() // expr.lt() * expr
....: reducing = True
....: return cur
....:
sage: travis_red(x^2*y + x*y^2 + y^2, [y^2 - 1, x*y - 1])
2*x + y^2

> Singular itself does this (see pp. 51-52 of "A Singular Introduction
> to Commutative Algebra") and the text even speaks of a non-unique "NF
> [normal form] w.r.t. a non-standard basis", illustrating how the
> results are different when the basis is not standard (aka Gröbner)
> and how there is a unique canonical result when the basis is
> standard.

Singular seems to have changed algorithm since. The example 1.6.13
you're referring to computes slightly different:

sage: A.<x,y,z> = PolynomialRing(QQ, order="degrevlex")
sage: f = x^2*y*z + x*y^2*z + y^2*z + z^3 + x*y
sage: f1 = x*y + y^2 - 1
sage: f2 = x*y
sage: f.reduce([f1, f2]) # same as in the example
y^2*z + z^3
sage: f.reduce([f2, f1]) # missing a xz monomial
y^2*z + z^3 - y^2 + 1


Luca

[Luca De Feo]

Hi Simon,

I hate to sound snarky, but...

> When reading `normal form` and `Groebner basis` in the same sentence,
> the meaning should be clear to anybody who took a course in commutative
> algebra. So, the question is: Whom should documentation be addressed to?
>
> I do *not* think that documentation should always be addressed to
> non-experts.

I am teaching a commutative algebra class to grad students. They and I
perfectly know what a Gröbner basis is. I have read many times Cox,
Little and O'Shea, and I hope my students have too.

Yet, none of us seems to be able to second guess what kind of normal
form is actually implemented by .reduce() (Singular's kNF, actually).
And from the answers to this thread, it seems to me that neither Sage
devs can.

So, maybe, we should at least modify the documentation so that
*experts* can understand it.

FTR, my question was raised by an exercise in my course book that is
several (about 6?) years old, whose goal was exactly to teach the
students that f.reduce([g, h]) can give different results from
f.reduce([h, g]) when [g, h] is not a Gröbner basis. At some point in
the past, Sage behaviour changed, and the two calls stopped giving
different results for the specific f,g,h in the exercise (at least I
think it used to give different results, my memory may be faulty).

Luca

[Simon King]

Hi Luca,

On 2017-10-18, Luca De Feo <de...@lix.polytechnique.fr> wrote:
> I hate to sound snarky, but...

No offence taken...

> Yet, none of us seems to be able to second guess what kind of normal
> form is actually implemented by .reduce() (Singular's kNF, actually).
> And from the answers to this thread, it seems to me that neither Sage
> devs can.

I agree that Sage devs don't (need to) know what *exact* algorithm is
used. But at least it is clear that the .reduce() method does polynomial
reductions and continues till no further reduction is possible; all this
for a fixed monomial ordering (determined by the polynomial ring), in
*some* order of execution that is considered an implementation detail
and is thus not necessarily backwards compatible or compatible with an
order of execution proposed in the literature.

Best regards,
Simon

[Luca De Feo]

> I agree that Sage devs don't (need to) know what *exact* algorithm is
> used. But at least it is clear that the .reduce() method does polynomial
> reductions and continues till no further reduction is possible; all this
> for a fixed monomial ordering (determined by the polynomial ring), in
> *some* order of execution that is considered an implementation detail
> and is thus not necessarily backwards compatible or compatible with an
> order of execution proposed in the literature.

Something along the lines of this paragraph would be a valuable
addition to the docstring of .reduce().

However, what's the use of such a function with
implementation-dependent outputs?

If I understand Martin's argument correctly, he is saying that
.reduce() *could* be used for a schoolbook "implementation of Gröbner

basis algorithms in Sage which call reduce, i.e. polynomial division
with remainders, on S-polynomials wrt to the current basis."

However for .reduce() to be useful for such applications, we would
need to have at least *some* documented guarantee on its output. You
say:

> But at least it is clear that the .reduce() method does polynomial
> reductions and continues till no further reduction is possible;

I agree with your analysis, but "no reduction is possible" is
ambiguous: different textbooks mean different things by it (some mean
Travis' algorithm, some others agree with Cox, Little and O'Shea).

I would formulate it more precisely as "no monomial in the output is
divisible by the leading monomial of any polynomial in I". This seems
to be true given the examples I have, it would be good to have a
confirmation (hidden somewhere in Singular's docs?)

Luca

[john_perry_usm]

I would formulate it more precisely as "no monomial in the output is
divisible by the leading monomial of any polynomial in I". This seems
to be true given the examples I have, it would be good to have a
confirmation (hidden somewhere in Singular's docs?)

Curiously, the definition of "normal form" in the Singular docs differs from that in the Singular text; see http://www.singular.uni-kl.de/Manual/latest/sing_845.htm#IDX415 where normal form is defined with a standard basis G. In the text G can be any list of polynomials.

It is possible in Singular to reduce by something that is not a Gröbner basis but Singular usually complains when that happens. So, for instance, here is a computation I just carried out in Singular:

> ring R = 43,(x,y),dp;
> ideal I = x2+y2-4,xy-1;
> reduce(x3,I);
// ** I is no standard basis
4x-y

The check that produces the warning is NOT a thorough check that G is a Gröbner basis; it merely inspects a flag to see if a Gröbner basis has been computed. You can circumvent the warning by setting the flag, so that even when it isn't a Gröbner basis there is no warning. I don't remember offhand how to set that flag.

Looking at the code, Sage first tries to see if it has cached a reduction already (at least, I think that's what I._groebner_strategy is about). If not, it invokes the Singular function kNF (line 4479 of sage/libs/singular/decl.pxd in Sage 8.0). Singular's source code is not well-known for documentation, and kNF is no exception: in my copy of the source there is more or less no documentation of what kNF is supposed to do. (My copy is old, but I know a developer pretty well & it's something we've talked about.) I worked with it some years back and as I recall kNF and the functions it calls together behave according to the description from the Singular text that I quoted earlier: "normal form" in that text is not the same as "normal form" in CLO. Figuring out exactly how they do this could be a challenge, because Singular can take one of many paths depending on the underlying data.

This is a long way of saying that I agree with the suggestion.

[Simon King]

Hi Luca,

On 2017-10-18, Luca De Feo <de...@lix.polytechnique.fr> wrote:

> However, what's the use of such a function with
> implementation-dependent outputs?
>
> If I understand Martin's argument correctly, he is saying that
> .reduce() *could* be used for a schoolbook "implementation of Gröbner
> basis algorithms in Sage which call reduce, i.e. polynomial division
> with remainders, on S-polynomials wrt to the current basis."
>
> However for .reduce() to be useful for such applications, we would
> need to have at least *some* documented guarantee on its output.

Why? The result of the schoolbook implementation of Buchberger's
algorithm will result in *a* Gröbner basis, regardless of the concrete
choice of execution order in the implementation of polynomial reduction;
and if a reduced Gröbner basis is computed, it is *unique*.

>> But at least it is clear that the .reduce() method does polynomial
>> reductions and continues till no further reduction is possible;
>
> I agree with your analysis, but "no reduction is possible" is
> ambiguous: different textbooks mean different things by it (some mean
> Travis' algorithm, some others agree with Cox, Little and O'Shea).
>
> I would formulate it more precisely as "no monomial in the output is
> divisible by the leading monomial of any polynomial in I". This seems
> to be true given the examples I have, it would be good to have a
> confirmation (hidden somewhere in Singular's docs?)

Well, reducing the leading term is enough to get a Gröbner basis,
but you need tail reductions to compute the reduced Gröbner basis.
And Singular provides both options; I am sure Sage can use that, too.

Best regards,
Simon

[Bill Hart]

According to Hans Schoenemann:

"Usually (i.e. in a ring with a well ordering and no additional flags)
reduce/kNF (p,I) computes p' (with p a polynomial and I a list of polynomials)
with p-p' is in the ideal generated by I and no monomial of p' is
divisible by any L(f) for f in I.
If I is a standard basis then p' is unique.
Otherwise, the result depends on the internally used algorithm.
In Singular's implementation p' does not depend on the order of the
polynomials in I because it starts with sorting I
(wrt. to the monomial ordering or the total degree)."

[Bill Hart]

Hans added that one should not rely on the behaviour that Singular reorders the list. This may change in a future version of Singular. He mentioned it only to explain the behaviour that is currently observed.

[Luca De Feo]

Thanks everyone for the discussion.

I opened https://trac.sagemath.org/ticket/24071 to improve the
docstring of reduce().

It's ready for review, let's continue the discussion there.

> In Singular's implementation p' does not depend on the order of the
> polynomials in I because it starts with sorting I
> (wrt. to the monomial ordering or the total degree)."

I had come to the same conclusion, but that's not a total ordering, so:

sage: A.<x,y> = QQ[]

sage: (x+y).reduce([x+y, x-y])
2*y
sage: (x+y).reduce([x-y, x+y])
0


Luca