Opened 12 years ago

Closed 12 years ago

#6482 closed defect (fixed)

[with patch, positive review] multivariate polynomial substitution has a design flaw

Reported by: was Owned by: malb
Priority: critical Milestone: sage-4.2
Component: commutative algebra Keywords:
Cc: Merged in: sage-4.2.alpha0
Authors: Martin Albrecht Reviewers: Mike Hansen
Report Upstream: Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Status badges

Description

On Wed, Jul 8, 2009 at 1:28 AM, Kwankyu<...> wrote:
>
> Hi,
>
> I was surprised to see
>
> sage: R.<x,y>=QQ[]
> sage: g=x+y
> sage: g.subs({x:x+1,y:x*y})
> x*y + x + y + 1
>
> So the order of substitution matters...unfortunately.
>
> sage: g.subs({x:x+1}).subs({y:x*y})
> x*y + x + 1
> sage: g.subs({y:x*y}).subs({x:x+1})
> x*y + x + y + 1
>
> So the order seems to be from right to left. This seems to me
> unnatural. Anyway this is undocumented. 

Actually, i guess it is documented.  However, I consider it a serious design flaw.
Many thanks for pointing this out!!

I consider this a serious design flaw for the following reasons:

 (1) it is too hard to understand the above behavior, since it depends on the hash values symbolic variables, which might possibly be system-dependent.

 (2) it is totally inconsistent with the behavior for symbolic expressions, where things are done right.

 (3) it is totally inconsistent with the behavior of *homomorphisms*, where things are also done right.

Here is a session to illustrate the above points:

# BAD
sage: R.<x,y>=QQ[]
sage: g=x+y
sage: g.subs({x:x+1,y:x*y})
x*y + x + y + 1

# BAD
sage: R.<x,y>=QQ[]
sage: g=x+y
sage: g.subs(x=x+1,y=x*y)
x*y + x + y + 1

# GOOD
sage: R.<x,y>=QQ[]
sage: phi = R.hom([x+1,x*y])
sage: g=x+y
sage: phi(g)
x*y + x + 1

# GOOD
sage: var('x,y')
sage: g = x+y
sage: g.subs({x:x+1,y:x*y})
x*y + x + 1

# GOOD
sage: var('x,y')
sage: g = x+y
sage: g.subs(x=x+1,y=x*y)
x*y + x + 1
        

> What should be done to this?

1. I suggest that for now you use Hom, as illustrated above, as a workaround.  

2. I think subs should be reimplemented using Hom ASAP.  Note that this could break existing code, so will have to be done carefully.    We can leave the old behavior in for speed, but as a non-default option.

3. Come up with a fast way to implement the new behavior. 

Attachments (1)

fix_mpoly_subs.patch (2.4 KB) - added by malb 12 years ago.

Download all attachments as: .zip

Change History (11)

comment:1 Changed 12 years ago by malb

the main use-case for which I wrote it is (which must be fast):

sage: R.<x,y>=QQ[]
sage: g=x+y
sage: %timeit g.subs({x:2,y:1})
10000 loops, best of 3: 62.9 µs per loop

The performance for elements in R is:

sage: %timeit g.subs({x:x+1,y:x*y})
10000 loops, best of 3: 153 µs per loop

However, to my surprise hom is faster:

sage: R.<x,y>=QQ[]
sage: phi = R.hom([x+1,x*y])
sage: g=x+y
sage: %timeit phi(g)
10000 loops, best of 3: 45.8 µs per loop

Is it because it caches or is really just better code?

comment:2 Changed 12 years ago by was

Is it because it caches or is really just better code?

Hom is implemented by Singular when the base ring is a singular ring.

sage: R.<x,y>=GF(next_prime(10^9))[]
sage: phi = R.hom([x+1,x*y])
sage: g=x+y
sage: print type(g)
sage: timeit('phi(g)')
<type 'sage.rings.polynomial.multi_polynomial_libsingular.MPolynomial_libsingular'>
625 loops, best of 3: 39.7 µs per loop
sage: R.<x,y>=GF(next_prime(10^10))[]
sage: phi = R.hom([x+1,x*y])
sage: g=x+y
sage: print type(g)
sage: timeit('phi(g)')
<class 'sage.rings.polynomial.multi_polynomial_element.MPolynomial_polydict'>
625 loops, best of 3: 305 µs per loop     

comment:3 Changed 12 years ago by malb

So it is _im_gens_ in MPolynomial_libsingular then (sorry for not checking myself earlier).

It seems the most straight-forward implementation possible including the comment #TODO: very slow in _im_gens_ wins for small examples, but for bigger ones we get:

sage: P.<a,b,c,d,e> = PolynomialRing(QQ)
sage: f = P.random_element(degree=3,terms=50)
sage: g = {a:b,b:c,c:d,d:e,e:a}
sage: %timeit f.subs(g)
10000 loops, best of 3: 138 µs per loop
sage: phi = P.hom([b,c,d,e,a])
sage: %timeit phi(f)
1000 loops, best of 3: 893 µs per loop

comment:4 Changed 12 years ago by malb

  • Summary changed from multivariate polynomial substitution has a design flaw to [with patch, needs review] multivariate polynomial substitution has a design flaw

Changed 12 years ago by malb

comment:5 Changed 12 years ago by malb

Performance

sage: P.<a,b,c,d,e> = PolynomialRing(QQ)
sage: f = P.random_element(degree=3,terms=50)
sage: g = {a:b,b:c,c:d,d:e,e:a}
sage: %timeit f.subs(g)
1000 loops, best of 3: 271 µs per loop
sage: phi = P.hom([b,c,d,e,a])
sage: %timeit phi(f)
1000 loops, best of 3: 939 µs per loop
sage: phi(f)
-a^2*b - 11/2*b^3 - 8*a^2*c + 1/51*b*c^2 + 1/2*c^3 - a^2*d + 1/3*a*b*d - 2/11*a*c*d + 195*b*c*d + 1/3*a*d^2 + 2*b*d^2 - c*d^2 - 2/3*a^2*e + 1/2*a*b*e - 203*b^2*e + 1/4*a*c*e + b*c*e - 5*c^2*e + 6*a*e^2 + b*e^2 - 1/3*c*e^2 - 5*d*e^2 + 3*e^3 + 1/3*a^2 - a*b - 7/48*a*c - 2*b*c - 53/2*c^2 - 1/3*a*d - 1/2*b*d + c*d - d^2 - a*e - 4*b*e - d*e + 13*e^2 - 2*a - 1/2*b - c + 9/2*d - 1/2
sage: f.sub
f.sub_m_mul_q  f.subs         f.substitute
sage: f.subs(g)
-a^2*b - 11/2*b^3 - 8*a^2*c + 1/51*b*c^2 + 1/2*c^3 - a^2*d + 1/3*a*b*d - 2/11*a*c*d + 195*b*c*d + 1/3*a*d^2 + 2*b*d^2 - c*d^2 - 2/3*a^2*e + 1/2*a*b*e - 203*b^2*e + 1/4*a*c*e + b*c*e - 5*c^2*e + 6*a*e^2 + b*e^2 - 1/3*c*e^2 - 5*d*e^2 + 3*e^3 + 1/3*a^2 - a*b - 7/48*a*c - 2*b*c - 53/2*c^2 - 1/3*a*d - 1/2*b*d + c*d - d^2 - a*e - 4*b*e - d*e + 13*e^2 - 2*a - 1/2*b - c + 9/2*d - 1/2

comment:6 Changed 12 years ago by jason

Note that since dictionaries do not preserve order, in python, there is no way to distinguish between {x:x+1,y:x*y} and {y:x*y, x:x+1} . There are classes out there that implement ordered dictionaries, but if the feature request depends on looping through a standard dictionary in the order that the dictionary was specified, that seems pretty impossible.

Note that this also means that there is no distinction between these calls: g(x=3, y=x) and g(y=x, x=3) .

comment:7 Changed 12 years ago by jason

  • Summary changed from [with patch, needs review] multivariate polynomial substitution has a design flaw to [with patch, needs work] multivariate polynomial substitution has a design flaw

To emphasize the point, here is the result after applying the patch:

sage: R.<x,y>=QQ[] 
sage: g=x+y 
sage: g.subs({x:x+1,y:x*y}) 
x*y + x + 1
sage: g.subs({y:x*y, x:x+1})
x*y + x + 1

Note that things are *not* in order in the second example.

I think this is probably hopeless as stated. If the substitutions were given as a list of tuples, then you could depend on the order. In other words, if you had something like g.subs([(y,x*y), (x,x+1)]) then you could say something about doing the substitutions in order. Or even if you did something like g.subs((y,x*y), (x,x+1)) you could do something in order, since *args is a list that preserves the order of the arguments.

comment:8 Changed 12 years ago by malb

It never even occurred to me that one could want to substitute in the order of the dictionary, so this is not what I tried to provide. 'ordered' to me only means 'by variables'. To me, the second example is in order but it might be a good idea to add a note to the docstring to clarify the behaviour.

comment:9 Changed 12 years ago by mhansen

  • Authors set to Martin Albrecht
  • Reviewers set to Mike Hansen
  • Summary changed from [with patch, needs work] multivariate polynomial substitution has a design flaw to [with patch, positive review] multivariate polynomial substitution has a design flaw

I think that this patch should get a positive review. The behavior after the patch is the correct behavior and is consistent with the rest of Sage.

I think Jason just misinterpreted what the ticket was for.

comment:10 Changed 12 years ago by mhansen

  • Merged in set to sage-4.2.alpha0
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.