users/dim-reduction/ssvd.page/ssvd.R

# standard SSVD
ssvd.svd <- function(x, k, p=25, qiter=0 ) { 

a <- as.matrix(x)
m <- nrow(a)
n <- ncol(a)
p <- min( min(m,n)-k,p)
r <- k+p

omega <- matrix ( rnorm(r*n), nrow=n, ncol=r)

y <- a %*% omega

q <- qr.Q(qr(y))

b<- t(q) %*% a

#power iterations
for ( i in 1:qiter ) { 
  y <- a %*% t(b)
  q <- qr.Q(qr(y))
  b <- t(q) %*% a
}

bbt <- b %*% t(b)

e <- eigen(bbt, symmetric=T)

res <- list()

res$svalues <- sqrt(e$values)[1:k]
uhat=e$vectors[1:k,1:k]

res$u <- (q %*% e$vectors)[,1:k]
res$v <- (t(b) %*% e$vectors %*% diag(1/e$values))[,1:k]

return(res)
}

#SSVD with Q=YR^-1 substitute.
# this is just a simulation, because it is suboptimal to verify the actual result
ssvd.svd1 <- function(x, k, p=25, qiter=0 ) { 

a <- as.matrix(x)
m <- nrow(a)
n <- ncol(a)
p <- min( min(m,n)-k,p)
r <- k+p

omega <- matrix ( rnorm(r*n), nrow=n, ncol=r)

# in reality we of course don't need to form and persist y
# but this is just verification
y <- a %*% omega

yty <- t(y) %*% y
R <- chol(yty, pivot = T)
q <- y %*% solve(R)

b<- t( q ) %*% a   

#power iterations
for ( i in 1:qiter ) { 
  y <- a %*% t(b)

  yty <- t(y) %*% y
  R <- chol(yty, pivot = T)
  q <- y %*% solve(R)
  b <- t(q) %*% a
}

bbt <- b %*% t(b)

e <- eigen(bbt, symmetric=T)

res <- list()

res$svalues <- sqrt(e$values)[1:k]
uhat=e$vectors[1:k,1:k]

res$u <- (q %*% e$vectors)[,1:k]
res$v <- (t(b) %*% e$vectors %*% diag(1/e$values))[,1:k]

return(res)
}


#############
## ssvd with pci options
ssvd.cpca <- function ( x, k, p=25, qiter=0, fixY=T ) { 

a <- as.matrix(x)
m <- nrow(a)
n <- ncol(a)
p <- min( min(m,n)-k,p)
r <- k+p


# compute median xi
xi<-colMeans(a)

omega <- matrix ( rnorm(r*n), nrow=n, ncol=r)

y <- a %*% omega

#fix y
if ( fixY ) { 
  #debug
  cat ("fixing Y...\n");

  s_o = t(omega) %*% cbind(xi)
  for (i in 1:r ) y[,i]<- y[,i]-s_o[i]
}


q <- qr.Q(qr(y))

b<- t(q) %*% a

# compute sum of q rows 
s_q <- cbind(colSums(q))

# compute B*xi
# of course in MR implementation 
# it will be collected as sums of ( B[,i] * xi[i] ) and reduced after.
s_b <- b %*% cbind(xi)


#power iterations
for ( i in 1:qiter ) { 

  # fix b 
  b <- b - s_q %*% rbind(xi) 

  y <- a %*% t(b)

  # fix y 
  if ( fixY )  
    for (i in 1:r ) y[,i]<- y[,i]-s_b[i]
  

  q <- qr.Q(qr(y))
  b <- t(q) %*% a

  # recompute s_{q}
  s_q <- cbind(colSums(q))

  #recompute s_{b}
  s_b <- b %*% cbind(xi)

}



#C is the outer product of S_q and S_b per doc
C <- s_q %*% t(s_b)

# fixing BB'
bbt <- b %*% t(b) -C -t(C) + sum(xi * xi)* (s_q %*% t(s_q))

e <- eigen(bbt, symmetric=T)

res <- list()

res$svalues <- sqrt(e$values)[1:k]
uhat=e$vectors[1:k,1:k]

res$u <- (q %*% e$vectors)[,1:k]

res$v <- (t(b- s_q %*% rbind(xi) ) %*% e$vectors %*% diag(1/e$values))[,1:k]

return(res)

}