Mathematica 9 is now available
Services & Resources / Wolfram Forums / MathGroup Archive
-----

MathGroup Archive 2011

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Algebraic substitution with PolynomialReduce

  • To: mathgroup at smc.vnet.net
  • Subject: [mg117119] Re: Algebraic substitution with PolynomialReduce
  • From: Daniel Lichtblau <danl at wolfram.com>
  • Date: Thu, 10 Mar 2011 06:09:21 -0500 (EST)

Guido Walter Pettinari wrote:
> Dear group,
> 
> in order to perform algebraic substitutions in expressions, I usually
> rely on replacementFunction, that is a very nice wrapper to
> PolynomialReduce made by Daniel Lichtblau and Andrzej Kozlowski (see
> links below).
> 
> However, I realized that replacementFunction is not enough when
> dealing with multiple substitutions.  Take this expression (which
> itself is a sub-expression of a more complicated one in my code):
> 
> terms = (2 kx1^2 kx2^2)/3 - (kx2^2 ky1^2)/3 + 2 kx1 kx2 ky1 ky2 - (
>  kx1^2 ky2^2)/3 + (2 ky1^2 ky2^2)/3 - (kx2^2 kz1^2)/3 - (
>  ky2^2 kz1^2)/3 + 2 kx1 kx2 kz1 kz2 + 2 ky1 ky2 kz1 kz2 - (
>  kx1^2 kz2^2)/3 - (ky1^2 kz2^2)/3 + (2 kz1^2 kz2^2)/3
> 
> If I define  k1 = {kx1, ky1, kz1}  and k2 = {kx2, ky2, kz2},  then the
> above is equal to the following combination of scalar products:
> 
> -(1/3 k1.k1 k2.k2 - (k1.k2)^2)
> 
> In order to explicitly reproduce the equality (and thus simplify the
> equation), I tried to apply replacementFunction several times, with no
> success.  So I turned to PolynomialReduce.  However, the following:
> 
> polylist = { k1.k1 k2.k2, (k1.k2)^2 };
> PolynomialReduce[ terms, polylist, Join[k1, k2] ]  // Last
> 
> did not return a zero remainder.
> 
> By reading old messages in the mailing list, I see that this may be
> due to the ordering of the variables in the third argument, i.e.
> Join[ k1,k2 ].  Using  GroebnerBasis, I managed to prove that 'terms'
> can be expressed in terms of the polynomials in 'polylist'.  In fact:
> 
> gb = GroebnerBasis[polylist, Join[k1, k2] ];
> PolynomialReduce[terms, gb, Join[k1, k2] ] // Last
> 
> yields zero.
> 
> The question is:  how can I obtain the coefficients of my expression
> ('term') with respect to my polynomials ('polylist'),  given that I
> already have the coefficients of 'terms' with respect to the Groebner
> basis of 'polylist'?
> 
> Thank you very much for your consideration!
> 
> Best wishes,
> 
> Guido W. Pettinari
> 
> 
> 
> REPLACEMENT FUNCTION LINKS:
> http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/thread/634a54e2e844837f/ddf3803cab4f3a5f?lnk=gst&q=replacementfunction#ddf3803cab4f3a5f
> 
> http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/thread/1397ca25b770b9af/59f63501668f7fbe?lnk=gst&q=replacementfunction#59f63501668f7fbe

Yeah. Okay. What you need is a conversion matrix. Not entirely trivial 
to get those.

This is discussed in a manuscript long ago abandoned, located here:

http://library.wolfram.com/infocenter/MathSource/7523/

I'll provide the needed code below, and skip minor details, like why 
anyone might expect it to do anything useful. If you know, you know, 
else just regard it as you would any other jinn.

(* Groebner basis of a set of vectors wrt POT order *)

moduleGroebnerBasis[polys_, vars_, cvars_, opts___] := Module[
   {newpols, rels, len = Length[cvars], gb, j, k, rul},
   rels = Flatten[
     Table[cvars[[j]]*cvars[[k]], {j, len}, {k, j, len}]];
   newpols = Join[polys, rels];
   gb = GroebnerBasis[newpols, Join[cvars, vars], opts];
   rul = Map[(# :> {}) &, rels];
   gb = Flatten[gb /. rul];
   Collect[gb, cvars]]

(* Finds a Groebner basis wrt lex order, and a conversion matrix. *)

conversionMatrix[polys_, vars_] := Module[
   {aa, coords, pmat, len = Length[polys], newpolys, mgb, gb, convmat,
    fvar, rvars},
   coords = Array[aa, len + 1];
   fvar = First[coords];
   rvars = Rest[coords];
   pmat = Transpose[Join[{polys}, IdentityMatrix[len]]];
   newpolys = pmat.coords;
   mgb = moduleGroebnerBasis[newpolys, vars, coords];
   gb = mgb /. Join[{fvar -> 1}, Thread[rvars -> 0]] /. 0 :> Sequence[];
   convmat = Select[mgb, ! FreeQ[#, fvar] &] /. fvar -> 0;
   {gb, convmat /. Thread[rvars -> Table[UnitVector[len, j], {j, len}]]}
   ]

Here is your example.

terms = (2 kx1^2 kx2^2)/3 - (kx2^2 ky1^2)/3 +
    2 kx1 kx2 ky1 ky2 - (kx1^2 ky2^2)/3 + (2 ky1^2 ky2^2)/
     3 - (kx2^2 kz1^2)/3 - (ky2^2 kz1^2)/3 + 2 kx1 kx2 kz1 kz2 +
    2 ky1 ky2 kz1 kz2 - (kx1^2 kz2^2)/3 - (ky1^2 kz2^2)/
     3 + (2 kz1^2 kz2^2)/3;
k1 = {kx1, ky1, kz1};
k2 = {kx2, ky2, kz2};
polylist = {k1.k1 k2.k2, (k1.k2)^2};
vars = Join[k1, k2];

gb1 = GroebnerBasis[polylist, vars];

In[82]:= {gb, convmat} = conversionMatrix[polylist, vars]

Out[82]= {{kx2^6*ky1^4 + 3*kx2^4*ky1^4*ky2^2 + 3*kx2^2*ky1^4*ky2^4 +
    ky1^4*ky2^6 +
        2*kx2^6*ky1^2*kz1^2 + 4*kx2^4*ky1^2*ky2^2*kz1^2 +
        2*kx2^2*ky1^2*ky2^4*kz1^2 + kx2^6*kz1^4 + kx2^4*ky2^2*kz1^4 +
        4*kx2^4*ky1^3*ky2*kz1*kz2 + 8*kx2^2*ky1^3*ky2^3*kz1*kz2 +
        4*ky1^3*ky2^5*kz1*kz2 + 4*kx2^4*ky1*ky2*kz1^3*kz2 +
        4*kx2^2*ky1*ky2^3*kz1^3*kz2 + kx2^4*ky1^4*kz2^2 +
        2*kx2^2*ky1^4*ky2^2*kz2^2 + ky1^4*ky2^4*kz2^2 +
        4*kx2^4*ky1^2*kz1^2*kz2^2 +
    10*kx2^2*ky1^2*ky2^2*kz1^2*kz2^2 +
        6*ky1^2*ky2^4*kz1^2*kz2^2 + 3*kx2^4*kz1^4*kz2^2 +
        2*kx2^2*ky2^2*kz1^4*kz2^2 + 4*kx2^2*ky1^3*ky2*kz1*kz2^3 +
        4*ky1^3*ky2^3*kz1*kz2^3 + 8*kx2^2*ky1*ky2*kz1^3*kz2^3 +
        4*ky1*ky2^3*kz1^3*kz2^3 + 2*kx2^2*ky1^2*kz1^2*kz2^4 +
        6*ky1^2*ky2^2*kz1^2*kz2^4 + 3*kx2^2*kz1^4*kz2^4 +
    ky2^2*kz1^4*kz2^4 +
        4*ky1*ky2*kz1^3*kz2^5 + kz1^4*kz2^6, kx2^6*ky1^3*ky2 +
        3*kx2^4*ky1^3*ky2^3 + 3*kx2^2*ky1^3*ky2^5 + ky1^3*ky2^7 +
        kx2^6*ky1*ky2*kz1^2 + 2*kx1*kx2^5*ky2^2*kz1^2 +
        4*kx2^4*ky1*ky2^3*kz1^2 + 2*kx1*kx2^3*ky2^4*kz1^2 +
        3*kx2^2*ky1*ky2^5*kz1^2 - kx2^6*ky1^2*kz1*kz2 +
        kx2^4*ky1^2*ky2^2*kz1*kz2 + 5*kx2^2*ky1^2*ky2^4*kz1*kz2 +
        3*ky1^2*ky2^6*kz1*kz2 - kx2^6*kz1^3*kz2 +
    2*kx2^4*ky2^2*kz1^3*kz2 +
        3*kx2^2*ky2^4*kz1^3*kz2 + kx2^4*ky1^3*ky2*kz2^2 +
        2*kx2^2*ky1^3*ky2^3*kz2^2 + ky1^3*ky2^5*kz2^2 +
        2*kx1*kx2^5*kz1^2*kz2^2 + 2*kx2^4*ky1*ky2*kz1^2*kz2^2 +
        4*kx1*kx2^3*ky2^2*kz1^2*kz2^2 + 7*kx2^2*ky1*ky2^3*kz1^2*kz2^2 +
        3*ky1*ky2^5*kz1^2*kz2^2 - kx2^4*ky1^2*kz1*kz2^3 +
        2*kx2^2*ky1^2*ky2^2*kz1*kz2^3 + 3*ky1^2*ky2^4*kz1*kz2^3 +
        5*kx2^2*ky2^2*kz1^3*kz2^3 + ky2^4*kz1^3*kz2^3 +
        2*kx1*kx2^3*kz1^2*kz2^4 + kx2^2*ky1*ky2*kz1^2*kz2^4 +
        3*ky1*ky2^3*kz1^2*kz2^4 + kx2^2*kz1^3*kz2^5 +
    ky2^2*kz1^3*kz2^5,
      (-kx2^4)*ky1^2 + 2*kx1*kx2^3*ky1*ky2 + 2*kx1*kx2*ky1*ky2^3 +
        ky1^2*ky2^4 - kx2^4*kz1^2 - kx2^2*ky2^2*kz1^2 +
    2*kx1*kx2^3*kz1*kz2 +
        2*kx2^2*ky1*ky2*kz1*kz2 + 2*kx1*kx2*ky2^2*kz1*kz2 +
        2*ky1*ky2^3*kz1*kz2 - kx2^2*ky1^2*kz2^2 +
    2*kx1*kx2*ky1*ky2*kz2^2 +
        ky1^2*ky2^2*kz2^2 + ky2^2*kz1^2*kz2^2 + 2*kx1*kx2*kz1*kz2^3 +
        2*ky1*ky2*kz1*kz2^3 + kz1^2*kz2^4, (-kx2^6)*ky1^3 -
        3*kx2^4*ky1^3*ky2^2 - 3*kx2^2*ky1^3*ky2^4 - ky1^3*ky2^6 -
        kx2^6*ky1*kz1^2 - 2*kx1*kx2^5*ky2*kz1^2 -
    4*kx2^4*ky1*ky2^2*kz1^2 -
        2*kx1*kx2^3*ky2^3*kz1^2 - 3*kx2^2*ky1*ky2^4*kz1^2 +
        2*kx1*kx2^5*ky1*kz1*kz2 - 4*kx2^2*ky1^2*ky2^3*kz1*kz2 -
        2*kx1*kx2*ky1*ky2^4*kz1*kz2 - 4*ky1^2*ky2^5*kz1*kz2 -
        2*kx2^4*ky2*kz1^3*kz2 - 2*kx2^2*ky2^3*kz1^3*kz2 -
    kx2^4*ky1^3*kz2^2 -
        2*kx2^2*ky1^3*ky2^2*kz2^2 - ky1^3*ky2^4*kz2^2 -
        4*kx1*kx2^3*ky2*kz1^2*kz2^2 - 7*kx2^2*ky1*ky2^2*kz1^2*kz2^2 -
        2*kx1*kx2*ky2^3*kz1^2*kz2^2 - 5*ky1*ky2^4*kz1^2*kz2^2 +
        2*kx1*kx2^3*ky1*kz1*kz2^3 - 2*kx1*kx2*ky1*ky2^2*kz1*kz2^3 -
        4*ky1^2*ky2^3*kz1*kz2^3 - 4*kx2^2*ky2*kz1^3*kz2^3 -
        2*ky2^3*kz1^3*kz2^3 + kx2^2*ky1*kz1^2*kz2^4 -
        2*kx1*kx2*ky2*kz1^2*kz2^4 - 5*ky1*ky2^2*kz1^2*kz2^4 -
    2*ky2*kz1^3*kz2^5,
      (-kx2^5)*ky1^3 - 4*kx2^3*ky1^3*ky2^2 +
    2*kx1*kx2^2*ky1^2*ky2^3 -
        3*kx2*ky1^3*ky2^4 + 2*kx1*ky1^2*ky2^5 - kx2^5*ky1*kz1^2 -
        2*kx1*kx2^4*ky2*kz1^2 - 5*kx2^3*ky1*ky2^2*kz1^2 -
        2*kx1*kx2^2*ky2^3*kz1^2 - 4*kx2*ky1*ky2^4*kz1^2 +
        2*kx1*kx2^4*ky1*kz1*kz2 - 2*kx2^3*ky1^2*ky2*kz1*kz2 +
        6*kx1*kx2^2*ky1*ky2^2*kz1*kz2 - 2*kx2*ky1^2*ky2^3*kz1*kz2 +
        4*kx1*ky1*ky2^4*kz1*kz2 - 4*kx2^3*ky2*kz1^3*kz2 -
        4*kx2*ky2^3*kz1^3*kz2 - kx2^3*ky1^3*kz2^2 -
    3*kx2*ky1^3*ky2^2*kz2^2 +
        2*kx1*ky1^2*ky2^3*kz2^2 - 3*kx2*ky1*ky2^2*kz1^2*kz2^2 +
        2*kx1*ky2^3*kz1^2*kz2^2 + 2*kx1*kx2^2*ky1*kz1*kz2^3 -
        2*kx2*ky1^2*ky2*kz1*kz2^3 + 4*kx1*ky1*ky2^2*kz1*kz2^3 -
        4*kx2*ky2*kz1^3*kz2^3 + kx2*ky1*kz1^2*kz2^4 +
    2*kx1*ky2*kz1^2*kz2^4,
      kx1*kx2^4*ky1^2 + 2*kx2^3*ky1^3*ky2 + 2*kx2*ky1^3*ky2^3 -
        kx1*ky1^2*ky2^4 + kx1*kx2^4*kz1^2 + 2*kx2^3*ky1*ky2*kz1^2 +
        kx1*kx2^2*ky2^2*kz1^2 + 2*kx2*ky1*ky2^3*kz1^2 +
    2*kx2^3*ky1^2*kz1*kz2 -
        2*kx1*kx2^2*ky1*ky2*kz1*kz2 + 2*kx2*ky1^2*ky2^2*kz1*kz2 -
        2*kx1*ky1*ky2^3*kz1*kz2 + 2*kx2^3*kz1^3*kz2 +
    2*kx2*ky2^2*kz1^3*kz2 +
        kx1*kx2^2*ky1^2*kz2^2 + 2*kx2*ky1^3*ky2*kz2^2 -
    kx1*ky1^2*ky2^2*kz2^2 +
        2*kx2*ky1*ky2*kz1^2*kz2^2 - kx1*ky2^2*kz1^2*kz2^2 +
        2*kx2*ky1^2*kz1*kz2^3 - 2*kx1*ky1*ky2*kz1*kz2^3 +
    2*kx2*kz1^3*kz2^3 -
        kx1*kz1^2*kz2^4,
   kx2^2*ky1^2 - 2*kx1*kx2*ky1*ky2 + kx1^2*ky2^2 +
        kx2^2*kz1^2 + ky2^2*kz1^2 - 2*kx1*kx2*kz1*kz2 -
    2*ky1*ky2*kz1*kz2 +
        kx1^2*kz2^2 + ky1^2*kz2^2, kx1^2*kx2^2 + 2*kx1*kx2*ky1*ky2 +
        ky1^2*ky2^2 + 2*kx1*kx2*kz1*kz2 + 2*ky1*ky2*kz1*kz2 +
    kz1^2*kz2^2},
    {{kx2^4*ky1^2 + 2*kx1*kx2^3*ky1*ky2 + 3*kx2^2*ky1^2*ky2^2 +
     kx2^4*kz1^2 +
          2*kx1*kx2^3*kz1*kz2 + 6*kx2^2*ky1*ky2*kz1*kz2 +
     3*kx2^2*kz1^2*kz2^2,
        (-kx2^4)*ky1^2 - 2*kx1*kx2^3*ky1*ky2 - 2*kx1*kx2*ky1*ky2^3 +
          ky1^2*ky2^4 - kx2^4*kz1^2 - kx2^2*ky2^2*kz1^2 -
     2*kx1*kx2^3*kz1*kz2 +
          2*kx2^2*ky1*ky2*kz1*kz2 - 2*kx1*kx2*ky2^2*kz1*kz2 +
          2*ky1*ky2^3*kz1*kz2 - kx2^2*ky1^2*kz2^2 -
     2*kx1*kx2*ky1*ky2*kz2^2 +
          ky1^2*ky2^2*kz2^2 + ky2^2*kz1^2*kz2^2 - 2*kx1*kx2*kz1*kz2^3 +
          2*ky1*ky2*kz1*kz2^3 + kz1^2*kz2^4},
      {kx2^4*ky1*ky2 + 2*kx1*kx2^3*ky2^2 + 3*kx2^2*ky1*ky2^3 -
     kx2^4*kz1*kz2 +
          3*kx2^2*ky2^2*kz1*kz2, (-kx2^4)*ky1*ky2 - 2*kx1*kx2^3*ky2^2 -
          2*kx1*kx2*ky2^4 + ky1*ky2^5 + kx2^4*kz1*kz2 +
     2*kx2^2*ky2^2*kz1*kz2 +
          ky2^4*kz1*kz2 - kx2^2*ky1*ky2*kz2^2 - 2*kx1*kx2*ky2^2*kz2^2 +
          ky1*ky2^3*kz2^2 + kx2^2*kz1*kz2^3 + ky2^2*kz1*kz2^3},
      {-kx2^2,
    kx2^2 + ky2^2 + kz2^2}, {(-kx2^4)*ky1 - 2*kx1*kx2^3*ky2 -
          3*kx2^2*ky1*ky2^2 - 2*kx2^2*ky2*kz1*kz2,
    kx2^4*ky1 + 2*kx1*kx2^3*ky2 +
          2*kx1*kx2*ky2^3 - ky1*ky2^4 - 2*kx2^2*ky2*kz1*kz2 -
     2*ky2^3*kz1*kz2 +
          kx2^2*ky1*kz2^2 + 2*kx1*kx2*ky2*kz2^2 - ky1*ky2^2*kz2^2 -
          2*ky2*kz1*kz2^3}, {(-kx2^3)*ky1 - 2*kx1*kx2^2*ky2 -
     4*kx2*ky1*ky2^2 -
          4*kx2*ky2*kz1*kz2,
    kx2^3*ky1 + 2*kx1*kx2^2*ky2 + kx2*ky1*ky2^2 +
          2*kx1*ky2^3 + kx2*ky1*kz2^2 + 2*kx1*ky2*kz2^2},
      {kx1*kx2^2 + 2*kx2*ky1*ky2 + 2*kx2*kz1*kz2, (-kx1)*kx2^2 -
     kx1*ky2^2 -
          kx1*kz2^2}, {1, -1}, {0, 1}}}

Quick check:

In[84]:= gb === gb1
Out[84]= True

As you did, we'll reduce terms with respect to the Groebner basis of 
polylist. We show the vector of reducing multipliers.

In[85]:= reduction = First[PolynomialReduce[terms, gb, vars]]
Out[85]= {0, 0, 0, 0, 0, 0, -(1/3), 2/3}

To get the reduction in terms of the original polylist, we compute 
reduction.mgb.

In[76]:= origreduction = Expand[reduction.convmat]
Out[76]= {-(1/3), 1}

Let's check this.

In[77]:= Expand[origreduction.polylist - terms]
Out[77]= 0

So we have now successfully rewritten 'terms' in terms of polylist.

I will remark that for production code purposes, one might wish to alter 
the above to use a faster term order. If that is a compelling need then 
I might be able to dredge up some old code for the occasion.


Daniel Lichtblau
Wolfram Research



  • Prev by Date: Re: Algebraic substitution with PolynomialReduce
  • Next by Date: Re: Execute on notebook creation
  • Previous by thread: Re: Algebraic substitution with PolynomialReduce
  • Next by thread: Re: Apply a rule to an expression only once