Difference between revisions of "Algorithms"

From Linear Mixed Models Toolbox
Jump to navigation Jump to search
Line 12: Line 12:
The direct solver requires the mixed model coefficient matrix to be build and all Kronecker products to be resolved. This can be quite memory demanding and should therefore be used carefully. The direct solver uses a [https://en.wikipedia.org/wiki/Cholesky_decomposition Cholesky decomposition] and [https://en.wikipedia.org/wiki/Triangular_matrix#Forward_and_back_substitution forward-backward-substitution] to solve the mixed model equation system, where especially the decomposition step can be very resource demanding and time consuming.
The direct solver requires the mixed model coefficient matrix to be build and all Kronecker products to be resolved. This can be quite memory demanding and should therefore be used carefully. The direct solver uses a [https://en.wikipedia.org/wiki/Cholesky_decomposition Cholesky decomposition] and [https://en.wikipedia.org/wiki/Triangular_matrix#Forward_and_back_substitution forward-backward-substitution] to solve the mixed model equation system, where especially the decomposition step can be very resource demanding and time consuming.


==Gibbs sampling of variance components==
==Variance component estimation==
===Single pass Gibbs sampling===
===Gibbs sampling====
====Single pass Gibbs sampling====
{{lmt}}'s single pass Gibbs sampling algorithm is described in <ref name="Sorensen2002" />. In short, all location parameters are drawn from their joint conditional posterior distribution. Note that this requires solving the mixed model equation system once per iteration which usually leads to a substantial increase in processing time.
{{lmt}}'s single pass Gibbs sampling algorithm is described in <ref name="Sorensen2002" />. In short, all location parameters are drawn from their joint conditional posterior distribution. Note that this requires solving the mixed model equation system once per iteration which usually leads to a substantial increase in processing time.


===Blocked Gibbs sampling===
====Blocked Gibbs sampling====
For random factors {{lmt}}'s blocked Gibbs sampler draws correlated location parameters within factor level from their joint conditional posterior distribution. Location parameters of fixed factors are drawn in scalar mode from their fully conditional posterior.
For random factors {{lmt}}'s blocked Gibbs sampler draws correlated location parameters within factor level from their joint conditional posterior distribution. Location parameters of fixed factors are drawn in scalar mode from their fully conditional posterior.
===Restricted Maximum Likelyhood===
====MC-EM-REML====
{{lmt}} provides a monte-carlo expectation-maximisation REML algorithms which uses the preconditioned gradient solver for solving the mixed model equations and a blocked Gibbs sampler to sample the necessary traces<ref name="Harville2004" />.


==Prediction error variance sampling==
 
==Prediction error variances==
===Gibbs Sampling===
Following the approach of Harville(1999)<ref name="Harville1999" /> {{lmt}} can sample for fixed factors the diagonal elements of the inverse of the mixed model coefficient matrix, for random factors the diagonal blocks of the inverse of the coefficient matrix where the block size is determined by the dimension of the related $$\Sigma$$ matrix. The blocks are the prediction error co-variance matrices of the factor levels of correlated sub-factors. When sampling prediction error variances {{lmt}} can run many Gibbs chains in parallel allowing to exploit multi-core hardware architecture. However, it is recommended to specify not more chains than the number of available '''real''' cores excluding hyper-threading technology.
Following the approach of Harville(1999)<ref name="Harville1999" /> {{lmt}} can sample for fixed factors the diagonal elements of the inverse of the mixed model coefficient matrix, for random factors the diagonal blocks of the inverse of the coefficient matrix where the block size is determined by the dimension of the related $$\Sigma$$ matrix. The blocks are the prediction error co-variance matrices of the factor levels of correlated sub-factors. When sampling prediction error variances {{lmt}} can run many Gibbs chains in parallel allowing to exploit multi-core hardware architecture. However, it is recommended to specify not more chains than the number of available '''real''' cores excluding hyper-threading technology.
===Solving===
{{lmt}} can obtain elements of the inverse of the coefficient matrix via solving the mixed model equations. This method is currently only supported for the diagonal prediction error co-variance blocks of random factors, where the block size is determined by the dimension of the related $$\Sigma$$ matrix. For this algorithm {{lmt}} can utilize either the [#Iterative Solver] or the [#Direct Solver].


==REML variance component estimation==
===MC-EM-REML===
{{lmt}} provides a monte-carlo expectation-maximisation REML algorithms which uses the preconditioned gradient solver for solving the mixed model equations and a blocked Gibbs sampler to sample the necessary traces<ref name="Harville2004" />.


==References==
==References==

Revision as of 04:26, 7 March 2021

Solving Linear Mixed model Equations

lmt supports two types of solver for solving MME's: a direct solver and an iterative solver

Iterative solver

The iterative solver uses the preconditioned conjugate gradient method and is lmt's default solver. It does not require the explicit construction of any mixed model equation, and is therefore less resource demanding than the direct solver. That is, many models which cannot be solved using the direct solver can still be solved using the iterative solver. Even for small models the iterative solver usually outperforms the direct solver in terms of total processing time.

The iterative solver has converged to a stable solution if $$log_e\left(\sqrt{\frac{(Cx-b)'(Cx-b)}{b'b}}\right)<t$$, where $$C$$ is the mixed-model coefficient matrix, $$x$$ is the solution vector, $$b$$ is the right-hand side and $$t$$ is the convergence threshold. The default convergence threshold is -18.42, which is equivalent to $$ \sqrt{\frac{\sum_{i=1}^n ((Ax)_i-b_i)^2}{n}/\frac{\sum_{i=1}^n b_i^2}{n}}<10^{-8} $$

Direct solver

The direct solver requires the mixed model coefficient matrix to be build and all Kronecker products to be resolved. This can be quite memory demanding and should therefore be used carefully. The direct solver uses a Cholesky decomposition and forward-backward-substitution to solve the mixed model equation system, where especially the decomposition step can be very resource demanding and time consuming.

Variance component estimation

Gibbs sampling=

Single pass Gibbs sampling

lmt's single pass Gibbs sampling algorithm is described in [1]. In short, all location parameters are drawn from their joint conditional posterior distribution. Note that this requires solving the mixed model equation system once per iteration which usually leads to a substantial increase in processing time.

Blocked Gibbs sampling

For random factors lmt's blocked Gibbs sampler draws correlated location parameters within factor level from their joint conditional posterior distribution. Location parameters of fixed factors are drawn in scalar mode from their fully conditional posterior.

Restricted Maximum Likelyhood

MC-EM-REML

lmt provides a monte-carlo expectation-maximisation REML algorithms which uses the preconditioned gradient solver for solving the mixed model equations and a blocked Gibbs sampler to sample the necessary traces[2].


Prediction error variances

Gibbs Sampling

Following the approach of Harville(1999)[3] lmt can sample for fixed factors the diagonal elements of the inverse of the mixed model coefficient matrix, for random factors the diagonal blocks of the inverse of the coefficient matrix where the block size is determined by the dimension of the related $$\Sigma$$ matrix. The blocks are the prediction error co-variance matrices of the factor levels of correlated sub-factors. When sampling prediction error variances lmt can run many Gibbs chains in parallel allowing to exploit multi-core hardware architecture. However, it is recommended to specify not more chains than the number of available real cores excluding hyper-threading technology.

Solving

lmt can obtain elements of the inverse of the coefficient matrix via solving the mixed model equations. This method is currently only supported for the diagonal prediction error co-variance blocks of random factors, where the block size is determined by the dimension of the related $$\Sigma$$ matrix. For this algorithm lmt can utilize either the [#Iterative Solver] or the [#Direct Solver].


References

  1. D. Sorensen and D. Gianola; Likelihood, Bayesian, and MCMC Methods in Quantitative Genetics; 2002; 584-588
  2. David A. Harville; Making REML computationally feasible for large data sets: use of the Gibbs sampler; Journal of Statistical Computation & Simulation; 2004
  3. David A. Harville; Use of the Gibbs sampler to invert large, possibly sparse, positive definite matrices; Linear Algebra and its Applications; 1999