Opened 9 years ago

Last modified 8 years ago

#9944 closed defect

categories for polynomial rings — at Version 22

Reported by: robertwb Owned by: nthiery
Priority: major Milestone: sage-4.7.1
Component: categories Keywords:
Cc: sage-combinat Merged in:
Authors: Robert Bradshaw, Simon King Reviewers: Nicolas M. Thiéry, Mike Hansen
Report Upstream: N/A Work issues: Speedup
Branch: Commit:
Dependencies: Stopgaps:

Change History (27)

Changed 9 years ago by robertwb

comment:1 Changed 9 years ago by robertwb

  • Status changed from new to needs_review

comment:2 Changed 9 years ago by nthiery

Hi Robert!

I have been through the patch, and it sounds good! I won't have the time to actually test it before some time, so please anyone beat me to it!

Just one micro question: does PolynomialRing? actually check that the ring is commutative?

Cheers

Changed 9 years ago by mhansen

comment:3 Changed 9 years ago by mhansen

  • Authors set to Robert Bradshaw
  • Reviewers set to Nicolas M. Thiéry, Mike Hansen

I went ahead and moved the functionality to it's own category since we want the mathematical information at the category level.  Could someone look over these changes?

comment:4 follow-up: Changed 9 years ago by robertwb

The first patch only concerned univarite polynomial rings, the logic is not all correct for multivariate polynomial rings (though on an orthogonal note, that could use some fixing up as well). It seems odd to have a category of univariate polynomial rings over a fixed basering, which is why I put the logic into the concrete object. I suppose the category should a be declared as a graded R-algebra as well (do we have join categories yet?).

I don't know if PolynomialRing? asserts its basering is commutative, but IIRC it's been assumed for a long time.

comment:5 Changed 9 years ago by robertwb

Apply only 9944-poly-cat.patch

Changed 9 years ago by robertwb

comment:6 Changed 9 years ago by robertwb

Apply 9944-poly-cat.patch and 9944-poly-cat-doctests.patch .

Changed 9 years ago by mraum

comment:7 Changed 9 years ago by mraum

  • Description modified (diff)

I would give this a positive review for Robert's idea and I would open a new ticket for the multivariate rings. I'll just send a mail to Mike whether he is ok with that or no.

comment:8 in reply to: ↑ 4 ; follow-up: Changed 9 years ago by nthiery

Replying to robertwb:

The first patch only concerned univarite polynomial rings, the logic is not all correct for multivariate polynomial rings (though on an orthogonal note, that could use some fixing up as well). It seems odd to have a category of univariate polynomial rings over a fixed basering, which is why I put the logic into the concrete object. I suppose the category should a be declared as a graded R-algebra as well (do we have join categories yet?).

Sorry for the very late answer. In MuPAD, we had a category for univariate polynomial rings: there are several possible implementations of such, and it's natural to factor out the generic code, together with the category inheritance logic, in a category.

And yes, we have join categories. See Category.join.

I let you see whether to create the UnivariatePolynomialRing? category in this ticket or in a later ticket.

comment:9 in reply to: ↑ 8 ; follow-up: Changed 9 years ago by SimonKing

Replying to nthiery:

Sorry for the very late answer. In MuPAD, we had a category for univariate polynomial rings: there are several possible implementations of such, and it's natural to factor out the generic code, together with the category inheritance logic, in a category.

Aparently there is a doctest failure. I fixed it, but unfortunately it went into my patch submitted for #9138. Therefore, "needs work".

Question: Do we really want a category of polynomial rings? Or do we want that (1) polynomial rings use the category framework (that's the purpose of my patch for #9138) and (2) the category to which a given polynomial ring belongs is a bit narrower than simply "category of rings"? I hope it is the latter.

My suggestion is that I submit a small patch fixing the doctests. Please tell whether my patch for #9138 improves the multivariate case. Then, perhaps it would be possible to give Roberts patches (+ doctest fix) a positive review, so that we can focus on #9138.

comment:10 Changed 9 years ago by SimonKing

  • Status changed from needs_review to needs_work
  • Work issues set to Fix one test

comment:11 Changed 9 years ago by SimonKing

At #9138, Jason Bandlow reported a slow-down, that is at least partially caused by the patches here. Do you have any idea what could be the reason?

comment:12 in reply to: ↑ 9 Changed 9 years ago by SimonKing

  • Description modified (diff)

Replying to SimonKing:

Aparently there is a doctest failure. I fixed it, but unfortunately it went into my patch submitted for #9138. Therefore, "needs work".

Strange: Although the patch bot did see that error in one run, I can not reproduce it (but I had to change that test in my patch for #9138, because it turns QQ['x'].category() into the join of the category of euclidean domains and commutative algebras over QQ.

The other issue, namely the performance loss, was studied on sage-devel.

Florent Hivert found that a long mro does not matter for Python, but it does matter if the classes inherit from a cdef class. That is the case for most classes in Sage (inheriting from SageObject), so, we should address the problem of a long mro.

Eventually, that should be fixed in Cython (and I think Florent reported it upstream). But for now, it seems to me we should think of a work-around.

Would it be acceptable coding practice to explicitly state in a derived class (say, MPolynomialRing_generic), that frequently used methods such as base or base_ring are the same as Parent.base or Parent.base_ring? David Roe stated that it might be dangerous to do so, at least if cpdef methods are involved.

comment:13 Changed 9 years ago by SimonKing

Concerning performance loss:

It seems to me that part of the reason is the fact that with the patchse, __init_extra__ is not executed when it should. Sometimes, the parent methods of a category provide a useful __init_extra__, for example the category of algebras.

I am not sure yet why that happens, but I think it would happen if Parent.__init__ is called without providing the category: Namely, doing self._init_category_(...) alone will not trigger the execution of __init_extra__.

comment:14 Changed 9 years ago by SimonKing

It seems I was right!

Namely, the whole ring stuff is (unfortunately) inherited from ParentWithGens, which inherits from ParentWithBase, which inherits from parent_old.Parent.

And parent_old.Parent inherits from "the one and only Parent" -- but forgets to call Parent.__init__!! Instead, it just does self._init_category_(...).

I'll change it and see what happens.

comment:15 Changed 9 years ago by SimonKing

Very bad things happen. As soon as parent.Parent.__init__ is called in parent_old.Parent.__init__, an infinite recursion occurs.

Changed 9 years ago by SimonKing

Call Parent.init during initialisation of a ring

comment:16 Changed 9 years ago by SimonKing

  • Description modified (diff)

I attached a small patch to solve part of the problem of the missing parent initialisation: I call Parent.__init__ and ParentWithGens.__init__ explicitly, during initialisation of a ring. In that way, the __init_extra__ methods are correctly picked up.

However, that does not solve the performance problem.

Question one: How can one come to speed?

Question two: Is my patch really trivial enough to be called a referee patch?

For the patchbot:

Apply 9944-poly-cat.patch 9944-poly-cat-doctests.patch trac-9944-poly-cat-review.patch trac9944_second_referee.patch

comment:17 follow-up: Changed 9 years ago by SimonKing

  • Status changed from needs_work to needs_info
  • Work issues Fix one test deleted

FWIW: The doc tests pass.

Here is another idea what the slow down may come from. It was pointed out by Nicolas that the mro from the polynomial ring to Parent does not become longer by initialising the category properly: The inheritance from category parent classes comes after Python inheritance. However, when all the parent classes of all super categories must be searched, it takes considerably longer before an AttributeError can be raised.

A similar issue has been studied at #10467. It seems to be important that an attribute error is raised as quickly as possible. That becomes difficult, if 60 parent classes need to be searched, before one eventually finds that the requested attribute does not exist.

But probably that question is out of the scope of this ticket.

So, what shall one do? Give it a positive review and accept the deceleration, or wait until someone has a model for improved attribute access?

comment:18 in reply to: ↑ 17 Changed 9 years ago by SimonKing

Replying to SimonKing:

A similar issue has been studied at #10467. It seems to be important that an attribute error is raised as quickly as possible. That becomes difficult, if 60 parent classes need to be searched, before one eventually finds that the requested attribute does not exist.

No, that is not the problem here! I was inserting print statements into getattr_from_other_class in order to find out what attributes are actually requested from the category when doing arithmetic. It turned out that, during the first computation of (2*x-1)^2+5, some attributes are requested. But when one repeats that computation, getattr_from_other_class is not involved.

But what else could be the reason?

comment:19 Changed 9 years ago by SimonKing

Perhaps it is a conversion map that is slower than necessary?

If you look at sage.categories.algebras.Algebras.ParentMethods.__init_extra__, you see that it tries to register a certain set morphism as a coercion from the base ring into the algebra (that obviously works only if the algebra is unital).

But aparently a different coercion is used -- a slower coercion!

Namely, together with my patch from #9138:

sage: R.<x> = ZZ[]
sage: R.category() # the __init_extra__ was supposed to be used.
Join of Category of unique factorization domains and Category of commutative algebras over Integer Ring
sage: c = R.convert_map_from(R.base_ring())
sage: c
Polynomial base injection morphism:
  From: Integer Ring
  To:   Univariate Polynomial Ring in x over Integer Ring

That is not what __init_extra__ attempted to register!

Let us compare:

sage: from sage.categories.morphism import SetMorphism
sage: H = R.base().Hom(R)
sage: f = SetMorphism(H,R.from_base_ring)
sage: timeit('c(100)',number=10^5)
100000 loops, best of 3: 8.13 µs per loop
sage: timeit('f(100)',number=10^5)
100000 loops, best of 3: 1.75 µs per loop

So, things could be considerably improved. Obvious questions: Will from_base always yield a faster approach than the base injection morphism? And can we enforce to use the faster coercion?

Aparently it is not so easy:

sage: AC = Algebras(ZZ).parent_class
sage: R._unset_coercions_used()
sage: AC.__init_extra__(R)
sage: R.convert_map_from(R.base_ring())
Polynomial base injection morphism:
  From: Integer Ring
  To:   Univariate Polynomial Ring in x over Integer Ring
sage: R._unset_coercions_used()
sage: f.register_as_coercion()
sage: R.convert_map_from(R.base_ring())
Polynomial base injection morphism:
  From: Integer Ring
  To:   Univariate Polynomial Ring in x over Integer Ring

Can you explain how to force the use of a particular map for coercion of the base ring?

comment:20 Changed 9 years ago by SimonKing

And aparently the univariate polynomial rings are special in their choice of a conversion from the base ring. Again with #9138

sage: R.<m> = ZZ[]
sage: R.convert_map_from(R.base_ring())
Polynomial base injection morphism:
  From: Integer Ring
  To:   Univariate Polynomial Ring in m over Integer Ring
sage: R.<x,y> = QQ['t'][]
sage: R.convert_map_from(R.base_ring())
Generic morphism:
  From: Univariate Polynomial Ring in t over Rational Field
  To:   Multivariate Polynomial Ring in x, y over Univariate Polynomial Ring in t over Rational Field

comment:21 Changed 9 years ago by SimonKing

  • Status changed from needs_info to needs_work
  • Work issues set to Speedup

I suggest to speed things up by modifying "Polynomial base injection morphism". Internally, it uses rather slow ways of creating a polynomial of degree zero. It is likely to be faster to do what R.from_base does: Take the One of the ring (if it exists!) and use its _lmul_ method (if it has _lmul_).

I also understand why Algebras(...).parent_class.__init_extra__(R) has no effect on the choice of a conversion map from R.base() to R: It calls R.one() in order to create a better coercion map; but R.one() will, internally, construct a conversion from the base to R. At that point, the "better" coercion is not defined, and thus the usual conversion is created and cached.

In other words, R.from_base will only be used for conversion if R does not obey the rules of the new coercion framework (e.g., if it has a custom __call__).

Since the polynomial base injection morphism is a specialised method, it should be possible to internally construct the One of R without invoking coercion. My plan is to combine it with trac9944_second_referee.patch and submit a patch that then certainly needs a reviewer.

comment:22 Changed 9 years ago by SimonKing

  • Authors changed from Robert Bradshaw to Robert Bradshaw, Simon King
  • Description modified (diff)
  • Status changed from needs_work to needs_review

Good news! Things are now faster than without the patches!

I found that one can considerably improve the conversion of an element of the base ring into a polynomial ring. Some polynomial rings used a generic conversion map, some used a polynomial base injection map -- and both were slow.

My inspiration came from Algebras.ParentMethods.__init_extra__: If R is a polynomial ring, then multiplication of a scalar with R.one() often is a very fast method to convert the scalar into R.

Problems:

  • We should not assume that any ring has a unit (ok, polynomial rings over a unital ring have...).
  • Calling R.one() will usually trigger the creation of a generic conversion - hence, it would be difficult to register it as conversion.
  • Not all flavours of polynomial elements have a _rmul_ (polynomial_element_generic has not).
  • Sometimes, other conversion maps are registered when one wants to register the polynomial base injection map.

So, I implemented _rmul_ and _lmul_ for polynomial_element_generic, try various ways (old and new coercion model) of creating a One bypassing conversion maps, and in one init method of polynomial rings I decided to re-initialise the conversion maps.

Timings

I tried to test as many cases (multi- versus univariate, libSingular and others, different base rings,...). Without all the patches, we have the following:

sage: R.<x> = ZZ[]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 23.4 µs per loop
sage: R.<x> = QQ[]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 24.6 µs per loop
sage: R.<x> = GF(3)[]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 87.9 µs per loop
sage: R.<x> = QQ['t'][]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 113 µs per loop
sage: R.<x,y> = ZZ[]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 13 µs per loop
sage: R.<x,y> = QQ[]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 16.6 µs per loop
sage: R.<x,y> = GF(3)[]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 10.8 µs per loop
sage: R.<x,y> = QQ['t'][]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 238 µs per loop
sage: R.<x,y> = Qp(3)[]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 511 µs per loop
sage: R.<x> = Qp(3)[]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 1.06 ms per loop

With the patches, I get

sage: R.<x> = ZZ[]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 8.97 µs per loop
sage: R.<x> = QQ[]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 8.3 µs per loop
sage: R.<x> = GF(3)[]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 70.3 µs per loop
sage: R.<x> = QQ['t'][]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 82.6 µs per loop
sage: R.<x,y> = ZZ[]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 12.6 µs per loop
sage: R.<x,y> = QQ[]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 16.4 µs per loop
sage: R.<x,y> = GF(3)[]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 10.5 µs per loop
sage: R.<x,y> = QQ['t'][]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 187 µs per loop
sage: R.<x,y> = Qp(3)[]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 503 µs per loop
sage: R.<x> = Qp(3)[]
sage: timeit('(2*x-1)^2+5', number=10^4)
10000 loops, best of 3: 1.08 ms per loop

So, there is no significant slow down at all, but a considerable speed up in most cases.

I suppose it can now be reviewed. I understood that the Robert's patches essentially have a positive review, except for the slow-down. So, would it suffice if some of you test my patch?

Apply 9944-poly-cat.patch 9944-poly-cat-doctests.patch trac-9944-poly-cat-review.patch trac9944_polynomial_speedup.patch

Note: See TracTickets for help on using tickets.