librla - Randomized Linear Algebra Library

A unified, multi-language library implementing randomized algorithms for low-rank matrix factorizations of both real and complex matrices. The library provides efficient sketching-based methods for large-scale matrix decompositions with consistent APIs across Python, MATLAB/Octave, and Julia.

Features

Core Algorithms

All algorithms support both tolerance mode (rtol < 1) for adaptive rank selection and rank mode (rtol >= 1) for fixed-rank approximations:

LinearOperator Support

All algorithms support LinearOperator abstraction for matrix-free computation:

• Dense matrices - Standard NumPy/MATLAB/Julia arrays
• Explicit LinearOperators - Matrix wrappers with unified interface
• Matrix-free LinearOperators - Function handles for A*x and A'*x operations only

Matrix-free operators enable sketching of implicit matrices (FFT, convolution, Toeplitz, circulant, etc.) without explicit storage.

Method Options for ID

The id_sketch and id_qrpiv functions support three methods for computing the interpolation matrix T:

Quick Start

Python

For Python, it is often a good idea to run in a virtual environment. A script named 'setup-venv.sh' is available for users. To utilize it, run the shell script in the terminal. To end the virtual environment, type 'deactivate' in the terminal.

The package librla can be used by either calling the .py file or using the following pip install command:

pip install librla

Here is an example using the library.

import numpy as np
from librla import orth_sketch, qr_sketch, svd_sketch, id_sketch

# Create a test matrix (Hilbert matrix is ill-conditioned)
A = 1.0 / (np.arange(1, 1001)[:, None] + np.arange(1, 501) - 1)

# Orthonormal basis with relative tolerance
Q, flag, diagR = orth_sketch(A, rtol=1e-6)

# Truncated QR factorization
Q, R, p = qr_sketch(A, rtol=1e-6)
# A[:, p] approx Q @ R

# Truncated SVD
U, s, Vh = svd_sketch(A, rtol=1e-6)
# A approx U @ np.diag(s) @ Vh

# Interpolative decomposition
k, piv, T = id_sketch(A, rtol=1e-6, method='lstsq')
# A[:, piv[k:]] approx A[:, piv[:k]] @ T

Usage Modes

Tolerance Mode (rtol < 1)

Adaptively determines rank to achieve specified relative tolerance:

# Python: Adaptive rank selection to achieve 10^-6 tolerance
Q, flag, diagR = orth_sketch(A, 1e-6)

Rank Mode (rtol >= 1)

Returns fixed-rank approximation:

# Python: Rank-20 approximation
U, s, Vh = svd_sketch(A, 20)

API Notes

svd_sketch Return Values

The svd_sketch function returns slightly different formats across languages:

Python returns V transposed/adjoint

Indexing

• Python: 0-based indexing (piv[:k] for skeleton columns)

Algorithm Details

Randomized Sketching (id_sketch)

• Uses random test matrix multiplication for fast column space approximation
• Geometric block growth for adaptive rank determination in tolerance mode
• Optional power iterations for improved accuracy (power_iter parameter)
• Stochastic (results vary slightly between runs)

Deterministic QR Pivoting (id_qrpiv)

• Uses LAPACK geqp3 column-pivoted QR factorization
• Deterministic and reproducible results
• Slower than randomized sketching but guaranteed behavior
• Same interface as id_sketch
• Useful for verification and when reproducibility is critical

Accuracy of Randomized Methods

Randomized sketching algorithms have inherent variance in reconstruction error. In rank mode (rtol >= 1), the reconstruction error of randomized methods is typically within a small factor of the optimal (deterministic) error. The validation tests use these thresholds:

This variance is expected for randomized algorithms and does not indicate a bug.

For matrices with slowly decaying singular values (small spectral gap), use power iterations to improve accuracy:

# Use power_iter=2 for matrices with slowly decaying singular values
U, s, Vh = svd_sketch(A, 20, power_iter=2)

LinearOperator Usage

Create matrix-free operators for implicit matrices:

Python

Python uses scipy.sparse.linalg.LinearOperator:

import numpy as np
from scipy.sparse.linalg import LinearOperator
from librla import orth_sketch

# Define matrix-vector products
def matvec(x):
    return np.fft.fft(x)

def rmatvec(x):
    return np.fft.ifft(x).real

# Create LinearOperator
n = 1000
A_op = LinearOperator(shape=(n, n), matvec=matvec, rmatvec=rmatvec, dtype=complex)

# Use with librla
Q, flag, diagR = orth_sketch(A_op, 20)  # Rank-20 approximation

Note: Matrix-free LinearOperators support both modes. In tolerance mode (rtol < 1) the deterministic fallback may materialize the operator as a dense matrix, one matvec per column (O(n) matvecs) — prefer rank mode when that cost is unacceptable.

Optional Parameters

All sketching functions support these optional parameters:

For id_sketch and id_qrpiv only:

Example with optional parameters:

# Python: Use power iterations for better accuracy
k, piv, T = id_sketch(A, 1e-6, power_iter=2, method='svd')

Requirements

Installation

Python

Requirements: Python 3.9+, NumPy >= 1.20, SciPy >= 1.7

Option 1: Add to Python path (recommended for testing)

import sys
sys.path.append('/path/to/distrib/python')
from librla import id_sketch, svd_sketch

Option 2: Add to PYTHONPATH environment variable

export PYTHONPATH="/path/to/distrib/python:$PYTHONPATH"

Option 3: Copy files to your project

cp distrib/python/librla.py your_project/

Troubleshooting:

• ImportError: No module named 'librla' - Verify path is correct and use absolute paths
• ImportError: No module named 'numpy' - Run pip install numpy scipy

Demos

Run the demos to see the library in action:

# Python
cd distrib/python/demo
python demo01_basic.py
python demo02_svd.py

The pip install librla wheel contains only the library module. To obtain the demos and tests after installing from PyPI, either clone the GitHub repository or download and unpack the source distribution:

pip download --no-binary :all: --no-deps librla
tar xf librla-*.tar.gz
cd librla-*/demo && python demo01_basic.py
cd ../test && python test_all.py

Examples

The demo suite provides examples for:

Tests

Tests for validating that the algorithms are working are:

Additional resources

Additional resources provided for users are:

• The folder '/distrib/compare' provides codes that compare 'librla' with other randomized linear algebra packages:
  • compare_id_scipy.py - Python: librla vs SciPy (ID)
  • compare_svd_torch.py - Python: librla vs PyTorch (SVD)
• The repo librla_appliations contains examples of using 'librla' for data compression. The examples are taken from common applications of interest and are taken from the randomized linear algebra literature.