rags2ridges: dot-armaRidgeP – R documentation

Become an expert in R — Interactive courses, Cheat Sheets, certificates and more!

dot-armaRidgeP

Core ridge precision estimators

Description

This is the interface to the C++ implementations of the ridge precision estimators. They perform core computations for many other functions.

Usage

.armaRidgeP(S, target, lambda, invert = 2L)

Arguments

`S`	A sample covariance `matrix`.
`target`	A `numeric` symmetric positive definite target `matrix` of the same size as `S`.
`lambda`	The ridge penalty. A `double` of length 1.
`invert`	An `integer`. Should the estimate be computed using inversion? Permitted values are `0L`, `1L`, and `2L` meaning "no", "yes", and "automatic" (default).
`alpha`	The diagonal value on the scalar target matrix. A `double` of length 1.

Details

The functions are R-interfaces to low-level C++ implementations of the ridge estimators in the reference below (cf. Lemma 1, Remark 6, Remark 7, and section 4.1 therein).

.armaRidgeP is simply a wrapper (on the C++ side) for .armaRidgePAnyTarget and .armaRidgePScalarTarget which are the estimators for arbitrary and scalar targets, respectively. The invert argument of the functions indicates if the computation uses matrix inversion or not.

Essentially, the functions compute

Ω(λ) = {[λ I + 1/4(S - λ T)^2 ]^{1/2} + 1/2(S - λ T)}^{-1},

if invert == 1 or

Ω(λ) = 1/{[λ I + 1/4(S - λ T)^2 ]^{1/2} - 1/2(S - λ T)},

if invert == 0 using appropriate eigenvalue decompositions. See the R implementations in the example section below.

Value

Returns a symmetric positive definite matrix of the same size as S.

Warning

The functions themselves perform no checks on the input. Correct input should be ensured by wrappers.

Author(s)

Anders Ellern Bilgrau, Carel F.W. Peeters <cf.peeters@vumc.nl>

References

van Wieringen, W.N. & Peeters, C.F.W. (2016). Ridge Estimation of Inverse Covariance Matrices from High-Dimensional Data, Computational Statistics & Data Analysis, vol. 103: 284-303. Also available as arXiv:1403.0904v3 [stat.ME].

Examples

S <- createS(n = 3, p = 4)
tgt <- diag(4)
rags2ridges:::.armaRidgeP(S, tgt, 1.2)

rags2ridges:::.armaRidgePAnyTarget(S, tgt, 1.2)
rags2ridges:::.armaRidgePScalarTarget(S, 1, 1.2)


################################################################################
# The C++ estimators essentially amount to the following functions.
################################################################################

rev_eig <- function(evalues, evectors) { # "Reverse" eigen decomposition
  evectors %*% diag(evalues) %*% t(evectors)
}

# R implementations

# As armaRidgePScalarTarget Inv
rRidgePScalarTarget <- function(S, a, l, invert = 2) {
  ED <- eigen(S, symmetric = TRUE)
  eigvals <- 0.5*(ED$values - l*a)
  sqroot <- sqrt(l + eigvals^2)

  if (l > 1e6 && (any(!is.finite(eigvals)) || any(!is.finite(sqroot)))) {
    return(diag(a, nrow(S)))
  }

  D_inv <- 1.0/(sqroot + eigvals)
  D_noinv <- (sqroot - eigvals)/l

  if (invert == 2) {   # Determine to invert or not
    if (l > 1) {  # Generally, use inversion for "small" lambda
      invert = 0;
    } else {
      invert <- ifelse(any(!is.finite(D_inv)), 0, 1)
    }
  }

  if (invert) {
    eigvals <- D_inv
  } else {
    eigvals <- D_noinv
  }
  return(rev_eig(eigvals, ED$vectors))
}

# As armaRidgePAnyTarget
rRidgePAnyTarget <- function(S, tgt, l, invert = 2) {
  ED <- eigen(S - l*tgt, symmetric = TRUE)
  eigvals <- 0.5*ED$values
  sqroot <- sqrt(l + eigvals^2)

  if (l > 1e6 && (any(!is.finite(eigvals)) || any(!is.finite(sqroot)))) {
    return(tgt)
  }

  D_inv <- 1.0/(sqroot + eigvals)
  D_noinv <- (sqroot - eigvals)/l

  if (invert == 2) {   # Determine to invert or not
    if (l > 1) {  # Generally, use inversion for "small" lambda
      invert = 0;
    } else {
      invert <- ifelse(any(!is.finite(D_inv)), 0, 1)
    }
  }

  if (invert == 1) {
    eigvals <- D_inv
  } else {
    eigvals <- D_noinv
  }
  return(rev_eig(eigvals, ED$vectors))
}

rRidgeP <- function(S, tgt, l, invert = 2) {
  if (l == Inf) {
    return(tgt)
  }
  a <- tgt[1,1]
  if (tgt == diag(a, nrow(tgt))) {
    rRidgePScalarTarget(S, tgt, l, invert)
  } else {
    rRidgePAnyTarget(S, tgt, l, invert)
  }

}

# Contrasted to the straight forward implementations:
sqrtm <- function(X) { # Matrix square root
  ed <- eigen(X, symmetric = TRUE)
  rev_eig(sqrt(ed$values), ed$vectors)
}

# Straight forward (Lemma 1)
ridgeP1 <- function(S, tgt, l) {
  solve(sqrtm( l*diag(nrow(S)) + 0.25*crossprod(S - l*tgt) ) + 0.5*(S - l*tgt))
}

# Straight forward  (Lemma 1 + remark 6 + 7)
ridgeP2 <- function(S, tgt, l) {
  1.0/l*(sqrtm(l*diag(nrow(S)) + 0.25*crossprod(S - l*tgt)) - 0.5*(S - l*tgt))
}

set.seed(1)
n <- 3
p <- 6
S <- covML(matrix(rnorm(p*n), n, p))
a <- 2.2
tgt <- diag(a, p)
l <- 1.21

(A <- ridgeP1(S, tgt, l))
(B <- ridgeP2(S, tgt, l))

(C  <- rags2ridges:::.armaRidgePAnyTarget(S, tgt, l))
(CR <-                   rRidgePAnyTarget(S, tgt, l))
(D  <- rags2ridges:::.armaRidgePScalarTarget(S, a, l))
(DR <-                   rRidgePScalarTarget(S, a, l))
(E  <- rags2ridges:::.armaRidgeP(S, tgt, l))

# Check
equal <- function(x, y) {isTRUE(all.equal(x, y))}
stopifnot(equal(A, B) & equal(A, C) & equal(A, D) & equal(A, E))
stopifnot(equal(C, CR) & equal(D, DR))

rags2ridges

Ridge Estimation of Precision Matrices from High-Dimensional Data

v2.2.4

GPL (>= 2)

Authors

Carel F.W. Peeters [cre, aut], Anders Ellern Bilgrau [aut], Wessel N. van Wieringen [aut]

Initial release