3
$\begingroup$

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]]))
*)
Improve this question
$\endgroup$
1
  • $\begingroup$ This appears to be a problem that periodically resurfaces. It looks like PolynomialReduce tends to be the answer... $\endgroup$
    – jjstankowicz
    Commented Jun 17, 2016 at 20:36

1 Answer 1

Reset to default
2
$\begingroup$

As usual, GroebnerBasis[] is very useful here:

det[aa_] := Det[Table[x[li1, li2], {li1, aa /@ {1, 2, 3}}, {li2, aa /@ {4, 5, 6}}]];

origExpr = det[a] det[b]/det[c] + det[d] det[e]/det[f] - det[g] det[h]/det[i];
oeet = origExpr // Expand // Together;

varsx = Cases[oeet, x[a_[i_], _], ∞];
vars = Union[Cases[oeet, x[a_[i_], _] :> a, ∞]];
rep = Table[d[ti] == det[ti], {ti, vars}];

Solve[First[GroebnerBasis[{w == oeet} ~Join~ rep, w, varsx]] == 0, w] // Simplify
   {{w -> (d[a] d[b])/d[c] + (d[d] d[e])/d[f] - (d[g] d[h])/d[i]}}
Improve this answer
$\endgroup$
7
  • $\begingroup$ This is a good first step, but there's a bit of a hiccup. Taking origExpr = x[a[1], x[a[2]]] det[a]*det[b]/det[c] + det[d]*det[e]/det[f] - det[g]*det[h]/det[i] requires I know in advance to exclude x[a[1], x[a[2]]] from varsx. Given oeet, I may not know a priori which x[_,_]'s can/can not be eliminated. $\endgroup$
    – jjstankowicz
    Commented Jun 15, 2016 at 22:30
  • $\begingroup$ Well, that's why the code just picks out all the x's (as you did) and feeds the entire list into GroebnerBasis[], leaving it to figure out how to perform the required eliminations. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Jun 15, 2016 at 23:08
  • $\begingroup$ What I meant to say is that the solution you propose: GroebnerBasis[{w == oeet}~Join~rep, w, varsx] returns {} when the original expression has simple additional x[_,_] dependence: origExpr = x[a[1], x[a[6]]] det[a]*det[b]/det[c] + det[d]*det[e]/det[f] - det[g]*det[h]/det[i]. If I happen know to remove x[a[1], x[a[6]]] from varsx by hand, then the result is as desired. But for a generic expression, I don't know which x[_,_]'s to exclude from the third argument of GroebnerBasis. Is there some way to address this issue as well? $\endgroup$
    – jjstankowicz
    Commented Jun 16, 2016 at 5:12
  • $\begingroup$ So, your expressions in general may have x objects that were not present in the det objects? If so, you really should have mentioned that to begin with in your question. $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Jun 16, 2016 at 5:15
  • $\begingroup$ Yes. And it actually doesn't affect the Simplify[oeet, rep] I originally included for the $2\times 2$ matrices, so I assumed it wouldn't affect other solutions either. Maybe a bad assumption. $\endgroup$
    – jjstankowicz
    Commented Jun 16, 2016 at 5:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.