R: Fast Computation of the Power of the Shrinkage Correlation...
powcor.shrink
R Documentation
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.
For instance, this function may be used for fast computation of the (inverse)
square root of the shrinkage correlation matrix (needed, e.g., for decorrelation).
crossprod.powcor.shrink efficiently computes R^{α} y without
actually computing the full matrix R^{α}.
a matrix, the number of rows of y must be the same as the number of columns of x
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)).
verbose
output status while computing (default: TRUE)
Details
This function employs a special matrix identity to speed up the computation of
the matrix power of the shrinkage correlation matrix (see Zuber and Strimmer 2009 for details).
Apart from a scaling factor the shrinkage correlation matrix computed
by cor.shrink takes
on the form
Z = I_p + V M V^T ,
where V M V^T is a multiple of the empirical correlation matrix.
Crucially, Z is a matrix of size p times p
whereas M is a potentially much smaller matrix of size m times m,
where m is the rank of the empirical correlation matrix.
In order to calculate the alpha-th power of Z
the function uses the identity
Z^α = I_p - V (I_m -(I_m + M)^α) V^T
requiring only the computation of the alpha-th power of the m by m matrix
I_m + M. This trick enables substantial computational savings especially when the number
of observations is much smaller than the number of variables.
Note that the above identity is related but not identical to the Woodbury matrix
identity for inversion of a matrix.
For α=-1 the above identity reduces to
Z^{-1} = I_p - V (I_m -(I_m + M)^{-1}) V^T ,
whereas the Woodbury matrix identity equals
Z^{-1} = I_p - V (I_m + M^{-1})^{-1} V^T .
Value
powcor.shrink returns a matrix of the same size as the correlation matrix R
crossprod.powcor.shrink returns a matrix of the same size as Ry.
Zuber, V., and K. Strimmer. 2009. Gene ranking and biomarker
discovery under correlation. Bioinformatics 25:2700-2707.
(http://arxiv.org/abs/0902.0751)
See Also
invcor.shrink, cor.shrink, mpower.
Examples
# load corpcor library
library("corpcor")
# generate data matrix
p = 500
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)
### fast computation of cross product
y = rnorm(p)
system.time(
(CP1 = crossprod(powcor.shrink(X, alpha=1.23, lambda=lambda), y))
)
system.time(
(CP2 = crossprod.powcor.shrink(X, y, alpha=1.23, lambda=lambda))
)
# no difference
sum((CP1-CP2)^2)