# R/irlba.R In irlba: Fast Truncated Singular Value Decomposition and Principal Components Analysis for Large Dense and Sparse Matrices

#### Documented in irlba

#' Find a few approximate singular values and corresponding
#' singular vectors of a matrix.
#'
#' The augmented implicitly restarted Lanczos bidiagonalization algorithm
#' (IRLBA) finds a few approximate largest (or, optionally, smallest) singular
#' values and corresponding
#' singular vectors of a sparse or dense matrix using a method of Baglama and
#' Reichel.  It is a fast and memory-efficient way to compute a partial SVD.
#'
#' @param A numeric real- or complex-valued matrix or real-valued sparse matrix.
#' @param nv number of right singular vectors to estimate.
#' @param nu number of left singular vectors to estimate (defaults to \code{nv}).
#' @param maxit maximum number of iterations.
#' @param work working subspace dimension, larger values can speed convergence at the cost of more memory use.
#' @param reorth if \code{TRUE}, apply full reorthogonalization to both SVD bases, otherwise
#'   only apply reorthogonalization to the right SVD basis vectors; the latter case is cheaper per
#'   iteration but, overall, may require more iterations for convergence. Automatically \code{TRUE}
#'   when \code{fastpath=TRUE} (see below).
#' @param tol convergence is determined when \eqn{\|A^TU - VS\| < tol\|A\|}{||A^T U - VS|| < tol*||A||},
#'   and when the maximum relative change in estimated singular values from one iteration to the
#'   next is less than \code{svtol = tol} (see \code{svtol} below),
#'   where the spectral norm ||A|| is approximated by the
#'   largest estimated singular value, and U, V, S are the matrices corresponding
#'   to the estimated left and right singular vectors, and diagonal matrix of
#'   estimated singular values, respectively.
#' @param v optional starting vector or output from a previous run of \code{irlba} used
#'   to restart the algorithm from where it left off (see the notes).
#' @param right_only logical value indicating return only the right singular vectors
#'  (\code{TRUE}) or both sets of vectors (\code{FALSE}). The right_only option can be
#'  cheaper to compute and use much less memory when \code{nrow(A) >> ncol(A)} but note
#'  that obtained solutions typically lose accuracy due to lack of re-orthogonalization in the
#'  algorithm and that \code{right_only = TRUE} sets \code{fastpath = FALSE} (only use this option
#'  for really large problems that run out of memory and when \code{nrow(A) >> ncol(A)}).
#'  Consider increasing the \code{work} option to improve accuracy with \code{right_only=TRUE}.
#' @param verbose logical value that when \code{TRUE} prints status messages during the computation.
#' @param scale optional column scaling vector whose values divide each column of \code{A};
#'   must be as long as the number of columns of \code{A} (see notes).
#' @param center optional column centering vector whose values are subtracted from each
#'   column of \code{A}; must be as long as the number of columns of \code{A} and may
#'   not be used together with the deflation options below (see notes).
#' @param shift optional shift value (square matrices only, see notes).
#' @param mult DEPRECATED optional custom matrix multiplication function (default is \code{\%*\%}, see notes).
#' @param fastpath try a fast C algorithm implementation if possible; set \code{fastpath=FALSE} to use the
#'     reference R implementation. See the notes for more details.
#' @param svtol additional stopping tolerance on maximum allowed absolute relative change across each
#' estimated singular value between iterations.
#' The default value of this parameter is to set it to \code{tol}. You can set \code{svtol=Inf} to
#' effectively disable this stopping criterion. Setting \code{svtol=Inf} allows the method to
#' terminate on the first Lanczos iteration if it finds an invariant subspace, but with less certainty
#' that the converged subspace is the desired one. (It may, for instance, miss some of the largest
#' singular values in difficult problems.)
#' @param smallest set \code{smallest=TRUE} to estimate the smallest singular values and associated
#' singular vectors. WARNING: this option is somewhat experimental, and may produce poor
#' estimates for ill-conditioned matrices.
#' @param ... optional additional arguments used to support experimental and deprecated features.
#'
#' @return
#' Returns a list with entries:
#' \describe{
#'   \item{d:}{ max(nu, nv) approximate singular values}
#'   \item{u:}{ nu approximate left singular vectors (only when right_only=FALSE)}
#'   \item{v:}{ nv approximate right singular vectors}
#'   \item{iter:}{ The number of Lanczos iterations carried out}
#'   \item{mprod:}{ The total number of matrix vector products carried out}
#' }
#'
#' @note
#' The syntax of \code{irlba} partially follows \code{svd}, with an important
#' exception. The usual R \code{svd} function always returns a complete set of
#' singular values, even if the number of singular vectors \code{nu} or \code{nv}
#' is set less than the maximum. The \code{irlba} function returns a number of
#' estimated singular values equal to the maximum of the number of specified
#' singular vectors \code{nu} and \code{nv}.
#'
#' Use the optional \code{scale} parameter to implicitly scale each column of
#' the matrix \code{A} by the values in the \code{scale} vector, computing the
#' truncated SVD of the column-scaled \code{sweep(A, 2, scale, FUN=/)}, or
#' equivalently, \code{A \%*\% diag(1 / scale)}, without explicitly forming the
#' scaled matrix. \code{scale} must be a non-zero vector of length equal
#' to the number of columns of \code{A}.
#'
#' Use the optional \code{center} parameter to implicitly subtract the values
#' in the \code{center} vector from each column of \code{A}, computing the
#' truncated SVD of \code{sweep(A, 2, center, FUN=-)},
#' without explicitly forming the centered matrix. \code{center}
#' must be a vector of length equal to the number of columns of \code{A}.
#' This option may be used to efficiently compute principal components without
#' explicitly forming the centered matrix (which can, importantly, preserve
#' sparsity in the matrix). See the examples.
#'
#' The optional \code{shift} scalar valued argument applies only to square matrices; use it
#' to estimate the partial svd of \code{A + diag(shift, nrow(A), nrow(A))}
#' (without explicitly forming the shifted matrix).
#'
#' (Deprecated) Specify an optional alternative matrix multiplication operator in the
#' \code{mult} parameter. \code{mult} must be a function of two arguments,
#' and must handle both cases where one argument is a vector and the other
#' a matrix. This option is deprecated and will be removed in a future version.
#' The new preferred method simply uses R itself to define a custom matrix class
#' with your user-defined matrix multiplication operator. See the examples.
#'
#' Use the \code{v} option to supply a starting vector for the iterative
#' method. A random vector is used by default (precede with \code{set.seed()}
#' for reproducibility). Optionally set \code{v} to
#' the output of a previous run of \code{irlba} to restart the method, adding
#' additional singular values/vectors without recomputing the solution
#' subspace. See the examples.
#'
#' The function may generate the following warnings:
#' \itemize{
#'   \item{"did not converge--results might be invalid!; try increasing work or maxit"
#'   means that the algorithm didn't
#'   converge -- this is potentially a serious problem and the returned results may not be valid. \code{irlba}
#'   reports a warning here instead of an error so that you can inspect whatever is returned. If this
#'   happens, carefully heed the warning and inspect the result. You may also try setting \code{fastpath=FALSE}.}
#'   \item{"You're computing a large percentage of total singular values, standard svd might work better!"
#'     \code{irlba} is designed to efficiently compute a few of the largest singular values and associated
#'      singular vectors of a matrix. The standard \code{svd} function will be more efficient for computing
#'      large numbers of singular values than \code{irlba}.}
#'    \item{"convergence criterion below machine epsilon" means that the product of \code{tol} and the
#'      largest estimated singular value is really small and the normal convergence criterion is only
#'      met up to round off error.}
#' }
#' The function might return an error for several reasons including a situation when the starting
#' vector \code{v} is near the null space of the matrix. In that case, try a different \code{v}.
#'
#' The \code{fastpath=TRUE} option only supports real-valued matrices and sparse matrices
#' of type \code{dgCMatrix} (for now). Other problems fall back to the reference
#' R implementation.
#'
#' @references
#' Baglama, James, and Lothar Reichel. "Augmented implicitly restarted Lanczos bidiagonalization methods." SIAM Journal on Scientific Computing 27.1 (2005): 19-42.
#'
#' @examples
#' set.seed(1)
#'
#' A <- matrix(runif(400), nrow=20)
#' S <- irlba(A, 3)
#' S$d #' #' # Compare with svd #' svd(A)$d[1:3]
#'
#' # Restart the algorithm to compute more singular values
#' # (starting with an existing solution S)
#' S1 <- irlba(A, 5, v=S)
#'
#' # Estimate smallest singular values
#' irlba(A, 3, smallest=TRUE)$d #' #' #Compare with #' tail(svd(A)$d, 3)
#'
#' # Principal components (see also prcomp_irlba)
#' P <- irlba(A, nv=1, center=colMeans(A))
#'
#' # Compare with prcomp and prcomp_irlba (might vary up to sign)
#' cbind(P$v, #' prcomp(A)$rotation[, 1],
#'       prcomp_irlba(A)$rotation[, 1]) #' #' # A custom matrix multiplication function that scales the columns of A #' # (cf the scale option). This function scales the columns of A to unit norm. #' col_scale <- sqrt(apply(A, 2, crossprod)) #' setClass("scaled_matrix", contains="matrix", slots=c(scale="numeric")) #' setMethod("%*%", signature(x="scaled_matrix", y="numeric"), #' function(x ,y) x@.Data %*% (y / x@scale)) #' setMethod("%*%", signature(x="numeric", y="scaled_matrix"), #' function(x ,y) (x %*% y@.Data) / y@scale) #' a <- new("scaled_matrix", A, scale=col_scale) #' irlba(a, 3)$d
#'
#' # Compare with:
#' svd(sweep(A, 2, col_scale, FUN=/))$d[1:3] #' #' #' @seealso \code{\link{svd}}, \code{\link{prcomp}}, \code{\link{partial_eigen}}, \code{\link{svdr}} #' @import Matrix #' @importFrom stats rnorm #' @importFrom methods slotNames #' @useDynLib irlba, .registration=TRUE, .fixes="C_" #' @export irlba <- function(A, # data matrix nv=5, nu=nv, # number of singular vectors to estimate maxit=1000, # maximum number of iterations work=nv + 7, # working subspace size reorth=TRUE, # TRUE=full reorthogonalization tol=1e-5, # stopping tolerance v=NULL, # optional starting vector or restart right_only=FALSE, # TRUE=only return V verbose=FALSE, # display status messages scale=NULL, # optional column scaling center=NULL, # optional column centering shift=NULL, # optional shift for square matrices mult=NULL, # optional custom matrix multiplication func. fastpath=TRUE, # use the faster C implementation if possible svtol=tol, # stopping tolerance percent change in estimated svs smallest=FALSE, # set to TRUE to estimate subspaces associated w/smallest singular values ...) # optional experimental or deprecated arguments { # --------------------------------------------------------------------- # Check input parameters # --------------------------------------------------------------------- ropts <- options(warn=1) # immediately show warnings mflag <- new.env() mflag$flag <- FALSE
on.exit(options(ropts))
interchange <- FALSE
eps <- .Machine$double.eps # hidden support for old, removed (previously deprecated) parameters # this is here as a convenience to keep old code working without change # also supports experimental features not yet promoted to the api mcall <- as.list(match.call()) random <- eval(mcall[["rng"]]) if (is.null(random)) random <- rnorm # default RNG # Maximum number of Ritz vectors to use in augmentation, may be less # depending on workspace size. maxritz <- eval(mcall[["maxritz"]]) # experimental if (is.null(maxritz)) maxritz <- 3 eps2 <- eval(mcall[["invariant_subspace_tolerance"]]) if (is.null(eps2)) eps2 <- eps ^ (4 / 5) du <- eval(mcall[["du"]]) # deprecated dv <- eval(mcall[["dv"]]) # deprecated ds <- eval(mcall[["ds"]]) # deprecated deflate <- is.null(du) + is.null(ds) + is.null(dv) if (is.logical(scale) && ! scale) scale <- NULL if (is.logical(shift) && ! shift) shift <- NULL if (is.logical(center) && ! center) center <- NULL if (smallest) fastpath <- FALSE # for now anyway if (any(dim(A) > 2 ^ 32 - 1)) fastpath <- FALSE # for now if (deflate == 3) { deflate <- FALSE } else if (deflate == 0) { deflate <- TRUE warning("The deflation options have been deprecated. Please modify your code to not use them.") if (length(ds) > 1) stop("deflation limited to one dimension") if (!is.null(dim(du))) du <- du[, 1] if (!is.null(dim(dv))) dv <- dv[, 1] } else stop("all three du ds dv parameters must be specified for deflation") if (!is.null(center)) { if (is.logical(center) && center) center <- colMeans(A) if (deflate) stop("the center parameter can't be specified together with deflation parameters") if (length(center) != ncol(A)) stop("center must be a vector of length ncol(A)") if (fastpath && ! right_only) du <- NULL else du <- 1 ds <- 1 dv <- center deflate <- TRUE } if ("integer" == typeof(A)) A <- A + 0.0 iscomplex <- is.complex(A) m <- nrow(A) n <- ncol(A) if (is.null(nu)) nu <- nv if (!is.null(mult) && deflate) stop("the mult parameter can't be specified together with deflation parameters") missingmult <- FALSE if (is.null(mult)) { missingmult <- TRUE mult <- %*% } k <- max(nu, nv) if (k <= 0) stop("max(nu, nv) must be positive") if (k > min(m - 1, n - 1)) stop("max(nu, nv) must be strictly less than min(nrow(A), ncol(A))") if (k >= 0.5 * min(m, n)) { warning("You're computing too large a percentage of total singular values, use a standard svd instead.") } if (work <= 1) stop("work must be greater than 1") if (tol < 0) stop("tol must be non-negative") if (maxit <= 0) stop("maxit must be positive") # work must be strictly larger than requested subspace dimension, except see right_only below if (work <= k && ! right_only) work <- k + 1 if (work >= min(n, m)) { work <- min(n, m) if (work <= k) { k <- work - 1 # the best we can do! Need to reduce output subspace dimension warning("Requested subspace dimension too large! Reduced to ", k) } } k_org <- k w_dim <- work if (right_only) { w_dim <- 1 fastpath <- FALSE } if (n > m && smallest) { # Interchange dimensions m,n so that dim(A'A) = min(m,n) when seeking the # smallest singular values; avoids finding zero-valued smallest singular values. interchange <- TRUE temp <- m m <- n n <- temp } if (verbose) { message("Working dimension size ", work) } # Check for tiny problem, use standard SVD in that case. Make definition of 'tiny' larger? if (min(m, n) < 6) { A <- as.matrix(A) # avoid need to define "+" and "/" for arbitrary matrix types. if (verbose) message("Tiny problem detected, using standard svd function.") if (!is.null(scale)) { A <- sweep(A, 2, scale, "/") dv <- dv / scale # scale the centering vector. } if (!is.null(shift)) A <- A + diag(shift, nrow(A), ncol(A)) if (deflate) { if (is.null(du)) du <- rep(1, nrow(A)) A <- A - (ds * du) %*% t(dv) } s <- svd(A) if (smallest) { return(list(d=tail(s$d, k), u=s$u[, tail(seq(ncol(s$u)), k), drop=FALSE],
v=s$v[, tail(seq(ncol(s$v), k)), drop=FALSE], iter=0, mprod=0))
}
return(list(d=s$d[1:k], u=s$u[, 1:nu, drop=FALSE],
v=s$v[, 1:nv, drop=FALSE], iter=0, mprod=0)) } # Try to use the fast C-language code path if (deflate) fastpath <- fastpath && is.null(du) # Only matrix, dgCMatrix supported by fastpath fastpath <- fastpath && (("Matrix" %in% attributes(class(A)) && ("dgCMatrix" %in% class(A))) || "matrix" %in% class(A)) if (fastpath && missingmult && !iscomplex && !right_only) { RESTART <- 0L RV <- RW <- RS <- NULL if (is.null(v)) { v <- random(n) if (verbose) message("Initializing starting vector v with samples from standard normal distribution. Use set.seed first for reproducibility.") } else if (is.list(v)) # restarted case { if (is.null(v$v) || is.null(v$d) || is.null(v$u)) stop("restart requires left and right singular vectors")
if (max(nu, nv) <= min(ncol(v$u), ncol(v$v))) return(v) # Nothing to do!
RESTART <- as.integer(length(v$d)) RND <- random(n) RND <- orthog(RND, v$v)
RV <- cbind(v$v, RND / norm2(RND)) RW <- v$u
RS <- v$d v <- NULL } SP <- ifelse(is.matrix(A), 0L, 1L) if (verbose) message("irlba: using fast C implementation") SCALE <- NULL SHIFT <- NULL CENTER <- NULL if (!is.null(scale)) { if (length(scale) != ncol(A)) stop("scale length must match number of matrix columns") SCALE <- as.double(scale) } if (!is.null(shift)) { if (length(shift) != 1) stop("shift length must be 1") SHIFT <- as.double(shift) } if (deflate) { if (length(center) != ncol(A)) stop("the centering vector length must match the number of matrix columns") CENTER <- as.double(center) } ans <- .Call(C_IRLB, A, as.integer(k), as.double(v), as.integer(work), as.integer(maxit), as.double(tol), as.double(eps2), as.integer(SP), as.integer(RESTART), RV, RW, RS, SCALE, SHIFT, CENTER, as.double(svtol)) if (ans[[6]] == 0 || ans[[6]] == -2) { names(ans) <- c("d", "u", "v", "iter", "mprod", "err") ans$u <- matrix(head(ans$u, m * nu), nrow=m, ncol=nu) ans$v <- matrix(head(ans$v, n * nv), nrow=n, ncol=nv) if (tol * ans$d[1] < eps) warning("convergence criterion below machine epsilon")
if (ans[[6]] == -2) warning("did not converge--results might be invalid!; try increasing work or maxit")
return(ans[-6])
}
errors <- c("invalid dimensions",
"didn't converge",
"out of memory",
"starting vector near the null space",
"linear dependency encountered")
erridx <- abs(ans[[6]])
if (erridx > 1)
warning("fast code path error ", errors[erridx], "; re-trying with fastpath=FALSE.", immediate.=TRUE)
}

# Allocate memory for W and F:
W <- matrix(0.0, m, w_dim)
F <- matrix(0.0, n, 1)
restart <- FALSE
if (is.list(v))
{
if (is.null(v$v) || is.null(v$d) || is.null(v$u)) stop("restart requires left and right singular vectors") if (max(nu, nv) <= min(ncol(v$u), ncol(v$v))) return(v) # Nothing to do! right_only <- FALSE W[, 1:ncol(v$u)] <- v$u d <- v$d
V <- matrix(0.0, n, work)
V[, 1:ncol(v$v)] <- v$v
restart <- TRUE
} else if (is.null(v))
{
# If starting matrix v is not given then set V to be an (n x 1) matrix of
# normally distributed random numbers.  In any case, allocate V appropriate to
# problem size:
V <- matrix(0.0, n, work)
V[, 1] <- random(n)
} else
{
# user-supplied starting subspace
V <- matrix(0.0, n, work)
V[1:length(v)] <- v
}

# ---------------------------------------------------------------------
# Initialize local variables
# ---------------------------------------------------------------------
B <- NULL                  # Bidiagonal matrix
Bsz <- NULL                # Size of B
eps23 <- eps ^ (2 / 3)     # Used for Smax/avoids using zero
iter <- 1                  # Man loop iteration count
mprod <- 0                 # Number of matrix-vector products
R_F <- NULL                # 2-norm of residual vector F
sqrteps <- sqrt(eps)       #
Smax <- 1                  # Max value of all computed singular values of
# B est. ||A||_2
Smin <- NULL               # Min value of all computed singular values of
# B est. cond(A)
lastsv <- c()              # estimated sv in last iteration

# Check for user-supplied restart condition
if (restart)
{
B <- cbind(diag(d), 0)
k <- length(d)

F <- random(n)
F <- orthog(F, V[, 1:k])
V[, k + 1] <- F / norm2(F)
}

# Change du to be non-NULL, for non-fastpath'able matrices with non-NULL scale.
if (deflate && is.null(du)) du <- 1

# ---------------------------------------------------------------------
# Main iteration
# ---------------------------------------------------------------------
while (iter <= maxit)
{
# ---------------------------------------------------------------------
# Compute the Lanczos bidiagonal decomposition:
# such that AV  = WB
# and       t(A)W = VB + Ft(E)
# This routine updates W, V, F, B, mprod
#
# Note on scale and center: These options are applied implicitly below
# for maximum computational efficiency. This complicates their application
# somewhat, but saves a few flops.
# ---------------------------------------------------------------------
j <- 1
#   Normalize starting vector:
if (iter == 1 && !restart)
{
V[, 1] <- V[, 1] / norm2(V[, 1])
}
else j <- k + 1
#   j_w is used here to support the right_only=TRUE case.
j_w <- ifelse(w_dim > 1, j, 1)

#   Compute W=AV
#   Optionally apply scale
VJ <- V[, j]
if (!is.null(scale))
{
VJ <- VJ / scale
}
if (interchange) avj <- mult(VJ, A)
else avj <- mult(A, VJ)

#   Handle non-ordinary arrays as products.
W[, j_w] <- as.vector(avj)
mprod <- mprod + 1

#   Optionally apply shift
if (!is.null(shift))
{
W[, j_w] <- W[, j_w] + shift * VJ
}

#   Optionally apply deflation
if (deflate)
{
W[, j_w] <- W[, j_w] - ds * drop(cross(dv, VJ)) * du
}

#   Orthogonalize W
if (iter != 1 && w_dim > 1 && reorth)
{
W[, j] <- orthog(W[, j, drop=FALSE], W[, 1:(j - 1), drop=FALSE])
}

S <- norm2(W[, j_w, drop=FALSE])
#   Check for linearly dependent vectors
if (is.na(S) || S < eps2 && j == 1) stop("starting vector near the null space")
if (is.na(S) || S < eps2)
{
if (verbose) message_once("invariant subspace found", flag=mflag)
W[, j_w] <- random(nrow(W))
if (w_dim > 1) W[, j] <- orthog(W[, j], W[, 1:(j - 1)])
W[, j_w] <- W[, j_w] / norm2(W[, j_w])
S <- 0
}
else W[, j_w] <- W[, j_w] / S

#   Lanczos process
while (j <= work)
{
j_w <- ifelse(w_dim > 1, j, 1)
if (iscomplex)
{
if (interchange) F <- Conj(t(drop(mult(A, Conj(drop(W[, j_w]))))))
else F <- Conj(t(drop(mult(Conj(drop(W[, j_w])), A))))
}
else
{
if (interchange) F <- t(drop(mult(A, drop(W[, j_w]))))
else F <- t(drop(mult(drop(W[, j_w]), A)))
}
#     Optionally apply shift, scale, deflate
if (!is.null(shift)) F <- F + shift * W[, j_w]
if (!is.null(scale)) F <- F / scale
if (deflate) {
sub <- sum(W[, j_w]) * dv
if (!is.null(scale)) sub <- sub / scale
F <- F - sub
}
mprod <- mprod + 1
F <- drop(F - S * V[, j])
#     Orthogonalize
F <- orthog(F, V[, 1:j, drop=FALSE])
if (j + 1 <= work)
{
R <- norm2(F)
#       Check for linear dependence
if (R < eps2)
{
if (verbose) message_once("invariant subspace found", flag=mflag)
F <- matrix(random(dim(V)[1]), dim(V)[1], 1)
F <- orthog(F, V[, 1:j, drop=FALSE])
V[, j + 1] <- F / norm2(F)
R <- 0
}
else V[, j + 1] <- F / R

#       Compute block diagonal matrix
if (is.null(B)) B <- cbind(S, R)
else            B <- rbind(cbind(B, 0), c(rep(0, ncol(B) - 1), S, R))

jp1_w <- ifelse(w_dim > 1, j + 1, 1)
w_old <- W[, j_w]

#       Optionally apply scale
VJP1 <- V[, j + 1]
if (!is.null(scale))
{
VJP1 <- VJP1 / scale
}
if (interchange) W[, jp1_w] <- drop(mult(drop(VJP1), A))
else W[, jp1_w] <- drop(mult(A, drop(VJP1)))
mprod <- mprod + 1

#       Optionally apply shift
if (!is.null(shift))
{
W[, jp1_w] <- W[, jp1_w] + shift * VJP1
}

#       Optionally apply deflation
if (deflate)
{
W[, jp1_w] <- W[, jp1_w] - ds * drop(cross(dv, VJP1)) * du
}

#       One step of the classical Gram-Schmidt process
W[, jp1_w] <- W[, jp1_w] - R * w_old

#       Full reorthogonalization of W
if (reorth && w_dim > 1) W[, j + 1] <- orthog(W[, j + 1], W[, 1:j])
S <- norm2(W[, jp1_w])
#       Check for linear dependence
if (S < eps2)
{
if (verbose) message_once("invariant subspace found", flag=mflag)
W[, jp1_w] <- random(nrow(W))
if (w_dim > 1) W[, j + 1] <- orthog(W[, j + 1], W[, 1:j])
W[, jp1_w] <- W[, jp1_w] / norm2(W[, jp1_w])
S <- 0
}
else W[, jp1_w] <- W[, jp1_w] / S
}
else
{
#       Add a last block to matrix B
B <- rbind(B, c(rep(0, j - 1), S))
}
j <- j + 1
}
# ---------------------------------------------------------------------
# (End of the Lanczos bidiagonalization part)
# ---------------------------------------------------------------------
Bsz <- nrow(B)
R_F <- norm2(F)
F <- F / R_F
#   Compute singular triplets of B, svd must return ordered singular
#   values from largest to smallest.
Bsvd <- svd(B)

#   Estimate ||A|| using the largest singular value over all iterations
#   and estimate the cond(A) using approximations to the largest and
#   smallest singular values. If a small singular value is less than sqrteps
#   require two-sided reorthogonalization.
if (iter == 1)
{
Smax <- Bsvd$d[1] Smin <- Bsvd$d[Bsz]
}
else
{
Smax <- max(Smax, Bsvd$d[1]) Smin <- min(Smin, Bsvd$d[Bsz])
}
Smax <- max(eps23, Smax)
if (! reorth && Smin / Smax < sqrteps)
{
warning("The matrix is ill-conditioned. Basis will be reorthogonalized.")
reorth <- TRUE
}
if (smallest)
{
jj <- seq(ncol(Bsvd$u), 1, by = -1) Bsvd$u <- Bsvd$u[, jj] Bsvd$d <- Bsvd$d[jj] Bsvd$v <- Bsvd$v[, jj] } # Compute the residuals R <- R_F * Bsvd$u[Bsz, , drop=FALSE]
#   Check for convergence
ct <- convtests(Bsz, tol, k_org, Bsvd, abs(R), k, Smax, lastsv, svtol, maxritz, work, S)
if (verbose)
{
message("iter= ", iter,
", mprod= ", mprod,
", sv[", k_org, "]=", sprintf("%.2e", Bsvd$d[k_org]), ", %change=", sprintf("%.2e", (Bsvd$d[k_org] - lastsv[k_org])/Bsvd$d[k_org]), ", k=", ct$k)
}
lastsv <- Bsvd$d k <- ct$k

#   If all desired singular values converged, then exit main loop
if (ct$converged) break if (iter >= maxit) break # Compute the starting vectors and first block of B if (smallest && (Smin / Smax > sqrteps)) { # Update the SVD of B to be the svd of [B ||F||E_m] Bsvd2.d <- Bsvd$d
Bsvd2.d <- diag(Bsvd2.d, nrow=length(Bsvd2.d))
Bsvd2 <- svd(cbind(Bsvd2.d, t(R)))
jj <- seq(ncol(Bsvd2$u), 1, by=-1) Bsvd2$u <- Bsvd2$u[, jj] Bsvd2$d <- Bsvd2$d[jj] Bsvd2$v <- Bsvd2$v[, jj] Bsvd$d <- Bsvd2$d Bsvd$u <- Bsvd$u %*% Bsvd2$u
Bsvd$v <- cbind(rbind(Bsvd$v, rep(0, Bsz)), c(rep(0, Bsz), 1)) %*% Bsvd2$v V_B_last <- Bsvd$v[Bsz + 1, 1:k]
s <- R_F * solve(B, cbind(c(rep(0, Bsz - 1), 1)))
Bsvd$v <- Bsvd$v[1:Bsz, , drop=FALSE] + s %*% Bsvd$v[Bsz + 1, ] qrv <- qr(cbind(rbind(Bsvd$v[, 1:k], 0), rbind(-s, 1)))
Bsvd$v <- qr.Q(qrv) R <- qr.R(qrv) V[, 1:(k + 1)] <- cbind(V, F) %*% Bsvd$v

#  Update and compute the k by k+1 part of B
UT <- t(R[1:(k + 1), 1:k] + R[, k + 1] %*% rbind(V_B_last))
B <- diag(Bsvd$d[1:k], nrow=k) %*% (UT * upper.tri(UT, diag=TRUE))[1:k, 1:(k+1)] } else { # using the Ritz vectors V[, 1:(k + dim(F)[2])] <- cbind(V[, 1:(dim(Bsvd$v)[1]), drop=FALSE] %*% Bsvd$v[, 1:k], F) B <- cbind(diag(Bsvd$d[1:k], nrow=k), R[1:k])
}

#   Update the left approximate singular vectors
if (w_dim > 1)
{
W[, 1:k] <- W[, 1:(dim(Bsvd$u)[1]), drop=FALSE] %*% Bsvd$u[, 1:k]
}

iter <- iter + 1
}
# ---------------------------------------------------------------------
# End of the main iteration loop
# Output results
# ---------------------------------------------------------------------
if (!ct$converged) warning("did not converge--results might be invalid!; try increasing maxit or work") d <- Bsvd$d[1:k_org]
if (!right_only)
{
u <- W[, 1:(dim(Bsvd$u)[1]), drop=FALSE] %*% Bsvd$u[, 1:k_org, drop=FALSE]
}
v <- V[, 1:(dim(Bsvd$v)[1]), drop=FALSE] %*% Bsvd$v[, 1:k_org, drop=FALSE]
if (smallest)
{
reverse <- seq(length(d), 1)
d <- d[reverse]
if (!right_only) u <- u[, reverse]
v <- v[, reverse]
}
if (tol * d[1] < eps) warning("convergence criterion below machine epsilon")
if (right_only)
return(list(d=d, v=v[, 1:nv, drop=FALSE], iter=iter, mprod=mprod))
return(list(d=d, u=u[, 1:nu, drop=FALSE],
v=v[, 1:nv, drop=FALSE], iter=iter, mprod=mprod))
}

## Try the irlba package in your browser

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

irlba documentation built on Oct. 4, 2022, 1:05 a.m.