Likelihood Geometry of the Squared Grassmannian
Abstract: We study projection determinantal point processes and their connection to the squared Grassmannian. We prove that the log-likelihood function of this statistical model has \((n−1)!/2\) critical points, all of which are real and positive, thereby settling a conjecture of Devriendt, Friedman, Reinke, and Sturmfels.
We provide here the code needed to reproduce the results in Example 4.1 and the Computation Experiments section of the paper. All code is written in \(\verb+Julia+\) and relies heavily on the numerical algebraic geometry package \(\verb+HomotopyContinuation.jl+\).
Once the requisite software is installed, files can be run in a terminal with \(\verb+ julia <filename>.jl+\).
Example 4.1
We claim that for \(d = 3, n = 6\) the parametric log-likelihood function has 17664 critical points and that typically 11904 of them are real.
This can be verified by running the code in realsolutions.jl.
The parameters \(d, n\) can be changed to count real solution solutions in other situations.
This file provides two functions. The function \(\verb+generate_start_solutions+\) uses the monodromy method implemented in \(\verb+HomotopyContinuation.jl+\) to compute critical points for random complex “data” vector. These start parameters and solutions are contained in a \(\verb+MonodromyResult+\) object. The function \(\verb+nreal_nmax_solutions+\) then uses a parameter homotopy to track the solutions as the parameters are moved to a vector of positive integers, representing real data; the function outputs the number of real solutions and local maxima.
We then compute the number of real solutions for a few random vectors with positive integer entries:
# Example: test a few random real vectors and see how many real critical points and
# local maxima the log-likelihood function has
d = 3
n = 6
println("d = ", d, ", n = ", n)
F, m = generate_start_solutions(n, d);
println("Number of complex critical points: ", nsolutions(m))
for _ in 1:5
nreal_sols, nlocal_max = nreal_nmax_solutions(F, m, rand(1:1000, binomial(n, d)))
println("Number of real critical points: ", nreal_sols,
", Number of local maxima: ", nlocal_max);
end
Sample output:
d = 3, n = 6
Number of complex critical points: 17664
Number of real critical points: 11904, Number of local maxima: 11904
Number of real critical points: 11904, Number of local maxima: 11904
Number of real critical points: 11904, Number of local maxima: 11904
Number of real critical points: 11904, Number of local maxima: 11904
Number of real critical points: 11904, Number of local maxima: 11904
We also compute a lower bound of 11904 on the number of sign vectors that can arise from Pluecker coordinates in \({\rm Gr}(3, 6)\). The code below will find all 11904 sign vectors with high probability, but it is possible that it might miss some. If this is the case, it can be run again or the number of trials can be increased.
using HomotopyContinuation, LinearAlgebra, Combinatorics
d = 3
n = 6
num_trials = 1000000
@var a[1:d, d+1:n]
A = [I a]
d_minors = [det(A[:, s]) for s in combinations(1:n, d)];
sgns = [sign.(evaluate.(d_minors, vec(a) => randn(d*(n-d)))) for _ in 1:num_trials];
# Checking the norm ensures that there are no zero entries in the sign vectors
println(length([v for v in unique(sgns) if norm(v) == sqrt(binomial(n, d))]))
Comptuational Experiments
The paper concludes with some computational experiments which show the complexity of computing maximum likelihood estimators for projection DPPs.
The runtimes are estimated using the \(\verb+Julia+\) package \(\verb+BenchmarkTools.jl+\).
The experiments can be recreated by running the code in runtime.jl for various values of \(n\).
(For these experiements, \(d\) is always 2.)
This file provides a function \(\verb+compute_MLE+\) which computes the maximum likelihood estimate given some data vector.
We then test the runtime of this function for small \(n\).
Sample output can be found below for \(n = 6\).
# Example: test the runtime for computing the MLE for small values of n
n = 6
println("n = ", n)
# Compute the runtime for random data
@benchmark compute_MLE(n, data) setup=(data=rand(1:1000, binomial(n, 2)))
Sample output:
n = 6
BenchmarkTools.Trial: 6 samples with 1 evaluation per sample.
Range (min … max): 838.224 ms … 1.019 s ┊ GC (min … max): 0.00% … 0.00%
Time (median): 954.079 ms ┊ GC (median): 0.00%
Time (mean ± σ): 953.124 ms ± 63.591 ms ┊ GC (mean ± σ): 0.00% ± 0.00%
▁ █ ▁ ▁ ▁
█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁█ ▁
838 ms Histogram: frequency by time 1.02 s <
Memory estimate: 46.40 MiB, allocs estimate: 1008546.