# Basis for multivariable polynomials

I have a bunch of two-variable polynomials and as part of a larger algorithm need to find a basis for them and express them in terms of this basis.

As an illustrative example, for one case my polynomials are:

• $$(1-x)(1-y)(y-x)$$,
• $$-y(1-x)(y-x)$$,
• $$y(1-x)(1-y)$$,
• $$-x(1-y)(y-x)$$,
• $$-x(1-x)(1-y)$$,
• $$xy(y-x)$$,
• $$-xy(1-y)$$,
• $$-xy(1-x)$$.

For this case it's easy to see by inspection that e.g. the sixth can be built from the seventh and eighth, and the first can be built from the third and fifth, so we can use the remaining six as a basis (and it's easy to check that this really is a basis e.g. using a matrix representation and row reduction).

However since I will be dealing with many different cases I would like a quick way to automate the process of finding a basis for such a set of polynomials and expressing them in terms of this basis.

• Do you need the basis vectors to come from this list, or would a GroebnerBasis be acceptable, because GroebnerBasis[ listOfPolynomials, {x, y}] returns {-y + y^2, x - y}. Feb 1 at 20:40
• @march basis needn't necessarily come from this list, but I don't understand how {-y+y^2, x-y} can possibly form a basis for this set of polynomials. I'm not entirely sure what a Groebner basis is but I am specifically after an additive basis in the usual sense for these polynomials (no need that they are specifically from the list).
– R.W
Feb 1 at 20:47
• So when you form linear combinations of these basis polynomials, do the coefficients need to be numbers (like integers or rationals, etc.), or are they also allowed to be polynomials? Groebner bases work for the latter (if I remember correctly). Feb 1 at 21:08
• @march They need to be constant numbers; they can't themselves be polynomial.
– R.W
Feb 1 at 23:27

You can use CoefficientRules to convert polynomial to a vector of coefficients. Below I use CoefficientRules on generic linear combination of original polynomials. Then everything is done by usual linear algebra tools, line RowReduce and NullSpace.

polys = {(1 - x) (1 - y) (y - x), -y (1 - x) (y - x),
y (1 - x) (1 - y), -x (1 - y) (y - x), -x (1 - x) (1 - y),
y (y - x), -y (1 - y), -y (1 - x)};
cs = Array[c, {Length@polys}];
crules = CoefficientRules[polys.cs, {x, y}];
vecs = Transpose@Outer[Coefficient, Last/@crules, cs];
vbasis = DeleteCases[RowReduce[vecs], {0 ..}];
coefs = Rest[#]/-First[#] &[First@NullSpace[Transpose@Prepend[vbasis, #]]] & /@ vecs;(*expansion coefficients of original polys via basis*)
pbasis = FromCoefficientRules[MapThread[Rule, {First/@ crules, #}], {x, y}] & /@ vbasis;(*polynomial basis*)

(*Check*)
Factor[coefs.pbasis - polys]

(*==> {0,0,0,0,0,0,0,0}*)

• Fantastic, thanks! This is pretty much what I had in mind but couldn't quite figure out the implementation.
– R.W
Feb 2 at 11:18
• Great! Could you maybe mark my answer as a solution? Feb 3 at 2:37

Not sure if this is what you want, but here is a way to get a vector space basis. We need to extract the different exponent vectors, create a matrix with each row representing a polynomial, augment with the actual polynomials, and reduce to row echelon form.

I use an internal function to obtain the exponent vectors and facilitate creating the matrix representation, but this can be done in other ways.

polylist = {(1 - x) (1 - y) (y - x), -y (1 - x) (y - x),
y (1 - x) (1 - y), -x (1 - y) (y - x), -x (1 - x) (1 - y),
y (y - x), -y (1 - y), -y (1 - x)};
vars = {x, y};
terms = GroebnerBasisDistributedTermsList[polylist, vars][[1]];
exponvecs = terms[[All, All, 1]];
exponbasis = Union[Flatten[exponvecs, 1]];
matrix0 = Map[# /. basisrules &, terms];
matrix =
Normal[SparseArray[
Flatten@MapIndexed[{#2[[1]], #1[[1]]} -> #1[[2]] &,
matrix0, {2}]]];
augmat = Join[matrix, Transpose[{Expand[polylist]}], 2]

(* Out[472]= {{1, -1, -1, 0, 1,
1, -1, -x + x^2 + y - x^2 y - y^2 + x y^2}, {0, -1, 0, 1, 1, 0, -1,
x y - x^2 y - y^2 + x y^2}, {1, -1, 0, -1, 1, 0, 0,
y - x y - y^2 + x y^2}, {0, 0, 0, -1, 1, 1, -1,
x^2 - x y - x^2 y + x y^2}, {0, 0, -1, 1, 0,
1, -1, -x + x^2 + x y - x^2 y}, {0, 1, 0, -1, 0, 0,
0, -x y + y^2}, {-1, 1, 0, 0, 0, 0, 0, -y + y^2}, {-1, 0, 0, 1, 0,
0, 0, -y + x y}} *)

redmat = RowReduce[augmat]

(* Out[462]= {{1, 0, 0, 0, 0, 0, -1, y - x^2 y}, {0, 1, 0, 0, 0,
0, -1, -x^2 y + y^2}, {0, 0, 1, 0, 0, 0, -1, x - x^2 y}, {0, 0, 0,
1, 0, 0, -1, x y - x^2 y}, {0, 0, 0, 0, 1,
0, -1, -x^2 y + x y^2}, {0, 0, 0, 0, 0, 1, -1, x^2 - x^2 y}, {0, 0,
0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0}} *)

basispolys = redmat[[All, -1]] /. 0 -> Nothing

(* Out[463]= {y - x^2 y, -x^2 y + y^2, x - x^2 y,
x y - x^2 y, -x^2 y + x y^2, x^2 - x^2 y} *)

• Thanks! It seems I should definitely learn about Groebner bases...
– R.W
Feb 2 at 11:17

I was going to tell you to use PolynomialReduce, but unfortunately, it's not so simple. In the code below, it expresses the first polynomial as the sum of the second, sixth, and negative seven, with the remainder $$-x+x^2$$:

poly = {(1 - x) (1 - y) (y - x), -y (1 - x) (y - x),
y (1 - x) (1 - y), -x (1 - y) (y - x), -x (1 - x) (1 - y),
y (y - x), -y (1 - y), -y (1 - x)}

PolynomialReduce[First@poly, Rest@poly, {x, y}]
(* {{1, 0, 0, 0, 1, -1, 0}, -x + x^2} *)

poly[[1]] - (poly[[2]] + poly[[6]] - poly[[7]]) // Simplify
(* (-1 + x)x *)


even though a simpler representation is possible, as you point out in your question:

poly[[1]] - poly[[3]] - poly[[5]] // Simplify
(* 0 *)


It seems that PolynomialReduce depends on the order of the polynomials in the second argument that is presented to it. The code below shows that out of all possible permutations of the remaining seven polynomials, PolynomialReduce can find a representation with a zero remainder only in about a fifth of the cases.

PolynomialReduce[First@poly, #, {x, y}] & /@
Permutations[Rest@poly] // CountsBy[Last]
(* <|-x + x^2 -> 3024, -x + y -> 1008, 0 -> 1008|> *)


You might consider this a bug - I suggest reporting this to Wolfram support. (Update: reported to Wolfram as [CASE:5002907]). PolynomialReduce behaves as designed according to Daniel's comment below.

• Thanks, and very interesting re the bug. So perhaps a naive fix for now (and not enormously problematic given the small sets I'll be working with -- certainly never more than eight polynomial) would be to iterate through the various permutations of the second argument.
– R.W
Feb 1 at 18:09
• Yes, or for bigger sets you can just select, say, a 1000 random permutations and hope for the best ;) Feb 1 at 18:21
• This is not a bug. PolynomialReduce` will not in general give a canonical result unless reducing against a Groebner basis. That was the whole idea to Buchberger's early work: solving the ideal membership problem in polynomial rings. Also see Scope section of documentation for this function. Feb 1 at 18:33
• Thanks for correcting me, @DanielLichtblau. Feb 1 at 18:35