Streaming POD vector space tutorial#
In this tutorial you will learn:
How to construct a POD vector space from snapshots that are loaded on demand through a
snapshot_loader, so that the full snapshot matrix is never held in memory at onceThe roles of
block_size,max_basis_dimension, and the basistruncaterHow to use the streaming vector space with no scaling and no orthogonalizer (the defaults)
How to use the streaming vector space with a scaler and an orthogonalizer, and why the two are paired
The math behind streaming POD#
Let the snapshot matrix be \(X \in \mathbb{R}^{N \times M}\), where each column is a state snapshot. For the tensor data used here, \(N = n_{\text{var}} \, n_x\) and \(M\) is the number of snapshots. Proper orthogonal decomposition (POD) seeks a rank-\(k\) orthonormal basis \(\Phi \in \mathbb{R}^{N \times k}\) that best approximates the snapshots,
whose solution is the leading \(k\) left singular vectors of \(X = U \Sigma V^\top\).
For large datasets, \(X\) may not fit in memory. Streaming POD never materializes all of \(X\): it reads one column block \(X_{[s:e]}\) at a time through a snapshot_loader callable, so the memory footprint scales with the block size, not with the total number of snapshots \(M\).
romtools uses a randomized two-pass algorithm. Given a maximum candidate-basis dimension \(r_{\max}\), draw a random matrix \(\Omega \in \mathbb{R}^{M \times r_{\max}}\):
Pass 1 accumulate the range sketch one block at a time,
then orthonormalize \(Y\) to obtain \(Q\) (with \(Q^\top Q = I\)) that approximately spans the range of \(X\).
Pass 2 form the small core matrix, again one block at a time,
compute its (inexpensive) SVD \(B = \tilde{U}\, \Sigma\, V^\top\) and form the candidate POD basis \(\Phi_{\mathrm{candidate}} = Q \tilde{U}\). The supplied truncater then selects the retained columns.
Because \(Q\) and \(\tilde{U}\) both have orthonormal columns, the candidate basis is orthonormal by construction. For fixed-rank truncation, choosing max_basis_dimension larger than the retained dimension provides headroom that improves the randomized range approximation. An energy-based truncater instead chooses the smallest available dimension that reaches its requested fraction of the complete snapshot energy; it raises an error if the maximum candidate dimension is too small. The snapshots are read twice (one pass each); a data-derived scaler (see below) adds one more pass to compute its scalings.
This tutorial runs in serial. For very large problems the same construction is row-distributed across MPI ranks by passing a communicator comm and an SVD functor svdFnc=vector_space.utils.SvdMethodOfSnapshots(comm); each rank’s loader then returns only its local rows for the requested snapshot range.
# First, let's import the relevant modules:
import inspect
import romtools
import numpy as np
from matplotlib import pyplot as plt
from romtools import vector_space
streaming_pod_parameters = inspect.signature(
vector_space.VectorSpaceFromStreamingPOD
).parameters
if 'max_basis_dimension' not in streaming_pod_parameters:
raise RuntimeError(
'This notebook imported an outdated romtools package from '
f'{romtools.__file__}. Reinstall this checkout and restart the kernel.'
)
Loading data through a snapshot_loader#
The streaming vector space does not take a snapshot array directly. Instead it takes a snapshot_loader: a callable loader(start, end) that returns the block of snapshots for the half-open range [start, end). The returned array must be three-dimensional with shape (n_var, n_dofs, end - start) (the last axis is the snapshot axis), and the leading axes must be the same for every call.
In a real large-scale problem the loader would read each block from disk so that the full matrix is never resident in memory. Any callable conforming to romtools.vector_space.utils.SnapshotLoader can provide those blocks. Here, for illustration, we simply slice the in-memory 1D Euler snapshots used throughout these tutorials.
# Load pre-computed snapshots of the 1D Euler equations (from pressio-demo-apps).
snapshots = np.load('snapshots.npz')['snapshots']
# The snapshots are in tensor form:
n_vars, nx, nt = snapshots.shape
print('snapshot tensor shape (n_vars, nx, nt):', snapshots.shape)
# A snapshot_loader returns the 3D block for the half-open range [start, end).
# In practice this would stream each block from disk; here we slice the
# in-memory array purely for illustration.
def my_snapshot_loader(start, end):
return snapshots[..., start:end]
snapshot tensor shape (n_vars, nx, nt): (3, 500, 400)
Case 1: streaming POD with no scaling and no orthogonalizer#
In the simplest case we pass the loader, block_size, max_basis_dimension, and a truncater. Here BasisSizeTruncater(basis_dimension) retains 20 modes from a 25-mode randomized candidate basis. The five-mode headroom improves the randomized range approximation. The scaler and orthogonalizer default to no-ops.
When no shifter is supplied, the shift (affine offset) vector is zero. Case 2 demonstrates a shifter initialized directly from streamed snapshot blocks.
block_size = 8 # maximum number of snapshots held in memory at once
basis_dimension = 20 # number of POD modes to retain
max_basis_dimension = 25 # size of the randomized candidate basis
truncater = vector_space.utils.BasisSizeTruncater(basis_dimension)
my_streaming_vector_space = vector_space.VectorSpaceFromStreamingPOD(
snapshot_loader=my_snapshot_loader,
block_size=block_size,
n_snapshots=nt,
max_basis_dimension=max_basis_dimension,
truncater=truncater)
# We can view the retained basis and the full candidate spectrum:
basis = my_streaming_vector_space.get_basis()
print('The dimension of the vector space is', my_streaming_vector_space.extents())
print('Candidate singular values:', my_streaming_vector_space.get_singular_values())
The dimension of the vector space is (3, 500, 20)
Candidate singular values: [809.24802331 99.93228461 42.50327706 24.91720389 17.46979565
13.78376717 10.96722205 8.77786738 7.20248157 6.72509262
5.79925724 4.86571602 4.61880349 4.27625551 3.68076115
3.47612479 2.90058272 2.72804175 2.51990746 2.45474475
2.16353006 2.02233188 1.8021061 1.63232862 1.54272927]
# We can look at the density component of the first basis vector:
plt.plot(basis[0, :, 0])
plt.xlabel(r'Index')
plt.ylabel(r'$\rho$')
plt.show()
# With no orthogonalizer, the streaming basis is already orthonormal in the
# Euclidean inner product (Phi^T Phi = I) by construction:
is_identity = np.einsum('ijk,ijl->kl', basis, basis)
assert np.allclose(is_identity, np.eye(basis.shape[-1]))
Case 2: streaming POD with shifting, scaling, and orthogonalization#
An average streaming shifter first computes the temporal mean in one loader pass and subtracts it from every subsequent snapshot block. The 1D Euler state also has variables of very different magnitudes (density, momentum, and energy), so a scaler non-dimensionalizes the shifted data before the SVD. Here we use a per-variable VariableScaler("variance").
Scaling changes the basis after the SVD (post_scale), which in general destroys the Euclidean orthonormality of the modes. We therefore pair the scaler with an orthogonalizer to restore \(\Phi^\top \Phi = I\). EuclideanL2Orthogonalizer re-orthonormalizes in the standard \(L^2\) inner product; if you need a physics-based inner product (e.g. weighting by cell volumes or quadrature weights \(w\)), use EuclideanVectorWeightedL2Orthogonalizer(w), which enforces \(\Phi^\top \mathrm{diag}(w)\, \Phi = I\).
Initialization follows the same order as standard POD: shift first, then scale. The average shifter and VariableScaler each perform one initialization pass before the two POD passes, for four passes total. A first-snapshot streaming shifter instead needs only loader(0, 1).
# Create an average shifter, a per-variable scaler, and an orthogonalizer:
my_shifter = vector_space.utils.create_streaming_average_shifter()
my_scaler = vector_space.utils.VariableScaler('variance')
my_orthogonalizer = vector_space.utils.EuclideanL2Orthogonalizer()
my_scaled_vector_space = vector_space.VectorSpaceFromStreamingPOD(
snapshot_loader=my_snapshot_loader,
block_size=block_size,
n_snapshots=nt,
max_basis_dimension=max_basis_dimension,
truncater=vector_space.utils.BasisSizeTruncater(basis_dimension),
shifter=my_shifter,
scaler=my_scaler,
orthogonalizer=my_orthogonalizer)
scaled_basis = my_scaled_vector_space.get_basis()
print('The dimension of the scaled vector space is', my_scaled_vector_space.extents())
The dimension of the scaled vector space is (3, 500, 20)
# Even though scaling was applied, the orthogonalizer guarantees the basis is
# orthonormal in the Euclidean inner product:
is_identity = np.einsum('ijk,ijl->kl', scaled_basis, scaled_basis)
assert np.allclose(is_identity, np.eye(scaled_basis.shape[-1]))
# Compare the density component of the first mode with and without scaling:
plt.plot(basis[0, :, 0], label='no scaling')
plt.plot(scaled_basis[0, :, 0], '--', label='variance scaling')
plt.xlabel(r'Index')
plt.ylabel(r'$\rho$')
plt.legend()
plt.show()
Summary#
VectorSpaceFromStreamingPODbuilds a POD basis while holding only one block of snapshots in memory at a time, loading data on demand through asnapshot_loader.block_sizetrades memory footprint against the number of loader calls;max_basis_dimensionsets the randomized candidate rank, and a truncater selects the retained modes. Leaving headroom above a fixed retained rank improves randomized accuracy.With no scaler and no orthogonalizer, the basis is orthonormal by construction.
Streaming average and first-snapshot shifters initialize their affine offsets directly from loader blocks; shifting occurs before scaling, as in standard POD.
A scaler (e.g.
VariableScaler) should be paired with an orthogonalizer to restore orthonormality after scaling.
See the API documentation for details: