# Ideal.variety() when working with symbolic rings

I have a system of multivariate polynomials I would like to solve (the equations come from chemical reaction networks). I have many symbolic constants, but I've managed to generate a Groebner basis with the following code (using a toy example with solutions {{x: x0, y:0},{x: 0,y: x0}):

x0 = var('x0') # In practice, there are many constants

S.<x,y> = PolynomialRing(SR,order='lex'); # in practice, there are 3-10 variables

I = S.ideal(
x**2 + y**2 - x0**2,
x + y - x0
)

print(I.dimension()) # dimension should always be zero
# >> 0

print(I.groebner_basis())
# >> [x + y - x0,
# >>  y^2 + (-x0)*y]


In practice, I get many resulting equations which are extremely long, and it is numerically challenging to sequentially find the roots of these polynomials (in downstream code, I will get values for the constant parameters and then need to return solutions to these equations on the fly). Using the I.variety() method would be great to symbolically isolate all possible solutions, but I get confusing internal errors when I try to evaluate it.

print(I.variety()) # should be equivalent to: I.variety(ring=SR)

#/Users/ksb/SageMath/local/lib/python2.7/site-packages/sage/rings/polynomial/multi_polynomial_ideal.pyc in _variety(T, V, v)
#   2340
#   2341             variable = f.variable(0)
#-> 2342             roots = f.univariate_polynomial().roots(ring=ring, multiplicities=False)
#   2343
#   2344             for root in roots:
#
#/Users/ksb/SageMath/local/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_element.pyx in sage.rings.polynomial.polynomial_element.Polynomial.roots (build/cythonized/sage/rings/polynomial/polynomial_element.c:68811)()
#   7714                     return l
#   7715                 else:
#-> 7716                     return [val for val,m in l]
#   7717             vname = 'do_not_use_this_name_in_a_polynomial_coefficient'
#   7718             var = SR(vname)

# TypeError: 'sage.symbolic.expression.Expression' object is not iterable


Is this expected behavior? When I try f.univariate_polynomial().roots(ring=ring, multiplicities=False) (setting multiplicities=True is what I.variety() does and leads to the error) I get [1/2*x0,2] for roots, which makes little sense (it seems like the roots should be [0,-x0]).

I would be happy to convert the MPolynomial_polydict (or the sage.rings.polynomial.polynomial_ring.PolynomialRing_field_with_category.element_class generated by univariate_polynomial()) into a symbolic expression that I could just set it equal to zero and use solve on it, but I don't know how to do this type conversion.

edit retag close merge delete

copying my comments to tmoteil here:

A less trivial problem looks like this:

R.<x0,x1,x2,x3,x4,x5,x6,x7,x8> = PolynomialRing(QQ)
F=R.fraction_field()
S.<x9,x10,x11> = PolynomialRing(F,order='lex');

I = S.ideal(
-x7*x10*x9 + x8*x2*x11**2 + x3*x0*x11 - x4*x10,
-x7*x10*x9 + x8*x2*x11**2 + 2*x5*x1*x11**2 - 2*x6*x9**2,
x9 + x10 + x11 - 1
)
print(I.dimension()) # 0
B = I.groebner_basis()
print('\n\n'.join(map(str,B)))

print(I.variety())


When I do this, I get an empty list as the result.

I suppose the problem may be that there are no solutions that are valid over all possible values of the constants? Here's an example of there existing a solution if we substitute 1.0 (or 0.5) for all the constants:

S.<x11,x10,x9> = PolynomialRing(QQ,order='lex');

I = S.ideal(
-x10*x9 + x11**2 + x11 - x10,
-x10*x9 + x11**2 + x11**2 - x9**2,
x9 + x10 + x11 - 1
)
print(I.dimension())
B = I.groebner_basis()
print('\n\n'.join(map(str,B)))

print(I.variety())


The solution is {1/3,1/3,1/3}

Secretly I know that the value for every constant will be a positive number. And, while there may be lots of possible (x9,x10,x11) triples that satisfy the equations, the only ones that are "physically meaningful" are those with all three values between 0 and 1 (inclusive), possibly irrational. The fact that I cannot find any physically meaningful solutions after taking the Groebner basis and doing the rootfinding myself is what alerted me to a problem.

I'm not sure if these physically-motivated assumptions of mine are at the root of my confusion.

Sort by » oldest newest most voted

Instead of using the symbolic ring, you can view x0 as an undeterminate in a polynomial ring (or in fact, a fracion field):

sage: R.<x0> = PolynomialRing(QQ)
sage: R
Univariate Polynomial Ring in x0 over Rational Field
sage: F = R.fraction_field()
sage: F
Fraction Field of Univariate Polynomial Ring in x0 over Rational Field
sage: S.<x,y> = PolynomialRing(F, order='lex')
sage: S
Multivariate Polynomial Ring in x, y over Fraction Field of Univariate Polynomial Ring in x0 over Rational Field
sage: I = S.ideal(
....:         x**2 + y**2 - x0**2,
....:         x + y - x0
....:         )
sage: I.variety()
verbose 0 (2088: multi_polynomial_ideal.py, variety) Warning: falling back to very slow toy implementation.
[{y: 0, x: x0}, {y: x0, x: 0}]


However, i am not sure if this will scale to your concrete examples (it could be nice to provide them), since the "field with constants" you are working on requires that we use a generic slow implementations for the computation of groebner bases.

more

I haven't yet adapted the very large chemical networks to this framework, but a less trivial problem looks like this:

R.<x0,x1,x2,x3,x4,x5,x6,x7,x8> = PolynomialRing(QQ)
F=R.fraction_field()
S.<x9,x10,x11> = PolynomialRing(F,order='lex');

I = S.ideal(
-x7*x10*x9 + x8*x2*x11**2 + x3*x0*x11 - x4*x10,
-x7*x10*x9 + x8*x2*x11**2 + 2*x5*x1*x11**2 - 2*x6*x9**2,
x9 + x10 + x11 - 1
)
print(I.dimension()) # 0
B = I.groebner_basis()
print('\n\n'.join(map(str,B)))

print(I.variety())


When I do this, I get an empty list as the result.

I suppose the problem may be that there are no solutions that are valid over all possible values of the constants? Here's an example of there existing a solution if we substitute 1.0 (or 0.5) for all the constants:

S.<x11,x10,x9> = PolynomialRing(QQ,order='lex');

I = S.ideal(
-x10*x9 + x11**2 + x11 - x10,
-x10*x9 + x11**2 + x11**2 - x9**2,
x9 + x10 + x11 - 1
)
print(I.dimension())
B = I.groebner_basis()
print('\n\n'.join(map(str,B)))

print(I.variety())


The solution is {1/3,1/3,1/3}

Secretly I know that the value for every constant will be a positive number. And, while there may be lots of possible (x9,x10,x11) triples that satisfy the equations, the only ones that are "physically meaningful" are those with all three values between 0 and 1 (inclusive), possibly irrational. The fact that I cannot find any physically meaningful solutions after taking the Groebner basis and doing the rootfinding myself is what alerted me to a problem.

I'm not sure if these physically-motivated assumptions of mine are at the root of my confusion.