R/randsccan.R

In neuroconductor/ANTsR: ANTs in R: Quantification Tools for Biomedical Images

#' randomized svd function.
#'
#' uses random matrix to estimate svd results
#'
#' @param A input matrix
#' @param k rank to use
#' @param seed for testing
#' @return randomized svd is output
#' @author Avants BB
#' @references N. Halko, P.G. Martinsson, J. Tropp "Finding structure with randomness: Stochastic algorithms for constructing approximate matrix decompositions" arXiv 0909.4061
#' @keywords svd
#' @examples
#'
#' A <- matrix(rnorm(3000), ncol=50 )
#' k=10
#' dr=randsvd(A,k)$d[1:k]
#' dt=svd(A)$d[1:k]
#' cor(dr,dt)
#' @export randsvd
randsvd <- function( A, k, seed=NA ) {
if ( !is.na(seed) ) set.seed(seed)
#  Let A be a matrix to be analyzed with n rows and m columns, and r be the ranks of a truncated SVD (if you choose r = min(n, m), then this is the original SVD).

# First a random Gaussian matrix G with m rows and r columns is sampled and computes Y = At G.
# Then apply the Gram-Schmidt process to Y so that each column of Y is ortho-normalized.
G = ( replicate( k , rnorm( nrow(A) ) ) )
Y = .gramschmidt( t(A) %*% G )  # this is a orthornomal basis in P (sparse?)

# Then we compute B = AY with n rows and r columns. Although the size of Y is much smaller than that of A, Y holds the informatin of A; that is AYYt = A. Intuitively, the row information is compresed by Y and can be decompressed by Yt
B = A %*% Y    # subject space

# Similarly, we take another random Gaussian matrix P with r rows and r columns, and compute Z = BP. As in the previous case, the columns of Z are ortho-normalized by the Gram-Schmidt process. ZZtt B = B.
P = ( replicate( k , rnorm( k ) ) )  # gauss mat
Z = .gramschmidt( B %*% P )  # small basis again

#  Then compute C = Zt B.
C = t(Z) %*% B

# Finally we compute SVD of C using the traditional SVD solver, and obtain C = USVt where U and V are orthonormal matrices, and S is the diagonal matrix whose entriesa are singular values. Since a matrix C is very small, this time is negligible.
csvd = svd( C )

# Now A is decomposed as A = AYYt = BYt = ZZtBYt = ZCYt = ZUSVtYt

# Both ZU and YV are othornormal, and ZU is the left singular vectors and YV is the right singular vector. S is the diagonal matrix with singular values.
normer<-function( x ) { sum( sqrt( x*x ) ) }
zu = Z %*% csvd$u
nn=apply( zu, MARGIN=2, FUN=normer )
zu=t( t(zu)/nn )
yv = Y %*% csvd$v   # could be made sparse and plugged back in above ...
nn=apply( yv, MARGIN=2, FUN=normer )
yv=t( t(yv)/nn )
return( list(
  d=csvd$d,
  u=zu,
  v=yv
  ) )

}









#' randomized whitening function.
#'
#' uses random matrix to whiten a matrix
#'
#' @param A input matrix
#' @param k rank to use
#' @param seed for testing
#' @return whitened matrix is output
#' @author Avants BB
#' @references N. Halko, P.G. Martinsson, J. Tropp "Finding structure with randomness: Stochastic algorithms for constructing approximate matrix decompositions" arXiv 0909.4061
#' @keywords whiten
#' @examples
#'
#' A <- matrix(rnorm(3000), ncol=50 )
#' k=10
#' Aw=randwhiten(A,k)
#' # (Aw) %*% t(Aw)
#' @export randwhiten
randwhiten <- function( A, k, seed=NA ) {
  if ( !is.na(seed) ) set.seed(seed)
  X <- t(A)
  n.comp=k
  n <- nrow(X)
  p <- ncol(X)
  X <- scale(X, scale = FALSE)
  X <- t(X)
#  V <- X %*% t(X)/n
  G = replicate( k , rnorm( ncol( t(X) ) ) )
  Y = .gramschmidt( X %*% ( t(X) %*% G ) ) # keep this
  rm( G )
  B = X %*% ( t(X)  %*% Y )
  P = replicate( k , rnorm( k ) )
  Z = .gramschmidt( B %*% P ) # keep this
  C = t(Z) %*% B
  rm( B )
  rm( P )
  csvd = svd( C )
  rm( C )
  normer<-function( x ) { sum( sqrt( x*x ) ) }
  zu = Z %*% csvd$u
  nn=apply( zu, MARGIN=2, FUN=normer )
  zu=t( t(zu)/nn )
  yv = Y %*% csvd$v   # could be made sparse and plugged back in above ...
  nn=apply( yv, MARGIN=2, FUN=normer )
  yv=t( t(yv)/nn )
  s = list(
    d=csvd$d,
    u=zu,
    v=yv
    )
  D <- diag(c(1/sqrt(s$d)))
  K <- D %*% t(s$u)
  K <- matrix(K[1:n.comp, ], n.comp, p)
  xw <- K %*% X # lowrank in row
}



#' randomized cca function.
#'
#' uses random matrix to estimate cca results
#'
#' @param x input matrix
#' @param y input matrix
#' @param k rank to use
#' @param seed for testing
#' @return outputs a list containing:
#' \itemize{
#'   \item{whitened: }{low-rank whitened joint matrix}
#'   \item{svd: }{low-rank svd of joint matrix}
#' }
#' @author Avants BB
#' @examples
#' set.seed(13)
#' x <- matrix(rnorm(3000), nrow=50 )
#' y <- x %*% matrix( rnorm( 60*100 ), nrow=60)
#' k=10
#' dr=randcca(x,y,k,1)
#'
#' @export randcca
randcca <- function( x, y, k, seed = NA ) {
  if ( !is.na(seed) ) set.seed(seed)
  xw = randwhiten( x , k )
  yw = randwhiten( y , k )
  G = replicate( k , rnorm( ncol(x) ) )
  Y = .gramschmidt( t(yw) %*% ( xw %*% G ) ) # keep this
  rm( G )
  B = t(xw) %*% ( yw  %*% Y )
  P = replicate( k , rnorm( k ) )
  Z = .gramschmidt( B %*% P ) # keep this
  C = t(Z) %*% B
  rm( B )
  rm( P )
  csvd = svd( C )
  rm( C )
  normer<-function( x ) { sum( sqrt( x*x ) ) }
  zu = Z %*% csvd$u
  nn=apply( zu, MARGIN=2, FUN=normer )
  zu=t( t(zu)/nn )
  yv = Y %*% csvd$v   # could be made sparse and plugged back in above ...
  nn=apply( yv, MARGIN=2, FUN=normer )
  yv=t( t(yv)/nn )
  s=list(
    d=csvd$d,
    u=zu,
    v=yv
    )
  D <- diag(c(1/sqrt(s$d)))
  K <- D %*% t(s$u)
  return(
    list( whitened=( K %*% t(xw) ) %*% yw, svd=s )
    )
}

.gramschmidt <- function( A, tol=1.e-8 ) {
  A.qr <- qr(A)
  return( qr.Q(A.qr) )
    stopifnot(is.numeric(A), is.matrix(A))
    m <- nrow(A)
    n <- ncol(A)
    if (m < n)
        stop("No. of rows of 'A' must be greater or equal no. of colums.")
    Q <- matrix(0, m, n)
    R <- matrix(0, n, n)
    for (k in 1:n) {
        Q[, k] <- A[, k]
        if (k > 1) {
            for (i in 1:(k - 1)) {
                R[i, k] <- t(Q[, i]) %*% Q[, k]
                Q[, k] <- Q[, k] - R[i, k] * Q[, i]
            }
        }
        R[k, k] <- sum(abs(Q[, k])^2)^(1/2.0) # Norm(Q[, k])
        if (abs(R[k, k]) <= tol)
            stop("Matrix 'A' does not have full rank.")
        Q[, k] <- Q[, k]/R[k, k]
    }
    Q # %*% R
}