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