# R/ssvd.R In s4vd: Biclustering via Sparse Singular Value Decomposition Incorporating Stability Selection

```ssvd <-
function(X, threu = 1, threv = 1, gamu = 0, gamv =0,  u0, v0,
merr = 1e-4, niter = 100)
{
if(missing(u0) || missing(v0))
{
##
## Use irlba efficient partial svd for large problems
if(min(dim(X)) < 500) S <- svd(X[],nu=1,nv=1)
else S <- irlba(X,nu=1,nv=1)
u0 <- S\$u
v0 <- S\$v
}
i <- NULL #to pass CRAN checks without NOTE
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;
## avoid explicitly forming the transpose of X:
z <- t(t(u0) %*% X)[,drop=FALSE]
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)
# Note: The 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...
#        }
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)  ## Note
}

Iv <- min(which(Bv==min(Bv)))
temp <- sort(c(0, abs(z*winv)))
lv <- temp[d+1-Iv]

temp2 <- 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" )
# Updating u
z <- (X %*% cbind(v1))[,drop=FALSE]      ## Note
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)
#        }
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)
}

Iu <- min(which(Bu==min(Bu)))
temp <- sort(c(0, abs(z*winu)))
lu <- temp[n+1-Iu]
temp2 <- 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))
}
```

## Try the s4vd package in your browser

Any scripts or data that you put into this service are public.

s4vd documentation built on May 2, 2019, 4:47 p.m.