5
$\begingroup$

Having a polynomial $f(x,y)$, I would like to compute the following quantity \begin{equation*} {\mathbb C}[X,Y,Z]/\langle f_{x}, f_{y}, f_{z} \rangle, \end{equation*} where $f_{x},f_{y},f_{z}$ are, respectively, derivatives of $f$ with respect to $x$, $y$ and $z$.

For $f=2x^2 + 4 y^3 + 5 z^5$, $\langle f_{x}, f_{y}, f_{z} \rangle =\{ x,y^2,z^4\} $. I have implemented this as

f = 2 x^2 + 4 y^3 + 5 z^5
vars = {x, y, z}
fp = D[f, {vars}]
monomialList[poly_, vars_] := 
 Times @@ (vars^#) & /@ CoefficientRules[poly, vars][[All, 1]]
Flatten[Table[monomialList[fp[[i]], {x, y, z}], {i, 1, 3, 1}]]

The next step is to obtain the basis of the quotient space, which should be $\{1,z,z^2,z^3,y,yz,yz^2,yz^3\}$. Do you have any suggestions on computing this set in Mathematica? As I should carry out these calculations for several polynomials, having a generic answer for this question is ideal.

Update: the length of this set is known as the Milnor number. Here is a screenshot of three different polynomials whose Milnor numbers are obtained using the Singular package. enter image description here

$\endgroup$
11
  • 2
    $\begingroup$ For harder cases you can crib from code for counting the solutions. Just toss the final line that takes the length of the normal set, and instead return the normal set itself. $\endgroup$ Commented Dec 24, 2023 at 17:24
  • $\begingroup$ @DanielLichtblau This is very interesting. I have test that on ``` poly = x z + y^2 Length[basisSet[poly]] solutionSetSize[poly, Variables[poly]]``` and I got an empty set and an $\infty$. I again get $\infty$ in the output for the example polynomial in my question. Do you know why? $\endgroup$
    – Shasa
    Commented Dec 24, 2023 at 17:32
  • $\begingroup$ @DanielLichtblau (+1) This is very interesting. If I apply it to the jacobian of my polynomial, it generates all the exponents of the set. I should figure out how to bring these exponents to polynomials. Thank you for sharing this exciting solution. $\endgroup$
    – Shasa
    Commented Dec 24, 2023 at 17:43
  • 1
    $\begingroup$ Sorry, I forgot it was the exponent vectors that I generate. You can do something like Map[vars^#&, exponvecs] (apologies for not formatting or testing, I’m traveling with just a phone). $\endgroup$ Commented Dec 24, 2023 at 18:31
  • $\begingroup$ Also I’ll look into the misbehaving example when I’m back at work. $\endgroup$ Commented Dec 24, 2023 at 18:32

3 Answers 3

5
$\begingroup$

Here is code for producing a monomial basis set for a polynomial algebra with a finite dimensional vector space. It is taken from the resource function CountPolynomialSolutions with very minor modification to return the basis set.

I will point out that this basis set is not unique. The one returned arises from a particular Groebner basis computation on the input polynomials. Also, in the case where the vector space generators are not finite it returns {1}. And if the system is inconsistent it returns {}.

Clear[PolynomialAlgebraGenerators]

headsTable[elist_List, maxlist_, nvars_Integer] /; 
  Length[elist] == nvars :=
 Module[{jj, indices, iterators},
  (* Easy case: 
  all exponent vectors surrounding the normal set are "pure," that is,
  unmixed. *)
  indices = Array[jj, {nvars}];
  iterators = Table[{indices[[j]], 0, maxlist[[j]] - 1}, {j, nvars}];
  Flatten[Table[indices, Evaluate[Apply[Sequence, iterators]]], 
   nvars - 1]]

getRoofTuples[elist_, height_] := 
 Module[{len = Length[elist], expvec, expvecs, tt, rest},
  (* Tricky but efficient method. 
  Similar to that of Corless reference. *)
  expvecs = tt[];
  Do[expvec = elist[[jj]];
   If[expvec[[height]] =!= 
      0 && (height === 1 || Max[Take[expvec, height - 1]] === 0),
    expvecs = tt[expvec, expvecs]], {jj, len}];
  expvecs = Apply[List, Flatten[expvecs, Infinity, tt]];
  rest = Complement[elist, expvecs];
  expvecs = Map[Drop[#, height - 1] &, expvecs];
  {expvecs, rest}]

isUnder[t1_, t2_] := Apply[And, Thread[t1 <= t2]]

hitRoof[tup_, roofs_] := 
 Catch[Module[{jj, len = Length[roofs]}, 
   For[jj = 1, jj <= len, jj++, 
    If[isUnder[roofs[[jj]], tup], Throw[True, "hR"]]];
   False], "hR"]

headsTable[elist_List, maxes_, nvars_Integer] := 
 Module[{tuples, ntups, tuple, ntup, mixedexpons, tt, roofs}, 
  mixedexpons = 
   Select[elist, (Max[Apply[Sequence, #]] != Apply[Plus, #]) &];
  tuples = Table[{j}, {j, 0, maxes[[nvars]] - 1}];
  Do[{roofs, mixedexpons} = getRoofTuples[mixedexpons, jj];
   ntups = tt[];
   Do[tuple = tuples[[kk]];
    ntups = tt[Prepend[tuple, 0], ntups];
    Do[ntup = Prepend[tuple, mm];
     If[hitRoof[ntup, roofs], Break[], ntups = tt[ntup, ntups]];
     , {mm, maxes[[jj]] - 1}]
    , {kk, Length[tuples]}];
   tuples = Flatten[ntups, Infinity, tt], {jj, nvars - 1, 1, -1}];
  Sort[Apply[List, tuples]]]

Options[PolynomialAlgebraGenerators] = {WorkingPrecision -> Automatic};

PolynomialAlgebraGenerators[polys_List, 
  opts : OptionsPattern[PolynomialAlgebraGenerators]] :=
 PolynomialAlgebraGenerators[polys, Variables[polys], opts]

PolynomialAlgebraGenerators[opoly_?(! ListQ[#] &), 
  var_?(! ListQ[#] &), 
  opts : OptionsPattern[PolynomialAlgebraGenerators]] := 
 PolynomialAlgebraGenerators[{opoly}, {var}, opts]

PolynomialAlgebraGenerators[opolys_, vars_List, 
   OptionsPattern[PolynomialAlgebraGenerators]] /; 
  AllTrue[Map[PolynomialQ[#, vars] &, opolys], TrueQ] := 
 Catch[Module[
   {polys = opolys, domain, wprec, params, gb, dtl, leadvecs, 
    pureterms, nv = Length[vars], orderedvars, maxexpons, expons, 
    monorder = (MonomialOrder -> DegreeReverseLexicographic)},
   wprec = OptionValue[WorkingPrecision];
   If[wprec === Automatic,
    If[Internal`EffectivePrecision[polys] === Infinity,
     domain = CoefficientDomain -> RationalFunctions
     ,(*else finite precision, 
     so substitute numeric values for any parameters*)
     wprec = Internal`EffectivePrecision[polys];
     domain = CoefficientDomain -> InexactNumbers[wprec];
     params = 
      DeleteCases[Complement[Variables[polys], vars], _?NumericQ];
     polys = 
      polys /. 
       Thread[params -> 
         RandomReal[{-10, 10}, Length[params], 
          WorkingPrecision -> wprec]];
     ]
    ,(*else not Automatic*)
    If[Internal`EffectivePrecision[polys] === Infinity && 
       wprec === Infinity,
      domain = CoefficientDomain -> RationalFunctions
      ,
      If[wprec === MachinePrecision, wprec = $MachinePrecision];
      domain = CoefficientDomain -> InexactNumbers[wprec];
      params = 
       DeleteCases[Complement[Variables[polys], vars], _?NumericQ];
      polys = 
       polys /. 
        Thread[params -> 
          RandomReal[{-10, 10}, Length[params], 
           WorkingPrecision -> wprec]];
      ];
    ];
   orderedvars = 
    Last[GroebnerBasis`DistributedTermsList[polys, vars, Sort -> True,
       domain, monorder]];
   gb = GroebnerBasis[polys, orderedvars, Sort -> False, domain, 
     monorder];
   If[Length[gb] == 1 && NumericQ[gb[[1]]], Throw[{}, "CPS"]];
   Quiet[
    dtl = First[
      GroebnerBasis`DistributedTermsList[gb, orderedvars, 
       Sort -> False, domain, monorder]]];
   leadvecs = Map[#[[1, 1]] &, dtl];
   pureterms = Select[leadvecs, Count[#, a_ /; a =!= 0] === 1 &];
   If[Length[pureterms] < nv, Throw[{1}, "CPS"]];
   maxexpons = Apply[Plus, pureterms];
   expons = headsTable[leadvecs, maxexpons, nv];
   Map[Inner[Power, orderedvars, #, Times] &, expons]
   ], "CPS"]

Here is the original example in the post.

f = 2 x^2 + 4 y^3 + 5 z^5;
vars = {x, y, z};
fp = D[f, {vars}];
gens = PolynomialAlgebraGenerators[fp, vars]

(* Out[81]= {1, z, z^2, z^3, y, y z, y z^2, y z^3} *)

Here is a case where apparently we disagree with the program Singular (assuming I transcribed the system correctly).

poly = x^7 + y^7 + (x - y)^2*x^2*y^2 + z^2;
derivs = Grad[poly, {x, y, z}];
gens = PolynomialAlgebraGenerators[Expand@derivs, {x, y, z}]
Length[gens]

(* Out[96]= {1, y, y^2, y^3, y^4, y^5, x, x y, x y^2, x y^3, x y^4, 
 x y^5, x^2, x^2 y, x^2 y^2, x^2 y^3, x^2 y^4, x^2 y^5, x^3, x^3 y, 
 x^3 y^2, x^3 y^3, x^3 y^4, x^3 y^5, x^4, x^4 y, x^4 y^2, x^4 y^3, 
 x^4 y^4, x^4 y^5, x^5, x^5 y, x^5 y^2, x^5 y^3, x^5 y^4, x^5 y^5}

Out[97]= 36 *)

This is in accord with CountPolynomialSolutions (as it better be).

ResourceFunction["CountPolynomialSolutions"][derivs == 0, {x, y, z}]

(* Out[73]= 36 *)

But as shown in the original post, Singular claims the count should be 28.

$\endgroup$
1
  • $\begingroup$ (+1)Thank you for sharing this solution and for adding further comparisons $\endgroup$
    – Shasa
    Commented Jan 1 at 7:54
4
$\begingroup$

I do't think @Domen's answer is quite correct. Take, for example,

basisSet[(x + y)^3]
(*    {1, y, x, x y}    *)

which I don't think is a correct answer: the argument $(x+y)^3$ is cubic in the variables and therefore the resulting basis polynomials should all be linear in the variables. But $x y$ is quadratic, hence something is wrong.

I'd propose a slightly different method that gives the same results as @Domen's in many cases, but gets the correlations between variables more correct:

basisSetR[f_] := Module[{F, p, b, e},
    (* the list of derivative polynomials *)
    F = D[f, {Variables[f]}];
    (* the LCM of these, decomposed into factors and their powers *)
    p = DeleteCases[FactorList[PolynomialLCM @@ F], {_?NumericQ, _}];
    (* make all combinations of these factors *)
    If[p == {}, {1},
        {b, e} = Transpose[p];
        (Times @@ (b^#)) & /@ Tuples[Range[0, #] & /@ (e - 1)]]]

The OP's example is solved correctly:

basisSetR[2 x^2 + 4 y^3 + 5 z^5]
(*    {1, z, z^2, z^3, y, y z, y z^2, y z^3}    *)

Let's try the above example again: it now gives only zeroth- and first-order basis polynomials,

basisSetR[(x + y)^3]
(*    {1, x + y}    *)

Some more edge cases:

basisSetR[x z + y^2]
(*    {1}    *)

basisSetR[1]
(*    {1}    *)

basisSetR[x]
(*    {1}    *)

basisSetR[x^2]
(*    {1}    *)
$\endgroup$
5
  • $\begingroup$ Thank you for your answer. I also think that Domen's answer should be improved. Daniel's suggested link looks better for more generic cases. I have tried your answer on poly = x^3 + x y^2 the output should be $\{1,x,y,x^2\}$ or $\{ 1,x,y,y^{2}\}$ [This is what I get from Daniel's suggetion]. From your code, I get $\{1 \}$. But your answer beats Daniel's suggestion when it comes to $poly=(x+y)^3$ as I get $\infty$ from Daniel's suggestion. $\endgroup$
    – Shasa
    Commented Dec 25, 2023 at 9:48
  • $\begingroup$ Hmm yes, my solution needs improvement. $\endgroup$
    – Roman
    Commented Dec 25, 2023 at 10:27
  • $\begingroup$ But first, I think I need some more clarity. What is the correct solution for $(x+y)^2$? Shouldn't the correct solution be $\{1, x-y\}$ instead of just $\{1\}$? $\endgroup$
    – Roman
    Commented Dec 25, 2023 at 10:50
  • $\begingroup$ The length of the final set is the Milnor number, which can become infinite. I have checked the Milnor number of $(x+y)^2$ in the SINGULAR package, and it gives $\infty$. I have tried two more examples in the screenshot added to the main question. $\endgroup$
    – Shasa
    Commented Dec 25, 2023 at 11:21
  • 2
    $\begingroup$ Hold on. If the polynomial is not squarefree then the gradient vanishes on a dimensional set, so infinity is correct there. $\endgroup$ Commented Dec 25, 2023 at 16:54
2
$\begingroup$
basisSet[f_] := Module[{vars, fp, maxExponents},
  vars = Variables[f];
  fp = D[f, {vars}]; 
  maxExponents = MapThread[Exponent[#1, #2] &, {fp, vars}];
  Times @@@ (vars^# & /@ Tuples[Range[0, # - 1] & /@ maxExponents])
  ]

f = 2 x^2 + 4 y^3 + 5 z^5;
basisSet[f]
(* {1, z, z^2, z^3, y, y z, y z^2, y z^3} *)
$\endgroup$
5
  • $\begingroup$ (+1) This is great. Thank you, Domen. $\endgroup$
    – Shasa
    Commented Dec 24, 2023 at 16:28
  • $\begingroup$ Please check for some edge cases ... For example, what should be the output of basisSet[x]? Is {} correct or should it be {1}? $\endgroup$
    – Domen
    Commented Dec 24, 2023 at 16:29
  • $\begingroup$ Sure! Thanks for pointing this out. $\endgroup$
    – Shasa
    Commented Dec 24, 2023 at 16:30
  • $\begingroup$ @Shasa, please unselect my answer as accepted – I don't think I can improve it in any simple way to make it work for all cases :) $\endgroup$
    – Domen
    Commented Dec 25, 2023 at 10:16
  • $\begingroup$ I have unaccepted that. Thank you, @Domen, anyway, for your effort. $\endgroup$
    – Shasa
    Commented Dec 25, 2023 at 10:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.