spaMM (version 3.1.2)

mat_sqrt: Computation of “square root” of symmetric positive definite matrix

Description

mat_sqrt is not usually directly called by users, but arguments may be passed to it through higher-level calls (see Examples). For given matrix C, it computes a factor L such that C = L * t(L), handling issues with nearly-singular matrices. The default behavior is to try Cholesky factorization, and use eigen if it fails. Matrix roots are not unique (for example, they are lower triangular for t(chol(.)), and symmetric for svd(.). As matrix roots are used to simulate samples under the fitted model (in particular in the parametric bootstrap implemented in fixedLRT), this implies that for given seed of random numbers, these samples will differ with these different methods (although their distribution should be identical).

designL.from.Corr is an older procedure with the same purpose. Set spaMM.options(mat_sqrt_fn="designL.from.Corr") to restore its use.

Usage

mat_sqrt(m = NULL, symSVD = NULL, try.chol = TRUE, condnum=1e12)

Arguments

m

The matrix which 'root' is to be computed. This argument is ignored if symSVD is provided.

symSVD

A list representing the symmetric singular value decomposition of the matrix which 'root' is to be computed. Must have elements $u, a matrix of eigenvectors, and $d, a vector of eigenvalues.

try.chol

If try.chol=TRUE, the Cholesky factorization will be tried.

condnum

(large) numeric value. In the case chol() was tried and failed, the matrix is regularized so that its (matrix 2-norm) condition number is reduced to condnum.

Value

For non-NULL m, its matrix root, with rows and columns labelled according to the columns of the original matrix. If eigen was used, the symmetric singular value decomposition (a list with members u (matrix of eigenvectors) and d (vector of eigenvalues)) is given as attribute.

Examples

# NOT RUN {
## try.chol argument passed to mat_sqrt 
## through the '...' argument of higher-level functions
## such as HLCor, corrHLfit, fixedLRT:
data("scotlip")
HLCor(cases~I(prop.ag/10) +adjacency(1|gridcode)+offset(log(expec)),
      ranPars=list(rho=0.174),adjMatrix=Nmatrix,family=poisson(),
      data=scotlip,try.chol=FALSE)
# }
# NOT RUN {
# }