Groebner examples

1414 days ago by hthall

Groebner basis calculations in Sage: three examples
# Pre-load zero forcing routines URL='https://raw.githubusercontent.com/jephianlin/mr_JG/py2/' files=['Zq_c.pyx','Zq.py','zero_forcing_64.pyx','zero_forcing_wavefront.pyx','minrank.py', 'inertia.py'] for f in files: print("Loading %s..."%f); load(URL+f) 
       

If you want to know what a Groebner basis actually is: "What is a Groebner Basis" (AMS Notices?)

Techniques to mention:

1) How to set up polynomial rings

2) Solving equations by changing term ordering or using elimination ideals

3) Reducing degrees of freedom by plugging in guesses

4) Breaking into cases with a complete primary decomposition

5) Generating a certificate via lifting

Under the hood of Sage
(1 + 1)^2 
       
preparse('(1 + 1)^2') 
       
sqrt 
       
1.parent() 
       
ZZ 
       
(2/3).parent() 
       
QQ 
       
preparse('2.3') 
       
(2.3).parent() 
       
RR 
       
CC 
       
QQbar 
       
AA 
       
(sqrt(2)).parent() 
       
x 
       
y 
       
var('y') 
       
y 
       
x.parent() 
       
y.parent() 
       

Question 1: What are all the points of intersection between the circle $$(x - 1)^2 + (y - 1)^2 = 1$$ and the parabola $$y = x^2?$$
# Define a polynomial ring R.<x, y> = QQ[] R 
       
# What does this shorthand actually do? preparse('R.<x, y> = QQ[]') 
       
# The variable x no longer belongs to the symbolic ring. x.parent() 
       
# Rewrite the equations as polynomial = 0, and create an ideal generated by those polynomials deal = ideal((x - 1)^2 + (y - 1)^2 - 1, y - x^2) 
       
deal 
       
# We can create a "polynomial sequence" (which looks like a list) to recover what the generating polynomials are deal.gens() 
       
# Most calculations using this ideal will trigger a Groebner basis computation. # Let's calculate it explicitly and have a look at it. groebner = deal.groebner_basis() groebner 
       

This Groebner basis is only slightly reduced from what we started with, and is not so helpful for solving the equations, since $x$ and $y$ are still mixed up together in both polynomials. (This is because the default term order for a polynomial ring in Sage, “degrevlex”, orders first by total degree, which is usually the fastest choice for calculations.)
# Does it at least have the expected dimension 0, meaning a finite set of isolated solutions? deal.dimension() 
       

Okay, so how do we actually solve it? Two approaches: Change the term order, or project down to an elimination ideal.
# If the example is small enough, you can use a degrevlex Groebner basis as a starting point for a more informative 'lex' Groebner basis. ideal(groebner).change_ring(R.change_ring(order='lex')).groebner_basis() 
       
# More generally (and usually more quickly), you can project the solution to a subspace # to find the restrictions on a smaller set of variables. xsolution = deal.elimination_ideal([y]) xsolution 
       
# Extracting the polynomial from the ideal xpol = xsolution.gens()[0] xpol 
       
# In this case, it does not factor over QQ. xpol.factor() 
       
# Although we got it by projecting away y, it still lives in the two-variable ring. xpol.parent() 
       
# We would like to know what all the solutions are . . . xpol.roots(ring=QQbar) 
       
# . . . but first we have to cast it to either a single-variable ring, or the symbolic ring. SR(xpol).roots(ring=QQbar) 
       

Question 2: Does $\mathcal{S}(C_4)$ contain a matrix $B$ with spectrum $(0, 0, 1, 1)$?

This is related to a question that came up during Shaun Fallat’s Combinatorial Matrix Theory session: For a tree $T$, $$q(T) \ge \mathrm{diam}(T) + 1.$$

For what other graphs is this true? Odd cycles? Even cycles?

It is not hard to show that $\mathrm{M}(C_k) = 2$ means that this inequality does hold for odd cycles $k = 3, 5, \dots$

So the first interesting case is $C_4$. For $B \in \mathcal{S}(C_4)$, can we have $q(B) = 2$?

Without loss of generality, can the spectrum of $B$ be $(0, 0, 1, 1)$?
G = graphs.CirculantGraph(4, [1]) G 
       
# Define a ring and define Python variables with its variable names P = PolynomialRing(QQ, ['a', 'b0', 'b1', 'c', 'd0', 'd1']) P.inject_variables() P 
       
a.parent() 
       
M = Matrix(P, [[a, 0], [b0, b1], [0, c], [d0, d1]]) show(M) 
       
B = M * M.transpose() show(B) 
       
G.adjacency_matrix() 
       
# One of the non-edge entries has a polynomial that is not identically zero. Graph(B) 
       
# So the first polynomial that we require to be 0 is that entry. deal = ideal(B[1, 3]) deal 
       
# If every eigenvalue is either 0 or 1, then since B is diagonalizable that gives the minimal polynomial. # We define a matrix all 16 of whose entries should be zero, and add those polynomials to the ideal. vanish = B * (B - 1) for ii in range(4): for jj in range(4): deal += ideal(vanish[ii, jj]) 
       
groebner = deal.groebner_basis() groebner 
       
# Hmm. Looks like a mess. The last polynomials in the list are the simplest ones. groebner[-1] 
       
groebner[-2] 
       
deal.dimension() 
       
# We can try plugging in arbitrary values for a few variables to reduce the degrees of freedom. bad_deal = deal + ideal(d0 - 3, d1 - 5) bad_deal.dimension() 
       
bad_deal.groebner_basis() 
       

Oops—that choice leads to imaginary values for $a$ and $c$. But maybe it is possible to break into simpler cases using a complete primary decomposition.
options = deal.complete_primary_decomposition() 
       
# options 
       
for option in options: for xx in option: print(list(xx.gens())) print 
       
options[2][1] 
       
pols = list(options[2][1].gens()) pols 
       
M2 = M(a=c, b1=-d0, b0=d1) M2 = M2(c=1/3, d0=2/3, d1=2/3) M2 = Matrix(QQ, M2) M2 
       
B2 = M2 * M2.transpose() B2 
       
Graph(B2) 
       
B2.eigenvalues() 
       

Question 3: For the graph $\mathrm{Circ} = $ CirculantGraph(11, [1, 3]), does $$\mathrm{M}(\mathrm{Circ}) = \mathrm{Z}(\mathrm{Circ})?$$
Circ = graphs.CirculantGraph(11, [1, 3]) Circ 
       
zs = zero_forcing_set_bruteforce(Circ) zs 
       

Useful fact: A zero forcing set of size $6$ guarantees a (non-principal) submatrix that is permutation equivalent to a necessarily invertible upper-triangular matrix of size $11 - 6 = 5$.
# Add two just to keep track of the diagonal entries, which might be either zero or nonzero. A = Circ.adjacency_matrix() + 2 A 
       
rowset = [6, 7, 8, 9, 10] colset = [3, 4, 5, 6, 7] A[rowset, colset] 
       
# Find a spanning tree of Circ that includes as many of the submatrix edges as possible (especially the diagonal entries that make it invertible) Tree = Circ.copy() Tree.name('Spanning tree of Circ') for xx in Tree.edges(): Tree.delete_edge(xx) Tree.add_edges([(6, 3), (7, 4), (8, 5), (9, 6), (10, 7), (6, 5), (8, 7), (0, 1), (1, 2), (2, 3)]) Tree 
       
Tree.is_tree() 
       

The question reduces to whether there exists any $B \in \mathcal{S}(\mathrm{Circ})$ whose rank is only $11 - \mathrm{Z}(\mathrm{Circ}) = 5$

Since we are only interested in rank rather than all of the spectrum, we can freely assume that all of the spanning tree edges have weight $1$.

For all other edges, and for the diagonal entries, we need variable names.
# Define names for the diagonal entries. dnames = ['d' + str(ii) for ii in range(11)] dnames 
       
# Define names for the off-diagonal non-spanning-tree edges. bnames = [] for xx in Circ.edges(): if xx not in Tree.edges(): bnames.append('b' + str(xx[0]) + '_' + str(xx[1])) bnames 
       
P = PolynomialRing(QQ, dnames + bnames) P.inject_variables() print P 
       
B = Matrix(P, Circ.adjacency_matrix()) for ii, xx in enumerate(P.gens()[:11]): B[ii, ii] = xx for xx in P.gens()[11:]: indices = str(xx)[1:].split('_') ii, jj = int(indices[0]), int(indices[1]) B[ii, jj] = xx B[jj, ii] = xx 
       
show(B) 
       

If $B$ is truly rank $5$, then any $6 \times 6$ submatrix must be singular. In fact, since our triangular submatrix is invertible, it suffices to check only those $6 \times 6$ which extend it by a single row and column.
# Collect a bunch of determinants giving polynomials that all must equal 0. deal = ideal(P) for ii in range(11): if ii in rowset: continue for jj in range(11): if jj in colset: continue deal += ideal(B[rowset + [ii], colset + [jj]].determinant()) 
       
# The moment of truth . . . deal.groebner_basis() 
       

We can just trust this answer, or we can ask Sage to provide us with a certificate that $1$ can be written as a combination of things all of which are supposed to equal $0$.
coeffs = P(1).lift(deal) coeffs 
       
# A routine to print out the linear combination of determinants that provides the certificate. just_started = True for ii, xx in enumerate(coeffs): if xx: if just_started: print '1 =', just_started = False else: print ' +', print '(', xx, ') * (', deal.gen(ii), ')' print