3

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]]))
*)
Community
  • 1
jjstankowicz
  • 687
  • 3
  • 14

1 Answers1

2

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]}}
J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
  • 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. – jjstankowicz Jun 15 '16 at 22:30
  • 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. – J. M.'s missing motivation Jun 15 '16 at 23:08
  • 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? – jjstankowicz Jun 16 '16 at 05:12
  • 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. – J. M.'s missing motivation Jun 16 '16 at 05:15
  • 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. – jjstankowicz Jun 16 '16 at 05:17
  • I may have misunderstood your question. The new factor x[a[1],a[6]] does appear in det[a]. My original augmented example with x[a[1],a[2]] does not. Either way the Groebner basis method doesn't work. A best case scenario would work even for x[_,_] that are in none of the det[_]s or any other functional dependence which does not appear in any of the det[_]s. – jjstankowicz Jun 16 '16 at 18:03
  • Just for reference, can you edit your question to include the example where the Gröbner approach goes bust? – J. M.'s missing motivation Jun 16 '16 at 18:08