Marginal Independence Models
ABSTRACT: We impose rank one constraints on marginalizations of a tensor, given by a simplicial complex. Following work of Kirkup and Sullivant, such marginal independence models can be made toric by a linear change of coordinates. We study their toric ideals, with emphasis on random graph models and independent set polytopes of matroids. We develop the numerical algebra of parameter estimation, using both Euclidean distance and maximum likelihood, and we present a comprehensive database of small models.
We consider \(n\) discrete random variables \(X_1, X_2, \ldots, X_n\) where \(X_i\) takes values in the finite set \([d_i] = \{ 1, 2, \ldots, d_i \}\). A joint distribution for these random variables is a tensor \(P = (p_{i_1 i_2 \cdots i_n})\) of format \(d_1 \times d_2 \times \cdots \times d_n\) whose entries are nonnegative real numbers that sum to \(1\).
Every subset \(\sigma \subseteq [n]\) defines a marginal distribution \(P_\sigma\) by summing out all dimensions of the tensor outside of \(\sigma\). This is a tensor of format \(\times_{i \in \sigma} d_i\) which represents the distribution of the subset of variables \(X_i\), \(i \in \sigma\). The random variables indexed by \(\sigma\) are completely independent if the marginal tensor \(P_\sigma\) has rank \(1\). Rank \(0\) is impossible due to the sum to one condition, so marginal independence is equivalent to the tensor having the minimal rank possible. This notion of complete independence is hereditary: if \(P_\sigma\) has rank \(1\), then all subsets of \(\sigma\) have the same property. Hence the set of completely independent subcollections of any tensor \(P\) is a simplicial complex. These complexes always include the singletons \(\{i\}\) whose marginal tensors are just non-zero vectors.
Let, conversely, \(\Sigma\) be any simplicial complex with vertex set \([n]\). We assume \(\{i\} \in \Sigma\) for all \(i \in [n]\). The marginal independence model \(\mathcal{M}_\Sigma\) is the set of distributions \(P \in \Delta_{D-1}\) such that the marginal tensor \(P_\sigma\) has rank \(1\) for all \(\sigma \in \Sigma\).
The rank \(1\) constraints translate to vanishing of all \(2 \times 2\) minors of flattenings of each marginal \(P_\sigma\). This gives the ideal of the marginal independence model in natural joint probability coordinates, i.e., the entries of \(P\). Its generators are quadratic polynomials. Under a linear change of coordinates, the ideal becomes toric, as proved in the paper. The new coordinates are called Möbius coordinates. This is implemented in our \(\verb|Macaulay2|\) package MarginalIndependenceModels.
Polytope, ideal, dimension and degree
The following example demonstrates the usage of this package at the binary \(3\)-cycle, i.e., the simplicial complex with facets \(\{ 12, 13, 23 \}\) imposing independence conditions on three binary random variables. The example also appears in S. Sullivant: Algebraic Statistics (2018) and G. Kirkup: Random variables with completely independent subcollections (2007).
needsPackage "MarginalIndependenceModels"
Sigma = {{1,2},{1,3},{2,3}}
I = marginalIndependenceIdeal sigma
It is enough to specify the simplicial complex Sigma in terms of its facets;
the set of all faces is computed on the fly if needed. The matrix defining this
toric ideal is given by the following polytope:
needsPackage "Polyhedra"
P = binaryPolytope Sigma
vertices P
To appreciate the linear change of coordinates, take a look at the map:
phi = moebius Sigma
The target and source rings are in natural probabilities and relevant Möbius coordinates, respectively.
target phi
source phi
Essential geometric information about this model can be computed and checked against the Kirkup’s paper:
codim I, degree I, betti mingens I
Caveat: The ambient polynomial ring of I uses only the relevant Möbius coordinates,
while the full model lives in a larger space. Hence there will be a discrepancy when computing
the model dimension in Macaulay2. The dimension of the model as a projective variety can
be obtained as 2^3 - 1 - codim I.
We refer to the documentation of the package for more information.
Minimal generators for the real projective plane model
Example 23 is now quick and easy to confirm:
RP = {{1,2,3},{1,2,4},{1,3,5},{1,4,6},{1,5,6},{2,3,6},{2,4,5},{2,5,6},{3,4,5},{3,4,6}}
betti mingens marginalIndependenceIdeal RP
The ternary 3-cycle
Our code only covers models on binary random variables. Checking the free resolution of the ternary model of the 3-cycle in Example 1 requires manual coding:
R = QQ[
-- Parameters for the Segre variety
z, a1, b1, c1, a2, b2, c2,
-- Relevant Moebius coordinates, S stands for "sum"
q11S, q12S, q21S, q22S, q1S1, q1S2, q2S1, q2S2,
qS11, qS12, qS21, qS22,
q1SS, q2SS, qS1S, qS2S, qSS1, qSS2, qSSS,
MonomialOrder => Eliminate 7
];
-- Parametrization of the relevant Moebius coordinates as in Section 3.
I = ideal(
q11S - z*a1*b1*(1-c1-c2),
q12S - z*a1*b2*(1-c1-c2),
q21S - z*a2*b1*(1-c1-c2),
q22S - z*a2*b2*(1-c1-c2),
-- etc.
);
-- Implicitize
M = eliminate({z, a1, b1, c1, a2, b2, c2}, I);
betti mingens M
betti res M
The entire file is attached to this page: Cycle3-3x3x3.m2.
Euclidean distance and maximum likelihood degrees
To compute aED, pED and ML degrees of our models, we used the \(\verb|Julia|\) package HomotopyContinuation.jl. We employ the parametrization of the model in Möbius coordinates to write the critical equations for the unconstrained optimization problem. In these coordinates, marginal independences correspond to factorizations of the relevant Möbius coordinates into the singletons: \(q_\sigma = \prod_{i \in \sigma} q_i\).
using HomotopyContinuation
@var q1 q2 q3 q12 q13 q23 q123 z
sample = randn(ComplexF64, 8)
p000,p001,p010,p011,
p100,p101,p110,p111 = sample
diffs = [
-p000 + z*(q123),
-p001 + z*(q12-q123),
-p010 + z*(q13-q123),
-p100 + z*(q23-q123),
-p101 + z*(q2-q12-q23+q123),
-p110 + z*(q3-q13-q23+q123),
-p011 + z*(q1-q12-q13+q123),
-p111 + z*(1-q1-q2+q12-q3+q13+q23-q123)
]
dist = sum([d^2 for d in diffs])
This sets up the parametrization and the objective function. The 3-cycle corresponds to the below factorization:
subslist = [ z=>1, q12=>q1*q2, q13=>q1*q3, q23=>q2*q3 ]
dist = subs(dist, subslist...)
HomotopyContinuation.jl can write out the critical equations and compute the number of complex solutions for the random complex data chosen above.
vars = variables(dist)
eqns = differentiate(dist, vars)
R = solve(eqns; show_progress = false,
tracker_options = TrackerOptions(
automatic_differentiation = 3,
parameters = :conservative
)
)
C = certify(eqns, R; show_progress = false)
println(ndistinct_certified(C))
To change the model, only the variable subslist has to be altered. To automate the
computation of degrees, the above program was embedded into a Perl program which
modifies the subslist for a given simplicial complex, runs the Julia program and
parses its output. The parametrizations are obtained via our \(\verb|Macaulay2|\) code and hardcoded
in each file:
Plain text database of complexes, models and degrees
Our findings are summarized in plain text database files. They list each simplicial complex and the data for its model in one line. This includes the facet description of the complex, its dimension, number of faces, facets and f-vector to make it easy to find a given complex by hand. We indicate whether the complex is a matroid or even the complex of forests of an undirected simple graph. For its toric ideal we list the dimension, codimension and degree, moreover its m-vector, i.e., for each degree the number of minimal generators. Finally, whenever available, we give the aED, pED and ML degree. The complexes are listed up to isomorphy and they always contain \(\{i\}\) for all \(i \in [n]\). For example:
complex cdim f-vector faces facets pure matroid graph pdim codim degree m-vector aED pED ML
[1,2,3] 0 (1,3) 4 3 Y Y Y 7 0 1 () 1 1 1
[12,3] 1 (1,3,1) 5 2 N N N 6 1 2 (1) 5 2 1
[12,13] 1 (1,3,2) 6 2 Y Y Y 5 2 3 (3) 5 2 9
[12,13,23] 1 (1,3,3) 7 3 Y Y Y 4 3 5 (5) 21 19 17
[123] 2 (1,3,3,1) 8 1 Y Y Y 3 4 6 (9) 17 6 1
Some degrees in these tables have a trailing “+” sign. As explained in the paper, these are certified lower bounds on the true degrees for which the homotopy continuation solver did not produce absolutely consistent results (but the relative error in each case was always below 5%).
Ideal data for matroids on \(n=7\)
Below are the logs of \(\verb|Macaulay2|\) computations for matroids (up to isomorphy) on \(n=7\) elements:
The file consists of records, one for each matroid, listing the facets of the simplicial complex (that is, the bases of the matroid), the dimension, codimension and degree of its binary ideal, as explained above, followed by its Betti table. The results were obtained using our \(\verb|MarginalIndependenceModels|\) package, as follows:
-- FF is initialized to the simplicial complex...
M = marginalIndependenceIdeal(FF);
print(FF);
<< "dim/codim/deg: " << dim(M) << "/" << codim(M) << "/" << degree(M) << endl;
G = mingens M;
<< "Betti mingens: " << betti G << endl;
-- Check for non-quadratic mingens:
H = select(first entries G, f -> sum degree(f) > 2);
if #H > 0 then (
<< "Non-quadratic generators:" << endl;
H / toString / print; << endl
);
<< endl;
In particular our Example 26, the Fano matroid, is the record beginning in line 497.
System information
These computations can in principle be performed on anyone’s laptop but for \(n \ge 5\) the degree computations may take a long time and the ideal computations may use a lot of memory. To compile the databases, we ran several tasks in parallel on a server with 88 cores and 1 TiB of memory, running Ubuntu 18.04.5 LTS, \(\verb|Macaulay2|\) v1.18, \(\verb|Julia|\) v1.7.0, HomotopyContinuation.jl v2.6.2, and Perl v5.26.
Colophon
Project page created: 18/12/2021
Project contributors: Tobias Boege, Sonja Petrović, Bernd Sturmfels
Corresponding author of this page: Tobias Boege, tobias.boege@mis.mpg.de
Software used: Macaulay2, Julia