GMMParameterEstimation.jl Documentation

Authors: Haley Colgate Kottler, Julia Lindberg, Jose Israel Rodriguez

GMMParameterEstimation.jl is a package for estimating the parameters of Gaussian k-mixture models using the method of moments. It can potentially find the parameters for arbitrary k with known or unknown mixing coefficients. However, since the number of possible solutions to the polynomial system that determines the first dimension parameters and mixing coefficients for $k>4$ is unknown, for the unknown mixing coefficient case with $k>4$ failure of the package to find the parameters might occur if an insufficient number of solutions to the system were found

Examples

The following code snippet will use the given moments to return an estimate of the parameters using the method of moments with unknown mixing coefficients and general covariance matrices.

using GMMParameterEstimation
d = 3
k = 2
first_moments = [1.0, -0.67, 2.44, -4.34, 17.4, -46.16, 201.67]
diagonal_moments = [-0.28 2.11 -2.46 15.29 -31.77; 0.4 4.25 3.88 54.75 59.10]
off_diag_system = Dict{Vector{Int64}, Float64}([2, 1, 0] => 1.8506, [1, 0, 1] => -0.329, [2, 0, 1] => 0.0291, [0, 2, 1] => 1.5869, [1, 1, 0] => -1.374, [0, 1, 1] => -0.333)
error_check, (mixing_coefficients, means, covariances) = estimate_parameters(d, k, first_moments, diagonal_moments, off_diag_system)

Inputs:

  1. The number of dimensions d

  2. The number of mixture components k

  3. Optional: A vector of mixing coefficients w with length k

  4. A list of the first $3k+1$ moments (including moment 0) of the first dimension as first_moments

  5. A matrix where row i contains the first $2k+1$ moments (not including moment 0) of the ith dimension as diagonal_moments

  6. Optional: A dictionary mapping the index of a mixed dimensional moment as a list of integers to the corresponding moment off_diag_system (See mixedMomentSystem for clarrification on which moments to include.)

Outputs:

  1. An indicator of success in finding the parameters error_check, 0 means no error, 1 means an error in the first dimension system with either finding real solutions or non-negative mixing coefficients or positive covariance, 2 means an error in finding real solutions or positive covariances in a higher dimension, and 3 means the resulting covariance matrices weren't positive definite.

  2. A tuple of the parameters (mixing_coefficients, means, covariances)

$\\~\\$

The following code snippet will generate the exact moments necessary for parameter recovery from the given parameters.

using GMMParameterEstimation

d=3
k=2
diagonal = false

means = [0.83 0.24 -1.53; 0.22 0.04 -0.71]
covariances = [0.8828527552401668 0.27735188899130847 1.6710529671002674; 2.257873093006253 -1.644707016523332 -0.533030022431624;;; 0.27735188899130847 1.2623673813995742 3.5270452552353238; -1.644707016523332 2.577324062116896 -0.5049891831614162;;; 1.6710529671002674 3.5270452552353238 16.696895556824817; -0.533030022431624 -0.5049891831614162 1.7733508773418585]
mixing_coefficients = [.3, .7]

if diagonal
    true_first, true_diag = diagonalPerfectMoments(d, k, w, true_means, true_covariances)
else
    true_first, true_diag, true_others = generalPerfectMoments(d, k, w, true_means, true_covariances)
end

$\\~\\$

Parameter estimation

The main functionality of this package stems from

GMMParameterEstimation.estimate_parameters — Function
estimate_parameters(d::Integer, k::Integer, first::Vector{Float64}, second::Matrix{Float64})

Compute an estimate for the parameters of a d-dimensional Gaussian k-mixture model from the moments.

For the unknown mixing coefficient diagonal covariance matrix case, first should be a list of moments 0 through 3k for the first dimension, and second should be a matrix of moments 1 through 2k+1 for the remaining dimensions.

estimate_parameters(d::Integer, k::Integer, first::Vector{Float64}, second::Matrix{Float64}, last::Dict{Vector{Int64}, Float64}; method = "recursive")

Compute an estimate for the parameters of a d-dimensional Gaussian k-mixture model from the moments.

For the unknown mixing coefficient general covariance matrix case, first should be a list of moments 0 through 3k for the first dimension, second should be a matrix of moments 1 through 2k+1 for the remaining dimensions, and last should be a dictionary of the indices as lists of integers and the corresponding moments.

estimate_parameters(k::Integer, shared_cov::Matrix{Float64}, first::Vector{Float64}, second::Matrix{Float64}; method = "recursive")

Compute an estimate for the means of a Gaussian k-mixture model with equal mixing coefficients and known shared covariances from the moments.

The shared covariance matrix shared_cov will determine the dimension. Then first should be a list of moments 0 through k for the first dimension, second should be a matrix of moments $m_{je_1+e_i}$ for j in 0 to k-1 and i in 2 to d as a matrix where d is the dimension, i varies across rows, and j varies down columns.

which computes the parameter recovery using Algorithm 1 from Estimating Gaussian Mixtures Using Sparse Polynomial Moment Systems. Note that the unknown mixing coefficient cases with $k\in\{2,3,4\}$ load a set of generic moments and the corresponding solutions to the first 1-D polynomial system from sys1_k2.jld2, sys1_k3.jld2, or sys1_k4.jld2 for a slight speedup. If k is not specified, k=1 will be assumed, and the resulting polynomial system will be solved explicitly and directly.

In one dimension, for a random variable $X$ with density $f$ define the $i$th moment as $m_i=E[X^i]=\int x^if(x)dx$. For a Gaussian mixture model, this results in a polynomial in the parameters. For a sample $\{y_1,y_2,\dots,y_N\}$, define the $i$th sample moment as $\overline{m_i}=\frac{1}{N}\sum_{j=1}^N y_j^i$. The sample moments approach the true moments as $N\rightarrow\infty$, so by setting the polynomials equal to the empirical moments, we can then solve the polynomial system to recover the parameters.

For a multivariate random variable $X$ with density $f_X$ define the moments as $m_{i_1,\dots,i_n} = E[X_1^{i_1}\cdots X_n^{i_n}] = \int\cdots\int x_1^{i_1}\cdots x_n^{i_n}f_X(x_1,\dots,x_n)dx_1\cdots dx_n$ and the empirical moments as $\overline{m}_{i_1,\dots,i_n} = \frac{1}{N}\sum_{j=1}^Ny_{j_1}^{i_1}\cdots y_{j_n}^{i_n}$. And again, by setting the polynomials equal to the empirical moments, we can then solve the system of polynomials to recover the parameters. However, choosing which moments becomes more complicated. If we know the mixing coefficients, we can use the first $2k+1$ moments of each dimension to find the means and the diagonal entries of the covariance matrices. If we do not know the mixing coefficients, we need the first $3k$ moments of the first dimension to also find the mixing coefficients. See mixedMomentSystem for which moments to include to fill in the off-diagonals of the covariance matrices if needed.

On a standard laptop we have successfully recovered parameters with unknown mixing coefficients for $k\leq 4$ and known mixing coefficients for $k\leq 5$, with $d\leq 10^5$ for the diagonal covariance case and $d\leq 50$ for the general covariance case. Higher k values or higher d values have led to issues with running out of RAM.

$\\~\\$

One potential difficulty in estimating the mixing coefficients is the resulting dependence on higher moments in the first dimension. In sample data, if another dimension leads to more accurate moments, using that dimension to recover mixing coefficients and then proceeding can address this difficulty.

Generate and sample from Gaussian Mixture Models

Generation

Note that the entries of the resulting covariance matrices are generated from a normal distribution centered at 0 with variance 1.

$\\~\\$

GMMParameterEstimation.generateGaussians — Function
generateGaussians(d::Integer, k::Integer, diagonal::Bool)

Generate means and covariances for k Gaussians with dimension d.

diagonal should be true for spherical case, and false for general covariance matrices.

The parameters are returned as a tuple, with weights in a 1D vector, means as a k x d array, and variances as a k x d x d array if diagonal is false or as a list of Diagonal{Float64, Vector{Float64}} if diagonal is true to save memory. Note that each entry of each parameter is generated from a normal distribution centered at 0 with variance 1.

$\\~\\$

GMMParameterEstimation.getSample — Function
getSample(numb::Integer, w::Vector{Float64}, means::Matrix{Float64}, covariances::Vector)

Generate a Gaussian mixture model sample with numb entries, mixing coefficients w, means means, and covariances covariances.

getSample(numb::Integer, w::Vector{Float64}, means::Matrix{Float64}, covariances::Array{Float64, 3})

Generate a Gaussian mixture model sample with numb entries, mixing coefficients w, means means, and covariances covariances.

This relies on the Distributions package.

$\\~\\$

Computing moments

GMMParameterEstimation.sampleMoments — Function
sampleMoments(sample::Matrix{Float64}, k; diagonal = false)

Use the sample to compute the moments necessary for parameter estimation using method of moments with general covariance matrices and mixing coefficients.

Returns moments 0 to 3k for the first dimension, moments 1 through 2k+1 for the other dimensions as a matrix, and a dictionary with indices and moments for the off-diagonal system if diagonal is false.

GMMParameterEstimation.diagonalPerfectMoments — Function
diagonalPerfectMoments(d, k, w, true_means, true_covariances)

Use the given parameters to compute the exact moments necessary for parameter estimation with diagonal covariance matrices.

Returns moments 0 to 3k for the first dimension, and moments 1 through 2k+1 for the other dimensions as a matrix.

GMMParameterEstimation.generalPerfectMoments — Function
generalPerfectMoments(d, k, w, true_means, true_covariances; method = "low")

Use the given parameters to compute the exact moments necessary for parameter estimation with general covariance matrices.

Returns moments 0 to 3k for the first dimension, moments 1 through 2k+1 for the other dimensions as a matrix, and a dictionary with indices and moments for the off-diagonal system.

GMMParameterEstimation.equalMixCovarianceKnown_moments — Function
equalMixCovarianceKnown_moments(k, mean, shared_cov)

Use the given parameters to compute the exact moments necessary for parameter estimation with equal mixing coefficients and shared known covariances.

Returns moments 0 to k for the first dimension, and moments math m\_{je\_1+e\_i} for j in 0 to k-1 and i in 2 to d as a matrix where d is the dimension, i varies across rows, and j varies down columns.

equalMixCovarianceKnown_moments(k, sample)

Use the given parameters to compute the sample moments necessary for parameter estimation with equal mixing coefficients and shared known covariances.

Returns moments 0 to k for the first dimension, and moments m{je1+e_i} for j in 0 to k-1 and i in 2 to d as a matrix where d is the dimension, i varies across rows, and j varies down columns.

These expect parameters to be given with weights in a 1D vector, means as a k x d array, and covariances as a k x d x d array for general covariance matrices or as a list of diagonal matrices for diagonal covariance matrices.

$\\~\\$

Build the polynomial systems

GMMParameterEstimation.build1DSystem — Function
build1DSystem(k::Integer, m::Integer)

Build the polynomial system for a mixture of 1D Gaussians where 'm'-1 is the highest desired moment and the mixing coefficients are unknown.

build1DSystem(k::Integer, m::Integer, a::Union{Vector{Float64}, Vector{Variable}})

Build the polynomial system for a mixture of 1D Gaussians where 'm'-1 is the highest desired moment, and a is a vector of the mixing coefficients.

This uses the usual recursive formula for moments of a univariate Gaussian in terms of the mean and variance, and then takes a convex combination with either variable mixing coefficients or the provided mixing coefficients.

$\\~\\$

The final step in our method of moments parameter recovery for non-diagonal covariance matrices is building and solving a system of $N:=\frac{k}{2}(d^2-d)$ linear equations in the same number of unknowns to fill in the off diagonal.

We present two options for moment selection. First, Lindberg et al., demonstrated that for fixed $i,j\in[d]$, $i\neq j$, the set of moment equations $\{m_{te_i+e_j}\}_{t=1}^k$ are linear in $\sigma_{\ell i j}$ for $\ell\in[k]$ and leads to a linear system that generically has a unique solution. Due to the symmetry of covariance matrices ($\sigma_{\ell i j}=\sigma_{\ell j i}$) we add the restriction that $i<j$. This is implemented and can be selected by setting style="k".

Second, further relying on the symmetry of the covariance matrices and their corresponding polynomials, we present a lower order system given by $\begin{cases} \{m_{te_i+e_j}, m_{e_i+te_j}\}_{t=1}^{\frac{k}{2}}\cup\{m_{te_i+e_j}\}_{t=1}^{\frac{k}{2}+1} & \text{ if k is even}\\ \{m_{te_i+e_j}, m_{e_i+te_j}\}_{t=1}^{\frac{k+1}{2}} & \text{ if k is odd}\end{cases}.$

We again assume $i,j\in[d]$ with $i<j$. This is the default setting and can be selected by setting style="low".

We also present two methods for generating the moment polynomials. First, we use tensor moments. Referring to Pereira et al. for a closed form method of generating the necessary moment polynomials, we generate the linear system using the already computed mixing coefficients, means, and diagonals of the covariances, and return it as a dictionary of index=>polynomial pairs that can then be matched with the corresponding moments. This can be accessed by selecting method="tensor" and uses the lower order system due to difficulties with higher order tensor moments.

Second, we use a recursive method, which is the default because it is faster, and can be selected via method="recursive".

Support Functions

GMMParameterEstimation.build1DSystem — Method
build1DSystem(k::Integer, m::Integer, a::Union{Vector{Float64}, Vector{Variable}})

Build the polynomial system for a mixture of 1D Gaussians where 'm'-1 is the highest desired moment, and a is a vector of the mixing coefficients.

GMMParameterEstimation.build1DSystem — Method
build1DSystem(k::Integer, m::Integer)

Build the polynomial system for a mixture of 1D Gaussians where 'm'-1 is the highest desired moment and the mixing coefficients are unknown.

GMMParameterEstimation.checkInputs — Method
checkInputs(d::Integer, k::Integer, first::Vector{Float64}, second::Matrix{Float64}, last::Dict{Vector{Int64}, Float64}, method)

Returns true if the inputs are the right format for estimate_parameters and an error otherwise.

GMMParameterEstimation.diagonalPerfectMoments — Method
diagonalPerfectMoments(d, k, w, true_means, true_covariances)

Use the given parameters to compute the exact moments necessary for parameter estimation with diagonal covariance matrices.

Returns moments 0 to 3k for the first dimension, and moments 1 through 2k+1 for the other dimensions as a matrix.

GMMParameterEstimation.equalMixCovarianceKnown_moments — Method
equalMixCovarianceKnown_moments(k, mean, shared_cov)

Use the given parameters to compute the exact moments necessary for parameter estimation with equal mixing coefficients and shared known covariances.

Returns moments 0 to k for the first dimension, and moments math m\_{je\_1+e\_i} for j in 0 to k-1 and i in 2 to d as a matrix where d is the dimension, i varies across rows, and j varies down columns.

GMMParameterEstimation.equalMixCovarianceKnown_moments — Method
equalMixCovarianceKnown_moments(k, sample)

Use the given parameters to compute the sample moments necessary for parameter estimation with equal mixing coefficients and shared known covariances.

Returns moments 0 to k for the first dimension, and moments m{je1+e_i} for j in 0 to k-1 and i in 2 to d as a matrix where d is the dimension, i varies across rows, and j varies down columns.

GMMParameterEstimation.estimate_parameters — Method
estimate_parameters(d::Integer, k::Integer, first::Vector{Float64}, second::Matrix{Float64}, last::Dict{Vector{Int64}, Float64}; method = "recursive")

Compute an estimate for the parameters of a d-dimensional Gaussian k-mixture model from the moments.

For the unknown mixing coefficient general covariance matrix case, first should be a list of moments 0 through 3k for the first dimension, second should be a matrix of moments 1 through 2k+1 for the remaining dimensions, and last should be a dictionary of the indices as lists of integers and the corresponding moments.

GMMParameterEstimation.estimate_parameters — Method
estimate_parameters(d::Integer, k::Integer, first::Vector{Float64}, second::Matrix{Float64})

Compute an estimate for the parameters of a d-dimensional Gaussian k-mixture model from the moments.

For the unknown mixing coefficient diagonal covariance matrix case, first should be a list of moments 0 through 3k for the first dimension, and second should be a matrix of moments 1 through 2k+1 for the remaining dimensions.

GMMParameterEstimation.estimate_parameters — Method
estimate_parameters(k::Integer, shared_cov::Matrix{Float64}, first::Vector{Float64}, second::Matrix{Float64}; method = "recursive")

Compute an estimate for the means of a Gaussian k-mixture model with equal mixing coefficients and known shared covariances from the moments.

The shared covariance matrix shared_cov will determine the dimension. Then first should be a list of moments 0 through k for the first dimension, second should be a matrix of moments $m_{je_1+e_i}$ for j in 0 to k-1 and i in 2 to d as a matrix where d is the dimension, i varies across rows, and j varies down columns.

GMMParameterEstimation.estimate_parameters_weights_known — Method
estimate_parameters_weights_known(d::Integer, k::Integer, w::Array{Float64}, first::Matrix{Float64}, last::Dict{Vector{Int64}, Float64}; method = "recursive")

Compute an estimate for the parameters of a d-dimensional Gaussian k-mixture model from the moments.

For the known mixing coefficient general covariance matrix case, w should be a vector of the mixing coefficients first should be a list of moments 0 through 3k for the first dimension, second should be a matrix of moments 1 through 2k+1 for the remaining dimensions, and last should be a dictionary of the indices as lists of integers and the corresponding moments.

GMMParameterEstimation.estimate_parameters_weights_known — Method
estimate_parameters_weights_known(d::Integer, k::Integer, w::Array{Float64}, moments::Matrix{Float64})

Compute an estimate for the parameters of a d-dimensional Gaussian k-mixture model from the moments.

For the known mixing coefficient diagonal covariance matrix case, w should be a vector of the mixing coefficients, and moments should be a matrix of moments 1 through 2k+1 for all dimensions.

GMMParameterEstimation.generalPerfectMoments — Method
generalPerfectMoments(d, k, w, true_means, true_covariances; method = "low")

Use the given parameters to compute the exact moments necessary for parameter estimation with general covariance matrices.

Returns moments 0 to 3k for the first dimension, moments 1 through 2k+1 for the other dimensions as a matrix, and a dictionary with indices and moments for the off-diagonal system.

GMMParameterEstimation.generateGaussians — Method
generateGaussians(d::Integer, k::Integer, diagonal::Bool)

Generate means and covariances for k Gaussians with dimension d.

diagonal should be true for spherical case, and false for general covariance matrices.

GMMParameterEstimation.getSample — Method
getSample(numb::Integer, w::Vector{Float64}, means::Matrix{Float64}, covariances::Array{Float64, 3})

Generate a Gaussian mixture model sample with numb entries, mixing coefficients w, means means, and covariances covariances.

GMMParameterEstimation.getSample — Method
getSample(numb::Integer, w::Vector{Float64}, means::Matrix{Float64}, covariances::Vector)

Generate a Gaussian mixture model sample with numb entries, mixing coefficients w, means means, and covariances covariances.

GMMParameterEstimation.kOrder_offDiagonalSolve — Method
kOrder_offDiagonalSolve(d, k, mixing_coefficients, means, covariances, last)

Use a recursive system to build and solve a system for the off-diagonal covariance entries using moments $m_{te_i+e_j}$ for t in [k] and i<j

GMMParameterEstimation.lowOrder_offDiagonalSolve — Method
lowOrder_offDiagonalSolve(d, k, mixing_coefficients, means, covariances, last)

Use a recursive system to build and solve a system for the off-diagonal covariance entries using minimal order moments $m_{te_i+e_j}$ and $m_{e_i+te_j}$

GMMParameterEstimation.sampleMoments — Method
sampleMoments(sample::Matrix{Float64}, k; diagonal = false)

Use the sample to compute the moments necessary for parameter estimation using method of moments with general covariance matrices and mixing coefficients.

Returns moments 0 to 3k for the first dimension, moments 1 through 2k+1 for the other dimensions as a matrix, and a dictionary with indices and moments for the off-diagonal system if diagonal is false.

GMMParameterEstimation.selectSol — Method
selectSol(k::Integer, solution::Result, polynomial::Expression, moment::Number)

Select a k mixture solution from solution accounting for polynomial and moment.

Sort out a k mixture statistically significant solutions from solution, and return the one closest to moment when polynomial is evaluated at those values.

GMMParameterEstimation.tensorMixedMomentSystem — Method
tensorMixedMomentSystem(d, k, mixing, ms, vs)

Build a linear system for finding the off-diagonal covariances entries using tensor moments.

For a d dimensional Gaussian k-mixture model with mixing coefficients mixing, means ms, and covariances vs where the diagonal entries have been filled in and the off diagonals are variables.

Index