Core Internals

These modules implement the two algorithmic building blocks that power tensorrsvd.ho_rsvd(). They are not part of the public API but are documented here for contributors and anyone who wants to use the primitives directly.

For the mathematical description of both components see Tensor-Free Approach and Randomized SVD in the Theory page.

Matricization operator

The key idea behind TensorRSVD is that a mode-m unfolding of the tensor can be exposed as a pylops.LinearOperator without ever forming the dense matrix. The operator evaluates the tensor callable on-the-fly whenever a matrix–vector (or matrix–matrix) product is requested.

class tensorrsvd.core.matricization.MatricizedTensorOperator(tensor_callable, tensor_shape, mode, dtype, backend)[source]

Bases: LinearOperator

Linear operator for the mode-m unfolding of a tensor defined by a callable.

Parameters:

tensor_callable (Callable) – Vectorized callable accepting k array arguments (coordinates in [0, 1]) and returning the tensor values at those coordinates. Must be JAX-traceable if backend == ‘jax’.
tensor_shape (tuple[int, ]) – Shape of the tensor (n_0, n_1, …, n_{k-1}).
mode (int) – Mode along which to unfold.
dtype (npt.DTypeLike) – Data type of the operator.
backend (str) – Array backend: ‘numpy’, ‘cupy’, or ‘jax’.

__init__(tensor_callable, tensor_shape, mode, dtype, backend)[source]

Parameters:

tensor_callable (Callable)
tensor_shape (tuple[int, ...])
mode (int)
dtype (DTypeLike)
backend (str)

Randomized SVD

The randomized SVD is implemented in two stages:

randomized_range_finder() builds an orthonormal basis \(Q\) for the approximate column space of the operator using a Gaussian sketch and optional power iterations.
rsvd_left() calls the range finder and then extracts the leading singular values and left singular vectors via a small dense SVD.

tensorrsvd.core.rsvd.randomized_range_finder(op, rank, num_oversamples, num_power_iterations, backend)[source]

Compute an orthonormal matrix Q whose range approximates the range of op, i.e., \(op \approx Q Q^\dagger op\)

Parameters:

op (LinearOperator) – The linear operator for which to find the range.
rank (int) – Target rank for the approximation.
num_oversamples (int) – Number of additional random vectors to sample (beyond the target rank) to improve accuracy.
num_power_iterations (int) – Number of power iterations for improved accuracy.
backend (str) – Backend to use: ‘numpy’, ‘jax’, or ‘cupy’.

Returns:

Q – Orthonormal matrix of shape (m, rank + num_oversamples) whose range approximates the range of op, where m is the number of rows of op.

Return type:

ArrayLike

References

[HMT2011]

tensorrsvd.core.rsvd.rsvd_left(op, rank, num_oversamples, num_power_iterations, backend)[source]

Compute the left singular vectors and values of op using randomized SVD.

Returns (U, S) such that op ≈ U @ diag(S) @ V^H (where V^H is not returned)

Parameters:

op (LinearOperator) – The linear operator to decompose.
rank (int) – Target rank for the approximation.
num_oversamples (int) – Number of additional random vectors to sample (beyond the target rank) to improve accuracy.
num_power_iterations (int) – Number of power iterations for improved accuracy.
backend (str) – Backend to use: ‘numpy’, ‘jax’, or ‘cupy’.

Returns:

U (ArrayLike) – Approximate left singular vectors of op, shape (m, rank).
S (ArrayLike) – Approximate singular values of op, shape (rank,).

Return type:

tuple[ArrayLike, ArrayLike]

References

Halko, N., Martinsson, P.-G., & Tropp, J. A. (2011). Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM Review, 53(2), 217–288. https://doi.org/10.1137/090771806