shrink_ooc.r

pvt.powscor_mod = function(x, alpha, lambda, w)
{
  #### determine correlation shrinkage intensity
  if (missing(lambda))
  {
    lambda = corpcor:::estimate.lambda(x, w, verbose=FALSE)
    lambda.estimated=TRUE
  }
  else
  {
    if (lambda < 0) lambda = 0
    if (lambda > 1) lambda = 1
    lambda.estimated=FALSE
  }
  
  n = nrow(x)
  w = corpcor:::pvt.check.w(w, n)
  xs = corpcor:::wt.scale(x, w, center=TRUE, scale=TRUE) # standardize data matrix

  # bias correction factor
  h1 = 1/(1-sum(w*w))   # for w=1/n this equals the usual h1=n/(n-1)
  p = ncol(xs)
  
  if (lambda == 1 | alpha == 0) # result in both cases is the identity matrix
  {
      # powr = diag(p)    # return identity matrix
      powr = diag(1.0, p, p)
  }
  else
  {
    # unbiased empirical estimator
    # for w=1/n  the following  would simplify to:  r = 1/(n-1)*crossprod(xs)
    #r0 = h1 * t(xs) %*% diag(w) %*% xs
    #r0 = h1 * t(xs) %*% sweep(xs, 1, w, "*") # sweep requires less memory
    xsw = sweep(xs, 1, sqrt(w), "*")
    
    h = hdfmat::hdfmat("/tmp/r0.h5", "crossprod", p, p, type="float")
    h$fill_crossprod(xsw)
    h$scale(h1*(1-lambda))
    h$fill_diag(1)
  }
  
  h
}

cor_shrink = function(x, lambda, w)
{
   # matrix power of shrinkage correlation
   powr = pvt.powscor_mod(x=x, alpha=1, lambda=lambda, w=w)
   powr
}



suppressMessages(library(corpcor))
suppressMessages(library(hdfmat))

m = 100
n = 1000
x = matrix(rnorm(m*n), m, n)

h = cor_shrink(x)
h$eigen(k=3)