#' @title Randomized QB Decomposition (rqb).
#
#' @description Compute the near-optimal QB decomposition of a rectangular matrix.
#
#' @details
#' The randomized QB decomposition factors a rectangular \eqn{(m,n)} matrix \eqn{A} as
#' \eqn{A = Q * B}. \eqn{Q} is an \eqn{(m,k)} matrix with orthogonal columns, and \eqn{B} a \eqn{(k,n)} matrix.
#' The target rank is assumed to be \eqn{k << min(m,n)}.
#'
#' \eqn{p} is an oversampling parameter to improve the approximation.
#' A value between 5 and 10 is recommended, and \eqn{p=10} is set by default.
#'
#' The parameter \eqn{q} specifies the number of power (subspace) iterations
#' to reduce the approximation error. This is recommended
#' if the the singular values decay slowly. In practice 1 or 2 iterations
#' achieve good results, however, computing power iterations increases the
#' computational time. The number of power iterations is set to \eqn{q=2} by default.
#'
#'
#' @param A array_like; \cr
#' real/complex \eqn{(m, n)} input matrix (or data frame).
#'
#' @param k integer, optional; \cr
#' target rank of the low-rank decomposition. It should satisfy \eqn{k << min(m,n)}.
#'
#' @param p integer, optional; \cr
#' oversampling parameter (default \eqn{p=10}).
#'
#' @param q integer, optional; \cr
#' number of power iterations (default \eqn{q=2}).
#'
#' @param sdist string \eqn{c( 'unif', 'normal', 'rademacher')}, optional; \cr
#' specifies the sampling distribution: \cr
#' \eqn{'unif'} : Uniform `[-1,1]`. \cr
#' \eqn{'normal}' (default) : Normal `~N(0,1)`. \cr
#' \eqn{'rademacher'} : Rademacher random variates. \cr
#'
#' @param rand bool, optional; \cr
#' If (\eqn{TRUE}), a probabilistic strategy is used, otherwise a deterministic algorithm is used.
#'
#'
#'@return \code{rqb} returns a list containing the following components:
#' \describe{
#'\item{Q}{ array_like; \cr
#' matrix with orthogonal columns; \eqn{(m, k)} dimensional array.
#'}
#'
#'\item{B}{ array_like; \cr
#' smaller matrix; \eqn{(k, n)} dimensional array.
#'}
#'}
#'
#'
#' @references
#' \itemize{
#' \item [1] N. B. Erichson, S. Voronin, S. L. Brunton and J. N. Kutz. 2019.
#' Randomized Matrix Decompositions Using {R}.
#' Journal of Statistical Software, 89(11), 1-48.
#' \url{http://doi.org/10.18637/jss.v089.i11}.
#'
#' \item [2] 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}).
#'}
#'
#' @author N. Benjamin Erichson, \email{erichson@berkeley.edu}
#'
#' @seealso \code{\link{svd}}
#'
#'
#' @export
rqb <- function(A, k=NULL, p=10, q=2, sdist="normal", rand = TRUE) UseMethod("rqb")
#' @export
rqb.default <- function(A, k=NULL, p=10, q=2, sdist="normal", rand = TRUE) {
A <- as.matrix(A)
# Create rqb object
rqbObj <- list()
#Dim of input matrix
m <- nrow(A)
n <- ncol(A)
#Set target rank
if(is.null(k)) k = n
if(k > n) k <- n
if(is.character(k)) stop("Target rank is not valid!")
if(k < 1) stop("Target rank is not valid!")
#Set oversampling parameter
l <- round(k) + round(p)
if(l > n) l <- n
if(l < 1) stop("Target rank is not valid!")
#Check if array is real or complex
if(is.complex(A)) {
isreal <- FALSE
} else {
isreal <- TRUE
}
if(rand == TRUE) {
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Generate a random sampling matrix O
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
O <- switch(sdist,
normal = matrix(stats::rnorm(l*n), n, l),
unif = matrix(stats::runif(l*n), n, l),
rademacher = matrix(sample(c(-1,1), (l*n), replace = TRUE, prob = c(0.5,0.5)), n, l),
stop("Selected sampling distribution is not supported!"))
if(isreal==FALSE) {
O <- O + switch(sdist,
normal = 1i * matrix(stats::rnorm(l*n), n, l),
unif = 1i * matrix(stats::runif(l*n), n, l),
rademacher = 1i * matrix(sample(c(-1,1), (l*n), replace = TRUE, prob = c(0.5,0.5)), n, l),
stop("Selected sampling distribution is not supported!"))
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Build sample matrix Y : Y = A * O
#Note: Y should approximate the range of A
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Y <- A %*% O
remove(O)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Orthogonalize Y using economic QR decomposition: Y=QR
#If q > 0 perfrom q subspace iterations
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if( q > 0 ) {
for( i in 1:q) {
Y <- qr.Q( qr(Y, complete = FALSE) , complete = FALSE )
Z <- crossprod_help(A , Y )
Z <- qr.Q( qr(Z, complete = FALSE) , complete = FALSE )
Y <- A %*% Z
}#End for
remove(Z)
}#End if
rqbObj$Q <- qr.Q( qr(Y, complete = FALSE) , complete = FALSE )
}else{
rqbObj$Q <- qr.Q( qr(A, complete = FALSE) , complete = FALSE )
}
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Project the data matrix a into a lower dimensional subspace
#B := Q.T * A
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
rqbObj$B <- crossprod_help(rqbObj$Q , A )
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#Return
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
class(rqbObj) <- "rqb"
return(rqbObj)
} # End rqb
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.