Quadrics vanishing on orbits in Pluecker coordinates

As in Example 6.9, we consider the following nets:

\[\begin{split}\mathcal{L}_1=\begin{pmatrix}x&0&0&0\\0&x&0&0\\0&0&y&0\\0&0&0&z\end{pmatrix},\mathcal{L}_2=\begin{pmatrix}x&y&0&0\\y&z&0&0\\0&0&x&y\\0&0&y&z\end{pmatrix},\mathcal{L}_3=\begin{pmatrix}0&0&x&y\\0&-2x&-y&z\\x&-y&-2z&0\\y&z&0&0\end{pmatrix}\end{split}\]

In Example 6.10, we claim that 85 quadrics (in dual Plücker coordinates) vanish on the orbit of \(\mathcal{L}_1\), 189 quadrics vanish on the orbit of \(\mathcal{L}_2\), which also vanish on the orbit of \(\mathcal{L}_3\), and that 918 quadrics vanish on the orbit of \(\mathcal{L}_3\). We also give examples of such quadrics. The goal of this code is to verify these claims.

Our first goal is to parametrize the congruence orbits of \(\mathcal{L}_1,\mathcal{L}_2,\mathcal{L}_3\) in the Grasmannian.

# We generate variables for the action of GL(4)
R = PolynomialRing(QQ,['x'+str(i)+str(j) for i in range(1,5) for j in range(1,5)])
R.inject_variables()

# In stead of nets of symmetric matrices, we consider nets of quadrics.
# We generate their variables.
S.<a,b,c,d> = PolynomialRing(R,4)

# The congruence action replaces a,b,c,d by general linear combinations of themselves.
# We generate these.
[A,B,C,D] = [sum([R('x'+str(i)+str(j))*[a,b,c,d][i-1] for i in range(1,5)]) for j in range(1,5)]

# Next, we generate a basis for the net of quadrics
forms = [A^2+B^2,C^2,D^2]          # L1 has basis a^2+b^2, c^2, d^2
#forms = [A^2+C^2,A*B+C*D,B^2+D^2]  # L2 has basis a^2+b^2, a*b+c*d, b^2+d^2
#forms = [A*C-B^2,A*D-B*C,B*D-C^2]  # L3 has basis a*c-b^2, a*d-b*c, b*d-c^2

# We generate the matrix seen in the example.
M = Matrix([form.coefficients() for form in forms])
# Note that SAGE orders the indices differently:
# 11, 12, 22, 13, 23, 33, 14, 24, 34, 44
# Also, off-diagonal matrix entries are half of the corresponding quadric entry.
tmp = Matrix(QQ,[[0 for t1 in range(10)] for t2 in range(10)])
tmp[0,0] = 1
tmp[1,1] = 1/2
tmp[3,2] = 1/2
tmp[6,3] = 1/2
tmp[2,4] = 1
tmp[4,5] = 1/2
tmp[7,6] = 1/2
tmp[5,7] = 1
tmp[8,8] = 1/2
tmp[9,9] = 1
M = M*tmp

# We generate all dual Pluecker indices and compute the parametrisation in dual Pluecker coordinates.
plueckerindex = [[i,j,k] for i in range(10) for j in range(i+1,10) for k in range(j+1,10)]
plueckercoords = [Matrix([[M[s,t] for t in index] for s in range(3)]).det() for index in plueckerindex]

Defining x11, x12, x13, x14, x21, x22, x23, x24, x31, x32, x33, x34, x41, x42, x43, x44

Next, we need to generate all degree-\(2\) monomials in the Plücker coordinates that are linearly independent modulo the Plücker relations.

# We first compute the pair of indices these monomials correspond to. There are 4950 such pairs.
standardpairs = [[s,t] for s in range(120) for t in range(s,120) if plueckerindex[s][0]<=plueckerindex[t][0] if plueckerindex[s][1]<=plueckerindex[t][1] if plueckerindex[s][2]<=plueckerindex[t][2]]
standardproducts = [plueckercoords[s]*plueckercoords[t] for [s,t] in standardpairs]

We now generated a list of \(4950\) homogeneous degree-\(12\) polynomials in \(16\) variables. We want to know the linear dependencies between them.

As this problem is too big to tackle directly, we consider in stead, for each polynomial, whether the polynomial is a linear combination of the previous ones, and if so how.

def Dothings(polys,exprs,i):
    # input:  a list of polynomials,
    #         a list of expressions of those polynomials as linear combinations of standardproducts
    #         an index
    # output: IF the i-th polynomial of standardproducts if a linear combination of the previous ones,
    #           we return [polys,exprs,expr] where expr is the list of coefficients of the dependence.
    #         ELSE we return [polys+[f],exprs+[expr],[]] where f is the reduced version of
    #           the i-th polynomial of standardproducts and expr expresses f as as a linear combination
    #           of standardproducts.
    #           NOTE: We sort polys such that the leading monomials are a descending chain.
    #                 This ensures that we can reduce f fast!

    # We take f the i-th polynomial of standardproducts.
    f = standardproducts[i]
    # It equals the i-th polynomial of standardproducts.
    expr = [0 for k in range(4950)]
    expr[i] = 1

    # We reduce f with elements of polys.
    # We change expr accordingly.
    for j in range(len(polys)):
        g = polys[j]
        if g.lm()==f.lm():
            [f,expr] = [g.lc()*f-f.lc()*g,[g.lc()*expr[k]-f.lc()*exprs[j][k] for k in range(4950)]]
            tmp = gcd(expr)
            f = f/tmp
            expr = [expr[k]/tmp for k in range(4950)]
            if f==0:
                return [polys,exprs,[expr]]
        elif (f+g).lm()==f.lm():
            return [polys[:j]+[f]+polys[j:],exprs[:j]+[expr]+exprs[j:],[]]
    return [polys+[f],exprs+[expr],[]]

# We add the first polynomial to polys (since it is not zero).
polys = [standardproducts[0]]
# The first polynomial equals the first polynomial.
exprs = [[1]+[0 for k in range(1,4950)]]
# We record the linear dependencies.
res = []

# For each polynomial in the list, we check if it is a linear combination of the others.
# If so, add how to res. If not, reduce the polynomial, add it to polys (sorted) and
# keep track of what the reduced polynomial is in terms of the orginal ones.
for i in range(1,4950):
    [polys,exprs,tmp] = Dothings(polys,exprs,i)
    res = res + tmp
    # This takes long... print how far we are and how many dependecies we found.
    #print (i,len(res))

# Save results
f = open('L1res.txt', 'w') # change name to match L2, L3.
f.write(str(res))
f.close()

def toString(eqn):
    # input:  a depencency
    # output: the dependency as a string
    tmp = gcd(eqn)
    eqn = [eqn[k]/tmp for k in range(4950)]
    string = ''
    for l in range(4950):
        if eqn[l]!=0:
            sign = (eqn[l]>0)
            value = str(abs(eqn[l]))
            if value=='1':
                value=''
            [a,b] = standardpairs[l]
            [i,j,k] = plueckerindex[a]
            [i2,j2,k2] = plueckerindex[b]
            if string=='':
                if sign==False:
                    string='-'
                string = string + value +'p_{'+str(i)+str(j)+str(k)+'}p_{'+str(i2)+str(j2)+str(k2)+'}'
            else:
                if sign==False:
                    string = string + '-'
                else:
                    string = string + '+'
                string = string + value +'p_{'+str(i)+str(j)+str(k)+'}p_{'+str(i2)+str(j2)+str(k2)+'}'
    return string

We verify the quadric vansihing on the orbit of \(\mathcal{L}_1\).

print toString(res[60])

p_{124}p_{789}-p_{125}p_{589}+p_{125}p_{679}-p_{126}p_{678}+p_{127}p_{569}-p_{128}p_{568}-p_{234}p_{579}+p_{234}p_{678}+p_{235}p_{479}+p_{235}p_{568}-p_{236}p_{478}-p_{236}p_{567}-p_{237}p_{459}+p_{238}p_{458}

# Import the relations for L1
with open('./L1res.txt', 'r') as f:
    L = f.readlines()
res1 = eval(L[0])
X = Matrix(res1).transpose()
# NOTE: the matrix is in echelon form if we transpose.

# Import the relations for L2
with open('./L2res.txt', 'r') as f:
    L = f.readlines()
res2 = eval(L[0])
Y = Matrix(res2).transpose()
# NOTE: the matrix is in echelon form if we transpose.

# Import the relations for L3
with open('./L3res.txt', 'r') as f:
    L = f.readlines()
res3 = eval(L[0])
Z = Matrix(res3).transpose()
# NOTE: the matrix is in echelon form if we transpose.

We find the numbers in the paper.

print len(res1)
print len(res2)
print len(res3)

85
189
918

# Turn a degree-2 monomial into our format
def stdpair(i,j,k,i2,j2,k2):
    res = [0 for tmp in range(4950)]
    a = plueckerindex.index([i,j,k])
    b = plueckerindex.index([i2,j2,k2])
    res[standardpairs.index([a,b])] = 1
    return Matrix(res).transpose()

We verify the quadric vansihing on the orbit of \(\mathcal{L}_2\).

eqn = stdpair(0,1,2,0,1,8)-stdpair(0,1,2,0,2,6)-stdpair(0,1,2,0,3,5)-2*stdpair(0,1,2,1,2,3)-stdpair(0,1,3,0,1,7)+2*stdpair(0,1,3,0,2,5)-stdpair(0,2,3,0,2,4)
# Check that the quadric is a linear combination of the ones we found.
block_matrix([[Y,eqn]]).rank() == Y.rank()

True

We verify the quadric vansihing on the orbit of \(\mathcal{L}_3\).

# Generate the quadric on L3 from the paper.
eqn = 8*stdpair(0,1,2,4,5,7)-4*stdpair(0,4,5,0,5,7)+stdpair(0,4,7,0,4,7)
# Check that the quadric is a linear combination of the ones we found.
block_matrix([[Z,eqn]]).rank() == Z.rank()

True

Finally, we check that the quadrics that vanish on the orbit of \(\mathcal{L}_2\) also vanish on the orbit of \(\mathcal{L}_3\).

# Check that the equations on L2 are among the equation on L3.
block_matrix([[Y,Z]]).rank() == Z.rank()

True