Smallest square root of matrix

Find a square root of a positive semi-definite matrix, having as few columns as possible. Uses either pivoted Cholesky decomposition or singular value decomposition to do this.

mroot(A,rank=NULL,method="chol")

Arguments

A

The positive semi-definite matrix, a square root of which is to be found.

rank

if the rank of the matrix A is known then it should be supplied. NULL or <1 imply that it should be estimated.

method

"chol" to use pivoted Cholesky decompositon, which is fast but tends to over-estimate rank. "svd" to use singular value decomposition, which is slower, but is the most accurate way to estimate rank.

Details

The function is primarily of use for turning penalized regression problems into ordinary regression problems. Given that A is positive semi-definite the SVD option actually uses a symmetric eigen routine, which gives the same result more efficiently.

Value

A matrix, \( {\bf B}\) with as many columns as the rank of \( {\bf A}\), and such that \( {\bf A} = {\bf BB}^\prime\).

Author

Simon N. Wood simon.wood@r-project.org

Examples

  require(mgcv)
  set.seed(0)
  a <- matrix(runif(24),6,4)
  A <- a%*%t(a) ## A is +ve semi-definite, rank 4
  B <- mroot(A) ## default pivoted choleski method
  tol <- 100*.Machine$double.eps
  chol.err <- max(abs(A-B%*%t(B)));chol.err
#> [1] 4.440892e-16
  if (chol.err>tol) warning("mroot (chol) suspect")
  B <- mroot(A,method="svd") ## svd method
  svd.err <- max(abs(A-B%*%t(B)));svd.err
#> [1] 4.440892e-15
  if (svd.err>tol) warning("mroot (svd) suspect")