# Coordinate Descent for Variance-Component Models
###### Abstract  
Variance-component models are an indispensable tool for statisticians wanting to capture both random and fixed model effects. They have applications in a wide range of scientific disciplines. While maximum likelihood estimation (MLE) is the most popular method for estimating the variance-component model parameters, it is numerically challenging for large data sets. In this article, we consider the class of coordinate descent (CD) algorithms for computing the MLE. We show that a basic implementation of coordinate descent is numerically costly to implement and does not easily satisfy the standard theoretical conditions for convergence. We instead propose two parameter-expanded versions of CD, called PX-CD and PXI-CD. These novel algorithms not only converge faster than existing competitors (MM and EM algorithms) but are also more amenable to convergence analysis. PX-CD and PXI-CD are particularly well-suited for large data sets--namely, as the scale of the model increases, the performance gap between the parameter-expanded CD algorithms and the current competitor methods increases. linear-mixed models; maximum likelihood estimation; numerical optimization +
Footnote †: journal:  
00 July 20223 September 20223 September 20223 September 20223 September 20230 July 20223 September 20230 July 20230 July 20230 September 20230 July  
To simplify the MLE estimate for \(\gamma\), one can adopt the restricted MLE (REML) method [4] to remove the mean effect in the likelihood expression by projecting \(y\) onto the null space of \(\mathbf{X}\). Let \(\
u=n-p\) and suppose we have the QR decomposition,  
\[\mathbf{X}= \begin{bmatrix}\mathbf{Q}_{[p]}&\mathbf{Q}_{[\
u]}\end{bmatrix} \begin{bmatrix}\mathbf{R}\\ \mathbf{0}_{\
u\times p}\end{bmatrix} ,\]  
where \(\mathbf{R}\) is an \(p\times p\) upper triangular matrix, \(\mathbf{0}_{[\
u\times p]}\) is an \(\
u\times p\) zero matrix, \(\mathbf{Q}_{[p]}\) is an \(n\times p\) matrix, \(\mathbf{Q}_{[\
u]}\) is \(n\times\
u\), and \(\mathbf{Q}_{[p]}\) and \(\mathbf{Q}_{[\
u]}\) both have orthogonal columns. If we take the Cholesky decomposition \(\mathbf{LL}^{\top}\) of the matrix \(\mathbf{Q}_{[\
u]}^{\top}\mathbf{V}_{0}\mathbf{Q}_{[\
u]}\), then the transformation, \(\mathbf{L}^{-1}\mathbf{Q}_{[\
u]}^{\top}:\mathbb{R}^{n}\mapsto\mathbb{R}^{\
u}\) removes the mean from the response, and we obtain \(\boldsymbol{y}^{\prime}:=\mathbf{L}^{-1}\mathbf{Q}_{[\
u]}^{\top}\boldsymbol{y }\sim\mathcal{N}(\mathbf{0},\mathbf{\Omega}^{\prime})\), where \(\mathbf{\Omega}^{\prime}:=\gamma_{0}\mathbf{I}+\sum_{i=1}^{m}\gamma_{i} \mathbf{V}_{i}^{\prime}\in\mathbb{R}^{\
u\times\
u}\) and \(\mathbf{V}_{i}^{\prime}:=\mathbf{L}^{-1}\mathbf{Q}_{[\
u]}^{\top}\mathbf{V}_ {i}\mathbf{Q}_{[\
u]}\mathbf{L}^{-\top}\). After this transformation, the restricted likelihood \(-\frac{1}{2}\ln\det\mathbf{\Omega}^{\prime}-\frac{1}{2}(\boldsymbol{y}^{ \prime})^{\top}(\mathbf{\Omega}^{\prime})^{-1}\boldsymbol{y}^{\prime}\) will not depend on \(\boldsymbol{\beta}\). Henceforth, and without loss of generality, we assume that such a transformation has been performed so that we can focus on minimizing an objective function of the form:  
\[-2L(\boldsymbol{\gamma})=\ln\det\mathbf{\Omega}+\boldsymbol{y}^{\top}\mathbf{ \Omega}^{-1}\boldsymbol{y}. \tag{2}\]  
There exists extensive literature on optimization methods for the log-likelihood expression (2), including Newton's Method [5], Fisher Scoring Method [6], the EM and MM Algorithms [7; 8; 9]. [PERSON]'s method is known to scale poorly as \(n\), or the number of variance components \(m+1\), increase due the cost of \(O(mn^{3})+O(m^{3})\) flops required to invert a Hessian matrix at each update. Both the EM and MM algorithms have simple updating steps; however, numerical experience shows that they are slow to identify the active set \(\{i:\hat{\gamma}_{i}=0\}\), where \(\hat{\gamma}\) is the MLE. One class of algorithms yet to be applied to this problem are coordinate-descent (CD) algorithms.