Computing the ML degree of the squared Grassmannian

In this file, we compute the ML degree of the squared Grassmannian. We do this computation using the birational parametrization of the (regular) Grassmannian, together with the \(2^{n-1}:1\) parametrization of the squared Grassmannian via the Grassmannian. 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)^2)_{I \in \binom{[n]}{d}}\) in the squared 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 squared Grassmannian \(\textrm{sGr}(d,n)\).

n = 5;
d = 3;

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);

The next cell defines the log-likelihood equation:

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

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

minors = [det(A[:,s]) for s in combinations(1:n, d)]
denominator = sum([minor^2 for minor in minors])
likelihood = sum(u[i] * log(minors[i]^2) for i in 1:binomial(n,d)) - sum(u[i] for i in 1:binomial(n,d))*log(denominator)

$$ u₂*log(a₃₋₁^2) + u₃*log(a₃₋₂^2) + u₄*log(a₂₋₁^2) + u₅*log(a₂₋₂^2) + u₇*log(a₁₋₁^2) + u₈*log(a₁₋₂^2) - log(1 + a₁₋₁^2 + a₁₋₂^2 + a₂₋₁^2 + a₂₋₂^2 + a₃₋₁^2 + a₃₋₂^2 + (a₁₋₁*a₂₋₂ - a₁₋₂*a₂₋₁)^2 + (-a₃₋₁*a₁₋₂ + a₃₋₂*a₁₋₁)^2 + (-a₃₋₁*a₂₋₂ + a₃₋₂*a₂₋₁)^2)*(u₁ + u₁₀ + u₂ + u₃ + u₄ + u₅ + u₆ + u₇ + u₈ + u₉) + log((a₁₋₁*a₂₋₂ - a₁₋₂*a₂₋₁)^2)*u₁₀ + log((-a₃₋₁*a₁₋₂ + a₃₋₂*a₁₋₁)^2)*u₉ + log((-a₃₋₁*a₂₋₂ + a₃₋₂*a₂₋₁)^2)*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 6
 6 variables: a₁₋₁, a₂₋₁, a₃₋₁, a₁₋₂, a₂₋₂, a₃₋₂
 10 parameters: u₁, u₂, u₃, u₄, u₅, u₆, u₇, u₈, u₉, u₁₀

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

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 manually compute a start pair. 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 squared Grassmannian.

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

Recall that because of the parameterization, the number of solutions is \(2^{n-1}\) times the ML degree of the squared Grassmannian. Thus the compute the ML degree of the squared Grassmannian, we must divide the output by \(2^{n-1}\).

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

Solutions found: 192    Time: 0:00:07
  tracked loops (queued):            1344 (0)
  solutions in current (last) loop:  0 (0)
  generated loops (no change):       7 (5)

nresults(solns) / 2^(n-1)

12.0