Likelihood Geometry of Reflexive Polytopes
ABSTRACT: We study the problem of maximum likelihood (ML) estimation for statistical models defined by reflexive polytopes. Our focus is on the maximum likelihood degree of these models as an algebraic measure of complexity of the corresponding optimization problem. We compute the ML degrees of all 4319 classes of three-dimensional reflexive polytopes, and observe some surprising behavior in terms of the presence of gaps between ML degrees and degrees of the associated toric varieties. We interpret these drops in the context of discriminants and prove formulas for the ML degree for families of reflexive polytopes, including the hypercube and its dual, the cross polytope, in arbitrary dimension. In particular, we determine a family of embeddings for the \(d\)-cube that implies ML degree one. Finally, we discuss generalized constructions of families of reflexive polytopes in terms of their ML degrees.
The section numbers (except for 0) refer to the sections in the paper.
0 Preliminaries
0.1 Required Software
The main computations mentioned in the article were performed using Macaulay2 (v1.20) and julia (v1.8.5). For a smooth running of all computations in Macaulay2 we recommend to ensure access to the following packages at the beginning.
The computations in Julia are based on HomotopyContinuation.jl (v2.8.2). In addition, Oscar.jl (v0.11.3) was used to access the polyDB .
Computations presented in Section 7 also used Mathematica. Mathematica (v12.3.1).
0.2 ML Degree Computations
The computations of most ML degrees presented in the paper are carried out in two steps.
Given \(A' \in \mathbb{Z}^{(d+1) \times n}\), the score equations \(A' \psi(\theta) = \frac{1}{u_+} A'u\) are determined using Macaulay2.
Given the score equations, the number of complex solutions for generic data \(u\) is computed numerically using Julia, especially HomotopyContinuation.jl.
This idea is based on https://www.juliahomotopycontinuation.org/guides/macaulay2/.
The code of this subsection can be downloaded here: subsection02.m2.
Example: Consider
To determine the score equations for \(A'\), the following code can be used.
A = matrix{{1,1,1,1},{2,1,1,0},{1,2,1,0},{1,1,2,0}};
-- determine score equations
MLeqs = {coeffRing=>QQ} >> o -> (expr) -> (
m = numgens source expr;
R := (o.coeffRing)[vars {53..(52+(numgens target expr))},u_1..u_m];
A := substitute(expr,R);
U = u_1..u_m;
U = toList U;
U = transpose matrix {for i in U list i_R};
N = sum for i from 1 to m list u_i;
X := for i in 0..(m-1) list R_(entries expr_i);
X = transpose matrix {X};
system := ideal (A*((N * X) - U));
return first entries gens system
)
system = MLeqs A;
In order to write the output into a suitable Julia file, that can be used to compute the number of non-singular solutions, one can use the following code.
1 -- write Julia file
2 polynomialjl = p -> (
3 replace("_1","[1]",replace("_2","[2]",replace("_3","[3]",
4 replace("_4","[4]",toExternalString p))))
5 )
6
7 writeSystem = method(Options=>{})
8
9 writeSystem (List, File) := o -> (P,f) -> (
10 variables = {"x1","x2","x3","x4","u[1:4]"};
11 f << "## This file was generated by Macaulay2 to compute the ML degree via homotopy continuation.\n\n";
12 f << "using HomotopyContinuation";
13 f << "\n\n";
14 f << "@var " | concatenate(mingle(variables, #variables:1));
15 eqnCommas:=(#P-1):", ";
16 f << "\n\n";
17 f << "F = System(\n";
18 f << " [\n";
19 f << " "|concatenate mingle(apply(P,e->polynomialjl e),
20 eqnCommas);
21 f << "\n";
22 f << " ];\n";
23 f << " parameters = u\n" ;
24 f << ")\n\n";
25 f << "## Solve generically:\n\n";
26 f << "u0 = rand(ComplexF64,4)\n";
27 f << "result_u0 = solve(F, target_parameters = u0)\n\n";
28 f << "## Solve for random parameter values\n\n";
29 f << "data = [rand((0:100000),4)]\n";
30 f << "data_points = solve(\n";
31 f << " F,\n";
32 f << " solutions(result_u0),\n";
33 f << " start_parameters = u0,\n";
34 f << " target_parameters = data,\n";
35 f << ")";
36 )
37
38 writeSystem (List, String) :=o -> (P,filename) -> (
39 f := openOut filename;
40 writeSystem(P,f);
41 close f;
42 )
43
44 -- create Julia file called example.jl
45 writeSystem(system, "example.jl")
Run the code in example.jl. The number of non-singular solutions corresponds to the ML degree of four.
Attention
The code above works for all \(A' \in \mathbb{Z}^{(d+1) \times n}\) with \(n \le 4\) and \(d \le 3\). For larger \(n\) or \(d\) the following adjustments are necessary:
In function \(\texttt{polynomialjl}\), the replace function must be nested at least \(n\) times. Replace \(\texttt{toExternalString p}\) with \(\texttt{replace("_5","[5]",toExternalString p)}\), etc.
In line [10] use variables \(\textup{x}k\) for all \(k \in [d+1]\) and \(\textup{u}[1:n]\).
In line [26] and [29] replace the \(4\) in the second argument with the specific \(n\).
3 Reflexive Polygons
The ML degrees of all reflexive polygons were computed using the Macaulay2 package AlgebraicOptimization.m2. The code is available for download here: MLdeg-2D.m2.
4 Reflexive Polyhedra
As presented in the paper, the following quantities were computed for each of the 4319 three-dimensional reflexive polytopes:
The ML degree
The degree of the toric variety
The f-vector
The number of lattice points contained in the polytope
The number of generators of the toric ideal
The results can be downloaded here: data.txt.
Fundamental to all computations is access to the Kreuzer and Skarke database, which is obtained with the following Macaulay2 code.
-- access to KS database
needsPackage "ReflexivePolytopesDB";
polytopes = kreuzerSkarkeDim3();
The polytopes in the KS database are numbered \(\texttt{0}\) to \(\texttt{4318}\). The following Macaulay2 code gives a matrix whose columns are given by the lattice points of polytope \(i\). For further computations it is necessary to consider a translation of the polytope such that all lattice points are nonnegative and to add the all-ones row vector.
-- sufficient statistic with all-ones vector and nonnegative entries
matrixA = i -> (
P = polytopes_i;
P = matrix P;
P = convexHull P;
L = latticePoints P;
L = matrix{L};
n = numgens source L;
d = numgens target L;
seqn = 0..n-1;
seqd = 0..d-1;
listn = toList seqn;
listd = toList seqd;
listd = apply(listd, i -> min apply(listn, j -> L_(i,j)));
T = matrix table(d,n,(i,j) -> -listd_i);
L = L + T;
E = transpose matrix vector(apply(listn, i -> 1));
L = E||L
)
-- example
matrixA 0 -- gives A' for Polytope 0
For all reflexive polyhedra, the ML degrees were computed as described in Subsection 0.2. The Macaulay2 code MLdeg-3D.m2 generates a file which contains the score equations for chosen polytopes. The systems of all reflexive polyhedra are available for download here: ScoreEquations.zip. The Julia file HC-3D.jl can access these files to compute the ML degrees.
Attention
Some adjustments are necessary. More detailed information is given in the files as comments.
The computations of the quantities mentioned in 2nd to 5th were carried out with computations-3D.m2.
Example 4.2
The ML degree of the scaled Polytope \(\texttt{0}\) was computed using ScaledPolytope0.jl.
Example 4.3
The ML degree of the scaled Polytope \(\texttt{132}\) was computed using ScaledPolytope132.jl.
Table 3 - Smooth Polytopes
The polyDB database of smooth reflexive polytopes can be accessed using Oscar. The vertices of these polytopes can be output with Smooth3D.jl. To check whether two smooth polytopes are isomorphic, the Macaulay2 package LatticePolytopes.m2 can be used. Table 3 in the paper can be reconstructed using Isomorphic.m2.
5 Hypercubes and Cross Polytopes
Example 5.13
The ML degree of the scaled \(2\)-cross polytope in Example 5.13 was computed using Example513.jl.
6 Reflexive Simplices
For each of the four families we computed the ML degrees as described in Subsection 0.2. The files are available for download here:
The matrices containing the vertices as columns were determined according to the constructions in the paper.
7 Constructing Reflexive Polytopes
The ML degrees in Table 4 and 5 were computed as described in Subsection 0.2. For each reflexive polygon, the following file contains a folder with the corresponding Macaulay2 and Julia files:
The common roots used in Theorem 7.9 were computed using Mathematica: CommonRoots.nb.
Project page created: 29/11/2023
Project contributors: Carlos Améndola, Janike Oldekop
Corresponding author of this page: Janike Oldekop, oldekop@math.tu-berlin.de
Software used: Macaulay2 (Version 1.20), Julia (version 1.8.5), Mathematica (version 12.3.1)
System setup used: MacBook Pro with macOS Venture 13.3.1, Processor 2.3 GHz Dual Core Intel Core i5, Memory 8 GB
License for code of this project page: MIT License (https://spdx.org/licenses/MIT.html)
License for all other content of this project page: CC BY 4.0 (https://creativecommons.org/licenses/by/4.0/)
Last updated 29/11/2023.