Arrangements and Likelihood
On this page, we provide Macaulay2 code to compute the main objects of the paper, and explain how to recreate all computations mentioned in examples throughout the text.
In order to run our code you should download the file ArrangementsLikelihood.m2 and save it in the directory from which you are running Macaulay2. Load the functions provided by typing load "ArrangementsLikelihood.m2".
We provide the following functions:
twoMatricesreturns a tuple \((A,B)\) of two matrices defining the kernel of the matrix \(Q\) as defined on page 2 in the paperteraoModulereturns the Terao module of the arrangementjacobianModulereturns the Jacobian syzygy module \(J(\mathcal{A})\) of the arrangementlogDerModulereturns the log-derivation module \(D(\mathcal{A})\) of the arrangementlikelihoodModulereturns the likelihood module \(M(\mathcal{A})\) of the arrangementpreLikelihoodIdealreturns the pre-likelihood ideal \(I_0(\mathcal{A})\) of the arrangementlikelihoodIdealreturns the likelihood ideal \(I(\mathcal{A})\) of the arrangement. This function might be slow as all potentially necessary saturations are computed. It might be faster to compute the saturations step-by-step and check for primality after each step, as described in Remark 2.11.isTamechecks whether an arrangement is tame. For arrangements with many variables this function is likely to be very slow.
Each of these functions takes as input a list of polynomials from a polynomial ring with no superfluous variables (i.e. every variable should appear in at least one of the polynomials provided in the list).
In the following we illustrate how to use our functions to recreate the examples from the paper.
Example 2.1 (Braid arrangement)
R = QQ[x,y,z,w]
A = {x-y, x-z, x-w, y-z, y-w, z-w}
twoMatrices A
We can check for gentleness of \(\mathcal{A}\) as follows:
I0 = preLikelihoodIdeal A
isPrime A
This returns true so \(\mathcal{A}\) is gentle. Freeness can be checked as follows:
D = logDerModule A
isFreeModule prune D
Moreover, we can check that \(\mathcal{A}\) is tame via isTame A.
Example 2.8 (3 lines and a cubic in the plane)
R = QQ[x,y,z]
A = {x, y, z, x^3 + y^3 + x*y*z}
I0 = preLikelihoodIdeal A
primaryDecomposition I0
To obtain the likelihood ideal we can saturate at any element in the ideal \(\langle x,y\rangle\), for example at \(x\):
I = saturate(I0, x)
isPrime I
This gives a prime ideal, therefore being the likelihood ideal, see Proposition 2.9.
Example 3.4 (Coin model)
The coin model is defined by the following parametrization which we first need to homogenize:
R = QQ[t1,t2,t3]
p1 = t1*(1-t2)^4 + (1-t1)*(1-t3)^4
p2 = 4*t1*t2*(1-t2)^3 + 4*(1-t1)*t3*(1-t3)^3
p3 = 6*t1*t2^2*(1-t2)^2 + 6*(1-t1)*t3^2*(1-t3)^2
p4 = 4*t1*t2^3*(1-t2) + 4*(1-t1)*t3^3*(1-t3)
p5 = t1*t2^4 + (1-t1)*t3^4
S = QQ[x1,x2,x3,x4]
K = frac S
L = frac R
phi = map(K, L, {x1/x4, x2/x4, x3/x4})
p = {p1,p2,p3,p4,p5}
A = p / (i -> numerator phi(i)) | {sub(x4, S)}
for f in A list #(terms f)
A#3
I0 = preLikelihoodIdeal A;
(flatten entries mingens I0) / degree
(flatten entries mingens I0)#0
(flatten entries mingens I0)#1
Example 3.5 (4 conics and 1 line)
The independence model on two binary random variables is given by the following (homogeneous) parametrization:
R = QQ[x_1..x_3]
A = {x_1*x_2, (x_3-x_1)*x_2, x_1*(x_3-x_2), (x_3-x_1)*(x_3-x_2), x_3}
I0 = preLikelihoodIdeal A
isPrime I0
We see that the arrangement is not gentle. We compute the ML degree as follows:
I = likelihoodIdeal A
multidegree I
The leading coefficient is 1, and this is the ML degree.
Example 3.5 (Random ternary quadrics)
The claims made in Example 3.5 about the arrangement defined by 4 random ternary quadrics and a single coordinate hyperplane can be verified as follows.
R = ZZ/10007[x_1, x_2, x_3]
A = toList(1..4) / (i -> random(2, R)) | {x_3}
I0 = preLikelihoodIdeal A
(flatten entries mingens I0) / (i -> degree i)
isPrime I0
multidegree I0
isTame A
Example 3.6 (Independence model)
The following code computes the likelihood ideal and the ML degree of the independence model in its standard parametrization. The pre-likelihood ideal is prime and hence equals the likelihood ideal.
R = QQ[a0,b0,a1,b1]
imgs = {a0*b0, a1*b0, a0*b1, a1*b1}
f = sum imgs
A = gens R | {f}
I0 = preLikelihoodIdeal A;
isPrime I0
multidegree I0
A minimal parametrization of the independence model is encoded in the following arrangement:
R = QQ[x,y,z]
A = {x,y,z, x*y + x*z + y*z + z^2}
I0 = preLikelihoodIdeal A;
isPrime I0
multidegree I0
Notice that the multidegrees of the likelihood ideals for the two parametrizations differ by a factor \(T_1\).
Example 5.1 (Octahedron)
With the help of the following code we show that the graphical arrangement of the octahedron is not gentle.
R = QQ[x_1..x_6]
L = {x_1-x_2, x_1-x_3, x_1-x_5, x_1-x_6, x_2-x_3, x_2 -x_4, x_2 - x_6, x_3-x_4, x_3-x_5, x_4-x_5, x_4-x_6, x_5-x_6}
I0 = preLikelihoodIdeal L
I1 = I0:sub(L#0, ring I0) -- takes about 2 mins
I1 == I0
isPrime I1 -- returns true after about 10 mins
This shows that the likelihood ideal is given by \(I_1\), not \(I_0\). The ideals \(I_1\) and \(I_0\) differ by a single generator, which we can compute via
f = ((I1)_*)#-1;
degree f
#(terms f)
Therefore, the prime decomposition of \(I_0\) is given by \(I_0 = I_1 \cap P\), where \(P\) can be computed via P = I0 : f.
Graphical arrangements up to 6 vertices (Theorem 5.2)
Theorem 5.2 states that the graph of the octahedron gives the only non-gentle graphical arrangements among all such arrangements for graphs with up to six vertices. We prove this via exhaustive computation which we present here.
We make use of the graph database provided by Brendan McKay at this website. Save the data for all connected graphs on six vertices in a file named graphs6c.txt. We use the Julia programming language to convert the data into arrangements readable by Macaulay2. Download the files Project.toml and Manifest.toml. Then type ] in a Julia terminal started from the directory where you saved the files to go into the package manager. Type instantiate to install all necessary dependencies. Then you can run the following code. It suffices to only consider non-chordal graphs, since chordal graphs give rise to free and hence gentle arrangements.
using Oscar
using Graphs
using GraphIO
graphs6 = collect(values(loadgraphs("graphs6c.txt", GraphIO.Graph6.Graph6Format())))
R, (x1,x2,x3,x4,x5,x6) = graded_polynomial_ring(QQ, 6)
function mycycles(G)
DG = SimpleDiGraph(Graphs.nv(G))
for e in Graphs.edges(G)
Graphs.add_edge!(DG, e.dst, e.src)
Graphs.add_edge!(DG, e.src, e.dst)
end
cycs = filter(k->length(k)>2,simplecycles(DG))
#return unique([sort(k) for k in cycs])
return cycs
end
function has_chord(cyc, G)
for e in Graphs.edges(G)
if e.src in cyc && e.dst in cyc
pos_src = findall(x->x==e.src, cyc)[1]
if pos_src == 1
if cyc[2] != e.dst && cyc[length(cyc)] != e.dst
return true
end
elseif pos_src == length(cyc)
if cyc[1] != e.dst && cyc[length(cyc)-1] != e.dst
return true
end
else
if cyc[pos_src-1] != e.dst && cyc[pos_src+1] != e.dst
return true
end
end
end
end
return false
end
function chordless_cycles(G)
cycles = mycycles(G)
return filter(c -> has_chord(c,G)==false,cycles)
end
function longest_cl_cycle(G)
mcs = chordless_cycles(G)
if length(mcs) == 0
return -1
end
return maximum([length(k) for k in mcs])
end
nc_graphs6 = filter(x-> longest_cl_cycle(x) > 3, graphs6) # only non-chordal graphs
p = sortperm(Graphs.ne.(nc_graphs6))
nc_graphs6 = nc_graphs6[p]
function getPolynomials(graph, ring)
polList = []
for edge in Graphs.edges(graph)
append!(polList, [gens(ring)[edge.src]-gens(ring)[edge.dst]])
end
return Vector{MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}}(polList)
end
function graphs_to_m2()
for (i, G) in enumerate(nc_graphs6)
P = getPolynomials(G, R)
open("GraphsLinearForms/$i.txt", "w") do file
print(file, "{")
l = length(P)
for (j, p) in enumerate(P)
if j == l
print(file, p)
else
print(file, p)
print(file, ",")
end
end
print(file, "}")
end
end
end
graphs_to_m2()
This creates a directory GraphsLinearForms with 54 text files, each containing the linear equations defining the arrangement corresponding to one of the 54 graphs. You can then run the following Macaulay2 code to check for gentleness in each instance. Note that we do a check after each saturation step to reduce the run time. Warning: it might take several days for the code to terminate.
R = QQ[x1,x2,x3,x4,x5,x6]
As = apply(toList(1..54), i -> value get concatenate {"GraphsLinearForms/", toString i, ".txt"})
for i from 0 to 53 do(
print(i);
A := As#i;
I0 := preLikelihoodIdeal A;
notGentle := false;
for f in A do(
I1 := saturate(I0, sub(f, ring I0));
if I1 != I0 then(
notGentle = true;
print("not gentle");
break;
);
);
if notGentle then(
continue;
);
S := transpose jacobian matrix {A};
Sr := sub (S, for i in gens ring S list i=>random(QQ));
n := rank Sr;
SL := sub(minors_n S, ring I0);
I1 := saturate(I0, SL);
if I1 != I0 then(
print("not gentle");
continue;
);
print("gentle");
)
Proposition 5.6 (Edelman-Reiner arrangement)
The following provides an example where the restriction \(B\) of a gentle arrangement \(A\) is not gentle. Moreover, the restriction \(B\) is not tame any more.
R = QQ[x_1..x_5]
A = apply(toList (set {-1,1})^**4 / splice / splice, i -> x_1 + (i#0)*x_2 + (i#1)*x_3 + (i#2)*x_4 + (i#3)*x_5) | (flatten entries vars R)
D = logDerModule A
isFreeModule prune D -- free
S = QQ[y_1..y_4]
B = apply(delete((0,0,0,0), toList (set {0,1})^**4 / splice / splice), i -> (i#0)*y_1 + (i#1)*y_2 + (i#2)*y_3 + (i#3)*y_4)
Drestr = logDerModule B
isFreeModule prune Drestr -- not free
Omega1 = Hom(Drestr, S)
res Omega1 -- pdim 2, not tame
I0 = preLikelihoodIdeal B
I1 = I0 : y_1
I0 == I1 -- false, B is not gentle
Example 6.3 (No 3-way interaction)
We first compute the likelihood ideal for the no 3-way interaction model in its most standard (over-)parametrization.
R = QQ[a00,a01,a10,a11, b00,b01,b10,b11, c00,c01,c10,c11]
imgs = {a00*b00*c00, a10*b10*c00, a01*b00*c10, a11*b10*c10, a00*b01*c01, a10*b11*c01, a01*b01*c11, a11*b11*c11}
f = sum imgs
A = gens R | {f}
I0 = preLikelihoodIdeal A;
I1 = I0 : a00;
I = I1 : sub(f^2, ring I0); -- this is the likelihood ideal
isPrime I -- true
The pre-likelihood ideal itself has 25 minimal primes as can be checked via minimalPrimes I0.
An alternative parametrization using fewer parameters is given as follows. This allows to compute the ML degree easily.
R = QQ[x_1..x_7]
imgs = {x_1^6, x_2*x_1^5, x_3*x_1^5, x_1^5*x_4, x_2*x_3*x_5*x_1^3, x_3*x_4*x_6*x_1^3, x_2*x_4*x_7*x_1^3, x_7*x_2*x_3*x_4*x_5*x_6}
f = sum imgs
A = gens R | {f}
I0 = preLikelihoodIdeal A;
I = I0 : sub(x_1*x_2*x_3*x_4*x_5, ring I0);
isPrime I -- true
multidegree I
Example 6.4 (CEGM model)
The following code computes the pre-likelihood ideal for the CEGM model.
R = QQ[x1,x2,x3,x4,x5];
F = {x1,x2,x3,x4,x5,x1-x2,x1-x3,x1-x5,x2-x5,x2-x4,x3-x4,x3-x5,x4-x5,
x1*x4-x2*x3,x1*x4-x2*x3-x1+x2+x3-x4};
I0 = preLikelihoodIdeal F;
Unfortunately, it is very challenging to compute a Gröbner bases of \(I_0\). We can compute a numerical irreducible decomposition using the HomotopyContinuation package by Paul Breiding and Sascha Timme in Julia. First we need to convert the Macaulay2 output into a format usable by HomotopyContinuation. Type
f = "equations.txt" << toExternalString flatten entries mingens I0;
close f
to save the defining equations of \(I_0\) in a file named equations.txt. Then, in a Julia terminal, execute the following commands.
using Oscar
using HomotopyContinuation
function polyOscarToHomCon(p, vars)
homconp = 0
for t in terms(p)
expo = [e for e in exponents(t)][1]
coef = [c for c in Oscar.coefficients(t)][1]
homconp += Rational(coef) * prod(map(^, vars, expo))
end
return homconp
end
function idealOscarToHomConSystem(I)
vars = gens(base_ring(I))
@var a[1:(size(vars)[1])]
return ([polyOscarToHomCon(p, a) for p in gens(I)], a)
end
R, (x1,x2,x3,x4,x5, s_1,s_2,s_3,s_4,s_5,s_6,s_7,s_8,s_9,s_10,s_11,s_12,s_13,s_14,s_15) = polynomial_ring(QQ, "y" => 1:20)
equs = []
open("equations.txt", "r") do file
append!(equs, eval(Meta.parse("[" * readline(file)[2:end-1] * "]")))
end
I0 = ideal(equs)
S = idealOscarToHomConSystem(I0)[1]
Finally, we can compute a numerical irreducible decomposition via
nid(S)
This finds 25 irreducible components and hence the arrangement cannot be gentle.
Numerical irreducible decomposition is an interesting experimental tool to check whether an arrangement is not gentle when symbolic techniques fail. However, one cannot conclude that if one finds only a single irreducible component that the arrangement is necessarily gentle. For example, numerical irreducible decomposition applied to the octahedron arrangement yields a single component, but the arrangement is not gentle.
Project page created: 31/10/2024.
Project contributors: Thomas Kahle, Lukas Kühne, Leonie Mühlherr, Bernd Sturmfels, Maximilian Wiesmann.
Corresponding author of this page: Maximilian Wiesmann, wiesmann@mis.mpg.de.
Code written by: Thomas Kahle, Lukas Kühne, Leonie Mühlherr, Bernd Sturmfels, Maximilian Wiesmann.
Software used: Macaulay2 (Version 1.20), Julia (Version 1.9.1).
System setup used: MacBook Pro with macOS Monterey 12.5, Processor Apple M1 Pro, Memory 16 GB LPDDR5.
License for code of this project page: MIT License (https://spdx.org/licenses/MIT.html).
License for all other content of this project page (text, images, …): CC BY 4.0 (https://creativecommons.org/licenses/by/4.0/).
Last updated 13/11/2024.