dot-armaRidgeP: Core ridge precision estimators

In rags2ridges: Ridge Estimation of Precision Matrices from High-Dimensional Data

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).

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

\hat{\mathbf{\Omega}}^{\mathrm{I}a}(\lambda_{a}) = \left\{\left[\lambda_{a}\mathbf{I}_{p} + \frac{1}{4}(\mathbf{S} - \lambda_{a}\mathbf{T})^{2}\right]^{1/2} + \frac{1}{2}(\mathbf{S} - \lambda_{a}\mathbf{T})\right\}^{-1},

if invert == 1 or

\hat{\mathbf{\Omega}}^{\mathrm{I}a}(\lambda_{a}) = \frac{1}{\lambda}\left\{\left[\lambda_{a}\mathbf{I}_{p} + \frac{1}{4}(\mathbf{S} - \lambda_{a}\mathbf{T})^{2}\right]^{1/2} - \frac{1}{2}(\mathbf{S} - \lambda_{a}\mathbf{T})\right\}

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 <[e-mail]>

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].

See Also

Used as a backbone in ridgeP, ridgeP.fused, etc.

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 documentation built on Oct. 14, 2023, 5:06 p.m.