Computing the ML degree of the positive Grassmannian

In this file, we compute the ML degree of the positive Grassmannian. We use a birational parametrization of the Grassmannian where a point in the Grassmannian \(\textrm{Gr}(d, n)\) is represented by a \(d \times n\) matrix \(A\) where the first \(d \times d\) block is the identity matrix and the remaining entries are filled with variables \(a_{ij}\). The matrix \(A\) parametrizes the point \((\det(A_I))_{I \in \binom{[n]}{d}}\) in the Grassmannian. For example, if \(n = 5, d = 2\), we have the matrix

\[\begin{split}A = \begin{pmatrix} 1 & 0 & a_{11} & a_{12} & a_{13} \\ 0 & 1 & a_{21} & a_{22} & a_{23} \end{pmatrix}.\end{split}\]

using HomotopyContinuation
using LinearAlgebra
using Combinatorics

When \(n\) and \(d\) are specified below, the remainder of the code computes the ML degree of the positive Grassmannian \(\textrm{Gr}_{\geq 0}(d,n)\).

n = 6;
d = 2;

Setting up the system

Here we define our variable matrix \(A\) as described above and our data vector \(u\).

@var a[1:d, 1:n - d]
@var u[1:binomial(n, d)]

A = [UniformScaling(1) a];
a = vec(a);

minors = [det(A[:,s]) for s in combinations(1:n, d)];
denominator = sum(minors);

Here we define the log-likelihood equation

\[L_u = \left (\sum_{I \in {[n] \choose d}} u_I \log(\det(A_I)) \right ) - \left( \sum_{I \in {[n] \choose d}} u_I \right) \log\left (\sum_{I \in {[n] \choose d}} \det(A_I) \right ),\]

where \(A_I\) is the submatrix of \(A\) given by selecting only the columns indexed by \(I\).

likelihood = sum(u[i] * log(minors[i]) for i in 1:binomial(n,d)) - sum(u[i] for i in 1:binomial(n,d))*log(denominator)

$$ u₂*log(a₂₋₁) + u₃*log(a₂₋₂) + u₄*log(a₂₋₃) + u₅*log(a₂₋₄) + u₆*log(-a₁₋₁) + u₇*log(-a₁₋₂) + u₈*log(-a₁₋₃) + u₉*log(-a₁₋₄) + log(a₁₋₁*a₂₋₃ - a₁₋₃*a₂₋₁)*u₁₁ + log(a₁₋₁*a₂₋₂ - a₁₋₂*a₂₋₁)*u₁₀ + log(a₁₋₂*a₂₋₃ - a₁₋₃*a₂₋₂)*u₁₃ + log(a₁₋₂*a₂₋₄ - a₁₋₄*a₂₋₂)*u₁₄ + log(a₁₋₁*a₂₋₄ - a₁₋₄*a₂₋₁)*u₁₂ + log(a₁₋₃*a₂₋₄ - a₁₋₄*a₂₋₃)*u₁₅ - log(1 - a₁₋₁ - a₁₋₂ - a₁₋₃ - a₁₋₄ + a₂₋₁ + a₂₋₂ + a₂₋₃ + a₂₋₄ + a₁₋₁*a₂₋₂ + a₁₋₁*a₂₋₃ + a₁₋₁*a₂₋₄ - a₁₋₂*a₂₋₁ + a₁₋₂*a₂₋₃ + a₁₋₂*a₂₋₄ - a₁₋₃*a₂₋₁ - a₁₋₃*a₂₋₂ + a₁₋₃*a₂₋₄ - a₁₋₄*a₂₋₁ - a₁₋₄*a₂₋₂ - a₁₋₄*a₂₋₃)*(u₁ + u₁₀ + u₁₁ + u₁₂ + u₁₃ + u₁₄ + u₁₅ + u₂ + u₃ + u₄ + u₅ + u₆ + u₇ + u₈ + u₉) $$

We want to find critical points of this function, i.e. points where all the derivatives vanish. Therefore we take the derivatives and put them in a \(\verb|HomotopyContinuation.jl|\) system with variables \(a\) and parameters \(u\).

F = System(differentiate(likelihood, a); parameters = u)

System of length 8
 8 variables: a₁₋₁, a₂₋₁, a₁₋₂, a₂₋₂, a₁₋₃, a₂₋₃, a₁₋₄, a₂₋₄
 15 parameters: u₁, u₂, u₃, u₄, u₅, u₆, u₇, u₈, u₉, u₁₀, u₁₁, u₁₂, u₁₃, u₁₄, u₁₅

 u₆/a₁₋₁ + a₂₋₃*u₁₁/(a₁₋₁*a₂₋₃ - a₁₋₃*a₂₋₁) + a₂₋₂*u₁₀/(a₁₋₁*a₂₋₂ - a₁₋₂*a₂₋₁) + a₂₋₄*u₁₂/(a₁₋₁*a₂₋₄ - a₁₋₄*a₂₋₁) - (u₁ + u₁₀ + u₁₁ + u₁₂ + u₁₃ + u₁₄ + u₁₅ + u₂ + u₃ + u₄ + u₅ + u₆ + u₇ + u₈ + u₉)*(-1 + a₂₋₂ + a₂₋₃ + a₂₋₄)/(1 - a₁₋₁ - a₁₋₂ - a₁₋₃ - a₁₋₄ + a₂₋₁ + a₂₋₂ + a₂₋₃ + a₂₋₄ + a₁₋₁*a₂₋₂ + a₁₋₁*a₂₋₃ + a₁₋₁*a₂₋₄ - a₁₋₂*a₂₋₁ + a₁₋₂*a₂₋₃ + a₁₋₂*a₂₋₄ - a₁₋₃*a₂₋₁ - a₁₋₃*a₂₋₂ + a₁₋₃*a₂₋₄ - a₁₋₄*a₂₋₁ - a₁₋₄*a₂₋₂ - a₁₋₄*a₂₋₃)
 u₂/a₂₋₁ - a₁₋₃*u₁₁/(a₁₋₁*a₂₋₃ - a₁₋₃*a₂₋₁) - a₁₋₂*u₁₀/(a₁₋₁*a₂₋₂ - a₁₋₂*a₂₋₁) - a₁₋₄*u₁₂/(a₁₋₁*a₂₋₄ - a₁₋₄*a₂₋₁) - (u₁ + u₁₀ + u₁₁ + u₁₂ + u₁₃ + u₁₄ + u₁₅ + u₂ + u₃ + u₄ + u₅ + u₆ + u₇ + u₈ + u₉)*(1 - a₁₋₂ - a₁₋₃ - a₁₋₄)/(1 - a₁₋₁ - a₁₋₂ - a₁₋₃ - a₁₋₄ + a₂₋₁ + a₂₋₂ + a₂₋₃ + a₂₋₄ + a₁₋₁*a₂₋₂ + a₁₋₁*a₂₋₃ + a₁₋₁*a₂₋₄ - a₁₋₂*a₂₋₁ + a₁₋₂*a₂₋₃ + a₁₋₂*a₂₋₄ - a₁₋₃*a₂₋₁ - a₁₋₃*a₂₋₂ + a₁₋₃*a₂₋₄ - a₁₋₄*a₂₋₁ - a₁₋₄*a₂₋₂ - a₁₋₄*a₂₋₃)
 u₇/a₁₋₂ - a₂₋₁*u₁₀/(a₁₋₁*a₂₋₂ - a₁₋₂*a₂₋₁) + a₂₋₃*u₁₃/(a₁₋₂*a₂₋₃ - a₁₋₃*a₂₋₂) + a₂₋₄*u₁₄/(a₁₋₂*a₂₋₄ - a₁₋₄*a₂₋₂) - (u₁ + u₁₀ + u₁₁ + u₁₂ + u₁₃ + u₁₄ + u₁₅ + u₂ + u₃ + u₄ + u₅ + u₆ + u₇ + u₈ + u₉)*(-1 - a₂₋₁ + a₂₋₃ + a₂₋₄)/(1 - a₁₋₁ - a₁₋₂ - a₁₋₃ - a₁₋₄ + a₂₋₁ + a₂₋₂ + a₂₋₃ + a₂₋₄ + a₁₋₁*a₂₋₂ + a₁₋₁*a₂₋₃ + a₁₋₁*a₂₋₄ - a₁₋₂*a₂₋₁ + a₁₋₂*a₂₋₃ + a₁₋₂*a₂₋₄ - a₁₋₃*a₂₋₁ - a₁₋₃*a₂₋₂ + a₁₋₃*a₂₋₄ - a₁₋₄*a₂₋₁ - a₁₋₄*a₂₋₂ - a₁₋₄*a₂₋₃)
 u₃/a₂₋₂ + a₁₋₁*u₁₀/(a₁₋₁*a₂₋₂ - a₁₋₂*a₂₋₁) - a₁₋₃*u₁₃/(a₁₋₂*a₂₋₃ - a₁₋₃*a₂₋₂) - a₁₋₄*u₁₄/(a₁₋₂*a₂₋₄ - a₁₋₄*a₂₋₂) - (u₁ + u₁₀ + u₁₁ + u₁₂ + u₁₃ + u₁₄ + u₁₅ + u₂ + u₃ + u₄ + u₅ + u₆ + u₇ + u₈ + u₉)*(1 + a₁₋₁ - a₁₋₃ - a₁₋₄)/(1 - a₁₋₁ - a₁₋₂ - a₁₋₃ - a₁₋₄ + a₂₋₁ + a₂₋₂ + a₂₋₃ + a₂₋₄ + a₁₋₁*a₂₋₂ + a₁₋₁*a₂₋₃ + a₁₋₁*a₂₋₄ - a₁₋₂*a₂₋₁ + a₁₋₂*a₂₋₃ + a₁₋₂*a₂₋₄ - a₁₋₃*a₂₋₁ - a₁₋₃*a₂₋₂ + a₁₋₃*a₂₋₄ - a₁₋₄*a₂₋₁ - a₁₋₄*a₂₋₂ - a₁₋₄*a₂₋₃)
 u₈/a₁₋₃ - a₂₋₁*u₁₁/(a₁₋₁*a₂₋₃ - a₁₋₃*a₂₋₁) - a₂₋₂*u₁₃/(a₁₋₂*a₂₋₃ - a₁₋₃*a₂₋₂) + a₂₋₄*u₁₅/(a₁₋₃*a₂₋₄ - a₁₋₄*a₂₋₃) - (u₁ + u₁₀ + u₁₁ + u₁₂ + u₁₃ + u₁₄ + u₁₅ + u₂ + u₃ + u₄ + u₅ + u₆ + u₇ + u₈ + u₉)*(-1 - a₂₋₁ - a₂₋₂ + a₂₋₄)/(1 - a₁₋₁ - a₁₋₂ - a₁₋₃ - a₁₋₄ + a₂₋₁ + a₂₋₂ + a₂₋₃ + a₂₋₄ + a₁₋₁*a₂₋₂ + a₁₋₁*a₂₋₃ + a₁₋₁*a₂₋₄ - a₁₋₂*a₂₋₁ + a₁₋₂*a₂₋₃ + a₁₋₂*a₂₋₄ - a₁₋₃*a₂₋₁ - a₁₋₃*a₂₋₂ + a₁₋₃*a₂₋₄ - a₁₋₄*a₂₋₁ - a₁₋₄*a₂₋₂ - a₁₋₄*a₂₋₃)
 u₄/a₂₋₃ + a₁₋₁*u₁₁/(a₁₋₁*a₂₋₃ - a₁₋₃*a₂₋₁) + a₁₋₂*u₁₃/(a₁₋₂*a₂₋₃ - a₁₋₃*a₂₋₂) - a₁₋₄*u₁₅/(a₁₋₃*a₂₋₄ - a₁₋₄*a₂₋₃) - (u₁ + u₁₀ + u₁₁ + u₁₂ + u₁₃ + u₁₄ + u₁₅ + u₂ + u₃ + u₄ + u₅ + u₆ + u₇ + u₈ + u₉)*(1 + a₁₋₁ + a₁₋₂ - a₁₋₄)/(1 - a₁₋₁ - a₁₋₂ - a₁₋₃ - a₁₋₄ + a₂₋₁ + a₂₋₂ + a₂₋₃ + a₂₋₄ + a₁₋₁*a₂₋₂ + a₁₋₁*a₂₋₃ + a₁₋₁*a₂₋₄ - a₁₋₂*a₂₋₁ + a₁₋₂*a₂₋₃ + a₁₋₂*a₂₋₄ - a₁₋₃*a₂₋₁ - a₁₋₃*a₂₋₂ + a₁₋₃*a₂₋₄ - a₁₋₄*a₂₋₁ - a₁₋₄*a₂₋₂ - a₁₋₄*a₂₋₃)
 u₉/a₁₋₄ - a₂₋₂*u₁₄/(a₁₋₂*a₂₋₄ - a₁₋₄*a₂₋₂) - a₂₋₁*u₁₂/(a₁₋₁*a₂₋₄ - a₁₋₄*a₂₋₁) - a₂₋₃*u₁₅/(a₁₋₃*a₂₋₄ - a₁₋₄*a₂₋₃) - (u₁ + u₁₀ + u₁₁ + u₁₂ + u₁₃ + u₁₄ + u₁₅ + u₂ + u₃ + u₄ + u₅ + u₆ + u₇ + u₈ + u₉)*(-1 - a₂₋₁ - a₂₋₂ - a₂₋₃)/(1 - a₁₋₁ - a₁₋₂ - a₁₋₃ - a₁₋₄ + a₂₋₁ + a₂₋₂ + a₂₋₃ + a₂₋₄ + a₁₋₁*a₂₋₂ + a₁₋₁*a₂₋₃ + a₁₋₁*a₂₋₄ - a₁₋₂*a₂₋₁ + a₁₋₂*a₂₋₃ + a₁₋₂*a₂₋₄ - a₁₋₃*a₂₋₁ - a₁₋₃*a₂₋₂ + a₁₋₃*a₂₋₄ - a₁₋₄*a₂₋₁ - a₁₋₄*a₂₋₂ - a₁₋₄*a₂₋₃)
 u₅/a₂₋₄ + a₁₋₂*u₁₄/(a₁₋₂*a₂₋₄ - a₁₋₄*a₂₋₂) + a₁₋₁*u₁₂/(a₁₋₁*a₂₋₄ - a₁₋₄*a₂₋₁) + a₁₋₃*u₁₅/(a₁₋₃*a₂₋₄ - a₁₋₄*a₂₋₃) - (u₁ + u₁₀ + u₁₁ + u₁₂ + u₁₃ + u₁₄ + u₁₅ + u₂ + u₃ + u₄ + u₅ + u₆ + u₇ + u₈ + u₉)*(1 + a₁₋₁ + a₁₋₂ + a₁₋₃)/(1 - a₁₋₁ - a₁₋₂ - a₁₋₃ - a₁₋₄ + a₂₋₁ + a₂₋₂ + a₂₋₃ + a₂₋₄ + a₁₋₁*a₂₋₂ + a₁₋₁*a₂₋₃ + a₁₋₁*a₂₋₄ - a₁₋₂*a₂₋₁ + a₁₋₂*a₂₋₃ + a₁₋₂*a₂₋₄ - a₁₋₃*a₂₋₁ - a₁₋₃*a₂₋₂ + a₁₋₃*a₂₋₄ - a₁₋₄*a₂₋₁ - a₁₋₄*a₂₋₂ - a₁₋₄*a₂₋₃)

Solving the system

We then compute a start system and then solve the system using monodromy. Depending on the system, it may take more tries than the default to compute a start pair, so we find a start pair manually. We can use threading here to speed up the monodromy which can be very time consuming. We then certify the solutions to check that we got distinct solutions. We now have a lower bound on the ML degree of this partiticular positive Grassmannian.

start_pair = find_start_pair(F; max_tries = 10000, atol = 0.0, rtol = 1e-12);

solns = monodromy_solve(F, start_pair[1], start_pair[2], threading=true)

MonodromyResult
===============
• return_code → :heuristic_stop
• 156 solutions
• 1092 tracked loops
• random_seed → 0xa039df43

certify(F, solns)

CertificationResult
===================
• 156 solution candidates given
• 156 certified solution intervals (0 real, 156 complex)
• 156 distinct certified solution intervals (0 real, 156 complex)