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: [mg117133] Re: Algebraic substitution with PolynomialReduce
  • From: Guido Walter Pettinari <coccoinomane at gmail.com>
  • Date: Thu, 10 Mar 2011 06:11:53 -0500 (EST)

Thank you Daniel!  This is exactly what I needed.

In the meanwhile I found a temporary solution, that is running PolynomialReduce on all permutations of 'vars' until a zero-remainder result was found, but yours is much faster and cleaner.

I am going to apply your solution to much tougher expression, and I'll let you know if I succeed.

Thank you again.

Cheers,

Guido


On 9 Mar 2011, at 16:26, Daniel Lichtblau wrote:

> 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_thr=
ead/thread/634a54e2e844837f/ddf3803cab4f3a5f?lnk==gst&q==replacementfunctio=
n#ddf3803cab4f3a5f
>> http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thr=
ead/thread/1397ca25b770b9af/59f63501668f7fbe?lnk==gst&q==replacementfunctio=
n#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 anyo=
ne 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 poly=
list. 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 reduct=
ion.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 m=
ight be able to dredge up some old code for the occasion.
>
>
> Daniel Lichtblau
> Wolfram Research
>


  • Prev by Date: Re: FindFit power law problem
  • Next by Date: Re: Algebraic substitution with PolynomialReduce
  • Previous by thread: Algebraic substitution with PolynomialReduce
  • Next by thread: Re: Algebraic substitution with PolynomialReduce