Edit on GitHub

Introduction

LinearCovarianceModels.jl is a package for computing Maximum Likelihood degrees of linear covariance models using numerical nonlinear algebra. In particular HomotopyContinuation.jl.

Introduction by Example

# Create a linear covariance model
julia> Σ = toeplitz(3)
3-dimensional LCModel:
 θ₁  θ₂  θ₃
 θ₂  θ₁  θ₂
 θ₃  θ₂  θ₁

# Compute a witness for the ML degree
julia> W = ml_degree_witness(Σ)
 MLDegreeWitness:
 • ML degree → 3
 • model dimension → 3
 • dual → false

# We offer the option to numerically verify the ML Degree
julia> verify(W)
 Compute additional witnesses for completeness...
 Found 10 additional witnesses
 Found 10 additional witnesses
 Compute trace...
 Norm of trace: 2.6521474798326718e-12
 true

# Consider the sample covariance matrix S
julia> S = [4/5 -9/5 -1/25
            -9/5 79/16 25/24
            -1/25 25/24 17/16];

# We use the ML degree witness set W to compute all critical points of the MLE
# problem.
julia> critical_points(W, S)
3-element Array{Tuple{Array{Float64,1},Float64,Symbol},1}:
 ([2.39038, -0.286009, 0.949965], -5.421751313919751, :local_maximum)
 ([2.52783, -0.215929, -1.45229], -5.346601549034418, :global_maximum)
 ([2.28596, -0.256394, 0.422321], -5.424161999175718, :saddle_point)  

# If we are just interested in the MLE, there is also a shorthand.
julia> mle(W, S)
3-element Array{Float64,1}:
  2.527832268219689  
 -0.21592947057775033
 -1.4522862659134732

Linear Covariance Models

The linear covariance models are wrapped in the LCModel type:

LinearCovarianceModels.LCModel — Type
LCModel(Σ::Matrix)

Create a linear covariance model from the parameterization Σ. This uses as input a matrix of polynomials created by the @var macro from HomotopyContinuation.jl.

Example

using HomotopyContinuation # load polynomials package

# use HomotopyContinuation to create variables θ₁, θ₂, θ₃.
@var θ[1:3]

# create our model as matrix of DynamicPolynomials
Σ = [θ[1] θ[2] θ[3]; θ[2] θ[1] θ[2]; θ[3] θ[2] θ[1]]

# create model
model = LCModel(Σ)
source

Default models

The following linear covariance models are provided by default

LinearCovarianceModels.generic_subspace — Function
generic_subspace(n::Integer, m::Integer); pos_def::Bool=true)

Generate a generic family of symmetric $n×n$ matrices living in an $m$-dimensional subspace. If pos_def is true then positive definite matrices are used as a basis.

source
LinearCovarianceModels.tree — Function
tree(n, id::String)

Get the covariance matrix corresponding to the tree with the given id on n leaves. Returns nothing if the tree was not found.

Example

julia> tree(4, "{{1, 2}, {3, 4}}") 4×4 Array{PolyVar{true},2}: t₁ t₅ t₇ t₇ t₅ t₂ t₇ t₇ t₇ t₇ t₃ t₆ t₇ t₇ t₆ t₄

source

ML Degree

LinearCovarianceModels.ml_degree_witness — Function
ml_degree_witness(Σ::LCModel; ml_degree=nothing, max_tries=5, dual=false, compile = false)

Compute a MLDegreeWitness for a given model Σ. If the ML degree is already known it can be provided to stop the computations early. The stopping criterion is based on a heuristic, max_tries indicates how many different parameters are tried a most until an agreement is found.

source
LinearCovarianceModels.verify — Function
verify(W::MLDegreeWitness; trace_tol=1e-5, options...)

Tries to verify that the computed ML degree witness is complete, i.e., that the ML degree is correct. This uses the verify_solution_completeness of HomotopyContinuation.jl. All caveats mentioned there apply. The options are also passed to verify_solution_completeness.

source

Compute MLE for specific instances

LinearCovarianceModels.mle — Function
mle(W::MLDegreeWitness, S::AbstractMatrix; only_positive_definite=true, only_positive=false)

Compute the MLE for the matrix S using the MLDegreeWitness W. Returns the parameters for the MLE covariance matrix or nothing if no solution was found satisfying the constraints (see options below).

Options

  • only_positive_definite: controls whether only positive definite

covariance matrices should be considered.

  • only_positive: controls whether only (entrywise) positive covariance matrices

should be considered.

source
LinearCovarianceModels.critical_points — Function
critical_points(W::MLDegreeWitness, S::AbstractMatrix;
        only_positive_definite=true, only_non_negative=false,
        options...)

Compute all critical points to the MLE problem of W for the given sample covariance matrix S. If only_positive_definite is true only positive definite solutions are considered. If only_non_negative is true only non-negative solutions are considered. The options are argument passed to the solve routine from HomotopyContinuation.jl.

source
LinearCovarianceModels.logl — Function
logl(M::LCModel, θ, S::AbstractMatrix)

Evaluate the log-likelihood $log(det(Σ⁻¹)) - tr(SΣ⁻¹)$ of the MLE problem.

source

Helper functions

LinearCovarianceModels.vec_to_sym — Function
vec_to_sym(v)

Converts a vector v to a symmetrix matrix by filling the lower triangular part columnwise.

Example

julia> v = [1,2,3, 4, 5, 6]; julia> vec_to_sym(v) 3×3 Array{Int64,2}: 1 2 3 2 4 5 3 5 6

source