# R/rid.R In Benli11/rSVD: Randomized Singular Value Decomposition

#### Documented in rid

#' @title  Randomized interpolative decomposition (ID).
#
#' @description Randomized interpolative decomposition.
#
#' @details
#' Algorithm for computing the ID of a rectangular \eqn{(m, n)} matrix \eqn{A}, with target rank
#' \eqn{k << min(m,n)}. The input matrix is factored as
#'
#' \deqn{A = C * Z}
#'
#' using the column pivoted QR decomposition. The factor matrix \eqn{C} is formed as a subset of
#' columns of \eqn{A}, also called the partial column skeleton.
#' If \code{mode='row'}, then the input matrix is factored as
#'
#' \deqn{A = Z * R}
#'
#' using the row pivoted QR decomposition. The factor matrix \eqn{R} is now formed as
#' a subset of rows of \eqn{A}, also called the partial row skeleton.
#' The factor matrix \eqn{Z} contains a \eqn{(k, k)} identity matrix as a submatrix,
#' and is well-conditioned.
#'
#' If \eqn{rand='TRUE'} a probabilistic strategy is used to compute the decomposition, otherwise a
#' deterministic algorithm is used.
#'
#'
#' @param A   array_like; \cr
#'            numeric \eqn{(m, n)} input matrix (or data frame). \cr
#'            If the data contain \eqn{NA}s na.omit is applied.
#'
#' @param k   integer, optional; \cr
#'            number of rows/columns to be selected.
#'            It is required that \eqn{k} is smaller or equal to \eqn{min(m,n)}.
#'
#' @param mode  string c('column', 'row'), optional; \cr
#'              columns or rows ID.
#'
#' @param p    integer, optional; \cr
#'             oversampling parameter (default \eqn{p=10}).
#'
#' @param q    integer, optional. \cr
#'             number of additional power iterations (default \eqn{q=0}).
#'
#' @param idx_only  bool, optional; \cr
#'              if (\eqn{TRUE}), the index set \code{idx} is returned, but not the matrix \code{C} or \code{R}.
#'              This is more memory efficient, when dealing with large-scale data.
#'
#' @param rand  bool, optional; \cr
#'              if (\eqn{TRUE}), a probabilistic strategy is used, otherwise a deterministic algorithm is used.
#'
#'
#'
#' @return \code{rid} returns a list containing the following components:
#'    \item{C}{ array_like; \cr
#'              column subset \eqn{C = A[,idx]}, if \code{mode='column'}; array with dimensions \eqn{(m, k)}.
#'    }
#'
#'    \item{R}{ array_like; \cr
#'              row subset \eqn{R = A[idx, ]}, if \code{mode='row'}; array with dimensions \eqn{(k, n)}.
#'    }
#'
#'    \item{Z}{ array_like; \cr
#'             well conditioned matrix; Depending on the selected mode, this is an
#'             array with dimensions \eqn{(k,n)} or \eqn{(m,k)}.
#'    }
#'
#'    \item{idx}{ array_like; \cr
#'               index set of the \eqn{k} selected columns or rows used to form \eqn{C} or \eqn{R}.
#'    }
#'
#'    \item{pivot}{ array_like; \cr
#'                information on the pivoting strategy used during the decomposition.
#'    }
#'
#'    \item{scores}{ array_like; \cr
#'                  scores of the columns or rows of the input matrix \eqn{A}.
#'    }
#'
#'    \item{scores.idx}{ array_like; \cr
#'                  scores of the \eqn{k} selected columns or rows in \eqn{C} or \eqn{R}.
#'    }
#'
#'
#' @references
#' \itemize{
#'   \item  [1] N. Halko, P. Martinsson, and J. Tropp.
#'          "Finding structure with randomness: probabilistic
#'          algorithms for constructing approximate matrix
#'          decompositions" (2009).
#'          (available at arXiv \url{http://arxiv.org/abs/0909.4061}).
#'   \item  [2] N. B. Erichson, S. Voronin, S. Brunton, J. N. Kutz.
#'          "Randomized matrix decompositions using R" (2016).
#'          (available at arXiv \url{http://arxiv.org/abs/1608.02148}).
#' }
#'
#' @author N. Benjamin Erichson, \email{[email protected]}
#'
#'
#'
#' @examples
#' \dontrun{
#' data("tiger")
#'
#' # Compute (column) randomized interpolative decompsition
#' # Note that the image needs to be transposed for correct plotting
#' out <- rid(t(tiger), k = 150)
#'
#' # Show selected columns
#' tiger.partial <- matrix(0, 1200, 1600)
#' tiger.partial[,out$idx] <- t(tiger)[,out$idx]
#' image(t(tiger.partial), col = gray((0:255)/255), useRaster = TRUE)
#'
#' # Reconstruct image
#' tiger.re <- t(out$C %*% out$Z)
#'
#' # Compute relative error
#' print(norm(tiger-tiger.re, 'F') / norm(tiger, 'F'))
#'
#' # Plot approximated image
#' image(tiger.re, col = gray((0:255)/255))
#' }

#' @export
rid <- function(A, k = NULL, mode = 'column', p = 10, q = 0, idx_only = FALSE, rand = TRUE) UseMethod("rid")

#' @export
rid.default <- function(A, k = NULL, mode = 'column', p = 10, q = 0, idx_only = FALSE, rand = TRUE) {

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Checks
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (any(is.na(A))) {
warning("Missing values are omitted: na.omit(A).")
A <- stats::na.omit(A)
}

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Init id object
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
idObj = list( C = NULL,
R = NULL,
Z = NULL,
idx = NULL,
pivot = NULL,
scores = NULL,
scores.idx = NULL,
mode = mode,
rand = rand)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Transpose input matrix if mode == 'row'
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if(mode == 'row') A = H(A)

m <- nrow(A)
n <- ncol(A)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Set target rank
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if(is.null(k)) k <- n
if(k > min(m,n)) k <- min(m,n)
if(k < 1) stop("Target rank is not valid!")

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Compute interpolative decompositon
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Pivoted QR decomposition
if(rand == TRUE) {
out_rqb <- rqb(A, k = k, p = p, q = q)
out <- qr(out_rqb$B, LAPACK=TRUE) } else { out <- qr(A, LAPACK=TRUE) } # Get index set idObj$idx <- out$pivot[1:k] # Get row set idObj$pivot <- out$pivot # Get row set ordered.pivot <- order(idObj$pivot)

# Get R =: S
S <- qr.R( out )

# Compute scores
idObj$scores <- abs(diag(S)) idObj$scores.idx <- idObj$scores[1:k] idObj$scores <- idObj$scores[ordered.pivot] # Compute Z if(k == n) { V = pinv(S[1:k, 1:k]) }else{ V = pinv(S[1:k, 1:k]) %*% S[1:k, (k+1):n] } idObj$Z = cbind(diag(k), V)
idObj$Z = matrix(idObj$Z[, ordered.pivot], nrow = k, ncol = n)

if(mode == 'row') {
idObj$Z = H(idObj$Z)

}

# Create column / row subset
if(idx_only == FALSE) {
if(mode == 'column') {
idObj$C = matrix(A[, idObj$idx], nrow = m, ncol = k)
colnames(idObj$C) <- colnames(A)[idObj$idx]
rownames(idObj$C) <- rownames(A) } if(mode == 'row') { idObj$R = H(matrix(A[, idObj$idx], nrow = m, ncol = k)) rownames(idObj$R) <- colnames(A)[idObj$idx] colnames(idObj$R) <- rownames(A)

}
}

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Return
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
class(idObj) <- "rid"
return( idObj )
}
`
Benli11/rSVD documentation built on Nov. 6, 2018, 10:46 p.m.