Skip to content
Snippets Groups Projects

Update of Code and Data in PositiveDelPezzoGeometry and GraphCurveMatroids

Merged Update of Code and Data in PositiveDelPezzoGeometry and GraphCurveMatroids
develop into master
Merged Alheydis Geiger requested to merge develop into master
Viewing commit 445f8ab8
Prev
Show latest version
1 file
+ 295
− 0
Compare changes
  • Side-by-side
  • Inline
source/LikelihoodReflexive/index.rst 0 → 100644
+ 295
− 0
==========================================
Likelihood Geometry of Reflexive Polytopes
==========================================
| This page contains the source code and explanations related to the computational
| results presented in the paper:
| Carlos Améndola, Janike Oldekop: Likelihood Geometry of Reflexive Polytopes
| ARXIV: `<https://arxiv.org/abs/2311.13572>`_
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 :math:`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 <http://www2.macaulay2.com/Macaulay2/>`_ (v1.20) and `julia <https://julialang.org/>`_ (v1.8.5). For a smooth running of all computations in Macaulay2 we recommend to ensure access to the following packages at the beginning.
****
|
|
| `ReflexivePolytopesDB.m2 <http://www2.macaulay2.com/Macaulay2/doc/Macaulay2-1.21/share/doc/Macaulay2/ReflexivePolytopesDB/html/index.html>`_
| `Polyhedra.m2 <https://faculty.math.illinois.edu/Macaulay2/doc/Macaulay2-1.20/share/doc/Macaulay2/Polyhedra/html/index.html>`_
| `FourTiTwo.m2 <http://www2.macaulay2.com/Macaulay2/doc/Macaulay2-1.20/share/doc/Macaulay2/FourTiTwo/html/index.html>`_
| `AlgebraicOptimization.m2 <https://github.com/Macaulay2/Workshop-2020-Cleveland/tree/ISSAC-AlgOpt/alg-stat/AlgebraicOptimization>`_
| `LatticePolytopes.m2 <http://www2.macaulay2.com/Macaulay2/doc/Macaulay2/share/doc/Macaulay2/LatticePolytopes/html/index.html>`_
The computations in Julia are based on `HomotopyContinuation.jl <https://www.juliahomotopycontinuation.org/>`_ (v2.8.2). In addition, `Oscar.jl <https://www.oscar-system.org>`_ (v0.11.3) was used to access the `polyDB <https://polydb.org/#collection=Polytopes.Lattice.SmoothReflexive>`_ .
Computations presented in Section 7 also used Mathematica. `Mathematica <https://www.wolfram.com/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.
1. Given :math:`A' \in \mathbb{Z}^{(d+1) \times n}`, the score equations :math:`A' \psi(\theta) = \frac{1}{u_+} A'u` are determined using Macaulay2.
2. Given the score equations, the number of complex solutions for generic data :math:`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: :download:`subsection02.m2 <subsection02.m2>`.
Example: Consider
.. math::
A' = \begin{bmatrix} 1 & 1 & 1 & 1 \\ 2 & 1 & 1 & 0 \\ 1 & 2 & 1 & 0 \\ 1 & 1 & 2 & 0 \end{bmatrix}.
To determine the score equations for :math:`A'`, the following code can be used.
.. code-block:: macaulay2
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.
.. code-block:: macaulay2
:linenos:
:emphasize-lines: 3,4,10,26,29
-- write Julia file
polynomialjl = p -> (
replace("_1","[1]",replace("_2","[2]",replace("_3","[3]",
replace("_4","[4]",toExternalString p))))
)
writeSystem = method(Options=>{})
writeSystem (List, File) := o -> (P,f) -> (
variables = {"x1","x2","x3","x4","u[1:4]"};
f << "## This file was generated by Macaulay2 to compute the ML degree via homotopy continuation.\n\n";
f << "using HomotopyContinuation";
f << "\n\n";
f << "@var " | concatenate(mingle(variables, #variables:1));
eqnCommas:=(#P-1):", ";
f << "\n\n";
f << "F = System(\n";
f << " [\n";
f << " "|concatenate mingle(apply(P,e->polynomialjl e),
eqnCommas);
f << "\n";
f << " ];\n";
f << " parameters = u\n" ;
f << ")\n\n";
f << "## Solve generically:\n\n";
f << "u0 = rand(ComplexF64,4)\n";
f << "result_u0 = solve(F, target_parameters = u0)\n\n";
f << "## Solve for random parameter values\n\n";
f << "data = [rand((0:100000),4)]\n";
f << "data_points = solve(\n";
f << " F,\n";
f << " solutions(result_u0),\n";
f << " start_parameters = u0,\n";
f << " target_parameters = data,\n";
f << ")";
)
writeSystem (List, String) :=o -> (P,filename) -> (
f := openOut filename;
writeSystem(P,f);
close f;
)
-- create Julia file called example.jl
writeSystem(system, "example.jl")
Run the code in :download:`example.jl <example.jl>`. The number of non-singular solutions corresponds to the ML degree of four.
.. attention::
The code above works for all :math:`A' \in \mathbb{Z}^{(d+1) \times n}` with :math:`n \le 4` and :math:`d \le 3`. For larger :math:`n` or :math:`d` the following adjustments are necessary:
In function :math:`\texttt{polynomialjl}`, the replace function must be nested at least :math:`n` times. Replace :math:`\texttt{toExternalString p}` with :math:`\texttt{replace("_5","[5]",toExternalString p)}`, etc.
In line [10] use variables :math:`\textup{x}k` for all :math:`k \in [d+1]` and :math:`\textup{u}[1:n]`.
In line [26] and [29] replace the :math:`4` in the second argument with the specific :math:`n`.
3 Reflexive Polygons
--------------------
The ML degrees of all reflexive polygons were computed using the Macaulay2 package `AlgebraicOptimization.m2 <https://github.com/Macaulay2/Workshop-2020-Cleveland/tree/ISSAC-AlgOpt/alg-stat/AlgebraicOptimization>`_. The code is available for download here: :download:`MLdeg-2D.m2 <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:
1. The ML degree
2. The degree of the toric variety
3. The f-vector
4. The number of lattice points contained in the polytope
5. The number of generators of the toric ideal
The results can be downloaded here: :download:`data.txt <data.txt>`.
.. At this point Julia (v1.6.2) and HomotopyContinuation.jl (v2.6.0) were used. With these versions, the computation of all ML degrees took about 15 hours.
Fundamental to all computations is access to the `Kreuzer and Skarke database <http://hep.itp.tuwien.ac.at/~kreuzer/CY/CYk3.html>`_, which is obtained with the following Macaulay2 code.
.. code-block:: macaulay2
-- access to KS database
needsPackage "ReflexivePolytopesDB";
polytopes = kreuzerSkarkeDim3();
The polytopes in the KS database are numbered :math:`\texttt{0}` to :math:`\texttt{4318}`. The following Macaulay2 code gives a matrix whose columns are given by the lattice points of polytope :math:`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.
.. code-block:: macaulay2
-- 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 :download:`MLdeg-3D.m2 <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: :download:`ScoreEquations.zip <ScoreEquations.zip>`. The Julia file :download:`HC-3D.jl <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 :download:`computations-3D.m2 <computations-3D.m2>`.
**Example 4.2**
The ML degree of the scaled Polytope :math:`\texttt{0}` was computed using :download:`ScaledPolytope0.jl <ScaledPolytope0.jl>`.
**Example 4.3**
The ML degree of the scaled Polytope :math:`\texttt{132}` was computed using :download:`ScaledPolytope132.jl <ScaledPolytope132.jl>`.
**Table 3 - Smooth Polytopes**
The `polyDB <https://polydb.org/#collection=Polytopes.Lattice.SmoothReflexive>`_ database of smooth reflexive polytopes can be accessed using `Oscar <https://www.oscar-system.org>`_. The vertices of these polytopes can be output with :download:`Smooth3D.jl <Smooth3D.jl>`. To check whether two smooth polytopes are isomorphic, the Macaulay2 package `LatticePolytopes.m2 <http://www2.macaulay2.com/Macaulay2/doc/Macaulay2/share/doc/Macaulay2/LatticePolytopes/html/index.html>`_ can be used. Table 3 in the paper can be reconstructed using :download:`Isomorphic.m2 <Isomorphic.m2>`.
5 Hypercubes and Cross Polytopes
--------------------------------
**Example 5.13**
The ML degree of the scaled :math:`2`-cross polytope in Example 5.13 was computed using :download:`Example513.jl <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:
****
|
|
| :download:`Q.zip <Q.zip>`
| :download:`R.zip <R.zip>`
| :download:`S.zip <S.zip>`
| :download:`T.zip <T.zip>`
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:
****
|
|
| :download:`Constructions.zip <Constructions.zip>`
The common roots used in Theorem 7.9 were computed using Mathematica: :download:`CommonRoots.nb <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.