shrink.r

pvt.powscor_mod = function(x, alpha, lambda, w)
{
  #### determine correlation shrinkage intensity
  if (missing(lambda))
  {
    if (comm.rank() == 0)
      lambda = corpcor:::estimate.lambda(x, w, verbose=FALSE)
    else
      lambda = NULL
    
    lambda = pbdMPI::bcast(lambda, rank.source=0)
    
    lambda.estimated=TRUE
  }
  else
  {
    if (lambda < 0) lambda = 0
    if (lambda > 1) lambda = 1
    lambda.estimated=FALSE
  }
  
  if (comm.rank() == 0)
  {
    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
      if (comm.rank() != 0)
        p = NULL
      
      p = pbdMPI::bcast(p, rank.source=0)
      powr = diag(1.0, p, p, type="ddmatrix")
  }
  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
    if (comm.rank() == 0)
      xsw = sweep(xs, 1, sqrt(w), "*")
    else
      xsw = h1 = NULL
    
    h1 = pbdMPI::bcast(h1, rank.source=0)
    
    d_xsw = as.ddmatrix(xsw)
    r0 = h1 * crossprod(d_xsw)
    
    # shrink off-diagonal elements
    powr = (1-lambda)*r0
    diag(powr) = 1
  }
  
  powr
}

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

ev_cor_shrink = function(x, lambda, w)
{
  co_mat = cor_shrink(x=x, lambda=lambda, w=w)
  ev = eigen(co_mat, symmetric=TRUE, only.values=TRUE)$values
  ev
}



suppressMessages(library(corpcor))
suppressMessages(library(pbdDMAT))
init.grid()

.pbd_env$BLDIM = c(16, 16)
.pbd_env$ICTXT = 0

comm.set.seed(1234, diff=TRUE)
m = 100
n = 1000
if (comm.rank() == 0){
  x = matrix(rnorm(m*n), m, n)
} else {
  x = NULL
}

ev = ev_cor_shrink(x)
comm.print(ev)

finalize()