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

Question

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: ###
[1] -2940.794
[1] -2926.269
[1] -2909.906
[1] -2904.717
[1] -2898.039
[1] -2890.579
[1] -2881.971
[1] -2879.668
[1] -2879.223
[1] -2879.206
### Monotone increase, great!

> EM(testdata, lambda = 0.4, testmixtures, testmeans, testcovariances)
[1] -3590.492
[1] -3568.838 # increase
[1] -3561.29  # increase
[1] -3566.702 # decrease!
[1] -3568.456 # decrease!
[1] -3567.888 # increase
[1] -3567.33  # increase
[1] -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.
    newEstimates <- list()
    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
    covariances <- newEstimates
    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))
}

Answer 1

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.