[med-svn] [r-cran-irlba] 03/05: New upstream version 2.2.1
Andreas Tille
tille at debian.org
Tue Oct 10 15:19:13 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-irlba.
commit 83dee92387cd47088327941677eed220f30b79df
Author: Andreas Tille <tille at debian.org>
Date: Tue Oct 10 17:17:01 2017 +0200
New upstream version 2.2.1
---
DESCRIPTION | 19 ++-
MD5 | 41 +++---
NAMESPACE | 3 +-
R/eigen.R | 3 +-
R/irlba.R | 301 +++++++++++++++++++++++++++---------------
R/prcomp.R | 4 +-
R/svdr.R | 118 +++++++++++++++++
R/utility.R | 42 +++---
README.md | 106 +++++++++++++--
build/vignette.rds | Bin 199 -> 197 bytes
demo/00Index | 1 -
demo/custom_matrix_multiply.R | 44 ------
inst/doc/irlba.Rnw | 2 +-
inst/doc/irlba.pdf | Bin 181643 -> 181657 bytes
man/irlba.Rd | 93 +++++++------
man/partial_eigen.Rd | 3 +-
man/svdr.Rd | 103 +++++++++++++++
src/irlb.c | 67 ++++++----
src/irlb.h | 16 ++-
src/utility.c | 21 +--
tests/edge.R | 127 ++++++++++--------
tests/svdr.R | 50 +++++++
tests/test.R | 53 ++++++++
vignettes/irlba.Rnw | 2 +-
24 files changed, 867 insertions(+), 352 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index d507e0d..35f8396 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,27 +1,26 @@
Package: irlba
Type: Package
-Title: Fast Truncated SVD, PCA and Symmetric Eigendecomposition for
- Large Dense and Sparse Matrices
-Version: 2.1.2
-Date: 2016-09-20
+Title: Fast Truncated Singular Value Decomposition and Principal
+ Components Analysis for Large Dense and Sparse Matrices
+Version: 2.2.1
+Date: 2017-05-16
Authors at R: c(
person("Jim", "Baglama", rol=c("aut", "cph"), email="jbaglama at uri.edu"),
person("Lothar", "Reichel", role=c("aut", "cph"), email="reichel at math.kent.edu"),
person("B. W.", "Lewis", role=c("aut","cre","cph"), email="blewis at illposed.net"))
-Description: Fast and memory efficient methods for truncated singular and
- eigenvalue decompositions and principal component analysis of large sparse or
- dense matrices.
+Description: Fast and memory efficient methods for truncated singular value
+ decomposition and principal components analysis of large sparse and dense matrices.
Depends: Matrix
LinkingTo: Matrix
Imports: stats, methods
License: GPL-3
BugReports: https://github.com/bwlewis/irlba/issues
-RoxygenNote: 5.0.0
+RoxygenNote: 5.0.1
NeedsCompilation: yes
-Packaged: 2016-09-21 13:57:09 UTC; blewis
+Packaged: 2017-05-16 20:24:46 UTC; blewis
Author: Jim Baglama [aut, cph],
Lothar Reichel [aut, cph],
B. W. Lewis [aut, cre, cph]
Maintainer: B. W. Lewis <blewis at illposed.net>
Repository: CRAN
-Date/Publication: 2016-09-21 19:27:40
+Date/Publication: 2017-05-17 07:28:01 UTC
diff --git a/MD5 b/MD5
index f6915e5..9100742 100644
--- a/MD5
+++ b/MD5
@@ -1,22 +1,23 @@
-1a4cc9a6d803b51099c7e8019203340d *DESCRIPTION
-a92cc95a84584b03f782a11a0f14ccac *NAMESPACE
-a5caf177747baf3939d8da5a7a65f20d *R/eigen.R
-e2e2bd1fc7d563008881ce2f7c02047e *R/irlba.R
-261b601f66e15a95762cfcb25ff7bc52 *R/prcomp.R
-a3aeb39b2c56f9626fcc500a2efbff3c *R/utility.R
-b370501f642fc4e79215983fd2851f20 *README.md
-e5726491a241e7d7cd30836a46ef2ee9 *build/vignette.rds
-f7a30b7455d7c846dd5e4baf3a7a2963 *demo/00Index
-fa21422f4cbda8c9083f67e47421f99c *demo/custom_matrix_multiply.R
-e9b9777c575e2b56f2c85aa758ba4779 *inst/doc/irlba.Rnw
-4122e6a46391d621cb23bf6635653e3b *inst/doc/irlba.pdf
-4285679ab7678c954d2c0ad0e63417b3 *man/irlba.Rd
-a30fddf6e2d613164a96320d0cbc22f0 *man/partial_eigen.Rd
+0bc045318da221e318dc9eca2f2c614a *DESCRIPTION
+2df3ba50b634767e4f69a4a9d9520166 *NAMESPACE
+46c031e01c7cb76685cb4ef0837be3b4 *R/eigen.R
+0b00fde03142679e5ff0b9d74d8c648f *R/irlba.R
+da63d5e8855861a4d27d79800bd4c8b0 *R/prcomp.R
+6a74c49b8e03319481133c72d3f4de08 *R/svdr.R
+20d6b2a573231dad638a6ed961cbcd40 *R/utility.R
+063e4582d2787b591a11d396b19189f7 *README.md
+adf849b0f1d127b206a67d3a451b64cf *build/vignette.rds
+751d1231241f7229377a7f3e066d1e60 *inst/doc/irlba.Rnw
+c9b5ff4a2272f0174c384d1f8c10045f *inst/doc/irlba.pdf
+3001af394110e4a3afb468f87dfa17b9 *man/irlba.Rd
+4c2aaeb543e9cd98eb6203fada578ce3 *man/partial_eigen.Rd
3e425d237604f5286d6db3c9d32435ba *man/prcomp_irlba.Rd
+8f38539535f6102529a42e0a40ab423f *man/svdr.Rd
619d13ce6900fca2139e2139408d39ad *src/Makevars
-fd6a9dd38f353ec10da475f55d662314 *src/irlb.c
-27b357884a2232d31d8d9f5390391d08 *src/irlb.h
-aba4af25c20becc25598407f487abac9 *src/utility.c
-98df9161ed239100dcbb95109a0cf0ac *tests/edge.R
-26d6d03bbcfa4c377d475d392f74da30 *tests/test.R
-e9b9777c575e2b56f2c85aa758ba4779 *vignettes/irlba.Rnw
+68674c958bdb72d325877006f974ce72 *src/irlb.c
+c3060f4abbe417e55fbed6f15af4ca62 *src/irlb.h
+07f335c271e72b7f0083588c9e214686 *src/utility.c
+7a09897e0456484e488df2f46591b797 *tests/edge.R
+9db2b0b6fd5ec7d6776e6b4bc30fbfe1 *tests/svdr.R
+e106ea2fb269fe4f63ed116b64f3c340 *tests/test.R
+751d1231241f7229377a7f3e066d1e60 *vignettes/irlba.Rnw
diff --git a/NAMESPACE b/NAMESPACE
index 1842492..42d4fe6 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -3,10 +3,11 @@
export(irlba)
export(partial_eigen)
export(prcomp_irlba)
+export(svdr)
import(Matrix)
importFrom(methods,slot)
importFrom(methods,slotNames)
importFrom(stats,prcomp)
importFrom(stats,rnorm)
importFrom(stats,sd)
-useDynLib(irlba)
+useDynLib(irlba, .registration=TRUE, .fixes="C_")
diff --git a/R/eigen.R b/R/eigen.R
index 287804f..0865e91 100644
--- a/R/eigen.R
+++ b/R/eigen.R
@@ -7,7 +7,8 @@
#' @param x numeric real-valued dense or sparse matrix.
#' @param n number of largest eigenvalues and corresponding eigenvectors to compute.
#' @param symmetric \code{TRUE} indicates \code{x} is a symmetric matrix (the default);
-#' specify \code{symmetric=FALSE} to compute the largest eigenvalues and corresponding eigenvectors of \code{t(x) \%*\% x} instead.
+#' specify \code{symmetric=FALSE} to compute the largest eigenvalues and corresponding
+#' eigenvectors of \code{t(x) \%*\% x} instead.
#' @param ... optional additional parameters passed to the \code{irlba} function.
#'
#' @return
diff --git a/R/irlba.R b/R/irlba.R
index 5371c6f..22fd0d8 100644
--- a/R/irlba.R
+++ b/R/irlba.R
@@ -1,8 +1,9 @@
-#' Find a few approximate largest singular values and corresponding
+#' 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 singular values and corresponding
+#' (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.
#'
@@ -13,9 +14,11 @@
#' @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. Always set to \code{TRUE}
+#' 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{\|AV - US\| < tol\|A\|}{||AV - US|| < tol*||A||},
+#' @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
@@ -26,19 +29,28 @@
#' (\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 \code{right_only = TRUE} sets \code{fastpath = FALSE} (only use this option
-#' for really large problems that run out of memory).
+#' for really large problems that run out of memory and when \code{nrow(A) >> ncol(A)}).
#' @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 du DEPRECATED optional subspace deflation vector (see notes).
-#' @param ds DEPRECATED optional subspace deflation scalar (see notes).
-#' @param dv DEPRECATED optional subspace deflation vector (see notes).
+#' not be used together with the deflation options below (see notes).
#' @param shift optional shift value (square matrices only, see notes).
-#' @param mult 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 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:
@@ -68,33 +80,25 @@
#' 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. This option may not be
-#' used together with the general rank 1 deflation options. \code{center}
+#' 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 deflation parameters are deprecated and will be removed in
-#' a future version. They could be used to compute the rank-one deflated
-#' SVD of \eqn{A - ds \cdot du dv^T}{A - ds*du \%*\% t(dv)}, where
-#' \eqn{du^T A - ds\cdot dv^T = 0}{t(du) \%*\% A - ds * t(dv) == 0}. For
-#' example, the triple \code{ds, du, dv} may be a known singular value
-#' and corresponding singular vectors. Or \code{ds=m} and \code{dv}
-#' and \code{du} represent a vector of column means of \code{A} and of ones,
-#' respectively, where \code{m} is the number of rows of \code{A}.
-#' This functionality can be effectively replaced with custom matrix
-#' product functions.
+#' 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).
#'
-#' Specify an optional alternative matrix multiplication operator in the
+#' (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. See the examples and
-#' \code{demo("custom_matrix_multiply", package="irlba")} for an
-#' alternative approach.
+#' 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()} to
+#' 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
@@ -138,6 +142,12 @@
#' # (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))
#'
@@ -148,6 +158,7 @@
#'
#' # 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.
+#' # This approach is deprecated (see below for a bettwe way to do this).
#' col_scale <- sqrt(apply(A, 2, crossprod))
#' mult <- function(x, y)
#' {
@@ -162,53 +173,73 @@
#' irlba(A, 3, mult=mult)$d
#'
#' # Compare with:
-#' irlba(A, 3, scale=col_scale)$d
-#'
-#' # Compare with:
#' svd(sweep(A, 2, col_scale, FUN=`/`))$d[1:3]
#'
-#' @seealso \code{\link{svd}}, \code{\link{prcomp}}, \code{\link{partial_eigen}}
+#' # Compare with the new recommended approach:
+#' setClass("scaled_matrix", contains="matrix", slots=c(scale="numeric"))
+#' setMethod("%*%", signature(x="scaled_matrix", y="numeric"),
+#' function(x ,y) x at .Data %*% (y / x at scale))
+#' setMethod("%*%", signature(x="numeric", y="scaled_matrix"),
+#' function(x ,y) (x %*% y at .Data) / y at scale)
+#' a <- new("scaled_matrix", A, scale=col_scale)
+#' irlba(a, 3)$d
+#'
+#' @seealso \code{\link{svd}}, \code{\link{prcomp}}, \code{\link{partial_eigen}}, \code{\link{svdr}}
#' @import Matrix
#' @importFrom stats rnorm
#' @importFrom methods slotNames
-#' @useDynLib irlba
+#' @useDynLib irlba, .registration=TRUE, .fixes="C_"
#' @export
irlba <-
-function (A, # data matrix
- nv=5, nu, # 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, # optional column scaling
- center, # optional column centering
- du, ds, dv, # optional general rank 1 deflation
- shift, # optional shift for square matrices
- mult, # optional custom matrix multiplication func.
- fastpath=TRUE) # use the faster C implementation if possible
+function(A, # data matrix
+ nv=5, nu=nv, # number of singular vectors to estimate
+ maxit=100, # 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 arguments (really just to support old removed args)
{
# ---------------------------------------------------------------------
# Check input parameters
# ---------------------------------------------------------------------
ropts <- options(warn=1) # immediately show warnings
on.exit(options(ropts)) # reset on exit
+ interchange <- FALSE
eps <- .Machine$double.eps
- deflate <- missing(du) + missing(ds) + missing(dv)
+ # hidden support for old, removed (previously deprecated) parameters
+ # this is here as a convenience to keep old code working without change
+ mcall <- as.list(match.call())
+ # Maximum number of Ritz vectors to use in augmentation, may be less
+ # depending on workspace size.
+ maxritz <- mcall[["maxritz"]]
+ if (is.null(maxritz)) maxritz <- 3
+ du <- mcall[["du"]]
+ dv <- mcall[["dv"]]
+ ds <- mcall[["ds"]]
+ deflate <- is.null(du) + is.null(ds) + is.null(dv)
+ if (smallest) fastpath <- FALSE
if (deflate == 3)
{
deflate <- FALSE
} else if (deflate == 0)
{
deflate <- TRUE
- warning("The deflation options are deprecated and will be removed in a future version.")
+ warning("The deflation options are deprecated and have been removed as formal arugments. 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 (!missing(center))
+ if (!is.null(center))
{
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)")
@@ -221,10 +252,10 @@ function (A, # data matrix
iscomplex <- is.complex(A)
m <- nrow(A)
n <- ncol(A)
- if (missing(nu)) nu <- nv
- if (!missing(mult) && deflate) stop("the mult parameter can't be specified together with deflation parameters")
+ 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 (missing(mult))
+ if (is.null(mult))
{
missingmult <- TRUE
mult <- `%*%`
@@ -234,7 +265,7 @@ function (A, # data matrix
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 a large percentage of total singular values, standard svd might work better!")
+ 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")
@@ -254,45 +285,51 @@ function (A, # data matrix
if (right_only)
{
w_dim <- 1
- work <- min(min(m, n), work + 20 ) # typically need this to help convergence
+ work <- min(min(m, n), work + 20) # typically need this to help convergence
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)
+ 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)
{
- if (verbose) message ("Tiny problem detected, using standard `svd` function.")
- if (!missing(scale)) A <- A / scale
- if (!missing(shift)) A <- A + diag(shift)
+ if (verbose) message("Tiny problem detected, using standard `svd` function.")
+ if (!is.null(scale)) A <- A / scale
+ 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 dgCMatrix supported by fastpath for now
- if ("Matrix" %in% attributes(class(A)) && ! ("dgCMatrix" %in% class(A)))
- {
- fastpath <- FALSE
- }
-# Check for custom class
- if ("matrix" %in% attributes(A)$.S3Class && ! ("matrix" %in% class(A)))
- {
- fastpath <- FALSE
- }
+# 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 <- 0
+ RESTART <- 0L
RV <- RW <- RS <- NULL
if (is.null(v))
{
@@ -316,12 +353,12 @@ function (A, # data matrix
SCALE <- NULL
SHIFT <- NULL
CENTER <- NULL
- if (!missing(scale))
+ if (!is.null(scale))
{
if (length(scale) != ncol(A)) stop("scale length must mactch number of matrix columns")
SCALE <- as.double(scale)
}
- if (!missing(shift))
+ if (!is.null(shift))
{
if (length(shift) != 1) stop("shift length must be 1")
SHIFT <- as.double(shift)
@@ -331,15 +368,15 @@ function (A, # data matrix
if (length(center) != ncol(A)) stop("the centering vector length must match the number of matrix columns")
CENTER <- as.double(center)
}
- ans <- .Call("IRLB", A, as.integer(k), as.double(v), as.integer(work),
- as.integer(maxit), as.double(tol), .Machine$double.eps, as.integer(SP),
- RESTART, RV, RW, RS, SCALE, SHIFT, CENTER, PACKAGE="irlba")
+ ans <- .Call(C_IRLB, A, as.integer(k), as.double(v), as.integer(work),
+ as.integer(maxit), as.double(tol), as.double(.Machine$double.eps), 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 (tol * ans$d[1] < eps) warning("convergence criterion below machine epsilon")
if (ans[[6]] == -2) warning("did not converge--results might be invlaid!; try increasing maxit or fastpath=FALSE")
return(ans[-6])
}
@@ -396,7 +433,7 @@ function (A, # data matrix
# B est. ||A||_2
Smin <- NULL # Min value of all computed singular values of
# B est. cond(A)
- SVTol <- tol # Tolerance for singular vector convergence
+ lastsv <- c() # estimated sv in last iteration
# Check for user-supplied restart condition
if (restart)
@@ -437,12 +474,13 @@ function (A, # data matrix
# Compute W=AV
# Optionally apply scale
VJ <- V[, j]
- if (!missing(scale))
+ if (!is.null(scale))
{
VJ <- VJ / scale
}
+ if (interchange) avj <- mult(VJ, A)
+ else avj <- mult(A, VJ)
# Handle sparse products.
- avj <- mult(A, VJ)
if ("Matrix" %in% attributes(class(avj)) && "x" %in% slotNames(avj))
{
if (length(avj at x) == nrow(W)) avj <- slot(avj, "x")
@@ -452,7 +490,7 @@ function (A, # data matrix
mprod <- mprod + 1
# Optionally apply shift
- if (!missing(shift))
+ if (!is.null(shift))
{
W[, j_w] <- W[, j_w] + shift * VJ
}
@@ -466,7 +504,7 @@ function (A, # data matrix
# Orthogonalize W
if (iter != 1 && w_dim > 1 && reorth)
{
- W[, j] <- orthog (W[, j, drop=FALSE], W[, 1:(j - 1), drop=FALSE])
+ W[, j] <- orthog(W[, j, drop=FALSE], W[, 1:(j - 1), drop=FALSE])
}
S <- norm2(W[, j_w, drop=FALSE])
@@ -474,6 +512,7 @@ function (A, # data matrix
if (is.na(S) || S < eps2 && j == 1) stop("starting vector near the null space")
if (is.na(S) || S < eps2)
{
+ if (verbose) message("invariant subspace found")
W[, j_w] <- rnorm(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])
@@ -487,12 +526,17 @@ function (A, # data matrix
j_w <- ifelse(w_dim > 1, j, 1)
if (iscomplex)
{
- F <- Conj(t(drop(mult(Conj(drop(W[, j_w])), A))))
+ 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)))
}
- else F <- t(drop(mult(drop(W[, j_w]), A)))
# Optionally apply shift and scale
- if (!missing(shift)) F <- F + shift * W[, j_w]
- if (!missing(scale)) F <- F / scale
+ if (!is.null(shift)) F <- F + shift * W[, j_w]
+ if (!is.null(scale)) F <- F / scale
mprod <- mprod + 1
F <- drop(F - S * V[, j])
# Orthogonalize
@@ -503,6 +547,7 @@ function (A, # data matrix
# Check for linear dependence
if (R < eps2)
{
+ if (verbose) message("invariant subspace found")
F <- matrix(rnorm(dim(V)[1]), dim(V)[1], 1)
F <- orthog(F, V[, 1:j, drop=FALSE])
V[, j + 1] <- F / norm2(F)
@@ -519,15 +564,16 @@ function (A, # data matrix
# Optionally apply scale
VJP1 <- V[, j + 1]
- if (!missing(scale))
+ if (!is.null(scale))
{
VJP1 <- VJP1 / scale
}
- W[, jp1_w] <- drop(mult(A, drop(VJP1)))
+ 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 (!missing(shift))
+ if (!is.null(shift))
{
W[, jp1_w] <- W[, jp1_w] + shift * VJP1
}
@@ -547,6 +593,7 @@ function (A, # data matrix
# Check for linear dependence
if (S < eps2)
{
+ if (verbose) message("invariant subspace found")
W[, jp1_w] <- rnorm(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])
@@ -561,18 +608,14 @@ function (A, # data matrix
}
j <- j + 1
}
- if (verbose)
- {
- message("Lanczos iter = ", iter, ", dim = ", j - 1, ", mprod = ", mprod)
- }
# ---------------------------------------------------------------------
# (End of the Lanczos bidiagonalization part)
# ---------------------------------------------------------------------
Bsz <- nrow(B)
R_F <- norm2(F)
F <- F / R_F
-# Compute singular triplets of B. Expect svd to return s.v.s in order
-# from largest to smallest.
+# 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
@@ -590,28 +633,71 @@ function (A, # data matrix
Smin <- min(Smin, Bsvd$d[Bsz])
}
Smax <- max(eps23, Smax)
- if (Smin / Smax < sqrteps && !reorth)
+ 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$u,
- Bsvd$d, Bsvd$v, abs(R), k, SVTol, Smax)
+ 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[1:k, 1:(k+1), drop=FALSE]
+# 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])
+ 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)
@@ -625,14 +711,21 @@ function (A, # data matrix
# End of the main iteration loop
# Output results
# ---------------------------------------------------------------------
- if (iter > maxit) warning("did not converge--results might be invalid!; try increasing maxit")
+ 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(tol * d[1] < eps) warning("convergence criterion below machine epsilon")
+ 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],
diff --git a/R/prcomp.R b/R/prcomp.R
index 33427fa..5a11207 100644
--- a/R/prcomp.R
+++ b/R/prcomp.R
@@ -61,7 +61,7 @@
#' @importFrom stats rnorm prcomp sd
#' @importFrom methods slotNames slot
#' @export
-prcomp_irlba <- function (x, n = 3, retx = TRUE, center = TRUE, scale. = FALSE, ...)
+prcomp_irlba <- function(x, n = 3, retx = TRUE, center = TRUE, scale. = FALSE, ...)
{
a <- names(as.list(match.call()))
if ("tol" %in% a)
@@ -88,7 +88,7 @@ control that algorithm's convergence tolerance. See `?prcomp_irlba` for help.")
ans$scale <- args$scale
if (retx)
{
- ans <- c(ans, list(x = s$d * s$u))
+ ans <- c(ans, list(x = sweep(s$u, 2, s$d, FUN=`*`)))
colnames(ans$x) <- paste("PC", seq(1, ncol(ans$rotation)), sep="")
}
class(ans) <- "prcomp"
diff --git a/R/svdr.R b/R/svdr.R
new file mode 100644
index 0000000..0be2d24
--- /dev/null
+++ b/R/svdr.R
@@ -0,0 +1,118 @@
+#' Find a few approximate largest singular values and corresponding
+#' singular vectors of a matrix.
+#'
+#' The randomized method for truncated SVD by P. G. Martinsson and colleagues
+#' finds a few approximate largest singular values and corresponding
+#' singular vectors of a sparse or dense matrix. It is a fast and
+#' memory-efficient way to compute a partial SVD, similar in performance
+#' for many problems to \code{\link{irlba}}. The \code{svdr} method
+#' is a block method and may produce more accurate estimations with
+#' less work for problems with clustered large singular values (see
+#' the examples). In other problems, \code{irlba} may exhibit faster
+#' convergence.
+#'
+#' @param x numeric real- or complex-valued matrix or real-valued sparse matrix.
+#' @param k dimension of subspace to estimate (number of approximate singular values to compute).
+#' @param it fixed number of algorithm iterations, larger values improve accuracy.
+#' @param extra number of extra vectors of dimension \code{ncol(x)}, larger values generally improve accuracy (with increased
+#' computational cost).
+#' @param center optional column centering vector whose values are implicitly subtracted from each
+#' column of \code{A} without explicitly forming the centered matrix (preserving sparsity).
+#' Optionally specify \code{center=TRUE} as shorthand for \code{center=colMeans(x)}.
+#' Use for efficient principal components computation.
+#' @param Q optional initial random matrix, defaults to a matrix of size \code{ncol(x)} by \code{k + extra} with
+#' entries sampled from a normal random distribution.
+#' @return
+#' Returns a list with entries:
+#' \describe{
+#' \item{d:}{ k approximate singular values}
+#' \item{u:}{ k approximate left singular vectors}
+#' \item{v:}{ k approximate right singular vectors}
+#' \item{mprod:}{ The total number of matrix vector products carried out}
+#' }
+#' @seealso \code{\link{irlba}}, \code{\link{svd}}
+#' @references
+#' Finding structure with randomness: Stochastic algorithms for constructing
+#' approximate matrix decompositions N. Halko, P. G. Martinsson, J. Tropp. Sep. 2009.
+#' @examples
+#' set.seed(1)
+#'
+#' A <- matrix(runif(400), nrow=20)
+#' svdr(A, 3)$d
+#'
+#' # Compare with svd
+#' svd(A)$d[1:3]
+#'
+#' # Compare with irlba
+#' irlba(A, 3)$d
+#'
+#' \dontrun{
+#' # A problem with clustered large singular values where svdr out-performs irlba.
+#' tprolate <- function(n, w=0.25)
+#' {
+#' a <- rep(0, n)
+#' a[1] <- 2 * w
+#' a[2:n] <- sin( 2 * pi * w * (1:(n-1)) ) / ( pi * (1:(n-1)) )
+#' toeplitz(a)
+#' }
+#'
+#' x <- tprolate(512)
+#' set.seed(1)
+#' tL <- system.time(L <- irlba(x, 20))
+#' tR <- system.time(R <- svdr(x, 20))
+#' S <- svd(x)
+#' plot(S$d)
+#' data.frame(time=c(tL[3], tR[3]),
+#' error=sqrt(c(crossprod(L$d - S$d[1:20]), crossprod(R$d - S$d[1:20]))),
+#' row.names=c("IRLBA", "Randomized SVD"))
+#'
+#' # But, here is a similar problem with clustered singular values where svdr
+#' # doesn't out-perform irlba as easily...clusters of singular values are,
+#' # in general, very hard to deal with!
+#' # (This example based on https://github.com/bwlewis/irlba/issues/16.)
+#' set.seed(1)
+#' s <- svd(matrix(rnorm(200 * 200), 200))
+#' x <- s$u %*% (c(exp(-(1:100)^0.3) * 1e-12 + 1, rep(0.5, 100)) * t(s$v))
+#' tL <- system.time(L <- irlba(x, 5))
+#' tR <- system.time(R <- svdr(x, 5))
+#' S <- svd(x)
+#' plot(S$d)
+#' data.frame(time=c(tL[3], tR[3]),
+#' error=sqrt(c(crossprod(L$d - S$d[1:5]), crossprod(R$d - S$d[1:5]))),
+#' row.names=c("IRLBA", "Randomized SVD"))
+#' }
+#' @export
+svdr <- function(x, k, it=3, extra=10, center=NULL, Q=NULL)
+{
+ n <- min(ncol(x), k + extra)
+ if (isTRUE(center)) center <- colMeans(x)
+ if (is.null(Q)) Q <- matrix(rnorm(ncol(x) * n), ncol(x))
+ for (j in 1:it)
+ {
+ if (is.null(center))
+ {
+ Q <- qr.Q(qr(x %*% Q))
+ Q <- qr.Q(qr(t(crossprod(Q, x)))) # a.k.a. Q <- qr.Q(qr(t(x) %*% Q)), but avoiding t(x)
+ } else
+ {
+ Q <- qr.Q(qr(x %*% Q - rep(1, nrow(x)) %*% crossprod(center, Q)))
+ Q <- qr.Q(qr(t(crossprod(Q, x) - tcrossprod(crossprod(Q, rep(1, nrow(x))), center))))
+ }
+ }
+ if (is.null(center))
+ {
+ Q <- qr.Q(qr(x %*% Q))
+ B <- t(Q) %*% x
+ } else
+ {
+ Q <- qr.Q(qr(x %*% Q - rep(1, nrow(x)) %*% crossprod(center, Q)))
+ B <- crossprod(Q, x) - tcrossprod(crossprod(Q, rep(1, nrow(x))), center)
+ }
+ s <- svd(B)
+ s$u <- Q %*% s$u
+ s$u <- s$u[, 1:k]
+ s$d <- s$d[1:k]
+ s$v <- s$v[, 1:k]
+ s$mprod <- 2 * it + 1
+ s
+}
diff --git a/R/utility.R b/R/utility.R
index e87e05d..dda9414 100644
--- a/R/utility.R
+++ b/R/utility.R
@@ -14,13 +14,13 @@ cross <- function(x, y)
}
# Euclidean norm
-norm2 <- function (x)
+norm2 <- function(x)
{
drop(sqrt(cross(x)))
}
# Orthogonalize vectors Y against vectors X.
-orthog <- function (Y, X)
+orthog <- function(Y, X)
{
dx2 <- dim(X)[2]
if (is.null(dx2)) dx2 <- 1
@@ -28,7 +28,7 @@ orthog <- function (Y, X)
if (is.null(dy2)) dy2 <- 1
if (dx2 < dy2) doty <- cross(X, Y)
else doty <- Conj(t(cross(Y, X)))
- return (Y - X %*% doty)
+ Y - X %*% doty
}
# Convergence tests
@@ -36,34 +36,30 @@ orthog <- function (Y, X)
# Bsz Number of rows of the bidiagonal matrix B
# tol
# k_org
-# U_B Left singular vectors of small matrix B
-# S_B Singular values of B
-# V_B Right singular vectors of B
+# Bsvd svd list of small matrix B
# residuals
# k
-# SVTol
# Smax
+# lastsv, svtol, work, S
#
# Output parameter list
# converged TRUE/FALSE
-# U_B Left singular vectors of small matrix B
-# S_B Singular values of B
-# V_B Right singular vectors of B
# k Number of singular vectors returned
-convtests <- function (Bsz, tol, k_org, U_B, S_B, V_B,
- residuals, k, SVTol, Smax)
+convtests <- function(Bsz, tol, k_org, Bsvd, residuals, k, Smax, lastsv, svtol, maxritz, work, S)
{
- len_res <- sum(residuals[1:k_org] < tol * Smax)
+# Converged singular triplets
+ subspace_converged <- residuals[1:k_org] < tol * Smax
+# Converged fixed point triplets
+ if (is.null(lastsv)) lastsv <- 0 * Bsvd$d
+ delta_converged <- (abs(Bsvd$d[1:k_org] - lastsv[1:k_org]) / Bsvd$d[1:k_org]) < svtol
+ len_res <- sum(subspace_converged & delta_converged) # both
if (is.na(len_res)) len_res <- 0
- if (len_res == k_org) {
- return (list(converged=TRUE, U_B=U_B[, 1:k_org, drop=FALSE],
- S_B=S_B[1:k_org, drop=FALSE],
- V_B=V_B[, 1:k_org, drop=FALSE], k=k))
- }
+ if (len_res >= k_org) return(list(converged=TRUE, k=k))
+ if (S == 0) return(list(converged=TRUE, k=k))
# Not converged yet...
-# Adjust k to include more vectors as the number of vectors converge.
- len_res <- sum(residuals[1:k_org] < SVTol * Smax)
- k <- max(k, k_org + len_res)
- if (k > Bsz - 3) k <- max(Bsz - 3, 1)
- return (list(converged=FALSE, U_B=U_B, S_B=S_B, V_B=V_B, k=k) )
+# Adjust k to include more vectors as the number of vectors converge, but not
+# too many (maxritz):
+ augment <- min(sum(subspace_converged), maxritz)
+ k <- min(max(k, k_org + augment), work - 1)
+ list(converged=FALSE, k=k)
}
diff --git a/README.md b/README.md
index 340d3da..836ec4d 100644
--- a/README.md
+++ b/README.md
@@ -6,28 +6,110 @@ for Augmented, <b>I</b>mplicitly <b>R</b>estarted <b>L</b>anczos
<b>B</b>idiagonalization <b>A</b>lgorithm. The package provides the following
functions (see help on each for details and examples).
-* `irlba` main partial SVD function
-* `prcomp_irlba` principal components function similar to the `prcomp` function in stats package for computing the first few principal components of large matrices
-* `partial_eigen` a very limited partial eigenvalue decomposition for symmetric matrices (see the [RSpectra](https://cran.r-project.org/package=RSpectra) package for more comprehensive truncated eigenvalue decomposition)
+* `irlba()` partial SVD function
+* `svdr()` alternate partial SVD function based on randomized SVD
+* `prcomp_irlba()` principal components function similar to the `prcomp` function in stats package for computing the first few principal components of large matrices
+* `partial_eigen()` a very limited partial eigenvalue decomposition for symmetric matrices (see the [RSpectra](https://cran.r-project.org/package=RSpectra) package for more comprehensive truncated eigenvalue decomposition)
-Help documentation for each function includes examples. Also see the package
-vignette, `vignette("irlba", package="irlba")`, and demo,
-`demo("custom_matrix_multiply", package="irlba")`.
+Help documentation for each function includes extensive documentation and
+examples. Also see the package vignette, `vignette("irlba", package="irlba")`.
-Version 2.1.2 of the package includes a convenience `prcomp`-like function for
-computing principal components and numerous bug fixes and numerical stability
-improvements in edge cases.
+## What's new in Version 2.2.1?
+We include stronger convergence criteria and a new argument `svtol`
+for that. The new approach helps guarantee more accurate solutions for some
+difficult problems. The tradeoff is that the default behavior is a little
+slower than before because at least two Lanczos iterations are always run. The
+new convergence behavior can be disabled with `svtol=Inf`.
+
+The new package version includes a new function, svdr()--another state of the
+art truncated SVD method based on the randomized SVD algorithm of Gunnar
+Martinsson and others. Both irlba() and svdr() work well. Svdr uses a block
+method and may exhibit better convergence in problems where the largest
+singular values are clustered. See the documentation and examples in the
+package. (Block versions of irlba exists, but are not yet implemented by this R
+package--something coming in the future.)
+
+We re-introduced a solver for estimating the smallest singular values of a
+matrix and associated singular vector spaces. The solver is based on the
+oringial Harmonic Ritz vector augmentation method of Baglama and Reichel.
+Beware that this method is somewhat experimental and may fail to converge, or
+may converge poorly, to estimated singular values for very ill-conditioned
+matrices. Along with block methods for irlba, this is an active area of
+work--feel free to contribute!
+
+
+## Deprecated features
+
+The `mult` argument is deprecated and will be removed in a future version. We
+now recommend simply defining a custom class with a custom multiplcation
+operator. The example below illustrates the old and new approaches.
+
+```{r}
+library(irlba)
+set.seed(1)
+A <- matrix(rnorm(100), 10)
+
+# ------------------ old way ----------------------------------------------
+# 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))
+mult <- function(x, y)
+ {
+ # check if x is a vector
+ if (is.vector(x))
+ {
+ return((x %*% y) / col_scale)
+ }
+ # else x is the matrix
+ x %*% (y / col_scale)
+ }
+irlba(A, 3, mult=mult)$d
+## [1] 1.820227 1.622988 1.067185
+
+# Compare with:
+irlba(A, 3, scale=col_scale)$d
+## [1] 1.820227 1.622988 1.067185
+
+# Compare with:
+svd(sweep(A, 2, col_scale, FUN=`/`))$d[1:3]
+## [1] 1.820227 1.622988 1.067185
+
+# ------------------ new way ----------------------------------------------
+setClass("scaled_matrix", contains="matrix", slots=c(scale="numeric"))
+setMethod("%*%", signature(x="scaled_matrix", y="numeric"), function(x ,y) x at .Data %*% (y / x at scale))
+setMethod("%*%", signature(x="numeric", y="scaled_matrix"), function(x ,y) (x %*% y at .Data) / y at scale)
+a <- new("scaled_matrix", A, scale=col_scale)
+
+irlba(a, 3)$d
+## [1] 1.820227 1.622988 1.067185
+```
+
+We have learned that using R's existing S4 system is simpler, easier, and more
+flexible than using custom arguments with idiosyncratic syntax and behavior.
+We've even used the new approach to implement distributed parallel matrix
+products for very large problems with amazingly little code.
+
+## Wishlist
+
+- Optional block implementation for some use cases
+- More Matrix classes supported in the fast code path
+- Help improving the solver for smallest singular values in tricky cases (basically, for ill-conditioned problems)
## References
* Augmented Implicitly Restarted Lanczos Bidiagonalization Methods, J. Baglama and L. Reichel, SIAM J. Sci. Comput. 2005. (http://www.math.uri.edu/~jbaglama/papers/paper14.pdf)
+* Finding structure with randomness: Stochastic algorithms for constructing approximate matrix decompositions N. Halko, P. G. Martinsson, J. Tropp. Sep. 2009.
+
## Status
<a href="https://travis-ci.org/bwlewis/irlba">
<img src="https://travis-ci.org/bwlewis/irlba.svg?branch=master" alt="Travis CI status"></img>
</a>
-[](https://codecov.io/github/bwlewis/irlba?branch=master)
-[](https://cran.r-project.org/package=irlba)
-
+<a href="https://codecov.io/gh/bwlewis/irlba">
+ <img src="https://codecov.io/gh/bwlewis/irlba/branch/master/graph/badge.svg" alt="Codecov" />
+</a>
+<a href="https://www.r-pkg.org/pkg/irlba">
+ <img src="http://cranlogs.r-pkg.org/badges/irlba" />
+</a>
diff --git a/build/vignette.rds b/build/vignette.rds
index 6bf9c7b..e919b8c 100644
Binary files a/build/vignette.rds and b/build/vignette.rds differ
diff --git a/demo/00Index b/demo/00Index
deleted file mode 100644
index b968cce..0000000
--- a/demo/00Index
+++ /dev/null
@@ -1 +0,0 @@
-custom_matrix_multiply Examples of custom matrix multiplication with irlba.
diff --git a/demo/custom_matrix_multiply.R b/demo/custom_matrix_multiply.R
deleted file mode 100644
index cee6f97..0000000
--- a/demo/custom_matrix_multiply.R
+++ /dev/null
@@ -1,44 +0,0 @@
-library(irlba)
-set.seed(1)
-A <- matrix(runif(400), nrow=20)
-
-# 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.
-# This approach uses the `mult` argument to irlba and is also found in the
-# examples.
-col_scale <- sqrt(apply(A, 2, crossprod))
-mult <- function(x, y)
- {
- # check if x is a vector
- if (is.vector(x))
- {
- return((x %*% y) / col_scale)
- }
- # else x is the matrix
- x %*% (y / col_scale)
- }
-print(irlba(A, 3, mult=mult)$d)
-
-# Compare with:
-print(svd(sweep(A, 2, col_scale, FUN=`/`))$d[1:3])
-
-
-# An alternative and maybe better approach simply uses R's native operator
-# overloading. You simply need to define a matrix x vector product method and a
-# vector x matrix product method.
-setClass("mymat", contains="matrix", S3methods=TRUE, slots=c(scale="numeric"))
-setMethod("%*%", signature(x="mymat", y="numeric"), function(x ,y)
- {
- x at .Data %*% (y / x at scale)
- })
-setMethod("%*%", signature(x="numeric", y="mymat"), function(x ,y)
- {
- (x %*% y at .Data) / y at scale
- })
-X <- new("mymat", A, scale=col_scale)
-
-# Compare with other approaches
-print(irlba(X, 3, fastpath=FALSE)$d)
-
-# See https://bwlewis.github.io/1000_genomes_examples/PCA_whole_genome.html
-# for a much more complex example using a custom class.
diff --git a/inst/doc/irlba.Rnw b/inst/doc/irlba.Rnw
index 546237b..5f11eb8 100644
--- a/inst/doc/irlba.Rnw
+++ b/inst/doc/irlba.Rnw
@@ -101,7 +101,7 @@ several easy ways to define custom matrix arithmetic that works with other
matrix classes including {\tt big.matrix} from the {\tt bigmemory} package and
others. The {\tt irlba} is both faster and more memory efficient than the
usual R {\tt svd} function for computing a few of the largest singular vectors
-and corresponding singular values of a matrix, and advantage of available
+and corresponding singular values of a matrix. It takes advantage of available
high-performance linear algebra libraries if R is compiled to use them. In
particular, the package uses the same BLAS and LAPACK libraries that R uses
(see
diff --git a/inst/doc/irlba.pdf b/inst/doc/irlba.pdf
index 26781a5..432a249 100644
Binary files a/inst/doc/irlba.pdf and b/inst/doc/irlba.pdf differ
diff --git a/man/irlba.Rd b/man/irlba.Rd
index 0c6096c..d62dff6 100644
--- a/man/irlba.Rd
+++ b/man/irlba.Rd
@@ -2,12 +2,13 @@
% Please edit documentation in R/irlba.R
\name{irlba}
\alias{irlba}
-\title{Find a few approximate largest singular values and corresponding
+\title{Find a few approximate singular values and corresponding
singular vectors of a matrix.}
\usage{
-irlba(A, nv = 5, nu, maxit = 1000, work = nv + 7, reorth = TRUE,
- tol = 1e-05, v = NULL, right_only = FALSE, verbose = FALSE, scale,
- center, du, ds, dv, shift, mult, fastpath = TRUE)
+irlba(A, nv = 5, nu = nv, maxit = 100, work = nv + 7, reorth = TRUE,
+ tol = 1e-05, v = NULL, right_only = FALSE, verbose = FALSE,
+ scale = NULL, center = NULL, shift = NULL, mult = NULL,
+ fastpath = TRUE, svtol = tol, smallest = FALSE, ...)
}
\arguments{
\item{A}{numeric real- or complex-valued matrix or real-valued sparse matrix.}
@@ -22,10 +23,12 @@ irlba(A, nv = 5, nu, maxit = 1000, work = nv + 7, reorth = TRUE,
\item{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. Always set to \code{TRUE}
+iteration but, overall, may require more iterations for convergence. Automatically \code{TRUE}
when \code{fastpath=TRUE} (see below).}
-\item{tol}{convergence is determined when \eqn{\|AV - US\| < tol\|A\|}{||AV - US|| < tol*||A||},
+\item{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
@@ -38,7 +41,7 @@ to restart the algorithm from where it left off (see the notes).}
(\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 \code{right_only = TRUE} sets \code{fastpath = FALSE} (only use this option
-for really large problems that run out of memory).}
+for really large problems that run out of memory and when \code{nrow(A) >> ncol(A)}).}
\item{verbose}{logical value that when \code{TRUE} prints status messages during the computation.}
@@ -46,19 +49,29 @@ for really large problems that run out of memory).}
must be as long as the number of columns of \code{A} (see notes).}
\item{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}
+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).}
-\item{du}{DEPRECATED optional subspace deflation vector (see notes).}
+\item{shift}{optional shift value (square matrices only, see notes).}
-\item{ds}{DEPRECATED optional subspace deflation scalar (see notes).}
+\item{mult}{DEPRECATED optional custom matrix multiplication function (default is \code{\%*\%}, see notes).}
-\item{dv}{DEPRECATED optional subspace deflation vector (see notes).}
+\item{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.}
-\item{shift}{optional shift value (square matrices only, see notes).}
+\item{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.)}
-\item{mult}{optional custom matrix multiplication function (default is \code{\%*\%}, see notes).}
+\item{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.}
-\item{fastpath}{try a fast C algorithm implementation if possible; set \code{fastpath=FALSE} to use the reference R implementation. See notes.}
+\item{...}{optional additional arguments used to support experimental and deprecated features.}
}
\value{
Returns a list with entries:
@@ -72,7 +85,8 @@ Returns a list with entries:
}
\description{
The augmented implicitly restarted Lanczos bidiagonalization algorithm
-(IRLBA) finds a few approximate largest singular values and corresponding
+(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.
}
@@ -94,33 +108,25 @@ 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. This option may not be
-used together with the general rank 1 deflation options. \code{center}
+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 deflation parameters are deprecated and will be removed in
-a future version. They could be used to compute the rank-one deflated
-SVD of \eqn{A - ds \cdot du dv^T}{A - ds*du \%*\% t(dv)}, where
-\eqn{du^T A - ds\cdot dv^T = 0}{t(du) \%*\% A - ds * t(dv) == 0}. For
-example, the triple \code{ds, du, dv} may be a known singular value
-and corresponding singular vectors. Or \code{ds=m} and \code{dv}
-and \code{du} represent a vector of column means of \code{A} and of ones,
-respectively, where \code{m} is the number of rows of \code{A}.
-This functionality can be effectively replaced with custom matrix
-product functions.
-
-Specify an optional alternative matrix multiplication operator in the
+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. See the examples and
-\code{demo("custom_matrix_multiply", package="irlba")} for an
-alternative approach.
+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()} to
+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
@@ -161,6 +167,12 @@ svd(A)$d[1:3]
# (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))
@@ -171,6 +183,7 @@ cbind(P$v,
# 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.
+# This approach is deprecated (see below for a bettwe way to do this).
col_scale <- sqrt(apply(A, 2, crossprod))
mult <- function(x, y)
{
@@ -185,16 +198,22 @@ mult <- function(x, y)
irlba(A, 3, mult=mult)$d
# Compare with:
-irlba(A, 3, scale=col_scale)$d
-
-# Compare with:
svd(sweep(A, 2, col_scale, FUN=`/`))$d[1:3]
+# Compare with the new recommended approach:
+setClass("scaled_matrix", contains="matrix", slots=c(scale="numeric"))
+setMethod("\%*\%", signature(x="scaled_matrix", y="numeric"),
+ function(x ,y) x at .Data \%*\% (y / x at scale))
+setMethod("\%*\%", signature(x="numeric", y="scaled_matrix"),
+ function(x ,y) (x \%*\% y at .Data) / y at scale)
+a <- new("scaled_matrix", A, scale=col_scale)
+irlba(a, 3)$d
+
}
\references{
Augmented Implicitly Restarted Lanczos Bidiagonalization Methods, J. Baglama and L. Reichel, SIAM J. Sci. Comput. 2005.
}
\seealso{
-\code{\link{svd}}, \code{\link{prcomp}}, \code{\link{partial_eigen}}
+\code{\link{svd}}, \code{\link{prcomp}}, \code{\link{partial_eigen}}, \code{\link{svdr}}
}
diff --git a/man/partial_eigen.Rd b/man/partial_eigen.Rd
index c21fba2..efc1685 100644
--- a/man/partial_eigen.Rd
+++ b/man/partial_eigen.Rd
@@ -12,7 +12,8 @@ partial_eigen(x, n = 5, symmetric = TRUE, ...)
\item{n}{number of largest eigenvalues and corresponding eigenvectors to compute.}
\item{symmetric}{\code{TRUE} indicates \code{x} is a symmetric matrix (the default);
-specify \code{symmetric=FALSE} to compute the largest eigenvalues and corresponding eigenvectors of \code{t(x) \%*\% x} instead.}
+specify \code{symmetric=FALSE} to compute the largest eigenvalues and corresponding
+eigenvectors of \code{t(x) \%*\% x} instead.}
\item{...}{optional additional parameters passed to the \code{irlba} function.}
}
diff --git a/man/svdr.Rd b/man/svdr.Rd
new file mode 100644
index 0000000..608a45b
--- /dev/null
+++ b/man/svdr.Rd
@@ -0,0 +1,103 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/svdr.R
+\name{svdr}
+\alias{svdr}
+\title{Find a few approximate largest singular values and corresponding
+singular vectors of a matrix.}
+\usage{
+svdr(x, k, it = 3, extra = 10, center = NULL, Q = NULL)
+}
+\arguments{
+\item{x}{numeric real- or complex-valued matrix or real-valued sparse matrix.}
+
+\item{k}{dimension of subspace to estimate (number of approximate singular values to compute).}
+
+\item{it}{fixed number of algorithm iterations, larger values improve accuracy.}
+
+\item{extra}{number of extra vectors of dimension \code{ncol(x)}, larger values generally improve accuracy (with increased
+computational cost).}
+
+\item{center}{optional column centering vector whose values are implicitly subtracted from each
+column of \code{A} without explicitly forming the centered matrix (preserving sparsity).
+Optionally specify \code{center=TRUE} as shorthand for \code{center=colMeans(x)}.
+Use for efficient principal components computation.}
+
+\item{Q}{optional initial random matrix, defaults to a matrix of size \code{ncol(x)} by \code{k + extra} with
+entries sampled from a normal random distribution.}
+}
+\value{
+Returns a list with entries:
+\describe{
+ \item{d:}{ k approximate singular values}
+ \item{u:}{ k approximate left singular vectors}
+ \item{v:}{ k approximate right singular vectors}
+ \item{mprod:}{ The total number of matrix vector products carried out}
+}
+}
+\description{
+The randomized method for truncated SVD by P. G. Martinsson and colleagues
+finds a few approximate largest singular values and corresponding
+singular vectors of a sparse or dense matrix. It is a fast and
+memory-efficient way to compute a partial SVD, similar in performance
+for many problems to \code{\link{irlba}}. The \code{svdr} method
+is a block method and may produce more accurate estimations with
+less work for problems with clustered large singular values (see
+the examples). In other problems, \code{irlba} may exhibit faster
+convergence.
+}
+\examples{
+set.seed(1)
+
+A <- matrix(runif(400), nrow=20)
+svdr(A, 3)$d
+
+# Compare with svd
+svd(A)$d[1:3]
+
+# Compare with irlba
+irlba(A, 3)$d
+
+\dontrun{
+# A problem with clustered large singular values where svdr out-performs irlba.
+tprolate <- function(n, w=0.25)
+{
+ a <- rep(0, n)
+ a[1] <- 2 * w
+ a[2:n] <- sin( 2 * pi * w * (1:(n-1)) ) / ( pi * (1:(n-1)) )
+ toeplitz(a)
+}
+
+x <- tprolate(512)
+set.seed(1)
+tL <- system.time(L <- irlba(x, 20))
+tR <- system.time(R <- svdr(x, 20))
+S <- svd(x)
+plot(S$d)
+data.frame(time=c(tL[3], tR[3]),
+ error=sqrt(c(crossprod(L$d - S$d[1:20]), crossprod(R$d - S$d[1:20]))),
+ row.names=c("IRLBA", "Randomized SVD"))
+
+# But, here is a similar problem with clustered singular values where svdr
+# doesn't out-perform irlba as easily...clusters of singular values are,
+# in general, very hard to deal with!
+# (This example based on https://github.com/bwlewis/irlba/issues/16.)
+set.seed(1)
+s <- svd(matrix(rnorm(200 * 200), 200))
+x <- s$u \%*\% (c(exp(-(1:100)^0.3) * 1e-12 + 1, rep(0.5, 100)) * t(s$v))
+tL <- system.time(L <- irlba(x, 5))
+tR <- system.time(R <- svdr(x, 5))
+S <- svd(x)
+plot(S$d)
+data.frame(time=c(tL[3], tR[3]),
+ error=sqrt(c(crossprod(L$d - S$d[1:5]), crossprod(R$d - S$d[1:5]))),
+ row.names=c("IRLBA", "Randomized SVD"))
+}
+}
+\references{
+Finding structure with randomness: Stochastic algorithms for constructing
+approximate matrix decompositions N. Halko, P. G. Martinsson, J. Tropp. Sep. 2009.
+}
+\seealso{
+\code{\link{irlba}}, \code{\link{svd}}
+}
+
diff --git a/src/irlb.c b/src/irlb.c
index 0bb265c..e269341 100644
--- a/src/irlb.c
+++ b/src/irlb.c
@@ -61,7 +61,8 @@ RNORM (int n)
return ans;
}
-/* irlb C implementation wrapper
+/* irlb C implementation wrapper for R
+ *
* X double precision input matrix
* NU integer number of singular values/vectors to compute must be > 3
* INIT double precision starting vector length(INIT) must equal ncol(X)
@@ -76,6 +77,7 @@ RNORM (int n)
* SCALE either NULL (no scaling) or a vector of length ncol(X)
* SHIFT either NULL (no shift) or a single double-precision number
* CENTER either NULL (no centering) or a vector of length ncol(X)
+ * SVTOL double tolerance max allowed per cent change in each estimated singular value
*
* Returns a list with 6 elements:
* 1. vector of estimated singular values
@@ -88,26 +90,27 @@ RNORM (int n)
SEXP
IRLB (SEXP X, SEXP NU, SEXP INIT, SEXP WORK, SEXP MAXIT, SEXP TOL, SEXP EPS,
SEXP MULT, SEXP RESTART, SEXP RV, SEXP RW, SEXP RS, SEXP SCALE,
- SEXP SHIFT, SEXP CENTER)
+ SEXP SHIFT, SEXP CENTER, SEXP SVTOL)
{
SEXP ANS, S, U, V;
double *V1, *U1, *W, *F, *B, *BU, *BV, *BS, *BW, *res, *T, *scale, *shift,
- *center;
+ *center, *SVRATIO;
int i, iter, mprod, ret;
R_xlen_t m, n;
int mult = INTEGER (MULT)[0];
- void *A;
+ void *AS = NULL;
+ double *A = NULL;
switch (mult)
{
case 1:
- A = (void *) AS_CHM_SP (X);
+ AS = (void *) AS_CHM_SP (X);
int *dims = INTEGER (GET_SLOT (X, install ("Dim")));
m = dims[0];
n = dims[1];
break;
default:
- A = (void *) REAL (X);
+ A = REAL (X);
m = nrows (X);
n = ncols (X);
}
@@ -115,6 +118,7 @@ IRLB (SEXP X, SEXP NU, SEXP INIT, SEXP WORK, SEXP MAXIT, SEXP TOL, SEXP EPS,
R_xlen_t work = INTEGER (WORK)[0];
int maxit = INTEGER (MAXIT)[0];
double tol = REAL (TOL)[0];
+ double svtol = REAL (SVTOL)[0];
int lwork = 7 * work * (1 + work);
int restart = INTEGER (RESTART)[0];
double eps = REAL (EPS)[0];
@@ -144,6 +148,7 @@ IRLB (SEXP X, SEXP NU, SEXP INIT, SEXP WORK, SEXP MAXIT, SEXP TOL, SEXP EPS,
{
center = REAL (CENTER);
}
+ SVRATIO = (double *) R_alloc (work, sizeof (double));
V1 = (double *) R_alloc (n * work, sizeof (double));
U1 = (double *) R_alloc (m * work, sizeof (double));
W = (double *) R_alloc (m * work, sizeof (double));
@@ -164,9 +169,9 @@ IRLB (SEXP X, SEXP NU, SEXP INIT, SEXP WORK, SEXP MAXIT, SEXP TOL, SEXP EPS,
B[i + work * i] = REAL (RS)[i];
}
ret =
- irlb (A, mult, m, n, nu, work, maxit, restart, tol, scale, shift, center,
+ irlb (A, AS, mult, m, n, nu, work, maxit, restart, tol, scale, shift, center,
REAL (S), REAL (U), REAL (V), &iter, &mprod, eps, lwork, V1, U1, W,
- F, B, BU, BV, BS, BW, res, T);
+ F, B, BU, BV, BS, BW, res, T, svtol, SVRATIO);
SET_VECTOR_ELT (ANS, 0, S);
SET_VECTOR_ELT (ANS, 1, U);
SET_VECTOR_ELT (ANS, 2, V);
@@ -188,8 +193,9 @@ IRLB (SEXP X, SEXP NU, SEXP INIT, SEXP WORK, SEXP MAXIT, SEXP TOL, SEXP EPS,
* all data must be allocated by caller, required sizes listed below
*/
int
-irlb (void *A, // Input data matrix
- int mult, // 0 -> A is double *, 1 -> A is cholmod
+irlb (double *A, // Input data matrix (double case)
+ void * AS, // input data matrix (sparse case)
+ int mult, // 0 -> use double *A, 1 -> use AS
int m, // data matrix number of rows, must be > 3.
int n, // data matrix number of columns, must be > 3.
int nu, // dimension of solution
@@ -201,7 +207,7 @@ irlb (void *A, // Input data matrix
double *shift, // optional shift (NULL for no shift)
double *center, // optional center (NULL for no center)
// output values
- double *s, // output singular vectors at least length nu
+ double *s, // output singular values at least length nu
double *U, // output left singular vectors m x work
double *V, // output right singular vectors n x work
int *ITER, // ouput number of Lanczos iterations
@@ -218,7 +224,9 @@ irlb (void *A, // Input data matrix
double *BS, // work
double *BW, // lwork
double *res, // work
- double *T) // lwork
+ double *T, // lwork
+ double svtol, // svtol limit
+ double *svratio) // convtest extra storage vector of length work
{
double d, S, R, alpha, beta, R_F, SS;
double *x;
@@ -266,7 +274,7 @@ irlb (void *A, // Input data matrix
switch (mult)
{
case 1:
- dsdmult ('n', m, n, (CHM_SP) A, x, W + j * m);
+ dsdmult ('n', m, n, AS, x, W + j * m);
break;
default:
alpha = 1;
@@ -307,7 +315,7 @@ irlb (void *A, // Input data matrix
switch (mult)
{
case 1:
- dsdmult ('t', m, n, (CHM_SP) A, W + j * m, F);
+ dsdmult ('t', m, n, AS, W + j * m, F);
break;
default:
alpha = 1.0;
@@ -363,7 +371,7 @@ irlb (void *A, // Input data matrix
switch (mult)
{
case 1:
- dsdmult ('n', m, n, (CHM_SP) A, x, W + (j + 1) * m);
+ dsdmult ('n', m, n, AS, x, W + (j + 1) * m);
break;
default:
alpha = 1.0;
@@ -428,19 +436,25 @@ irlb (void *A, // Input data matrix
/* Force termination after encountering linear dependence */
if (R_F < 2 * eps)
R_F = 0;
+
Smax = 0;
for (jj = 0; jj < j; ++jj)
- if (BS[jj] > Smax)
- Smax = BS[jj];
+ {
+ if (BS[jj] > Smax)
+ Smax = BS[jj];
+ svratio[jj] = fabs (svratio[jj] - BS[jj]) / BS[jj];
+ }
for (kk = 0; kk < j; ++kk)
res[kk] = R_F * BU[kk * work + (j - 1)];
/* Update k to be the number of converged singular values. */
- convtests (j, nu, tol, Smax, res, &k, &converged);
+ convtests (j, nu, tol, svtol, Smax, svratio, res, &k, &converged, S);
if (converged == 1)
{
iter++;
break;
}
+ for (jj = 0; jj < j; ++jj)
+ svratio[jj] = BS[jj];
alpha = 1;
beta = 0;
@@ -493,14 +507,21 @@ irlba_R_cholmod_error (int status, const char *file, int line,
warning ("Cholmod warning '%s' at file:%s, line %d", message, file, line);
}
+static const R_CallMethodDef CallEntries[] = {
+ {"IRLB", (DL_FUNC) & IRLB, 16},
+ {NULL, NULL, 0}
+};
+
#ifdef HAVE_VISIBILITY_ATTRIBUTE
__attribute__ ((visibility ("default")))
#endif
- void R_init_irlba (DllInfo * dll)
+void
+R_init_irlba (DllInfo * dll)
{
+ R_registerRoutines (dll, NULL, CallEntries, NULL, NULL);
+ R_useDynamicSymbols (dll, 0);
M_R_cholmod_start (&chol_c);
chol_c.final_ll = 1; /* LL' form of simplicial factorization */
-
/* need own error handler, that resets final_ll (after *_defaults()) : */
chol_c.error_handler = irlba_R_cholmod_error;
}
@@ -513,7 +534,7 @@ R_unload_irlba (DllInfo * dll)
void
-dsdmult (char transpose, int m, int n, void *a, double *b, double *c)
+dsdmult (char transpose, int m, int n, void * a, double *b, double *c)
{
DL_FUNC sdmult = R_GetCCallable ("Matrix", "cholmod_sdmult");
int t = transpose == 't' ? 1 : 0;
@@ -539,8 +560,6 @@ dsdmult (char transpose, int m, int n, void *a, double *b, double *c)
chc.x = (void *) c;
chc.z = (void *) NULL;
- double one[] = { 1, 0 }, zero[] =
- {
- 0, 0};
+ double one[] = { 1, 0 }, zero[] = { 0, 0};
sdmult (cha, t, one, zero, &chb, &chc, &chol_c);
}
diff --git a/src/irlb.h b/src/irlb.h
index 6a2a630..7f963d7 100644
--- a/src/irlb.h
+++ b/src/irlb.h
@@ -11,11 +11,14 @@ void
convtests (int Bsz, // Number of rows of bidiagonal matrix B
int n, // requested number of singular values
double tol, // convergence tolerance
+ double svtol, // max change each singular value tolerane
double Smax, // largest singular value of B
+ double *svratio, // vector of relative singular value ratio compared to last iteration
double *residuals, // vector of residual values
int *k, // number of estimated singular values (INPUT)
// adjusted subspace size (OUTPUT)
- int *converged); // 0 = FALSE, 1 = TRUE
+ int *converged, // 0 = FALSE, 1 = TRUE
+ double S); // If S == 0 then invariant subspace found.
/*
* Simple cholmod double-precision sparse matrix times dense vector multiplication interface
@@ -31,10 +34,11 @@ dsdmult(char transpose, // 't' -> op(a) = t(a), non-transposed a otherwise
double *b, // double precision dense vector
double *c); // output
-/* IRLB method for sparse or dense double-precision valued matrices */
+/* IRLB function for sparse or dense double-precision valued matrices */
int
-irlb(void *A, // Input data matrix
- int mult, // 0->A is double *, 1-> A is sparse double *
+irlb(double *A, // input data matrix (dense case)
+ void *AS, // input data matrix (sparse case)
+ int mult, // 0 -> A is double *, 1-> A is sparse double *
int m, // data matrix number of rows
int n, // data matrix number of columns
int nu, // dimension of solution
@@ -62,4 +66,6 @@ irlb(void *A, // Input data matrix
double *BS, // working storage work
double *BW, // working storage lwork * lwork
double *res, // working storage work
- double *T); // working storage lwork
+ double *T, // working storage lwork
+ double svtol, // svtol tolerance on maximum ratio change per singular value per iteration
+ double *SVRATIO); // working storage nu
diff --git a/src/utility.c b/src/utility.c
index a359c8f..a2d6385 100644
--- a/src/utility.c
+++ b/src/utility.c
@@ -43,7 +43,7 @@ orthog (double *X, double *Y, double *T, int xm, int xn, int yn)
{
double a = 1, b = 1;
int inc = 1;
- memset(T, 0, xn * yn * sizeof(double));
+ memset (T, 0, xn * yn * sizeof (double));
// T = t(X) * Y
F77_NAME (dgemv) ("t", &xm, &xn, &a, X, &xm, Y, &inc, &b, T, &inc);
// Y = Y - X * T
@@ -56,28 +56,31 @@ orthog (double *X, double *Y, double *T, int xm, int xn, int yn)
/*
* Convergence tests
* Input parameters
- * Bsz Number of rows of the bidiagonal matrix B (scalar)
+ * Bsz number of rows of the bidiagonal matrix B (scalar)
* tol convergence tolerance (scalar)
- * n Requested number of singular values
- * Smax Largest singular value of B
+ * svtol max change in each singular value tolerance (scalar)
+ * n requested number of singular values
+ * Smax largest singular value of B
+ * svratio vector of abs(current - previous) / current singular value ratios
* residuals vector of residual values
* k number of estimated signular values (scalar)
+ * S check for invariant subspace when S == 0
*
* Output
* converged 0 = FALSE, 1 = TRUE (all converged)
- * k Adjusted subspace size.
+ * k adjusted subspace size.
*/
void
-convtests (int Bsz, int n, double tol, double Smax,
- double *residuals, int *k, int *converged)
+convtests (int Bsz, int n, double tol, double svtol, double Smax,
+ double *svratio, double *residuals, int *k, int *converged, double S)
{
int j, Len_res = 0;
for (j = 0; j < Bsz; ++j)
{
- if (fabs(residuals[j]) < tol * Smax)
+ if ((fabs (residuals[j]) < tol * Smax) && (svratio[j] < svtol))
Len_res++;
}
- if (Len_res >= n)
+ if (Len_res >= n || S == 0)
{
*converged = 1;
return;
diff --git a/tests/edge.R b/tests/edge.R
index fc1bcba..19da719 100644
--- a/tests/edge.R
+++ b/tests/edge.R
@@ -1,66 +1,81 @@
# Tests for a few edge cases
-
require("irlba")
-# Dense matrix
-set.seed(1)
-A <- matrix(rnorm(16), 4)
-L <- irlba(A, nu=1, nv=1, tol=1e-9, fastpath=FALSE)
-L1 <- irlba(A, nu=1, nv=1, tol=1e-9, fastpath=TRUE)
-S <- svd(A, nu=1, nv=1)
-if (!isTRUE(all.equal(L$d, S$d[1])))
-{
- stop("Failed tiny reference example ")
-}
-if (!isTRUE(all.equal(L1$d, S$d[1])))
+loc <- ""
+
+test <- function()
{
- stop("Failed tiny fastpath example")
-}
+ on.exit(message("Error occured in: ", loc))
+ # Dense matrix
+ loc <<- "dense"
+ set.seed(1)
+ A <- matrix(rnorm(16), 4)
+ L <- irlba(A, nu=1, nv=1, tol=1e-9, fastpath=FALSE)
+ L1 <- irlba(A, nu=1, nv=1, tol=1e-9, fastpath=TRUE)
+ S <- svd(A, nu=1, nv=1)
+ if (!isTRUE(all.equal(L$d, S$d[1])))
+ {
+ stop("Failed tiny reference example ")
+ }
+ if (!isTRUE(all.equal(L1$d, S$d[1])))
+ {
+ stop("Failed tiny fastpath example")
+ }
-# Tickle misc. checks
-set.seed(1)
-A <- matrix(rnorm(100), 10)
-L <- tryCatch(irlba(A, nv=3, tol=1e-9, fastpath=FALSE, work=2, v=rep(0, nrow(A))), error=function(e) "NULLSPACE")
-S <- svd(A)
-L <- irlba(A, nv=3, tol=1e-9, fastpath=FALSE, work=2, v=S$v[, 1])
-A <- S$u %*% diag(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1e-12)) %*% t(S$v)
-L <- irlba(A, nv=3, tol=1e-9, fastpath=FALSE, work=2, reorth=FALSE)
+ # Tickle misc. checks
+ loc <<- "misc"
+ set.seed(1)
+ A <- matrix(rnorm(100), 10)
+ L <- tryCatch(irlba(A, nv=3, tol=1e-9, fastpath=FALSE, work=2, v=rep(0, nrow(A))), error=function(e) "NULLSPACE")
+ S <- svd(A)
+ L <- irlba(A, nv=3, tol=1e-9, fastpath=FALSE, work=2, v=S$v[, 1])
+ A <- S$u %*% diag(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1e-12)) %*% t(S$v)
+ L <- irlba(A, nv=3, tol=1e-9, fastpath=FALSE, work=2, reorth=FALSE)
-# Convergence
-A <- S$u %*% (c(1e-5, rep(1e-9, 9)) * t(S$v))
-for (tol in 10 ^ -(7:12))
-{
- L <- irlba(A, 3, tol=tol)
- converged <- svd(A %*% L$v - L$u %*% diag(L$d))$d[1] < tol * L$d[1]
- stopifnot(converged)
-}
+ # Convergence
+ loc <<- "convergence"
+ A <- S$u %*% (c(1e-5, rep(1e-9, 9)) * t(S$v))
+ for (tol in 10 ^ - (7:12))
+ {
+ L <- irlba(A, 3, tol=tol, svtol=Inf)
+ converged <- svd(A %*% L$v - L$u %*% diag(L$d))$d[1] < tol * L$d[1]
+ stopifnot(converged)
+ }
-# Sparse but not dgCMatrix (issue #6)
-A <- Matrix(matrix(rnorm(100), 10))
-L <- irlba(A, nv=1)
-S <- svd(A, nu=1, nv=1)
-if (!isTRUE(all.equal(L$d, S$d[1])))
-{
- stop("Failed general sparse matrix example ")
-}
+ # Sparse but not dgCMatrix (issue #6)
+ loc <<- "misc sparse"
+ A <- Matrix(matrix(rnorm(100), 10))
+ L <- irlba(A, nv=1)
+ S <- svd(A, nu=1, nv=1)
+ if (!isTRUE(all.equal(L$d, S$d[1])))
+ {
+ stop("Failed general sparse matrix example ")
+ }
-A <- Matrix(sample(c(FALSE, TRUE), 100, replace=TRUE), 10, 10)
-L <- irlba(A, nv=1)
-S <- svd(A, nu=1, nv=1)
-if (!isTRUE(all.equal(L$d, S$d[1])))
-{
- stop("Failed logical sparse matrix example ")
-}
+ A <- Matrix(sample(c(FALSE, TRUE), 100, replace=TRUE), 10, 10)
+ L <- irlba(A, nv=1)
+ S <- svd(A, nu=1, nv=1)
+ if (!isTRUE(all.equal(L$d, S$d[1])))
+ {
+ stop("Failed logical sparse matrix example ")
+ }
-# Test for issue #7, a really dumb bug.
-mx <- matrix(sample(1:10, 10 * 100, replace=TRUE), nrow=10)
-S <- irlba(mx, nv=2, verbose=TRUE, center=colMeans(mx), right_only=TRUE)
+ # Test for issue #7, a really dumb bug.
+ loc <<- "issue 7"
+ mx <- matrix(sample(1:10, 10 * 100, replace=TRUE), nrow=10)
+ S <- irlba(mx, nv=2, verbose=TRUE, center=colMeans(mx), right_only=TRUE)
-# test for issue #9
-set.seed(2)
-s1 <- irlba(diag(c(1,2,3,4,5,0,0,0,0)), 4)
-set.seed(2)
-s2 <- irlba(diag(c(1,2,3,4,5,0,0,0,0)), 4, fastpath=FALSE)
-if (!isTRUE(all.equal(s1$d, s2$d)))
-{
- stop("Failed fastpath invariant subspace detection")
+ # test for issue #9
+ loc <<- "issue 9"
+ set.seed(2)
+ s1 <- irlba(diag(c(1, 2, 3, 4, 5, 0, 0, 0, 0)), 4, fastpath=FALSE)
+ set.seed(2)
+ s2 <- irlba(diag(c(1, 2, 3, 4, 5, 0, 0, 0, 0)), 4, fastpath=TRUE)
+ stopifnot(all.equal(s1$d, s2$d))
+ # Repeat this test with different seed
+ set.seed(3)
+ s2 <- irlba(diag(c(1, 2, 3, 4, 5, 0, 0, 0, 0)), 4, fastpath=TRUE)
+ stopifnot(all.equal(s1$d, s2$d))
+ on.exit()
}
+
+test()
diff --git a/tests/svdr.R b/tests/svdr.R
new file mode 100644
index 0000000..7feecdb
--- /dev/null
+++ b/tests/svdr.R
@@ -0,0 +1,50 @@
+# Tests for svdr
+require("irlba")
+loc <- ""
+
+test <- function()
+{
+ on.exit(message("Error occured in: ", loc))
+ # Dense matrix
+ loc <<- "svdr dense"
+ set.seed(1)
+ A <- matrix(rnorm(16), 4)
+ L <- svdr(A, 1)
+ S <- svd(A, nu=1, nv=1)
+ stopifnot(isTRUE(all.equal(L$d, S$d[1])))
+
+ loc <<- "svdr dense m > n"
+ A <- matrix(rnorm(50 * 40), 50)
+ L <- svdr(A, 5, 10, extra=15)
+ S <- svd(A, nu=5, nv=5)
+ stopifnot(isTRUE(all.equal(L$d, S$d[1:5])))
+
+ loc <<- "svdr dense m < n"
+ A <- matrix(rnorm(50 * 40), 40)
+ L <- svdr(A, 5, 10, extra=15)
+ S <- svd(A, nu=5, nv=5)
+ stopifnot(isTRUE(all.equal(L$d, S$d[1:5])))
+
+ # Sparse but not dgCMatrix (issue #6)
+ loc <<- "svdr misc sparse"
+ A <- Matrix(matrix(rnorm(100), 10))
+ L <- svdr(A, 1)
+ S <- svd(A, nu=1, nv=1)
+ stopifnot(isTRUE(all.equal(L$d, S$d[1])))
+
+ loc <<- "svdr logical sparse"
+ A <- Matrix(sample(c(FALSE, TRUE), 100, replace=TRUE), 10, 10)
+ L <- svdr(A, 1)
+ S <- svd(A, nu=1, nv=1)
+ stopifnot(isTRUE(all.equal(L$d, S$d[1])))
+
+ loc <<- "svdr center only, sparse"
+ A <- Matrix(matrix(rnorm(100), 10))
+ m <- colMeans(A)
+ L <- svdr(A, 3, center=m)
+ S <- svd(scale(A, center=TRUE, scale=FALSE))
+ stopifnot(isTRUE(all.equal(L$d, S$d[1:3])))
+
+ on.exit()
+}
+test()
diff --git a/tests/test.R b/tests/test.R
index fa932e6..ecaabdf 100644
--- a/tests/test.R
+++ b/tests/test.R
@@ -19,6 +19,13 @@ for (FAST in c(FALSE, TRUE))
stop("Failed restart", " fastpath=", FAST)
}
+ # unequal nu, nv
+ L <- irlba(A, nv=2, nu=3, fastpath=FAST)
+ if (!isTRUE(ncol(L$v) == 2 && ncol(L$u) == 3))
+ {
+ stop("Failed unequal nu,nv", " fastpath=", FAST)
+ }
+
# Scaling and centering, dense
s <- sqrt(apply(A, 2, crossprod))
m <- colMeans(A)
@@ -126,4 +133,50 @@ for (FAST in c(FALSE, TRUE))
{
stop("Failed nonsquare test", " fastpath=", FAST)
}
+
+# This pathological example was provided by Giuseppe Rodriguez, http://bugs.unica.it/~gppe/
+# The singular values cluster at 1 and 0, making it hard to converge to a truncated
+# subspace containing the largest few singular values (they are all very close).
+# Or, for that matter, the smallest.
+#
+# Reference:
+# J. M. Varah. The Prolate matrix. Linear Algebra and Appl.,
+# 187:269-278, 1993.
+# Michela Redivo-Zaglia, University of Padova, Italy
+# Email: Michela.RedivoZaglia at unipd.it
+# Giuseppe Rodriguez, University of Cagliari, Italy
+# Email: rodriguez at unica.it
+ tprolate <- function(n, w=0.25)
+ {
+ a <- rep(0, n)
+ a[1] <- 2 * w
+ a[2:n] <- sin(2 * pi * w * (1:(n-1))) / (pi * (1:(n-1)))
+ toeplitz(a)
+ }
+
+ x <- tprolate(512)
+ set.seed(1)
+ l <- irlba(x, nv=20, fastpath=FAST)
+ if (isTRUE(max(abs(l$d - 1)) > 1e-3))
+ {
+ stop("Failed tprolate test fastpath=", FAST)
+ }
+}
+
+# smallest=TRUE, m > n (fastpath always FALSE in this case)
+x <- matrix(rnorm(5000), 100)
+set.seed(1)
+L <- irlba(x, nv=5, smallest=TRUE)
+if (!isTRUE(all.equal(L$d, tail(svd(x)$d, 5))))
+{
+ stop("Failed smallest svd test")
+}
+
+# smallest=TRUE, n > m
+x <- matrix(rnorm(5000), 50)
+set.seed(1)
+L <- irlba(x, nv=5, smallest=TRUE)
+if (!isTRUE(all.equal(L$d, tail(svd(x)$d, 5))))
+{
+ stop("Failed smallest svd test")
}
diff --git a/vignettes/irlba.Rnw b/vignettes/irlba.Rnw
index 546237b..5f11eb8 100644
--- a/vignettes/irlba.Rnw
+++ b/vignettes/irlba.Rnw
@@ -101,7 +101,7 @@ several easy ways to define custom matrix arithmetic that works with other
matrix classes including {\tt big.matrix} from the {\tt bigmemory} package and
others. The {\tt irlba} is both faster and more memory efficient than the
usual R {\tt svd} function for computing a few of the largest singular vectors
-and corresponding singular values of a matrix, and advantage of available
+and corresponding singular values of a matrix. It takes advantage of available
high-performance linear algebra libraries if R is compiled to use them. In
particular, the package uses the same BLAS and LAPACK libraries that R uses
(see
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/r-cran-irlba.git
More information about the debian-med-commit
mailing list