Hyper-reduction#

romtools provides discrete empirical interpolation utilities for approximating high-dimensional model terms from a small set of sampled entries. The examples below apply DEIM and QDEIM to the complete semi-discrete right-hand side of the H2-air flame model rather than to only the reaction term. This keeps the example close to the operation required by a hyper-reduced projection ROM: construct a basis for the full RHS, evaluate only selected entries, and reconstruct the full tensor.

Tensor convention and multistate systems#

The DEIM and QDEIM class APIs use the standard romtools tensor layout. Function snapshots have shape (n_vars, n_dofs, n_snapshots) and tensor-form function bases have shape (n_vars, n_dofs, n_basis). Users therefore do not need to flatten multistate data before constructing a hyper-reduction operator.

Per-state basis construction#

The default basis_mode="per_state" is intended for multistate PDE systems. A separate POD basis is constructed for each state variable,

\[F_i \approx U_i A_i, \qquad i=1,\ldots,n_{vars},\]

and DEIM or QDEIM is applied independently to each \(U_i\). The resulting state-specific spatial sample sets \(S_i\) are combined into one sample mesh,

\[S = \bigcup_i S_i.\]

Every state variable is evaluated on this union mesh. Each state is then reconstructed with all points in the union,

\[f_i \approx U_i (U_i[S,:])^{\dagger} f_i[S].\]

Thus the union acts as an oversampled DEIM system for each state rather than throwing away points selected by another state. sample_indices contains the union of spatial sample points and a sampled RHS has shape (n_vars, n_sample_points). state_sample_indices exposes the individual pre-union DEIM/QDEIM selections for diagnostics.

The per-state formulation avoids combining variables with different physical units in a single POD/SVD. It is also invariant to independent constant scaling of the state variables, and an energy-based truncater can naturally retain a different number of modes for each state. The retained dimensions are available through basis_sizes.

Global coupled basis#

A coupled basis remains available with basis_mode="global". In this mode, the first two tensor axes are collapsed internally using explicit C ordering,

\[F \in \mathbb{R}^{(n_{vars}n_{dofs})\times n_{snapshots}},\]

and one POD/DEIM problem is constructed for the full multistate vector. This can exploit correlations between state variables, but the variables should first be nondimensionalized or weighted consistently so that the Euclidean inner product used by the SVD is physically meaningful.

With a global basis, sample_all_states=True remains the default. A scalar DEIM/QDEIM row selected for any state is converted to its spatial index and all states are sampled at that point. Independent flattened state-DOF sampling is available only in global mode with sample_all_states=False.

All tensor-to-matrix and matrix-to-tensor transformations in the global path use explicit order="C" so that the state-major ordering is unambiguous.

DEIM#

The shortest path through the public API is:

from romtools.hyper_reduction import DEIM
from romtools.vector_space.utils.truncater import BasisSizeTruncater

# rhs_training_snapshots.shape == (n_vars, n_dofs, n_snapshots)
deim = DEIM.from_snapshots(
    rhs_training_snapshots,
    truncater=BasisSizeTruncater(10),
    basis_mode="per_state",
)

# full_rhs.shape == (n_vars, n_dofs)
sampled_rhs = full_rhs[:, deim.sample_indices]
reconstructed_rhs = deim.reconstruct(sampled_rhs)

DEIM.from_snapshots performs a separate POD for each state by default, applies DEIM to each retained basis, unions the selected spatial points, and constructs an oversampled reconstruction operator for each state. With BasisSizeTruncater(10) each state retains ten modes. With an energy-based truncater, the individual states may retain different numbers of modes.

In the H2-air flame example, the four RHS fields are kept as the first tensor axis while the two spatial dimensions are collapsed into the n_dofs axis using C ordering. The training snapshots come from one flame trajectory, while reconstruction is measured on a second trajectory at different parameter values. This avoids measuring only interpolation of the data used to construct the basis.

Run the example with

python examples/h2_air_flame_deim/example.py

Full implementation#

  1import argparse
  2import os
  3import sys
  4
  5import numpy as np
  6
  7from romtools.hyper_reduction import DEIM
  8from romtools.vector_space.utils.truncater import BasisSizeTruncater
  9
 10
 11EXAMPLE_DIR = os.path.abspath(os.path.dirname(__file__))
 12MODELS_DIR = os.path.abspath(os.path.join(EXAMPLE_DIR, "..", "models"))
 13if MODELS_DIR not in sys.path:
 14    sys.path.insert(0, MODELS_DIR)
 15
 16from h2_air_flame import H2AirFlame  # noqa: E402
 17
 18
 19def _rhs_snapshots(model, states, parameters):
 20    snapshots = [
 21        np.reshape(
 22            model.rhs(state, *parameters),
 23            (len(model.field_names), model.nx * model.ny),
 24            order="C",
 25        )
 26        for state in states
 27    ]
 28    return np.stack(snapshots, axis=2)
 29
 30
 31def main(smoke=False):
 32    if smoke:
 33        model = H2AirFlame(
 34            nx=10,
 35            ny=7,
 36            dt=1.0e-4,
 37            t_end=5.0e-4,
 38            snapshot_stride=1,
 39        )
 40        basis_dimension = 4
 41    else:
 42        model = H2AirFlame(
 43            nx=18,
 44            ny=10,
 45            dt=1.0e-4,
 46            t_end=2.0e-3,
 47            snapshot_stride=1,
 48        )
 49        basis_dimension = 10
 50
 51    training_parameters = (2.0, 8.0, 40.0, 7.0)
 52    test_parameters = (2.2, 8.4, 36.0, 8.0)
 53
 54    training_states, _ = model.solve(*training_parameters)
 55    rhs_training_snapshots = _rhs_snapshots(
 56        model, training_states, training_parameters
 57    )
 58
 59    # Use a separate POD/DEIM basis for each field. The state-specific sample
 60    # points are unioned into one spatial sample mesh and all fields are
 61    # evaluated at every selected point.
 62    deim = DEIM.from_snapshots(
 63        rhs_training_snapshots,
 64        truncater=BasisSizeTruncater(basis_dimension),
 65        basis_mode="per_state",
 66    )
 67
 68    test_states, _ = model.solve(*test_parameters)
 69    rhs_test_snapshots = _rhs_snapshots(model, test_states, test_parameters)
 70
 71    relative_errors = []
 72    for snapshot_index in range(rhs_test_snapshots.shape[2]):
 73        rhs = rhs_test_snapshots[:, :, snapshot_index]
 74        sampled_rhs = rhs[:, deim.sample_indices]
 75        reconstructed_rhs = deim.reconstruct(sampled_rhs)
 76        denominator = max(np.linalg.norm(rhs), np.finfo(float).eps)
 77        relative_errors.append(
 78            np.linalg.norm(reconstructed_rhs - rhs) / denominator
 79        )
 80
 81    full_rhs_dimension = rhs_test_snapshots.shape[0] * rhs_test_snapshots.shape[1]
 82    print(f"Full RHS dimension: {full_rhs_dimension}")
 83    print(f"Per-state basis sizes: {deim.basis_sizes}")
 84    print(f"DEIM spatial sample points: {deim.sample_indices.size}")
 85    print(f"Sampled state entries: {rhs_test_snapshots.shape[0] * deim.sample_indices.size}")
 86    print(f"Mean relative RHS reconstruction error: {np.mean(relative_errors):.3e}")
 87    print(f"Max relative RHS reconstruction error: {np.max(relative_errors):.3e}")
 88
 89
 90if __name__ == "__main__":
 91    parser = argparse.ArgumentParser(
 92        description="Reconstruct the H2-air flame RHS with DEIM."
 93    )
 94    parser.add_argument(
 95        "--smoke",
 96        action="store_true",
 97        help="Run a reduced configuration for CI validation.",
 98    )
 99    args = parser.parse_args()
100    main(smoke=args.smoke)

QDEIM#

QDEIM uses the same tensor-native and per-state API. Replacing DEIM with QDEIM changes the sample selection performed on each state basis to pivoted QR:

from romtools.hyper_reduction import QDEIM
from romtools.vector_space.utils.truncater import BasisSizeTruncater

qdeim = QDEIM.from_snapshots(
    rhs_training_snapshots,
    truncater=BasisSizeTruncater(10),
    basis_mode="per_state",
)

sampled_rhs = full_rhs[:, qdeim.sample_indices]
reconstructed_rhs = qdeim.reconstruct(sampled_rhs)

The same spatial sample_indices and tensor reconstruct interface can therefore be used when comparing DEIM and QDEIM in an application.

Run the example with

python examples/h2_air_flame_qdeim/example.py

Full implementation#

 1import argparse
 2import os
 3import sys
 4
 5import numpy as np
 6
 7from romtools.hyper_reduction import QDEIM
 8from romtools.vector_space.utils.truncater import BasisSizeTruncater
 9
10
11EXAMPLE_DIR = os.path.abspath(os.path.dirname(__file__))
12MODELS_DIR = os.path.abspath(os.path.join(EXAMPLE_DIR, "..", "models"))
13if MODELS_DIR not in sys.path:
14    sys.path.insert(0, MODELS_DIR)
15
16from h2_air_flame import H2AirFlame  # noqa: E402
17
18
19def _rhs_snapshots(model, states, parameters):
20    snapshots = [
21        np.reshape(
22            model.rhs(state, *parameters),
23            (len(model.field_names), model.nx * model.ny),
24            order="C",
25        )
26        for state in states
27    ]
28    return np.stack(snapshots, axis=2)
29
30
31def main(smoke=False):
32    if smoke:
33        model = H2AirFlame(
34            nx=10,
35            ny=7,
36            dt=1.0e-4,
37            t_end=5.0e-4,
38            snapshot_stride=1,
39        )
40        basis_dimension = 4
41    else:
42        model = H2AirFlame(
43            nx=18,
44            ny=10,
45            dt=1.0e-4,
46            t_end=2.0e-3,
47            snapshot_stride=1,
48        )
49        basis_dimension = 10
50
51    training_parameters = (2.0, 8.0, 40.0, 7.0)
52    test_parameters = (2.2, 8.4, 36.0, 8.0)
53
54    training_states, _ = model.solve(*training_parameters)
55    rhs_training_snapshots = _rhs_snapshots(
56        model, training_states, training_parameters
57    )
58
59    # Use a separate POD/QDEIM basis for each field. Pivoted QR is applied to
60    # each basis independently, then the selected spatial points are unioned.
61    qdeim = QDEIM.from_snapshots(
62        rhs_training_snapshots,
63        truncater=BasisSizeTruncater(basis_dimension),
64        basis_mode="per_state",
65    )
66
67    test_states, _ = model.solve(*test_parameters)
68    rhs_test_snapshots = _rhs_snapshots(model, test_states, test_parameters)
69
70    relative_errors = []
71    for snapshot_index in range(rhs_test_snapshots.shape[2]):
72        rhs = rhs_test_snapshots[:, :, snapshot_index]
73        sampled_rhs = rhs[:, qdeim.sample_indices]
74        reconstructed_rhs = qdeim.reconstruct(sampled_rhs)
75        denominator = max(np.linalg.norm(rhs), np.finfo(float).eps)
76        relative_errors.append(
77            np.linalg.norm(reconstructed_rhs - rhs) / denominator
78        )
79
80    full_rhs_dimension = rhs_test_snapshots.shape[0] * rhs_test_snapshots.shape[1]
81    print(f"Full RHS dimension: {full_rhs_dimension}")
82    print(f"Per-state basis sizes: {qdeim.basis_sizes}")
83    print(f"QDEIM spatial sample points: {qdeim.sample_indices.size}")
84    print(f"Sampled state entries: {rhs_test_snapshots.shape[0] * qdeim.sample_indices.size}")
85    print(f"Mean relative RHS reconstruction error: {np.mean(relative_errors):.3e}")
86    print(f"Max relative RHS reconstruction error: {np.max(relative_errors):.3e}")
87
88
89if __name__ == "__main__":
90    parser = argparse.ArgumentParser(
91        description="Reconstruct the H2-air flame RHS with QDEIM."
92    )
93    parser.add_argument(
94        "--smoke",
95        action="store_true",
96        help="Run a reduced configuration for CI validation.",
97    )
98    args = parser.parse_args()
99    main(smoke=args.smoke)