Distributed SVD#

DistributedSvd computes a thin singular value decomposition of a two-dimensional matrix distributed by rows across MPI ranks. It uses NumPy for local QR and SVD operations and mpi4py for communication. The global input matrix is never gathered on any rank.

API#

class romtools.linalg.DistributedSvd(comm=None)[source]

Bases: object

NumPy-compatible thin SVD callable for a row-distributed matrix.

Bind an MPI communicator once, then pass this object anywhere a callable compatible with numpy.linalg.svd() is expected, including romtools.vector_space.VectorSpaceFromPOD.

Parameters:

comm (MPI_Comm) – Communicator over which matrix rows are distributed. If None, supported calls delegate to NumPy in serial.

Example

>>> distributed_svd = DistributedSvd(comm)
>>> U_local, s, Vh = distributed_svd(
...     A_local, full_matrices=False, compute_uv=True, hermitian=False
... )
__call__(a, full_matrices=True, compute_uv=True, hermitian=False)[source]

Compute a thin SVD collectively over the bound communicator.

Parameters:

  • a (np.ndarray) – The calling rank’s local matrix rows.

  • full_matrices (bool) – Must be False when compute_uv=True.

  • compute_uv (bool) – Return singular vectors when True; otherwise return only singular values.

  • hermitian (bool) – Must be False.

Returns:

(U_local, s, Vh) when compute_uv=True; otherwise s.

Result distribution#

For a global matrix \(A \in \mathbb{C}^{m \times n}\), let \(k = \min(m,n)\). Each rank supplies its local rows A_local and calls the operation collectively. The result has the following distribution:

  • U_local has shape (A_local.shape[0], k) and the same row distribution as the input.

  • s has shape (k,) and is replicated on all ranks.

  • Vh has shape (k, n) and is replicated on all ranks.

Uneven partitions and ranks with zero local rows are supported. Real and complex numeric dtypes, wide matrices, tall matrices, rank-deficient matrices, and globally empty dimensions are supported. Input arrays and the communicator are not modified.

Supported NumPy options#

The callable has the same argument names and defaults as numpy.linalg.svd. The following behavior is intentional:

  • full_matrices=False, compute_uv=True, hermitian=False returns (U_local, s, Vh).

  • compute_uv=False, hermitian=False returns only the replicated singular values. As in NumPy, full_matrices does not affect this result.

  • full_matrices=True with compute_uv=True raises NotImplementedError because constructing the global full left null space is not supported.

  • hermitian=True raises NotImplementedError.

  • Inputs must be two-dimensional, and every rank must use the same column count, dtype, and options. Invalid distributed metadata raises ValueError collectively.

Algorithm#

The implementation uses a two-level Tall-Skinny QR (TSQR) factorization:

  1. Each rank computes A_local = Q_local @ R_local.

  2. Rank zero gathers only the local R_local factors and performs a second reduced QR factorization of their vertical stack.

  3. Rank zero computes the SVD of the final reduced factor.

  4. The singular values and right singular vectors are broadcast, while the reduced left transformations are scattered and multiplied by each Q_local.

This avoids forming \(A^H A\), which would square the condition number.

MPI example#

The following program can be run with mpiexec -n 4 python distributed_pod.py:

import numpy as np
from mpi4py import MPI

from romtools.linalg import DistributedSvd
from romtools.vector_space import VectorSpaceFromPOD

comm = MPI.COMM_WORLD
rank = comm.Get_rank()

# Every rank owns a potentially different number of matrix rows.
local_rows = rank + 2
rng = np.random.default_rng(1000 + rank)
A_local = rng.normal(size=(local_rows, 6))

distributed_svd = DistributedSvd(comm)
U_local, s, Vh = distributed_svd(
   A_local,
   full_matrices=False,
   compute_uv=True,
   hermitian=False,
)

# VectorSpaceFromPOD accepts the same communicator-bound callable. With one
# variable, the tensor's spatial axis is distributed across ranks.
local_snapshots = A_local.reshape(1, local_rows, 6)
pod_space = VectorSpaceFromPOD(
   local_snapshots,
   svdFnc=distributed_svd,
)
local_basis = pod_space.get_basis()