Expectation-Maximization

Expectation-Maximization (EM) algorithm is an iterative method for finding the maximum likelihood and maximum a posteriori estimates of parameters in models that typically depend on hidden variables.

While serving as a clustering technique, EM is also used in non-linear dimensionality reduction, missing value problems, and other areas.

Details

Given a set \(X\) of \(n\) feature vectors \(x_1 = (x_{11}, \ldots, x_{1p}), \ldots, x_n = (x_{n1}, \ldots, x_{np})\) of dimension \(p\), the problem is to find a maximum-likelihood estimate of the parameters of the underlying distribution when the data is incomplete or has missing values.

Expectation-Maximization (EM) Algorithm in the General Form

Let \(X\) be the observed data which has log-likelihood \(l(\theta; X)\) depending on the parameters \(\theta\). Let \(X^m\) be the latent or missing data, so that \(T=(X, X^m)\) is the complete data with log-likelihood \(l_0(\theta; X)\). The algorithm for solving the problem in its general form is the following EM algorithm ([Dempster77], [Hastie2009]):

  1. Choose initial values of the parameters \(\theta^{(0)}\).

  2. Expectation step: in the \(j\)-th step, compute \(Q(\theta', \theta^{(j)}) = E (l_0(\theta'; T) | X, \theta^{(j)})\) as a function of the dummy argument \(\theta'\).

  3. Maximization step: in the \(j\)-th step, calculate the new estimate \(\theta^{(j+1)}\) by maximizing \(Q(\theta', \theta^{(j)})\) over \(\theta'\).

  4. Repeat steps 2 and 3 until convergence.

EM algorithm for the Gaussian Mixture Model

Gaussian Mixture Model (GMM) is a mixture of k p-dimensional multivariate Gaussian distributions represented as

\[F\left(x{|\alpha }_{1},\ldots,{\alpha }_{k};{\theta }_{1},\ldots,{\theta }_{k}\right)=\sum _{i=1}^{k}{\alpha }_{i}{\int }_{-\infty }^{x}pd\left(y|{\theta }_{i}\right),\]

where \(\sum _{i=1}^{k}{\alpha_i} = 1\) and \(\alpha_i \geq 0\).

The \(pd(x|\theta_i)\) is the probability density function with parameters \(\theta_i = (m_i, \Sigma_i)\), where \(m_i\) the vector of means, and \(\Sigma_i\) is the variance-covariance matrix. The probability density function for a \(p\)-dimensional multivariate Gaussian distribution is defined as follows:

\[pd\left(x|{\theta }_{i}\right)=\frac{\mathrm{exp}\left(-\frac{1}{2}{\left(x-{m}_{i}\right)}^{T}\sum _{i}^{-1}\left(x-{m}_{i}\right)\right)}{\sqrt{{\left(2\pi \right)}^{p}|\sum _{i}|}}.\]

Let \(x_{ij} = I\{x_i \text{belongs to j mixture component}\}\) be the indicator function and \(\theta = (\alpha_1, \ldots, \alpha_k; \theta_1, \ldots, \theta_k)\).

Computation

The EM algorithm for GMM includes the following steps:

Define the weights as follows:

\[{w}_{ij}=\frac{pd\left({x}_{i}|{z}_{ij},{\theta }_{j}\right){\alpha }_{j}}{\sum _{r=1}^{k}pd\left({x}_{i}|{z}_{ir},{\theta }_{r}\right){\alpha }_{r}}\]

for \(i = 1, \ldots, n\) and \(j=1, \ldots, k\).

  1. Choose initial values of the parameters: \({\theta }^{\left(0\right)}=\left({\alpha }_{1}^{\left(0\right)},...,{\alpha }_{k}^{\left(0\right)};{\theta }_{1}^{\left(0\right)},...,{\theta }_{k}^{\left(0\right)}\right)\)

  2. Expectation step: in the \(j\)-th step, compute the matrix \(W = {(w_{ij})}_{nxk}\) with the weights \(w_{ij}\)

  3. Maximization step: in the \(j\)-th step, for all \(r=1, \ldots, k\) compute:

    1. The mixture weights \({\alpha }_{r}^{\left(j+1\right)}=\frac{{n}_{r}}{n}\), where \({n}_{r}=\sum _{i=1}^{n}{w}_{ir}\) is the “amount” of the feature vectors that are assigned to the \(r\)-th mixture component

    2. Mean estimates \({m}_{r}^{\left(j+1\right)}=\frac{1}{{n}_{r}}\sum _{i=1}^{n}{w}_{ir}{x}_{i}\)

    3. Covariance estimate \(\sum _{r}^{(j+1)}=({\sigma }_{r,hg}^{(j+1)})\) of size \(p \times p\) with \(\sigma_{r,hg}^{(j+1)}=\frac{1}{n_r}\sum_{l=1}^{n}{w}_{lr}(x_{lh}-m_{r,h}^{(j+1)})(x_{lg}-m_{r,g}^{(j+1)})\)

  4. Repeat steps 2 and 3 until any of these conditions is met:

    • \(|\log({\theta }^{(j+1)}-{\theta }^{(j)})|<\epsilon\), where the likelihood function is:

      \(\log(\theta)=\sum_{i=1}^{n}\log(\sum _{j=1}^{k}{pd(x}_{i}|{z}_{j},{\theta }_{j}){\alpha }_{j})\)

    • The number of iterations exceeds the predefined level.

Initialization

The EM algorithm for GMM requires initialized vector of weights, vectors of means, and variance-covariance [Biernacki2003, Maitra2009].

The EM initialization algorithm for GMM includes the following steps:

  1. Perform nTrials starts of the EM algorithm with nIterations iterations and start values:

    • Initial means - \(k\) different random observations from the input data set

    • Initial weights - the values of \(1/k\)

    • Initial covariance matrices - the covariance of the input data

  2. Regard the result of the best EM algorithm in terms of the likelihood function values as the result of initialization

Initialization

The EM algorithm for GMM requires initialized vector of weights, vectors of means, and variance-covariance. Skip the initialization step if you already calculated initial weights, means, and covariance matrices.

Batch Processing

Algorithm Input

The EM for GMM initialization algorithm accepts the input described below. Pass the Input ID as a parameter to the methods that provide input for your algorithm.

Input ID

Input

data

Pointer to the \(n \times p\) numeric table with the data to which the EM initialization algorithm is applied. The input can be an object of any class derived from NumericTable.

Algorithm Parameters

The EM for GMM initialization algorithm has the following parameters:

Parameter

Default Value

Description

algorithmFPType

float

The floating-point type that the algorithm uses for intermediate computations. Can be float or double.

method

defaultDense

Performance-oriented computation method, the only method supported by the algorithm.

nComponents

Not applicable

The number of components in the Gaussian Mixture Model, a required parameter.

nTrials

\(20\)

The number of starts of the EM algorithm.

nIterations

\(10\)

The maximal number of iterations in each start of the EM algorithm.

accuracyThreshold

1.0e-04

The threshold for termination of the algorithm.

covarianceStorage

full

Covariance matrix storage scheme in the Gaussian Mixture Model:

  • full - covariance matrices are stored as numeric tables of size \(p \times p\). All elements of the matrix are updated during the processing.

  • diagonal - covariance matrices are stored as numeric tables of size \(1 \times p\). Only diagonal elements of the matrix are updated during the processing, and the rest are assumed to be zero.

engine

SharePtr< engines:: mt19937:: Batch>()

Pointer to the random number generator engine that is used internally to get the initial means in each EM start.

Algorithm Output

The EM for GMM initialization algorithm calculates the results described below. Pass the Result ID as a parameter to the methods that access the results of your algorithm.

Result ID

Result

weights

Pointer to the \(1 \times k\) numeric table with mixture weights.

Note

By default, this result is an object of the HomogenNumericTable class, but you can define the result as an object of any class derived from NumericTable except PackedTriangularMatrix, PackedSymmetricMatrix, and CSRNumericTable.

means

Pointer to the \(k \times p\) numeric table with each row containing the estimate of the means for the \(i\)-th mixture component, where \(i=0, 1, \ldots, k-1\).

Note

By default, this result is an object of the HomogenNumericTable class, but you can define the result as an object of any class derived from NumericTable except PackedTriangularMatrix, PackedSymmetricMatrix, and CSRNumericTable.

covariances

Pointer to the DataCollection object that contains \(k\) numeric tables, each with the \(p \times p\) variance-covariance matrix for the \(i\)-th mixture component of size:

  • \(p \times p\) - for the full covariance matrix storage scheme

  • \(1 \times p\) - for the diagonal covariance matrix storage scheme

Note

By default, the collection contains objects of the HomogenNumericTable class, but you can define them as objects of any class derived from NumericTable except PackedTriangularMatrix and CSRNumericTable.

Computation

Batch Processing

Algorithm Input

The EM for GMM algorithm accepts the input described below. Pass the Input ID as a parameter to the methods that provide input for your algorithm.

Input ID

Input

data

Pointer to the \(n \times p\) numeric table with the data to which the EM algorithm is applied. The input can be an object of any class derived from NumericTable.

inputWeights

Pointer to the \(1 \times k\) numeric table with initial mixture weights. This input can be an object of any class derived from NumericTable.

inputMeans

Pointer to a \(k \times p\) numeric table. Each row in this table contains the initial value of the means for the \(i\)-th mixture component, where \(i = 0, 1, \ldots, k-1\). This input can be an object of any class derived from NumericTable.

inputCovariances

Pointer to the DataCollection object that contains \(k\) numeric tables, each with the \(p \times p\) variance-covariance matrix for the \(i\)-th mixture component of size:

  • \(p \times p\) - for the full covariance matrix storage scheme

  • \(1 \times p\) - for the diagonal covariance matrix storage scheme

The collection can contain objects of any class derived from NumericTable.

inputValues

Pointer to the result of the EM for GMM initialization algorithm. The result of initialization contains weights, means, and a collection of covariances. You can use this input to set the initial values for the EM for GMM algorithm instead of explicitly specifying the weights, means, and covariance collection.

Algorithm Parameters

The EM for GMM algorithm has the following parameters:

Parameter

Default Value

Description

algorithmFPType

float

The floating-point type that the algorithm uses for intermediate computations. Can be float or double.

method

defaultDense

Performance-oriented computation method, the only method supported by the algorithm.

nComponents

Not applicable

The number of components in the Gaussian Mixture Model, a required parameter.

maxIterations

\(10\)

The maximal number of iterations in the algorithm.

accuracyThreshold

1.0e-04

The threshold for termination of the algorithm.

covariance

Pointer to an object of the BatchIface class

Pointer to the algorithm that computes the covariance matrix.

Note

By default, the respective oneDAL algorithm is used, implemented in the class derived from BatchIface.

regularizationFactor

\(0.01\)

Factor for covariance regularization in the case of ill-conditional data.

covarianceStorage

full

Covariance matrix storage scheme in the Gaussian Mixture Model:

  • full - covariance matrices are stored as numeric tables of size \(p \times p\). All elements of the matrix are updated during the processing.

  • diagonal - covariance matrices are stored as numeric tables of size \(1 \times p\). Only diagonal elements of the matrix are updated during the processing, and the rest are assumed to be zero.

Algorithm Output

The EM for GMM algorithm calculates the results described below. Pass the Result ID as a parameter to the methods that access the results of your algorithm.

Result ID

Result

weights

Pointer to the \(1 \times k\) numeric table with mixture weights.

Note

By default, this result is an object of the HomogenNumericTable class, but you can define the result as an object of any class derived from NumericTable except PackedTriangularMatrix, PackedSymmetricMatrix, and CSRNumericTable.

means

Pointer to the \(k \times p\) numeric table with each row containing the estimate of the means for the \(i\)-th mixture component, where \(i=0, 1, \ldots, k-1\).

Note

By default, this result is an object of the HomogenNumericTable class, but you can define the result as an object of any class derived from NumericTable except PackedTriangularMatrix, PackedSymmetricMatrix, and CSRNumericTable.

covariances

Pointer to the DataCollection object that contains \(k\) numeric tables, each with the \(p \times p\) variance-covariance matrix for the \(i\)-th mixture component of size:

  • \(p \times p\) - for the full covariance matrix storage scheme

  • \(1 \times p\) - for the diagonal covariance matrix storage scheme

Note

By default, the collection contains objects of the HomogenNumericTable class, but you can define them as objects of any class derived from NumericTable except PackedTriangularMatrix and CSRNumericTable.

goalFunction

Pointer to the \(1 \times 1\) numeric table with the value of the logarithm of the likelihood function after the last iteration.

Note

By default, this result is an object of the HomogenNumericTable class.

nIterations

Pointer to the \(1 \times 1\) numeric table with the number of iterations computed after completion of the algorithm.

Note

By default, this result is an object of the HomogenNumericTable class.

Examples

Batch Processing:

Note

There is no support for Java on GPU.

Batch Processing:

Batch Processing:

Performance Considerations

To get the best overall performance of the expectation-maximization algorithm at the initialization and computation stages:

  • If input data is homogeneous, provide the input data and store results in homogeneous numeric tables of the same type as specified in the algorithmFPType class template parameter.

  • If input data is non-homogeneous, use AOS layout rather than SOA layout.

Optimization Notice

Intel’s compilers may or may not optimize to the same degree for non-Intel microprocessors for optimizations that are not unique to Intel microprocessors. These optimizations include SSE2, SSE3, and SSSE3 instruction sets and other optimizations. Intel does not guarantee the availability, functionality, or effectiveness of any optimization on microprocessors not manufactured by Intel. Microprocessor-dependent optimizations in this product are intended for use with Intel microprocessors. Certain optimizations not specific to Intel microarchitecture are reserved for Intel microprocessors. Please refer to the applicable product User and Reference Guides for more information regarding the specific instruction sets covered by this notice.

Notice revision #20110804