R/P4_ssvd.R

In Paradigm4/s4vdp4: Performance enhancements and SciDB matrix support for the s4vd biclustering package.

`P4_ssvd` =
function(X, threu = 1, threv = 1, gamu = 0, gamv =0,  u0, v0,
                   merr = 1e-4, niter = 100)
{
    if(missing(u0) || missing(v0))
    {
      if("scidb" %in% class(X)  &&
         length(grep("tsvd",scidb:::.scidbenv$ops$name)) > 0)
      {
        S = scidb:::tsvd(X, nu=1)
        u0 = S$u[]
        v0 = S$v[]
      } else
      {
        if(min(dim(X)) < 1000) S = svd(X[],nu=1,nv=1)
        else S = irlba(X,nu=1,nv=1, matmul=s4vdp4:::matvec)
        u0 = S$u[]
        v0 = S$v[]
      }
    }
    n = dim(X)[1]
    d = dim(X)[2]
    stop = FALSE
    ud = 1
    vd = 1
    iter = 0
    SST = sum(X^2)
    while (ud > merr || vd > merr) {
        iter = iter+1
        cat("iter: ",iter,"\n")
        cat("v: ",length(which(v0!=0)),"\n")
        cat("u: ",length(which(u0!=0)),"\n")
        # Updating v
##        z =  t(X)%*% u0;
## We avoid explicitly forming the transpose of X:
        z = t((t(u0) %*% X)[])  ## P4 alt
        winv = abs(z)^gamv
        sigsq = abs(SST - sum(z^2))/(n*d-d)

        tv = sort(c(0, abs(z*winv)))
        rv = sum(tv>0)
        Bv = rep(1,d+1)*Inf
    
#        for (i in 1:rv){
#            lvc  =  tv[d+1-i]
#            temp1 = which(winv!=0)
#            
#            temp2 = s4vd:::thresh(z[temp1], type = threv, delta = lvc/winv[temp1])
#            vc = rep(0,d)
#            vc[temp1] = temp2
##            Bv[i] = sum((X - u0%*%t(vc))^2)/sigsq + i*log(n*d)
#            Bv[i] = sum((z - vc)^2)/sigsq + i*log(n*d)  ## P4
# Note: The P4 modified Bv[i] differs from Bv[i] by a constant, which
# doesn't change the BIC optimization problem. This observation also
# applies to Bu[i] below...
#        }
# The for loop is parallel:
        temp1 = which(winv!=0)
#        Bv[1:rv] = foreach(i=1:rv,.combine=c,.multicombine=TRUE) %dopar% {
#          lvc  =  tv[d+1-i]
#          delta = lvc/winv[temp1]
#          temp2 = sign(z[temp1]) * (abs(z[temp1]) >= delta) * (abs(z) - delta)
#          vc = rep(0,d)
#          vc[temp1] = temp2
#          sum((z - vc)^2)/sigsq + i*log(n*d)  ## P4
#        }
        fbv = function(i) {
          lvc  =  tv[d+1-i]
          delta = lvc/winv[temp1]
          temp2 = sign(z[temp1]) * (abs(z[temp1]) >= delta) * (abs(z) - delta)
          vc = rep(0,d)
          vc[temp1] = temp2
          sum((z - vc)^2)/sigsq + i*log(n*d)  ## P4
        }
        for(i in 1:rv) Bv[i] = fbv(i)


        Iv = min(which(Bv==min(Bv)))
        temp = sort(c(0, abs(z*winv)))
        lv = temp[d+1-Iv]
    
        temp2 = s4vd:::thresh(z[temp1],type = threv, delta = lv/winv[temp1])
        v1 = rep(0,d)
        v1[temp1] = temp2
        v1 = v1/sqrt(sum(v1^2)) #v_new
    
        cat("v1", length(which(v1!=0)) ,"\n" )
        #str(X)
        #str(v1)
        # Updating u
        z = (X %*% cbind(v1))[,drop=FALSE]      # P4
        winu = abs(z)^gamu
        sigsq = abs(SST - sum(z^2))/(n*d-n)

        tu = sort(c(0,abs(z*winu)))
        ru = sum(tu>0)
        Bu = rep(1,n+1)*Inf

#        for (i in 2:ru){
#            luc  =  tu[n+1-i]
#            temp1 = which(winu!=0)
#            temp2 = s4vd:::thresh(z[temp1], type = threu, delta = luc/winu[temp1])
#            uc = rep(0,n)
#            uc[temp1] = temp2
##            Bu[i] = sum((X - uc%*%t(v1))^2)/sigsq + i*log(n*d)
#            Bu[i] = sum((z - uc)^2)/sigsq + i*log(n*d)  ## P4
#        }
        temp1 = which(winu!=0)
#        Bu[1:ru] = foreach(i=1:ru,.combine=c) %dopar% {
#          luc  =  tu[n+1-i]
#          delta = luc/winu[temp1]
#          temp2 = sign(z[temp1])*(abs(z[temp1]) >= delta)*(abs(z[temp1])-delta)
#          uc = rep(0,n)
#          uc[temp1] = temp2
#          sum((z - uc)^2)/sigsq + i*log(n*d)
#        }

#        fbu = function(i) {
#          luc  =  tu[n+1-i]
#          delta = luc/winu[temp1]
#          temp2 = sign(z[temp1])*(abs(z[temp1]) >= delta)*(abs(z[temp1])-delta)
#          uc = rep(0,n)
#          uc[temp1] = temp2
#          sum((z - uc)^2)/sigsq + i*log(n*d)
#        }
#        for(i in 1:ru) Bu[i] = fbu(i)

        zt = z[temp1]
        at = abs(zt)      # vector of doubles
        st = sign(zt)     # vector of doubles
        wt = winu[temp1]  # vector of doubles
        lnd = log(n*d)    # scalar double
        lzt = as.integer(length(zt))
        lz = as.integer(length(z))
        for(i in 1:ru) Bu[i] = .Call("bcssvd", lzt, lz, as.integer(i), lnd, sigsq, temp1, zt, at, st, wt, tu, z,PACKAGE="s4vdp4")
 

        Iu = min(which(Bu==min(Bu)))
        temp = sort(c(0, abs(z*winu)))
        lu = temp[n+1-Iu]
        temp2 = s4vd:::thresh(z[temp1],delta = lu/winu[temp1])
        u1 = rep(0,n)
        u1[temp1] = temp2
        u1 = u1/sqrt(sum(u1^2))

        ud = sqrt(sum((u0-u1)^2))
        vd = sqrt(sum((v0-v1)^2))

        if (iter > niter){
          print("Fail to converge! Increase the niter!")
          stop = TRUE
          break
        }
        
        u0 = u1
        v0 = v1
    }
    return(list(u = u1, v = v1, iter = iter,stop=stop))
}