Skip to main content

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

Popular posts from this blog

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...