powcor.shrink {corpcor} | R Documentation |
The function powcor.shrink
efficiently computes the alpha
-th power
of the shrinkage correlation matrix produced by cor.shrink
.
powcor.shrink(x, alpha, lambda, w, collapse=FALSE, verbose=TRUE)
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.shrink .
For lambda=0 the empirical correlations are recovered. |
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) |
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^α = I_p - V (I_m -(I_m + M)^α) 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' .
powcor.shrink
returns a matrix containing the output from cor.shrink
taken to the power alpha
.
Korbinian Strimmer (http://strimmerlab.org).
invcor.shrink
, cor.shrink
, mpower
.
## Not run: # 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) ## End(Not run)