Basis for multivariable polynomials

Question

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.

Answer 1

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}*)

Answer 2

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 = GroebnerBasis`DistributedTermsList[polylist, vars][[1]];
exponvecs = terms[[All, All, 1]];
exponbasis = Union[Flatten[exponvecs, 1]];
basisrules = Thread[exponbasis -> Range[Length[exponbasis]]];
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} *)

Answer 3

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.