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
>