# Coordinate Descent for Variance-Component Models
## 3 Parameter-Expanded CD
### Univariate Minimization via Newton's Method  
Unlike the basic CD Algorithm 1, for which each coordinate update costs \(O(n^{3})\), here we show that a coordinate update for the PX-CD algorithm costs only \(O(jn^{2})\) for some constant \(j\) where typically \(j\ll n\).  
The function \(h_{k}(x)\) can be minimized via the second-order [PERSON]'s method, which numerically finds the root of \(h_{k}^{\prime}(x)\). The basic algorithm starts with an initial guess \(x_{0}\) of the root, and then  
\[x_{n+1}=x_{n}-h_{k}^{\prime}(x_{n})[h_{k}^{\prime\prime}(x_{n})]^{-1}, \tag{14}\]  
are successive better approximations. The algorithm can be terminated once successive iterates are sufficiently close together. The first and second derivatives of \(h\) are given as  
\[h_{k}^{\prime}(x) =-\mathbf{y}^{\top}(\mathbf{\Omega}_{-k}+x\mathbf{V}_{k})^{-1}\mathbf{V}_{k}(\mathbf{ \Omega}_{-k}+x\mathbf{V}_{k})^{-1}\mathbf{y}+\text{tr}\Big{(}\mathbf{C}^{-1}\mathbf{V}_{k} \Big{)}, \tag{15}\] \[h_{k}^{\prime\prime}(x) =2\mathbf{y}^{\top}(\mathbf{\Omega}_{-k}+x\mathbf{V}_{k})^{-1}\mathbf{V}_{k}(\bm {\Omega}_{-k}+x\mathbf{V}_{k})^{-1}\mathbf{V}_{k}(\mathbf{\Omega}_{-k}+x\mathbf{V}_{k})^{-1} \mathbf{y}, \tag{16}\]  
where we used differentiation of a matrix inverse, which implies that  
\[\frac{\partial(\mathbf{\Omega}_{-k}+x\mathbf{V}_{k})^{-1}}{\partial x}=-(\mathbf{\Omega}_ {-k}+x\mathbf{V}_{k})^{-1}\frac{\partial(\mathbf{\Omega}_{-k}+x\mathbf{V}_{k})}{\partial x }(\mathbf{\Omega}_{-k}+x\mathbf{V}_{k})^{-1}.\]  
Similar to the basic CD implementation, computing the algebraic expression in (11) via GEV is expensive. Evaluating (15) and (16) by explicitly calculating \((\mathbf{\Omega}_{-k}+x\mathbf{V}_{k})^{-1}\) is also expensive for large \(n\) and is of time complexity \(O(n^{3})\). Instead, we utilize the conjugate gradient (CG) algorithm [13] to efficiently solve linear systems. At each iteration of [PERSON]'s method, we approximately solve,  
\[(\mathbf{\Omega}_{-k}+x\mathbf{V}_{k})\mathbf{b}=\mathbf{y},\qquad(\mathbf{\Omega}_{-k}+x\mathbf{V}_{ k})\mathbf{c}=\mathbf{V}_{k}\mathbf{b}, \tag{17}\]  
and store the solution in \(\mathbf{b}\) and \(\mathbf{c}\), respectively, via CG algorithm. Generally, \(\|\mathbf{b}-(\mathbf{\Omega}_{-k}+x\mathbf{V}_{k})^{-1}\mathbf{y}\|\) and \(\|\mathbf{c}-(\mathbf{\Omega}_{-k}+x\mathbf{V}_{k})^{-1}\mathbf{V}_{k}\|\) can be made small with \(l\ll n\) iterations, where each iteration requires a matrix-vector-multiplication operation with a \(n\times n\) matrix. The CG algorithm has complexity \(O(ln^{2})\) and can be easily implemented with standard Linear Algebra packages. With the stored approximate solutions, we evaluate the first and second derivatives as  
\[h_{k}^{\prime}(x)=-\mathbf{b}^{\top}\mathbf{V}_{k}\mathbf{b}+\text{tr}\Big{(}\mathbf{C}^{ -1}\mathbf{V}_{k}\Big{)},\qquad h_{k}^{\prime\prime}(x)=2\mathbf{b}^{\top}\mathbf{V}_{k} \mathbf{c}. \tag{18}\]  
Before initiating [PERSON]'s method, we can check if \(k\) is in the active constraint set \(\{k:\hat{\gamma}_{k}=0\}\). Following from Lemma 1, if \(h_{k}^{\prime}(0)\geq 0\), then \(h_{k}(x)\) is non-decreasing on \([0,\infty)\). Then, \(h_{k}(0)\) is the global minimum for \(x\in[0,\infty)\), and we let \(\gamma_{i}^{(t+1)}=0\) if we are in cycle \(t+1\) of PX-CD. If \(h_{k}^{\prime}(0)<0\), we initiate [PERSON]'s method at the current value of the variance component, \(x_{0}=\gamma_{k}^{(t)}\). If \(\dim(\text{span}\{\mathbf{V}_{1},\dots,\mathbf{V}_{m}\})<n\), we require \(\gamma_{0}>0\) so that \(\mathbf{\Omega}\) is invertible.  
In this case, \(k=0\) cannot be in the active constraint set and we immediately initiate [PERSON]'s method at the starting point \(x_{0}=\gamma_{0}^{(t)}\). In rare cases \(h_{k}^{\prime}(x)\) is sufficiently flat at \(x_{n}\) and (14) may significantly overstep the location of the minimizer and return an approximation \(x_{n+1}<0\). In this case, we dampen the step size until \(x_{n+1}>0\).