low_rank_toolbox.cssp

Column Subset Selection Problem (CSSP) Algorithms.

This submodule implements various algorithms for selecting representative columns from a matrix.

Algorithms

DEIM : Discrete Empirical Interpolation Method QDEIM : QR-based Discrete Empirical Interpolation Method sQDEIM : Strong QDEIM ARP : Adaptive Randomized Pivoting gpode : Greedy POD with Energy constraint gpodr : Greedy POD with Residual constraint Osinsky : Osinsky’s algorithm oversampling_sQDEIM : Oversampling strong QDEIM

Author: Benjamin Carrel, Paul Scherrer Institute, 2025

low_rank_toolbox.cssp.ARP(U, seed=None, return_projector=False, return_inverse=False, **extra_args)[source]

Implements the Adaptive Randomized Pivoting (ARP) algorithm for Column Subset Selection Problem (CSSP).

Reference: ADAPTIVE RANDOMIZED PIVOTING FOR COLUMN SUBSET SELECTION, DEIM, AND LOW-RANK APPROXIMATION by Cortinovis and Kressner.

Note: The algorithm is similar to Osinsky’s algorithm for the CSSP. The randomization step allows for better error bounds (in expectation). (See Algorithm 2.1)

Parameters:

U (ndarray) – An (n x r) matrix with orthonormal columns.
seed (int, optional) – Random seed for reproducibility.
return_projector (bool, optional) – If True, return also the matrix U @ inv(U[S, :])
return_inverse (bool, optional) – If True, return also the inverse matrix inv(U[S, :])

Return type:

ndarray | tuple

Returns:

J (ndarray) – Selection of r row indices.
P_U (ndarray (n x r) (optional)) – Matrix U @ inv(U[J, :]) where U[J, :] is the (r x r) submatrix. Only returned if return_projector=True.
inv_U (ndarray (r x r) (optional)) – Matrix inv(U[J, :]). Only returned if return_inverse=True (requires return_projector=True).

low_rank_toolbox.cssp.DEIM(U, return_projector=False, return_inverse=False, **extra_args)[source]

DEIM - Discrete empirical interpolation method

Construct a matrix P = [e_p1, e_p2, …, e_pk] where the indexes pi are obtained via the DEIM procedure.

Reference: Nonlinear Model Reduction via Discrete Empirical Interpolation Saifon Chaturantabut and Danny C. Sorensen SIAM Journal on Scientific Computing 2010 32:5, 2737-2764

Parameters:

U (ndarray) – Orthonormal matrix of size n x k
return_projector (bool) – If True, return also the matrix U @ inv(U[S, :])
return_inverse (bool) – If True, return also the inverse matrix inv(U[S, :])
extra_args (dict) –

Additional arguments:

solve_kwargs: dict
Additional arguments for the solve function

Return type:

list | tuple

Returns:

p (list) – List of indexes selected by DEIM
P_U (ndarray (n x k) (optional)) – Projector matrix U @ inv(U[S, :])
inv_U (ndarray (k x k) (optional)) – Inverse matrix inv(U[S, :])

low_rank_toolbox.cssp.Osinsky(U, return_projector=False, return_inverse=False, **extra_args)[source]

Osinsky’s quasi optimal column subset selection algorithm.

Reference:: “Close to optimal column approximations with a single SVD.” by A.I. Osinsky, 2023.

Parameters:

U (numpy.ndarray) – Orthonormal real matrix of shape (n, r) defining a row space approximation
return_projector (bool) – If True, return also the matrix U @ inv(U[S, :])
return_inverse (bool) – If True, return also the inverse matrix inv(U[S, :])
extra_args (dict) –

Additional arguments:

solve_kwargsdict
Additional arguments for the solve function

Return type:

ndarray | tuple

Returns:

J (ndarray) – Selected column indices (r elements)
P_U (ndarray (n x r) (optional)) – Matrix U @ inv(U[J, :]) where U[J, :] is the (r x r) submatrix. Only returned if return_projector=True.
inv_U (ndarray (r x r) (optional)) – Matrix inv(U[J, :]). Only returned if return_inverse=True (requires return_projector=True).

low_rank_toolbox.cssp.QDEIM(U, return_projector=False, return_inverse=False, **extra_args)[source]

QDEIM - QR based DEIM of U (size n x k)

Reference: A new selection operator for the discrete empirical interpolation method - improved a priori error bound and extensions. Zlatko Drmač and Serkan Gugercin. SIAM Journal on Scientific Computing, 38(2), A631-A648.

Original Matlab code from Zlatko Drmač

Parameters:

U (ndarray) – Orthonormal matrix of size n x k
return_projector (bool) – If True, return also the matrix U @ inv(U[S, :])
return_inverse (bool) – If True, return also the matrix inv(U[S, :])
extra_args (dict) –

Additional arguments:

qr_kwargs: dict
Additional arguments for the QR factorization

solve_kwargs: dict
Additional arguments for the solve function

Return type:

ndarray | tuple

Returns:

p (list) – Selection of m row indices with guaranteed upper bound: norm(inv(U[S,:])) <= sqrt(n-k+1) * O(2^m).
P_U (ndarray (n x k) (optional)) – Matrix U @ inv(U[S, :])
inv_U (ndarray (k x k) (optional)) – Matrix inv(U[S, :])

low_rank_toolbox.cssp.gpode(U, oversampling_size, return_projector=False, return_inverse=False, **extra_args)[source]

Gappy POD+E - Greedy algorithm for the selection of m rows of U Minimize the norm of the pseudoinverse of U[S, :] Total cost is O(k^2 + m^2) + QR of U.T.conj() of cost O(nk^2)

Reference: Stability of discrete empirical interpolation and gappy proper orthogonal decomposition with randomized and deterministic sampling points Benjamin Peherstorfer, Zlatko Drmač, and Serkan Gugercin SIAM Journal on Scientific Computing, 42(5), A2837-A2864.

Parameters:

U (ndarray) – Orthonormal matrix of size n x k
oversampling_size (int) – Oversampling size p such that m = k + p
return_projector (bool) – If True, return also the matrix U @ pinv(U[p, :])
return_inverse (bool) – If True, return also the inverse of U[p, :]
extra_args (dict) –

Additional arguments:

qr_kwargs: dict
Additional arguments for the QR factorization

lstsq_kwargs: dict
Additional arguments for the lstsq function

Return type:

tuple | list

Returns:

p (list) – Selection of m row indices
P_U (ndarray (n x k) (optional)) – Matrix U @ pinv(U[p, :])
inv_U (ndarray (k x k) (optional)) – Inverse of U[p, :]

low_rank_toolbox.cssp.gpodr(U, oversampling_size=None, tol=None, max_iter=None, return_projector=False, return_inverse=False, **extra_args)[source]

Gappy POD+R - QDEIM and randomized oversampling.

When tol is given, then p is ignored and the algorithm will select the number of rows to satisfy the condition: sigma_{min}(U[p, :])^{-1} <= tol

Reference: Stability of discrete empirical interpolation and gappy proper orthogonal decomposition with randomized and deterministic sampling points Benjamin Peherstorfer, Zlatko Drmač, and Serkan Gugercin SIAM Journal on Scientific Computing, 42(5), A2837-A2864.

Parameters:

U (ndarray) – Orthonormal matrix of size n x k
oversampling_size (int) – Oversampling size
tol (float) – Tolerance for the quantity sigma_{min}(U[p, :])^{-1}
return_projector (bool) – If True, return also the matrix U @ pinv(U[p, :])
return_inverse (bool) – If True, return also the inverse of U[p, :]
extra_args (dict) –

Additional arguments:

qr_kwargs: dict
Additional arguments for the QR factorization

lstsq_kwargs: dict
Additional arguments for the lstsq function
max_iter (int | None)

Returns:

p (list) – Selection of m row indices
P_U (ndarray (n x k) (optional)) – Matrix U @ pinv(U[p, :])
inv_U (ndarray (k x k) (optional)) – Inverse of U[p, :]

low_rank_toolbox.cssp.oversampling_sQDEIM(U, oversampling_size, tol=None, return_projection=False, return_inverse=False)[source]

Oversampling sQDEIM - Oversampled version of sQDEIM

Reference: ACCURACY AND STABILITY OF CUR DECOMPOSITIONS WITH OVERSAMPLING (Taejun Park and Yuji Nakatsukasa)

Parameters:

U (ndarray) – Orthonormal matrix of size n x k
oversampling_size (int) – Oversampling size p < k such that m = k + p
tol (float) – Tolerance for the strong rank-revealing QR factorization If None, use the rank-revealing QR factorization with eta=2
return_projection (bool) – If True, return also the matrix U @ pseudoinv(U[S, :])
return_inverse (bool) – If True, return also the inverse of U[S, :]

Return type:

ndarray | Tuple[ndarray, ndarray] | Tuple[ndarray, ndarray, ndarray]

Returns:

p (list) – Selection of m = k + oversampling_size row indices.
P_U (ndarray (n x m) (optional)) – Matrix U @ pseudoinv(U[p, :]) where U[p, :] is the (m x k) submatrix. Only returned if return_projection=True.
inv_U (ndarray (k x m) (optional)) – Matrix U.T @ P_U, representing the pseudoinverse relationship. Only returned if return_inverse=True (requires return_projection=True).

low_rank_toolbox.cssp.sQDEIM(U, eta=2, return_projector=False, return_inverse=False, **extra_args)[source]

sQDEIM - Strong RRQR based DEIM of U (size n x k)

Key advantage: the selection of the indexes is guaranteed to satisfy the condition: sigma_{min}(U[p, :])^{-1} <= sqrt(1 + eta * r (n-k))

By default, eta = 2

Parameters:

U (ndarray) – Orthonormal matrix of size n x k
eta (float) – Bounding factor for R_11^{-1} R_12, must be >= 1
mode (str) – Specifies the truncation criterion. Must be ‘rank’ or ‘tol’.
param (int or float) –
The parameter for the chosen mode.
- If mode is ‘rank’, param is the desired rank k.
- If mode is ‘tol’, param is the error tolerance.
return_projector (bool) – If True, return also the matrix U @ inv(U[S, :])
return_inverse (bool) – If True, return also the matrix inv(U[S, :])

Return type:

ndarray | tuple

Returns:

p (list) – Selection of m row indices with guaranteed upper bound: norm(inv(U[S,:])) <= sqrt(n-k+1) * O(2^m).
P_U (ndarray (n x k) (optional)) – Matrix U @ inv(U[S, :])
inv_U (ndarray (k x k) (optional)) – Matrix inv(U[S, :])

Modules

`arp`	Adaptive Randomized Pivoting (ARP) algorithm.
`deim`	Discrete Empirical Interpolation Method (DEIM).
`gpode`(U, oversampling_size[, ...])	Gappy POD+E - Greedy algorithm for the selection of m rows of U Minimize the norm of the pseudoinverse of U[S, :] Total cost is O(k^2 + m^2) + QR of U.T.conj() of cost O(nk^2)
`gpodr`(U[, oversampling_size, tol, max_iter, ...])	Gappy POD+R - QDEIM and randomized oversampling.
`osinsky`	Osinsky's quasi-optimal column subset selection algorithm.
`oversampling_sqdeim`	Oversampling Strong QDEIM algorithm.
`qdeim`	QR-based Discrete Empirical Interpolation Method (QDEIM).
`sqdeim`	Strong QDEIM - Strong Rank-Revealing QR based DEIM.
`utils`	Utility Functions for Column Subset Selection.