# 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(Σ)
```

`LinearCovarianceModels.model`

— Function`model(W::MLDegreeWitness)`

Obtain the model corresponding to the `MLDegreeWitness`

`W`

.

`HomotopyContinuation.dim`

— Method`dim(M::LCModel)`

Returns the dimension of the model.

### 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.

`LinearCovarianceModels.generic_diagonal`

— Function`generic_diagonal(n::Integer, m::Integer)`

Generate a generic family of $n×n$ diagonal matrices living in an $m$-dimensional subspace.

`LinearCovarianceModels.toeplitz`

— Function`toeplitz(n::Integer)`

Returns a symmetric `n×n`

toeplitz matrix.

`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₄`

`LinearCovarianceModels.trees`

— Function`trees(n)`

Return all trees with `n`

leaves as a tuple (id, tree).

## 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.

`LinearCovarianceModels.MLDegreeWitness`

— Type`MLDegreeWitness`

Data structure holding an MLE model. This also holds a set of solutions for a generic instance, which we call a witness.

`LinearCovarianceModels.ml_degree`

— Function`ml_degree(W::MLDegreeWitness)`

Returns the ML degree.

`HomotopyContinuation.ModelKit.parameters`

— Method`parameters(W::MLDegreeWitness)`

Obtain the parameters of the `MLDegreeWitness`

`W`

.

`HomotopyContinuation.solutions`

— Method`solutions(W::MLDegreeWitness)`

Obtain the witness solutions corresponding to the `MLDegreeWitness`

`W`

with given parameters.

`LinearCovarianceModels.is_dual`

— Function`is_dual(W::MLDegreeWitness)`

Indicates whether `W`

is a witness for the dual MLE.

`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`

.

## 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.

`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`

.

`LinearCovarianceModels.covariance_matrix`

— Function`covariance_matrix(M::LCModel, θ)`

Compute the covariance matrix corresponding to the value of `θ`

and the given model `M`

.

`LinearCovarianceModels.logl`

— Function`logl(M::LCModel, θ, S::AbstractMatrix)`

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

`LinearCovarianceModels.gradient_logl`

— Function`gradient_logl(M::LCModel, θ, S::AbstractMatrix)`

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

`LinearCovarianceModels.hessian_logl`

— Function`hessian_logl(M::LCModel, θ, S::AbstractMatrix)`

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

`LinearCovarianceModels.classify_point`

— Function`classify_point(M::LCModel, θ, S::AbstractMatrix)`

Classify the critical point `θ`

of the log-likelihood function.

## Helper functions

`LinearCovarianceModels.sym_to_vec`

— Function`sym_to_vec(S)`

Converts a symmetric matrix `S`

to a vector by filling the vector with lower triangular part iterating columnwise.

`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`