svd.HDF5Matrix: Singular Value Decomposition of an HDF5Matrix

In BigDataStatMeth: Scalable Statistical Computing with HDF5-Backed Matrices

Singular Value Decomposition of an HDF5Matrix

Description

Block-wise SVD entirely on disk. The matrix x is decomposed into x = u %*% diag(d) %*% t(v).

Usage

## S3 method for class 'HDF5Matrix'
svd(
  x,
  nu = min(dim(x)),
  nv = min(dim(x)),
  center = TRUE,
  scale = TRUE,
  k = 2L,
  q = 1L,
  method = "auto",
  rankthreshold = 0,
  overwrite = FALSE,
  threads = -1L,
  ...
)

Arguments

x	An HDF5Matrix object.
nu	Number of left singular vectors to compute (default = min(dim(x))).
nv	Number of right singular vectors to compute (default = min(dim(x))).
center	Logical. Center columns before decomposition (default TRUE).
scale	Logical. Scale columns before decomposition (default TRUE).
k	Number of local SVDs per incremental level (default 2).
q	Number of incremental levels (default 1).
method	Computation method: "auto" (default), "blocks", or "full".
rankthreshold	Numeric in [0, 0.1]. Rank approximation threshold (default 0).
overwrite	Logical. Overwrite existing SVD results (default FALSE).
threads	Integer. OpenMP threads (-1 = auto-detect). The effective number may be lower than requested: the system ceiling (OMP_NUM_THREADS, OMP_THREAD_LIMIT and, by default, half of the detected cores) is applied silently by the OpenMP runtime. Requesting more than that is not an error, but it now emits a warning stating how many threads can actually be used.
...	Ignored (S3 compatibility).

Details

Singular values d are loaded into a plain numeric vector (they are always small: at most min(nrow(x), ncol(x)) values). u and v are returned as HDF5Matrix objects.

Constant columns and scale = TRUE. A column with zero variance cannot be rescaled to unit variance, so the call stops with an error reporting how many columns are affected and a few of their positions – the same contract as stats::prcomp(x, scale. = TRUE). Previously the division by zero propagated silently and every singular value came back as 0. Either pass scale = FALSE or remove the constant columns first. scale.HDF5Matrix is deliberately unaffected: like base::scale() it still returns NaN for such a column.

Value

Named list with:

d

Numeric vector of non-negative singular values, decreasing.

u

HDF5Matrix of left singular vectors, nrow(x) x nu.

v

HDF5Matrix of right singular vectors, ncol(x) x nv.

The list carries attributes describing what was actually computed:

method

"full" (exact LAPACK) or "blocks" (hierarchical block algorithm).

exact

TRUE when the exact LAPACK algorithm was used, i.e. method == "full" and no per-block truncation. It reports the algorithm, not a measured accuracy: the block path at full rank is not exact by construction, though in our tests it agreed with LAPACK to ~1e-15 (see below).

nev

Per-block truncation rank actually applied (0 = none). Requesting fewer vectors below the boundary does not truncate anything: the exact path runs and the reduction happens afterwards in R, losslessly.

truncated

TRUE when per-block truncation was applied, which is the dominant source of error (see below).

elements

nrow(x) * ncol(x), the quantity compared against the boundary.

auto_threshold

The boundary itself, see svd_auto_threshold.

blocking

The k and q actually used.

rank

Number of singular values returned.

Exact versus approximate decomposition

This is the most important behavioural difference from base::svd(). Two independent things can make the result approximate, and they are not equally dangerous.

1. The algorithm (minor). method = "full" reads the whole matrix into RAM and computes an exact LAPACK SVD, needing nrow * ncol * 8 bytes. method = "blocks" computes a hierarchical block SVD entirely on disk. method = "auto" (the default) picks "full" while nrow * ncol < svd_auto_threshold() (currently 53,687,091 elements, about 410 MB of doubles) and "blocks" at or above it. When the full rank is requested, the block path is essentially exact: relative errors of order 1e-15 across the whole spectrum, independent of k and q, in our measurements. Crossing the boundary is, on its own, not a numerical event.

2. Per-block rank truncation (major). Asking for fewer singular triplets than min(dim(x)) – nu or nv below full rank, or ncomponents/rank. in prcomp.HDF5Matrix – switches on truncation of every local block SVD to that rank before the blocks are merged. That is what actually costs accuracy, and it compounds with blocking depth. On a 1600 x 200 matrix with a decaying spectrum, relative error in the singular values:

requested rank	k, q	leading	trailing
200 (full)	2, 1	3e-16	7e-15
200 (full)	4, 3	1e-15	4e-15
50	2, 1	7e-06	8e-02
50	2, 3	5e-05	1e-01
10	2, 1	2e-03	1e-01
10	2, 3	1e-02	2e-01

Note that the error is not confined to the tail: at rank 10 even the leading singular value was wrong by 1-3\ of accurate components, ask for a generously oversampled rank and discard the extra ones – in the table above, computing 50 and keeping 10 is about 250 times more accurate in the leading value than computing 10 directly – or request the full rank when you can afford it.

Detecting it. Both switches are silent by default, so the result records them: check attr(res, "exact"), attr(res, "truncated") and attr(res, "method"). svd() also emits a message() when "auto" falls through to the block path and when per-block truncation is applied (silence with suppressMessages()).

There is no cheap residual diagnostic: verifying a decomposition means reconstructing u %*% diag(d) %*% t(v) and comparing it with x, another full pass over the data. The practical checks are to recompute with a larger requested rank, or with method = "full" on a subset, and see whether the leading components move.

There is no cheap residual diagnostic: verifying the decomposition means reconstructing u %*% diag(d) %*% t(v) and comparing it with x, which costs another full pass over the data. When accuracy matters, the practical check is to recompute a small trailing portion with method = "full" on a subset, or to compare across two values of q.

Examples

tmp <- tempfile(fileext = ".h5")
X   <- hdf5_create_matrix(tmp, "data/M", data = matrix(rnorm(500), 50, 10))
res <- svd(X)
length(res$d)   # 10  (min(50,10))
dim(res$u)      # 50 x 10
dim(res$v)      # 10 x 10
X$close()
res$u$close(); res$v$close()
unlink(tmp)

BigDataStatMeth documentation built on Sept. 15, 2026, 1:08 a.m.