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

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]

plotting - Filling between two spheres in SphericalPlot3D

Manipulate[ SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, Mesh -> None, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], {n, 0, 1}] I cant' seem to be able to make a filling between two spheres. I've already tried the obvious Filling -> {1 -> {2}} but Mathematica doesn't seem to like that option. Is there any easy way around this or ... Answer There is no built-in filling in SphericalPlot3D . One option is to use ParametricPlot3D to draw the surfaces between the two shells: Manipulate[ Show[SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], ParametricPlot3D[{ r {Sin[t] Cos[1.5 Pi], Sin[t] Sin[1.5 Pi], Cos[t]}, r {Sin[t] Cos[0 Pi], Sin[t] Sin[0 Pi], Cos[t]}}, {r, 1, 2 - n}, {t, 0, Pi}, PlotStyle -> Yellow, Mesh -> {2, 15}]], {n, 0, 1}]

plotting - Mathematica: 3D plot based on combined 2D graphs

I have several sigmoidal fits to 3 different datasets, with mean fit predictions plus the 95% confidence limits (not symmetrical around the mean) and the actual data. I would now like to show these different 2D plots projected in 3D as in but then using proper perspective. In the link here they give some solutions to combine the plots using isometric perspective, but I would like to use proper 3 point perspective. Any thoughts? Also any way to show the mean points per time point for each series plus or minus the standard error on the mean would be cool too, either using points+vertical bars, or using spheres plus tubes. Below are some test data and the fit function I am using. Note that I am working on a logit(proportion) scale and that the final vertical scale is Log10(percentage). (* some test data *) data = Table[Null, {i, 4}]; data[[1]] = {{1, -5.8}, {2, -5.4}, {3, -0.8}, {4, -0.2}, {5, 4.6}, {1, -6.4}, {2, -5.6}, {3, -0.7}, {4, 0.04}, {5, 1.0}, {1, -6.8}, {2, -4.7}, {3, -1....