Mixed-model Log-likelihood Evaluation Via a Blocked Cholesky Factorization

TL;DR

This paper introduces a computationally efficient method for evaluating the log-likelihood of linear mixed-effects models using a blocked Cholesky factorization, implemented in Julia, enabling analysis of very large datasets.

Contribution

The paper presents a novel derivation and implementation of a blocked Cholesky approach that improves efficiency in mixed-model likelihood evaluation, especially for large-scale data.

Findings

Achieved significant computational speed-up in fitting large mixed models.

Successfully applied the method to a dataset with over 32 million ratings.

Reduced fill-in during Cholesky factorization through optimized ordering and matrix-specific algorithms.

Abstract

Bates et al. (2015) described the evaluation of the profiled log-likelihood of a linear mixed-effects model by updating a sparse, symmetric positive-definite matrix and computing its Cholesky factor, as implemented in the lme4 package for R. Here we present enhancements to the derivation and theoretical presentation of the result and to its implementation using a blocked Cholesky factorization in the MixedModels$.$jl package for Julia (Bezanson et al., 2017). The gain in computational efficiency is primarily due to three factors: (1) the new derivation allows us to compute the penalized residual sum of squares without computing the conditional estimates of the fixed-effects parameters and the conditional modes of the random effects at each optimization step, (2) the blocked Cholesky representation and careful ordering of the random effects terms reduces the amount of "fill-in" that occurs during the Cholesky factorization, and (3) the multiple dispatch feature of the Julia language allows us to use specialized algorithms for different kinds of matrices instead of relying on generic algorithms during the Cholesky factorization. To show the effectiveness of the blocked Cholesky approach we use it to fit a linear mixed model to over 32 million ratings of movies in the MovieLens ml-32m (Harper and Konstan, 2016) data set. The model incorporates random effects for over 200,000 movies and over 80,000 participants. Further enhancements to these computational methods are suggested.