# Gaussian Mixture model - Penalized log-likelihood in EM algorithm not monotone increasing

I am working on a multivariate Gaussian Mixture Model in R. The goal is to do regularized clustering on the data, where each component represents a cluster.

I wrote an EM algorithm to maximize a penalized expected log-likelihood function given by $$\ell(\Theta)=\sum_{i=1}^n \log{(\pi_k)} + \log{(f_k(x_i | \mu_k, \Sigma_k))}- \frac{n}{2}\lambda \sum^K_{k=1} \pi_k ||\Omega_k||$$

with $$K$$ components, $$f_k$$ giving the probability density, $$\pi_k$$ mixing proportions, $$\mu_k$$ means, $$\Sigma_k$$ variance-covariance matrices, and $$\Sigma_k^{-1}=\Omega_k$$ precision matrices, also known as concentration matrices. $$\frac{n}{2}\lambda \sum \pi_k ||\Omega_k||$$ is the penalty term, with tuning parameter $$\lambda$$ and $$||\Omega_k|| = \sum_{i,j} |\omega^k_{i,j}|$$.

Given some data set, I am interested in estimating the parameters that maximize this log-likelihood for some value of the tuning parameter $$\lambda$$. To obtain penalized estimates of the precision matrices $$\Omega_k$$, I am using the R package glasso, which is a quick estimator of sparse precision matrices.

I implemented this in the Expectation-Maximization algorithm, which is supposed to lead to a monotone increase of the log-likelihood until convergence, as follows:

1. Initialize parameters $$\pi_k^{(0)}, \mu_k^{(0)}, \Sigma_k^{(0)}$$ randomly.
2. E-step: compute posterior probabilities $$\tau_{ik}^{(t)}$$ given the (new) parameters.
3. M-step (1): update parameters $$\pi_k^{(t+1)}$$ and $$\mu_k^{(t+1)}$$, and calculate sample covariance matrix $$S_k^{(t)} = \frac{\sum^n_{i=1} \tau_{ik}^{(t)}\left( x_i - \mu_k^{(t+1)}\right)\left( x_i - \mu_k^{(t+1)}\right)^T}{\sum^n_{i=1} \tau_{ik}^{(t)}}$$.
4. M-step (2): provide $$S_k^{(t)}$$ to glasso to obtain $$\Omega_k^{(t+1)}$$. Invert $$\Omega_k^{(t+1)}$$ to obtain the penalized estimate of $$\Sigma_k^{(t+1)}$$.
5. Repeat steps 2 - 4 until convergence.

The problem: my implementation of the EM produces increasing log-likelihoods when $$\lambda = 0$$, but does not monotonically increase when $$\lambda > 0$$. Usually, the value increases for a number of iterations, then decreases for some, and then goes back up by the end.

The question: is it supposed to do this? If not, what is the error in my implementation?

Example:

> EM(testdata, lambda = 0, testmixtures, testmeans, testcovariances)
### OUTPUT: ###
 -2940.794
 -2926.269
 -2909.906
 -2904.717
 -2898.039
 -2890.579
 -2881.971
 -2879.668
 -2879.223
 -2879.206
### Monotone increase, great!

> EM(testdata, lambda = 0.4, testmixtures, testmeans, testcovariances)
 -3590.492
 -3568.838 # increase
 -3561.29  # increase
 -3566.702 # decrease!
 -3568.456 # decrease!
 -3567.888 # increase
 -3567.33  # increase
 -3567.234 # increase


Context: I have been testing with a dataset of size $$N = 100$$ with 18 variables. I am certain that my unpenalized log-likelihood function is calculating properly, as I compared it with mvnormalmixEM from the R package mixtools. I cannot for the life of me find any errors in the calculations of the updates. I suspect the problem comes from applying the penalty wrong (e.g. double penalization somewhere).

Below you will find my code, I will try to keep it as concise as possible. I'm not sure how to properly simulate a suitable dataset for this and I'm hesitant to upload my own data here, so apologies for that. If there is any other information that would be helpful in solving this, please ask. Thanks in advance for even looking at this!

library("glasso")
library("mvtnorm")

# Initialize parameters in order to start the EM algorithm
init <- function(data, nclust){
mixtures <- c()
means <- list()
covariances <- list()

init.assign <- sample(1:nclust, nrow(data), replace=T)

for(k in 1:nclust){
mixtures[k] <- table[names(table)==k] / nrow(data)
means[[k]] <- as.vector(colMeans(data[which(init.assign %in% k),]))
covariances[[k]] <- cov(data[which(init.assign %in% k),])
}
return(list(mixtures = mixtures, means = means, covariances = covariances))
}

# Calculate posterior probabilities tau_ik
soft.assign <- function(data, mixtures, means, covariances){
tauTable <- matrix(NA, nrow = nrow(data), ncol = length(mixtures))

for (c in 1:ncol(tauTable)){
tauTable[,c] <- dmvnorm(data, mean=as.numeric(means[[c]]),
sigma=covariances[[c]])
tauTable[,c] <- mixtures[c] * tauTable[,c]
}
tauTable <- tauTable / rowSums(tauTable)
return(tauTable)
}

# Calculate (penalized) log-likelihood given model parameters
LL <- function(data, lambda, mixtures, means, covariances, penalize = F){
# Calculate log-likelihood for mixture model
LL <- sum(log(apply(sapply(lapply(1:length(mixtures),
function(i) mixtures[i] * dmvnorm(data, means[[i]],
covariances[[i]])), cbind), 1, sum)))

# Calculate and apply penalty through parameter lambda
if(penalize){
O <- lapply(covariances, solve)
sum <- sum(sapply(1:length(mixtures), function(i) mixtures[i] *
sum(abs(O[[i]]))))

penalty <- (nrow(data)/2) * lambda * sum
pLL <- LL - penalty
return(pLL)
}
else{ return(LL) }
}

# Main EM function. Returns parameters at convergence.
EM <- function(data, lambda, mixtures, means, covariances,
iterations = 100, = 10^-4){
data <- as.matrix(data)

K <- length(mixtures)
P <- ncol(data)
N <- nrow(data)

deltaLL <- 1

prevLL <- LL(data, lambda, mixtures, means, covariances, penalize=T)
LLpath <- c(prevLL)

mixtures <- mixtures
means <- means
covariances <- covariances
precisions <- list()

increment <- 1
while(increment < iterations && deltaLL > threshold){

# Expectation step: calculate posterior probabilities tau_{ik}
tauTable <- soft.assign(data, mixtures, means, covariances)

# Maximization step: update parameters mixtures, means, covariances, precisions
# Calculating mixtures update.
newMixtures <- colMeans(tauTable)

# Calculating means update.
newMeans <- lapply(1:K, function(j)
sapply(1:18, function(i)
apply(tauTable * data[, i], 2, sum))
[j,]/sum(tauTable[, j]
)
)

# Calculating sample covariances update.
sampleCov <- lapply(1:K, function(j) matrix(
apply(
sapply(1:N, function(i)
tauTable[i, j] * (data[i, ] - newMeans[[j]]) %*%
t(data[i, ] - newMeans[[j]])
), 1, sum
)
, P, P)/sum(tauTable[,j]))

# Calculating updates for covariance estimates & precision matrix estimates.
newPrecisions <- list()
for (k in 1:K){
gl <- glasso(s = sampleCov[[k]], rho = lambda, penalize.diagonal = F)
newEstimates[[k]] <- gl$$w newPrecisions[[k]] <- gl$$wi
}

mixtures <- newMixtures
means <- newMeans
precisions <- newPrecisions

# Calculate model penalized log-likelihood and compare with previous

pLL <- LL(data, lambda, mixtures, means, covariances, penalize=T)
print(pLL)
deltaLL <- abs((pLL / prevLL) - 1)
prevLL <- pLL
LLpath <- c(LLpath, pLL)

increment <- increment + 1

}
return(list(
LL = prevLL,
LLpath = LLpath,
soft.assign = tauTable,
mixtures = mixtures,
means = means,
covariances = covariances,
precisions = precisions))
}

• Is the log likelihood + penalty term monotone increasing? – jbowman Feb 4 '20 at 14:51
• In the EM proof it is only guaranteed to increase log likelihood. – user2939212 Feb 4 '20 at 23:56
• @jbowman It is not. I'm not even 100% sure that it should be. When I set the penalty term to 0, it is strictly monotone increasing, otherwise, it has ups and downs. This entire question is basically me trying to figure out why that is. (: – VHazeleger Feb 6 '20 at 8:29

For posterity and anyone else working on something similar, I found out what my error was: in the glasso call, I set penalize.diagonal == FALSE, which caused the penalty to not apply fully. Setting that parameter to TRUE everywhere has fixed the issue.