corpcor (version 1.5.0)

powcor.shrink: Fast Computation of the Power of the Shrinkage Correlation Matrix

Description

The function powcor.shrink efficiently computes the alpha-th power of the shrinkage correlation matrix produced by cor.shrink.

Usage

powcor.shrink(x, alpha, lambda, w, collapse=FALSE, verbose=TRUE)

Arguments

x

a data matrix

alpha

exponent

lambda

the correlation shrinkage intensity (range 0-1). If lambda is not specified (the default) it is estimated using an analytic formula from Sch"afer and Strimmer (2005) - see cor.shri

w

optional: weights for each data point - if not specified uniform weights are assumed (w = rep(1/n, n) with n = nrow(x)).

collapse

return vector instead of matrix if estimated or specified lambda equals 1.

verbose

output status while computing (default: TRUE)

Value

• powcor.shrink returns a matrix containing the output from cor.shrink taken to the power alpha.

Details

This function employs the following identity to compute the matrix power of the correlation matrix.

Apart from a scaling factor the shrinkage correlation matrix computed by cor.shrink takes on the form

$$Z = I_p + V M V' .$$

Note that Z is of size p times p whereas M is a matrix of size m times m, where m is the rank of the matrix Z. In order to calculate the alpha-th power of Z the identity

$$Z^\alpha = I_p - V (I_m -(I_m + M)^\alpha) V'$$

requires only the computation of the alpha-th power of the matrix $I_m + M$. This trick enables substantial computational savings when the number of samples is much smaller than the number of observations. Note that the above identity is related but not identical to the Woodbury matrix identity for inversion of a matrix. For $alpha=-1$ the above identity reduces to

$$Z^{-1} = I_p - V (I_m -(I_m + M)^{-1}) V' ,$$

whereas the Woodbury matrix identity equals

$$Z^{-1} = I_p - V (I_m + M^{-1})^{-1} V' .$$

See Also

invcor.shrink, cor.shrink, mpower.

Examples

# load corpcor library
library("corpcor")

# generate data matrix
p = 2000
n = 10
X = matrix(rnorm(n*p), nrow = n, ncol = p)

lambda = 0.23  # some arbitrary lambda

### computing the inverse ###
# slow
system.time(
  (W1 = solve(cor.shrink(X, lambda=lambda)))
)

# very fast
system.time(
  (W2 = powcor.shrink(X, alpha=-1, lambda=lambda))
)

# no difference
sum((W1-W2)^2)

### computing the square root ###

system.time(
  (W1 = mpower(cor.shrink(X, lambda=lambda), alpha=0.5))
)

# very fast
system.time(
  (W2 = powcor.shrink(X, alpha=0.5, lambda=lambda))
)

# no difference
sum((W1-W2)^2)


### computing an arbitrary power (alpha=1.23) ###

system.time(
  (W1 = mpower(cor.shrink(X, lambda=lambda), alpha=1.23))
)

# very fast
system.time(
  (W2 = powcor.shrink(X, alpha=1.23, lambda=lambda))
)

# no difference
sum((W1-W2)^2)