Demo : bases of multipolynomials

This is a jupyter notebook demo for the package multipolynomial_bases. The package is installed by default on CoCalc. If you want to run it on your local sage install, you need to install the package as explained here

The documentation is available here

# Run this cell to load the librairy
from multipolynomial_bases import *

This package offers an implementation of multivariate polynomials seen as a multi bases algebra (it is not a classical implementation for symbolic computation on multivariate polynomials).

The monomial basis

The main algebra is created by:

A.<x> = MultivariatePolynomialAlgebra(QQ); A

The Multivariate polynomial algebra on x over Rational Field

As you can see, we don’t specify the number of variables. The name \(x\) does not refer to a single variable but to the infinite alphabet \(x_1, x_2, x_3, \dots\). Or from an implementation point of view: to the monomial bases of the algebra.

x

The Multivariate polynomial algebra on x over Rational Field on the monomial basis

A basis element of the monomial basis is indexed by a vector. So an element of the algebra is just a formal sum of vectors.

x[2,2,1] + x[3,4,2]

x[2, 2, 1] + x[3, 4, 2]

In the world of classical multivariate polynomials, the values of the vectors are the exponents of the monomial.

x[2,2,1].to_expr()

x1^2*x2^2*x3

(x[2,2,1] + x[3,4,2]).to_expr()

x1^3*x2^4*x3^2 + x1^2*x2^2*x3

When you multiply two monomials, we make the sum of the vectors.

x[2,2,1]*x[3,1,4]

x[5, 3, 5]

(x[2,2,1] + x[3,4,2]) * x[1,1,2]

x[3, 3, 3] + x[4, 5, 4]

The number of variables is adjusted depending on the size of the vectors.

x[1] + x[0,1] + x[0,0,1]

x[1, 0, 0] + x[0, 1, 0] + x[0, 0, 1]

x[1] == x[1,0,0]

True

x[0,1].to_expr()

x2

A.var(6)

x[0, 0, 0, 0, 0, 1]

A.var(6).to_expr()

x6

In the monomial basis, the polynomial are always expanded, if you want to factorize, you have to move to the usual multivariate polynomials.

p1 = (x[1] - x[0,1])**2
p1

x[2, 0] + x[0, 2] - 2*x[1, 1]

K.<x1,x2,x3> = QQ[]
p2 = K(p1.to_expr())
p2

x1^2 - 2*x1*x2 + x2^2

p2.factor()

(x1 - x2)^2

On the other hand, this implementation allows for easy computation of divided differences. The following polynomial corresponds to

\(x_1^2 x_2^4 x_3 + 2 x_1 x_2 + x_4^2\)

p = x[2,4,1] + 2*x[1,1] + x[0,0,0,4]

By applying \(\delta_2\) you get

p.divided_difference(2)

x[2, 3, 1, 0] + 2*x[1, 0, 0, 0] + x[2, 2, 2, 0] + x[2, 1, 3, 0]

Which is

\(x_1^2x_2^3x_3 + 2x_1 + x_1^2x_2^2x_3^2 + x_1^2x_2x_3^3\)

Schubert basis

Expanding in the monomials is just one way of working with polynomials. The purpose of this program is to be able to use other bases (as for symmetric functions). The Schubert basis is one of the bases of the multivariate polynomial algebra.

Y = A.schubert_basis()
Y

The Multivariate polynomial algebra on x over Rational Field on the Schubert basis of type A

Y.an_element()

2*Y[1, 0, 0] + Y[2, 2, 3] + Y[0, 0, 0] + 3*Y[0, 1, 0]

Each element is index by a vector, you can create a Schubert polynomial by entering the vector.

pol = Y[2,1,2]
pol

Y[2, 1, 2]

If you want to compute the expansion of this polynomial in a sum of monomials, you just convert it.

x(pol)

x[2, 2, 1] + x[3, 1, 1] + x[2, 1, 2]

Similarly, you can write any sum of monomials into a sum of Schubert polynomials.

Y(x[1,1,2] + x[2,3])

-Y[3, 2, 0] - Y[1, 2, 1] + Y[2, 3, 0] + Y[1, 1, 2]

But what are those vectors?

You may be used to index the Schubert polynomials with permutations. The vectors we use are called Lehmer codes, they are in direct correspondence with permutations. Let \(\sigma = \sigma_1 \sigma_2 \dots \sigma_n\) the one line notation of a permutation of size \(n\). The corresponding Lehmer code \(v = v_1 \dots v_n\) is given by \(v_i = \# \lbrace j > i ; \sigma_j < \sigma_i \rbrace\). In particular, the sum of \(v\) is the number of inversions of the permutation.

Sage knows how to compute Lehmer codes.

Permutation([5,2,1,4,3]).to_lehmer_code()

[4, 1, 0, 1, 0]

Permutation([3,2,5,4,1]).to_lehmer_code()

[2, 1, 2, 1, 0]

import sage.combinat.permutation as permutation
permutation.from_lehmer_code([4,1,0,1,0])

[5, 2, 1, 4, 3]

permutation.from_lehmer_code([2,1,2,1,0])

[3, 2, 5, 4, 1]

Why using Lehmer codes instead of permutations?

This package has been written following the notations of Alain Lascoux who had good reasons to find vectors (Lehmer codes) a better indexing set than permutations. Indeed, this way, you don’t restrict yourself to one specific symmetric group. Especially, when you multiply two Schubert polynomials indexed by permutations of \(S_n\) you might end up in \(S_m\) with \(m > n\). On the other end, the size of the vector won’t change.

Y[2,1,0] * Y[0,1]

Y[3, 1, 0] + Y[2, 2, 0]

The code \(\left[ 2,1,0 \right]\) corresponds to \(\sigma = 321\) and \(\left[0,1 \right]\) is \(132\). But then \(\left[ 3,1,0 \right]\) is \(4213\) and \(\left[ 2,2,0 \right]\) is \(3412\).

Actually, the number of variables used to expand the polynomial in the monomial basis is given by the position of the last non zero value.

x(Y[2,1])

x[2, 1]

Y[2,1].to_expr()

x1^2*x2

x(Y[0,1])

x[0, 1] + x[1, 0]

Y[0,1].to_expr()

x1 + x2

In particular, you have:

Y[2,1] == Y[2,1,0]

True

Y[2,1] == Y[2,1,0,0]

True

Also, Lehmer codes are a very natural indexing set for applying divided differences. If \(v_{i} > v_{i+1}\) the divided difference \(\delta_i\) does \(v_i = v_{i+1}\) and \(v_{i+1} = v_i - 1\). In this example \(\left[1,3,2 \right]\) becomes \(\left[1,2,2\right]\).

Y[1,3,2].divided_difference(2)

Y[1, 2, 2]

When \(v_i \leq v_{i+1}\) the result is 0.

Y[1,3,2].divided_difference(1)

0

To expand Schubert polynomials into the monomial, we use divided differences. You can check that the Lehmer code of the maximum permutation \(\omega_n\) is \(\left[n-1, n-2, \dots 0 \right]\) and the corresponding Schubert polynomial is given by \(x\left[n-1, n-2, \dots 0 \right]\).

x(Y[4,3,2,1,0])

x[4, 3, 2, 1, 0]

This is actually true for any vectors where values are weakly decreasing.

x(Y[4,4,2])

x[4, 4, 2]

x(Y[6,2])

x[6, 2]

The other polynomial are obtained recursively by applying divided differences.

\(Y[4,2,2,1] = \delta_2(Y[4,3,2,1])\)

\(Y[2,3,2,1] = \delta_1(Y[4,2,2,1])\)

x(Y[4,3,2,1]).divided_difference(2).divided_difference(1)

x[2, 3, 2, 1] + x[3, 2, 2, 1]

x(Y[2,3,2,1])

x[2, 3, 2, 1] + x[3, 2, 2, 1]

To apply multiple divided differences, you can also write

x(Y[4,3,2,1]).apply_reduced_word([2,1])

x[2, 3, 2, 1] + x[3, 2, 2, 1]

Remark: when you apply the maximal reduced word to the the Schubert polynomial indexed by the maximal permutation, you obtain \(Y[0,0,\dots,0] = 1\).

max_reduced_word = Permutation([5,4,3,2,1]).reduced_word()

Y[4,3,2,1,0].apply_reduced_word(max_reduced_word)

Y[0, 0, 0, 0, 0]

In general, if you apply the maximal reduced word of size \(n\) on a Shubert polynomial indexed by any decreasing vector of size \(n\), you get a Schubert polynomial indexed by a weakly increasing vector.

s = Y[6,4,3,1,0].apply_reduced_word(max_reduced_word)
s

Y[0, 0, 1, 1, 2]

This is actually a symmetric polynomial and more precisely, this is the expansion of a schur function indexed by the partition given by the vector.

x(s)

x[0, 2, 1, 1, 0] + x[0, 1, 0, 2, 1] + 3*x[1, 1, 1, 0, 1] + x[2, 0, 0, 1, 1] + x[1, 0, 2, 1, 0] + x[0, 1, 2, 1, 0] + x[0, 0, 2, 1, 1] + x[2, 1, 1, 0, 0] + x[0, 1, 2, 0, 1] + x[0, 0, 1, 2, 1] + x[1, 0, 1, 0, 2] + 3*x[1, 1, 0, 1, 1] + x[2, 1, 0, 1, 0] + x[1, 2, 0, 0, 1] + x[0, 1, 1, 2, 0] + 3*x[1, 1, 1, 1, 0] + x[0, 2, 0, 1, 1] + x[0, 1, 1, 0, 2] + x[1, 1, 0, 2, 0] + x[0, 0, 1, 1, 2] + x[0, 1, 0, 1, 2] + x[1, 0, 0, 1, 2] + x[1, 0, 0, 2, 1] + x[2, 0, 1, 1, 0] + x[1, 1, 0, 0, 2] + x[1, 1, 2, 0, 0] + 3*x[1, 0, 1, 1, 1] + x[2, 1, 0, 0, 1] + x[1, 0, 1, 2, 0] + x[1, 2, 1, 0, 0] + x[0, 2, 1, 0, 1] + x[2, 0, 1, 0, 1] + x[1, 2, 0, 1, 0] + 3*x[0, 1, 1, 1, 1] + x[1, 0, 2, 0, 1]

S = SymmetricFunctions(QQ)
schur = S.schur()

p1 = schur[2,1,1].expand(5)
p1

x0^2*x1*x2 + x0*x1^2*x2 + x0*x1*x2^2 + x0^2*x1*x3 + x0*x1^2*x3 + x0^2*x2*x3 + 3*x0*x1*x2*x3 + x1^2*x2*x3 + x0*x2^2*x3 + x1*x2^2*x3 + x0*x1*x3^2 + x0*x2*x3^2 + x1*x2*x3^2 + x0^2*x1*x4 + x0*x1^2*x4 + x0^2*x2*x4 + 3*x0*x1*x2*x4 + x1^2*x2*x4 + x0*x2^2*x4 + x1*x2^2*x4 + x0^2*x3*x4 + 3*x0*x1*x3*x4 + x1^2*x3*x4 + 3*x0*x2*x3*x4 + 3*x1*x2*x3*x4 + x2^2*x3*x4 + x0*x3^2*x4 + x1*x3^2*x4 + x2*x3^2*x4 + x0*x1*x4^2 + x0*x2*x4^2 + x1*x2*x4^2 + x0*x3*x4^2 + x1*x3*x4^2 + x2*x3*x4^2

# let's make both our expansions live in the same world
var("x0 x1 x2 x3 x4")
K.<x0,x1,x2,x3,x4> = QQ[]
p1 = K(p1)
p1

x0^2*x1*x2 + x0*x1^2*x2 + x0*x1*x2^2 + x0^2*x1*x3 + x0*x1^2*x3 + x0^2*x2*x3 + 3*x0*x1*x2*x3 + x1^2*x2*x3 + x0*x2^2*x3 + x1*x2^2*x3 + x0*x1*x3^2 + x0*x2*x3^2 + x1*x2*x3^2 + x0^2*x1*x4 + x0*x1^2*x4 + x0^2*x2*x4 + 3*x0*x1*x2*x4 + x1^2*x2*x4 + x0*x2^2*x4 + x1*x2^2*x4 + x0^2*x3*x4 + 3*x0*x1*x3*x4 + x1^2*x3*x4 + 3*x0*x2*x3*x4 + 3*x1*x2*x3*x4 + x2^2*x3*x4 + x0*x3^2*x4 + x1*x3^2*x4 + x2*x3^2*x4 + x0*x1*x4^2 + x0*x2*x4^2 + x1*x2*x4^2 + x0*x3*x4^2 + x1*x3*x4^2 + x2*x3*x4^2

p2 = K(x(s).to_expr(alphabet=[x0,x1,x2,x3,x4]))

p2

x0^2*x1*x2 + x0*x1^2*x2 + x0*x1*x2^2 + x0^2*x1*x3 + x0*x1^2*x3 + x0^2*x2*x3 + 3*x0*x1*x2*x3 + x1^2*x2*x3 + x0*x2^2*x3 + x1*x2^2*x3 + x0*x1*x3^2 + x0*x2*x3^2 + x1*x2*x3^2 + x0^2*x1*x4 + x0*x1^2*x4 + x0^2*x2*x4 + 3*x0*x1*x2*x4 + x1^2*x2*x4 + x0*x2^2*x4 + x1*x2^2*x4 + x0^2*x3*x4 + 3*x0*x1*x3*x4 + x1^2*x3*x4 + 3*x0*x2*x3*x4 + 3*x1*x2*x3*x4 + x2^2*x3*x4 + x0*x3^2*x4 + x1*x3^2*x4 + x2*x3^2*x4 + x0*x1*x4^2 + x0*x2*x4^2 + x1*x2*x4^2 + x0*x3*x4^2 + x1*x3*x4^2 + x2*x3*x4^2

p1 == p2

True

So, now I’d like to manipulate your polynomials and play with them, how do I do that?

Many operations are implemented, you can check the documentation to know more. Also, it is always possible to just to get list of index vectors and coefficients and just do whatever you want with those.

p = Y[1,2,1] + 2*Y[3,2] - Y[2,2,3]

list(p)

[([3, 2, 0], 2), ([1, 2, 1], 1), ([2, 2, 3], -1)]

# getting the set of permutations along with multiplicity that appear
import sage.combinat.permutation as permutation
p = p.change_nb_variables(6) # to have enough 0 at the end of the Lehmer codes so that they correspond to permutations of S6
{permutation.from_lehmer_code(k): c for k,c in p}

{[2, 4, 3, 1, 5, 6]: 1, [3, 4, 6, 1, 2, 5]: -1, [4, 3, 1, 2, 5, 6]: 2}