algebraic manipulation - Replacing patterns (determinants) in a rational multivariate function and/or getting `replacementFunction` to work on a multi-linear expression
Given a rational multivariate function known to contain determinants, what is the best way to rewrite the function in terms of the determinants?
A simple example:
det[aa_] := Det[Table[x[li1, li2], {li1, aa /@ {1, 2}}, {li2, aa /@ {3, 4}}]];
origExpr = det[a]*det[b]/det[c] + det[d]*det[e]/det[f] - det[g]*det[h]/det[i];
oeet = origExpr //Expand // Together;
(* I paste the long, crappy expression at the end of this quesion *)
That the entries of the matrix be x[i_,j_]
isn't necessarily important, that just happens to be the problem I'm considering -- the determinants happen to be Gram determinants.
Question: Given the output of oeet
, what is the best way to rediscover origExpr
?
Edit 2: See "Edit 2" near the bottom of this post for an updated discussion of this problem.
One solution: This works for simple-enough examples:
vars = Union[Cases[oeet, x[a_[i_], _] :> a, Infinity]];(* Determine the variables *)
rep = Table[d[ti] == det[ti], {ti, vars}];
Simplify[oeet, rep]
(* (d[a] d[b])/d[c] + (d[d] d[e])/d[f] - (d[g] d[h])/d[i] *)
This doesn't do so well for $3\times 3$ determinants though:
det[aa_] := Det[Table[x[li1, li2], {li1, aa /@ {1, 2, 3}}, {li2, aa /@ {4, 5, 6}}]];
The analogous Simplify[oeet, rep]
runs for over a minute before I kill it. The expressions that I'm actually manipulating are rational in the x[_,_]
's, and are not completely rational in in the determinants, and I think that adds complexity. (I couldn't come up with a non-trivial alteration to the minimal working example here.)
Additional attempts: I also tried to use the replacementFunction
for example here (from @DanielLichtblau, I think?), but it appears not to even alter oeet
:
aVars =Select[Union[Cases[oeet, x[__], Infinity]],! FreeQ[#, a, Infinity] &]
(* {x[a[1], a[3]], x[a[1], a[4]], x[a[2], a[3]], x[a[2], a[4]]} *)
replacementFunction[oeet, detA - det[a], aVars] - oeet
(* 0 (* Should be detA-depdendent *) *)
Assuming I'm somewhow misusing the function, I tried PolynomialReduce
which does a fairly good job on the $2\times 2$ matrices, but misses the denominators.
xVars = Union[Cases[oeet, x[_, _], Infinity]]
redEqns = (safe[#] - det[#]) & /@ vars
PolynomialReduce[oeet, redEqns, xVars][[2]]
(*
(safe[c] safe[f] safe[g] safe[h])/((x[c[1], c[4]] x[c[2], c[3]] -
x[c[1], c[3]] x[c[2], c[4]]) (x[f[1], f[4]] x[f[2], f[3]] -
x[f[1], f[3]] x[f[2], f[4]]) (x[i[1], i[4]] x[i[2], i[3]] -
x[i[1], i[3]] x[i[2], i[4]])) - (safe[c] safe[d] safe[e]
safe[i])/((x[c[1], c[4]] x[c[2], c[3]] -
x[c[1], c[3]] x[c[2], c[4]]) (x[f[1], f[4]] x[f[2], f[3]] -
x[f[1], f[3]] x[f[2], f[4]]) (x[i[1], i[4]] x[i[2], i[3]] -
x[i[1], i[3]] x[i[2], i[4]])) - (safe[a] safe[b] safe[f]
safe[i])/((x[c[1], c[4]] x[c[2], c[3]] -
x[c[1], c[3]] x[c[2], c[4]]) (x[f[1], f[4]] x[f[2], f[3]] -
x[f[1], f[3]] x[f[2], f[4]]) (x[i[1], i[4]] x[i[2], i[3]] -
x[i[1], i[3]] x[i[2], i[4]]))
*)
It again takes a while on $3\times 3$ matrices.
Edit: Here is a simple example where the Groebner basis approach suggested in the first solution doesn't work:
det[aa_] :=
Det[Table[
x[li1, li2], {li1, aa /@ {1, 2, 3}}, {li2, aa /@ {4, 5, 6}}]];
origExpr = x[a[1], x[a[6]]] det[a] det[b]/det[c]; (* <- Alteration here *)
oeet = origExpr // Expand // Together;
varsx = Cases[oeet, x[a_[i_], _], \[Infinity]];
vars = Union[Cases[oeet, x[a_[i_], _] :> a, \[Infinity]]];
rep = Table[d[ti] == det[ti], {ti, vars}];
GroebnerBasis[{w == oeet}~Join~rep, w, varsx]
(* {} *)
It does work if I know to remove x[a[1],a[6]]
from the "Elimination" argument of GroebnerBasis
:
GroebnerBasis[{w == oeet}~Join~rep, w,DeleteCases[varsx, x[a[1], a[6]]]]
(* {w d[c] - d[a] d[b] x[a[1], a[6]]} *)
but given a general polynomial in the x[_,_]
's, I won't know which x[_,_]
's to remove from the elimination argument of GroebnerBasis
. I guess in general I'm looking for something that can implement polynomial changes of variables in a multivariate rational function.
Edit 2: I'm now pretty sure the best way to do this is with the replacementFunction
I already referenced. The function is effectively an iterative PolynomialReduce
. It doesn't seem to like exactly multi-linear functions, such as determinants, however. Here's a version of replacementFunction
that has an optional fourth argument that when set to ON
prints output to track what's happening. This link has the barebones version of the function.
replacementFunction // ClearAll;
replacementFunction[expr_, rep_, vars_, TS_: 0] :=
Module[
{num = Numerator[expr], den = Denominator[expr], hed = Head[expr],
base, expon, out, tsp}
,
tsp[x_] := If[TS === ON, Print[x];];
If[
PolynomialQ[num, vars] && PolynomialQ[den, vars] && ! NumberQ[den]
,
tsp["T1 - A rational function"];
tsp[expr];
replacementFunction[num, rep, vars, TS]/replacementFunction[den, rep, vars, TS]
,
tsp["F1 - Not a rational function"];
tsp[expr];
If[
hed === Power && Length[expr] == 2
,
tsp["T2 - A power function"];
tsp[expr];
base = replacementFunction[expr[[1]], rep, vars, TS];
expon = replacementFunction[expr[[2]], rep, vars, TS];
out = PolynomialReduce[base^expon, rep, vars];
tsp["===T2out==="];
tsp[out // Flatten // TableForm];
out[[2]]
,
tsp["F2 - Not a power function"];
tsp[expr];
If[
Head[hed] === Symbol && MemberQ[Attributes[Evaluate[hed]], NumericFunction]
,
tsp["T3 - A numeric function"];
tsp[expr];
Map[replacementFunction[#, rep, vars, TS] &, expr]
,
tsp["F3 - Not a numeric function"];
tsp["***Reduce***"];
tsp["Divide ", expr];
tsp["by ", rep];
out = PolynomialReduce[expr, rep, vars];
tsp["===out==="];
tsp[out // Flatten // TableForm];
out[[2]]
], TS
]
]
];
This doesn't work on a simple determinant:
expr = (a[1, 1] a[2, 2] - a[1, 2] a[2, 1]);
replacementFunction[expr, expr + d, Variables[expr]]
(* Actually running will print out what the function is doing *)
(* -a[1, 2] a[2, 1] + a[1, 1] a[2, 2] (*Output*)*)
It appears to work when the expression is not multi-linear:
exprA = (a[1, 1] a[2, 2] - a[1, 2] a[2, 1]);
exprB = (b[1, 1] b[2, 2] - b[1, 2] b[2, 1]);
replacementFunction[(exprA + exprB)^2, {exprA - dA, exprB - dB}, Join[Variables[exprA],Variables[exprB]]]
(* dA^2 + 2 dA dB + dB^2 *)
But again breaks when the expression is multi-linear (like a determinant is):
replacementFunction[exprA exprB, {exprA - dA, exprB - dB},Join[Variables[exprA],Variables[exprB]]]
(* (-a[1, 2] a[2, 1] + a[1, 1] a[2, 2]) (-b[1, 2] b[2, 1] + b[1, 1] b[2, 2]) *)
Updated question: How do I get replacementFunction
to work on multi-linear expressions such as determinants.
It appears to already solve the more general problem of changing variables in a polynomial. From what I can tell, replacementFunction
jumps through every part at every level of the expression and performs PolynomialReduce
where it can. I can't currently see why it doesn't catch the multi-linear terms.
Below is the crappy output of oeet
.
(* The output from oeet = origExpr //Expand // Together:
(x[c[1], c[4]] x[c[2], c[3]] x[f[1], f[4]] x[f[2], f[3]] x[g[1],
g[4]] x[g[2], g[3]] x[h[1], h[4]] x[h[2], h[3]] -
x[c[1], c[3]] x[c[2], c[4]] x[f[1], f[4]] x[f[2], f[3]] x[g[1],
g[4]] x[g[2], g[3]] x[h[1], h[4]] x[h[2], h[3]] -
x[c[1], c[4]] x[c[2], c[3]] x[f[1], f[3]] x[f[2], f[4]] x[g[1],
g[4]] x[g[2], g[3]] x[h[1], h[4]] x[h[2], h[3]] +
x[c[1], c[3]] x[c[2], c[4]] x[f[1], f[3]] x[f[2], f[4]] x[g[1],
g[4]] x[g[2], g[3]] x[h[1], h[4]] x[h[2], h[3]] -
x[c[1], c[4]] x[c[2], c[3]] x[f[1], f[4]] x[f[2], f[3]] x[g[1],
g[3]] x[g[2], g[4]] x[h[1], h[4]] x[h[2], h[3]] +
x[c[1], c[3]] x[c[2], c[4]] x[f[1], f[4]] x[f[2], f[3]] x[g[1],
g[3]] x[g[2], g[4]] x[h[1], h[4]] x[h[2], h[3]] +
x[c[1], c[4]] x[c[2], c[3]] x[f[1], f[3]] x[f[2], f[4]] x[g[1],
g[3]] x[g[2], g[4]] x[h[1], h[4]] x[h[2], h[3]] -
x[c[1], c[3]] x[c[2], c[4]] x[f[1], f[3]] x[f[2], f[4]] x[g[1],
g[3]] x[g[2], g[4]] x[h[1], h[4]] x[h[2], h[3]] -
x[c[1], c[4]] x[c[2], c[3]] x[f[1], f[4]] x[f[2], f[3]] x[g[1],
g[4]] x[g[2], g[3]] x[h[1], h[3]] x[h[2], h[4]] +
x[c[1], c[3]] x[c[2], c[4]] x[f[1], f[4]] x[f[2], f[3]] x[g[1],
g[4]] x[g[2], g[3]] x[h[1], h[3]] x[h[2], h[4]] +
x[c[1], c[4]] x[c[2], c[3]] x[f[1], f[3]] x[f[2], f[4]] x[g[1],
g[4]] x[g[2], g[3]] x[h[1], h[3]] x[h[2], h[4]] -
x[c[1], c[3]] x[c[2], c[4]] x[f[1], f[3]] x[f[2], f[4]] x[g[1],
g[4]] x[g[2], g[3]] x[h[1], h[3]] x[h[2], h[4]] +
x[c[1], c[4]] x[c[2], c[3]] x[f[1], f[4]] x[f[2], f[3]] x[g[1],
g[3]] x[g[2], g[4]] x[h[1], h[3]] x[h[2], h[4]] -
x[c[1], c[3]] x[c[2], c[4]] x[f[1], f[4]] x[f[2], f[3]] x[g[1],
g[3]] x[g[2], g[4]] x[h[1], h[3]] x[h[2], h[4]] -
x[c[1], c[4]] x[c[2], c[3]] x[f[1], f[3]] x[f[2], f[4]] x[g[1],
g[3]] x[g[2], g[4]] x[h[1], h[3]] x[h[2], h[4]] +
x[c[1], c[3]] x[c[2], c[4]] x[f[1], f[3]] x[f[2], f[4]] x[g[1],
g[3]] x[g[2], g[4]] x[h[1], h[3]] x[h[2], h[4]] -
x[c[1], c[4]] x[c[2], c[3]] x[d[1], d[4]] x[d[2], d[3]] x[e[1],
e[4]] x[e[2], e[3]] x[i[1], i[4]] x[i[2], i[3]] +
x[c[1], c[3]] x[c[2], c[4]] x[d[1], d[4]] x[d[2], d[3]] x[e[1],
e[4]] x[e[2], e[3]] x[i[1], i[4]] x[i[2], i[3]] +
x[c[1], c[4]] x[c[2], c[3]] x[d[1], d[3]] x[d[2], d[4]] x[e[1],
e[4]] x[e[2], e[3]] x[i[1], i[4]] x[i[2], i[3]] -
x[c[1], c[3]] x[c[2], c[4]] x[d[1], d[3]] x[d[2], d[4]] x[e[1],
e[4]] x[e[2], e[3]] x[i[1], i[4]] x[i[2], i[3]] +
x[c[1], c[4]] x[c[2], c[3]] x[d[1], d[4]] x[d[2], d[3]] x[e[1],
e[3]] x[e[2], e[4]] x[i[1], i[4]] x[i[2], i[3]] -
x[c[1], c[3]] x[c[2], c[4]] x[d[1], d[4]] x[d[2], d[3]] x[e[1],
e[3]] x[e[2], e[4]] x[i[1], i[4]] x[i[2], i[3]] -
x[c[1], c[4]] x[c[2], c[3]] x[d[1], d[3]] x[d[2], d[4]] x[e[1],
e[3]] x[e[2], e[4]] x[i[1], i[4]] x[i[2], i[3]] +
x[c[1], c[3]] x[c[2], c[4]] x[d[1], d[3]] x[d[2], d[4]] x[e[1],
e[3]] x[e[2], e[4]] x[i[1], i[4]] x[i[2], i[3]] -
x[a[1], a[4]] x[a[2], a[3]] x[b[1], b[4]] x[b[2], b[3]] x[f[1],
f[4]] x[f[2], f[3]] x[i[1], i[4]] x[i[2], i[3]] +
x[a[1], a[3]] x[a[2], a[4]] x[b[1], b[4]] x[b[2], b[3]] x[f[1],
f[4]] x[f[2], f[3]] x[i[1], i[4]] x[i[2], i[3]] +
x[a[1], a[4]] x[a[2], a[3]] x[b[1], b[3]] x[b[2], b[4]] x[f[1],
f[4]] x[f[2], f[3]] x[i[1], i[4]] x[i[2], i[3]] -
x[a[1], a[3]] x[a[2], a[4]] x[b[1], b[3]] x[b[2], b[4]] x[f[1],
f[4]] x[f[2], f[3]] x[i[1], i[4]] x[i[2], i[3]] +
x[a[1], a[4]] x[a[2], a[3]] x[b[1], b[4]] x[b[2], b[3]] x[f[1],
f[3]] x[f[2], f[4]] x[i[1], i[4]] x[i[2], i[3]] -
x[a[1], a[3]] x[a[2], a[4]] x[b[1], b[4]] x[b[2], b[3]] x[f[1],
f[3]] x[f[2], f[4]] x[i[1], i[4]] x[i[2], i[3]] -
x[a[1], a[4]] x[a[2], a[3]] x[b[1], b[3]] x[b[2], b[4]] x[f[1],
f[3]] x[f[2], f[4]] x[i[1], i[4]] x[i[2], i[3]] +
x[a[1], a[3]] x[a[2], a[4]] x[b[1], b[3]] x[b[2], b[4]] x[f[1],
f[3]] x[f[2], f[4]] x[i[1], i[4]] x[i[2], i[3]] +
x[c[1], c[4]] x[c[2], c[3]] x[d[1], d[4]] x[d[2], d[3]] x[e[1],
e[4]] x[e[2], e[3]] x[i[1], i[3]] x[i[2], i[4]] -
x[c[1], c[3]] x[c[2], c[4]] x[d[1], d[4]] x[d[2], d[3]] x[e[1],
e[4]] x[e[2], e[3]] x[i[1], i[3]] x[i[2], i[4]] -
x[c[1], c[4]] x[c[2], c[3]] x[d[1], d[3]] x[d[2], d[4]] x[e[1],
e[4]] x[e[2], e[3]] x[i[1], i[3]] x[i[2], i[4]] +
x[c[1], c[3]] x[c[2], c[4]] x[d[1], d[3]] x[d[2], d[4]] x[e[1],
e[4]] x[e[2], e[3]] x[i[1], i[3]] x[i[2], i[4]] -
x[c[1], c[4]] x[c[2], c[3]] x[d[1], d[4]] x[d[2], d[3]] x[e[1],
e[3]] x[e[2], e[4]] x[i[1], i[3]] x[i[2], i[4]] +
x[c[1], c[3]] x[c[2], c[4]] x[d[1], d[4]] x[d[2], d[3]] x[e[1],
e[3]] x[e[2], e[4]] x[i[1], i[3]] x[i[2], i[4]] +
x[c[1], c[4]] x[c[2], c[3]] x[d[1], d[3]] x[d[2], d[4]] x[e[1],
e[3]] x[e[2], e[4]] x[i[1], i[3]] x[i[2], i[4]] -
x[c[1], c[3]] x[c[2], c[4]] x[d[1], d[3]] x[d[2], d[4]] x[e[1],
e[3]] x[e[2], e[4]] x[i[1], i[3]] x[i[2], i[4]] +
x[a[1], a[4]] x[a[2], a[3]] x[b[1], b[4]] x[b[2], b[3]] x[f[1],
f[4]] x[f[2], f[3]] x[i[1], i[3]] x[i[2], i[4]] -
x[a[1], a[3]] x[a[2], a[4]] x[b[1], b[4]] x[b[2], b[3]] x[f[1],
f[4]] x[f[2], f[3]] x[i[1], i[3]] x[i[2], i[4]] -
x[a[1], a[4]] x[a[2], a[3]] x[b[1], b[3]] x[b[2], b[4]] x[f[1],
f[4]] x[f[2], f[3]] x[i[1], i[3]] x[i[2], i[4]] +
x[a[1], a[3]] x[a[2], a[4]] x[b[1], b[3]] x[b[2], b[4]] x[f[1],
f[4]] x[f[2], f[3]] x[i[1], i[3]] x[i[2], i[4]] -
x[a[1], a[4]] x[a[2], a[3]] x[b[1], b[4]] x[b[2], b[3]] x[f[1],
f[3]] x[f[2], f[4]] x[i[1], i[3]] x[i[2], i[4]] +
x[a[1], a[3]] x[a[2], a[4]] x[b[1], b[4]] x[b[2], b[3]] x[f[1],
f[3]] x[f[2], f[4]] x[i[1], i[3]] x[i[2], i[4]] +
x[a[1], a[4]] x[a[2], a[3]] x[b[1], b[3]] x[b[2], b[4]] x[f[1],
f[3]] x[f[2], f[4]] x[i[1], i[3]] x[i[2], i[4]] -
x[a[1], a[3]] x[a[2], a[4]] x[b[1], b[3]] x[b[2], b[4]] x[f[1],
f[3]] x[f[2], f[4]] x[i[1], i[3]] x[i[2],
i[4]])/((x[c[1], c[4]] x[c[2], c[3]] -
x[c[1], c[3]] x[c[2], c[4]]) (x[f[1], f[4]] x[f[2], f[3]] -
x[f[1], f[3]] x[f[2], f[4]]) (x[i[1], i[4]] x[i[2], i[3]] -
x[i[1], i[3]] x[i[2], i[4]]))
*)
Comments
Post a Comment