2.1. Hyper-reduction and Galerkin projection at a glance

Let us consider the simple PDE problem:

(2.1) $$ {\partial}_tu\left(x,t\right)=f\left(u\left(x,t\right)\right), $$

where $ u\left(x,t\right) $ is a scalar dependent on space and time, while $ f:u\to f(u) $ is a nonlinear function, possibly involving space derivatives of solution $ u $ . The general procedure is the following:

1. 1. Finite element discretization. Introducing a basis of local functions $ \Theta (x)=\left({\theta}_1(x),\cdots, {\theta}_N(x)\right)\in {\unicode{x211D}}^{1\times N} $ where $ N $ is the total number of degrees of freedom. Introducing the unknown $ u\left(x,t\right)=\Theta (x)\mathbf{u}(t) $ , where $ \mathbf{u}(t)\in {\unicode{x211D}}^{N\times 1} $ , the finite element methodology allows one to write the weak-form:

   (2.2) $$ \mathbf{M}{\partial}_t\mathbf{u}=\mathbf{F}\left(\mathbf{u}\right), $$
   where
   (2.3) $$ \mathbf{F}\left(\mathbf{u}\right)=\sum \limits_{\mathrm{\ell}}\Theta {\left({x}_{\mathrm{\ell}}\right)}^Tf\left(\Theta \left({x}_{\mathrm{\ell}}\right)\mathbf{u}\right){\delta}_{\mathrm{\ell}},\hskip1.48em \mathbf{M}=\sum \limits_{\mathrm{\ell}}\Theta {\left({x}_{\mathrm{\ell}}\right)}^T\Theta \left({x}_{\mathrm{\ell}}\right){\delta}_{\mathrm{\ell}} $$
   with $ \left({x}_{\mathrm{\ell}},{\delta}_{\mathrm{\ell}}\right) $ being $ \tilde{m}\sim N $ quadrature points in the mesh.

2. 2. Generation of snapshots. The finite element code allows to generate $ P $ snapshots stacked in $ \mathbf{U}=\left({\mathbf{u}}_1,\cdots, {\mathbf{u}}_P\right)\in {\unicode{x211D}}^{N\times P} $ . These snapshots are rewritten as $ {\mathbf{u}}_p=\overline{\mathbf{u}}+{\mathbf{u}}_p^{\prime } $ , where $ \overline{\mathbf{u}} $ is an arbitrary flow-field (lifting function) containing the non-homogeneous boundary conditions (for example, the time-average of the snapshots $ {\mathbf{u}}_p $ or a fixed point of $ f\left(\cdot \right)=0 $ ). The homogeneous snapshots are then stacked in $ {\mathbf{U}}^{\prime }=\left({\mathbf{u}}_1^{\prime}\cdots {\mathbf{u}}_P^{\prime}\right)\in {\unicode{x211D}}^{N\times P} $ .

3. 3. Generation of POD basis. In a classical way, we generate POD modes by solving:

   (2.4) $$ \left({{\mathbf{U}}^{\prime}}^T\mathbf{M}{\mathbf{U}}^{\prime}\right)\mathbf{T}=\mathbf{T}\Lambda, $$
   where $ \mathbf{T} $ are the chronos modes and $ \Lambda $ is a diagonal matrix containing positive eigenvalues ranked in decreasing order. We can then select the leading $ r\le P $ eigenvectors/eigenvalues and compute $ r $ POD modes $ \Phi \in {\unicode{x211D}}^{N\times r} $ , with $ r\ll N $ :
   (2.5) $$ \Phi ={\mathbf{U}}^{\prime}\mathbf{T}{\Lambda}^{-1/2},\hskip0.48em {\Phi}^T\mathbf{M}\Phi =\mathbf{I}. $$

4. 4. Projection of the finite element solver on the POD basis:

   (2.6) $$ {\partial}_t\mathbf{a}=\sum \limits_{\mathrm{\ell}}\Phi {\left({x}_{\mathrm{\ell}}\right)}^Tf\left(\overline{u}\left({x}_{\mathrm{\ell}}\right)+\Phi \left({x}_{\mathrm{\ell}}\right)\mathbf{a}\right){\delta}_{\mathrm{\ell}}, $$
   where $ \mathbf{u}(t)=\Phi \mathbf{a}(t) $ in which $ \mathbf{a}(t)\in {\unicode{x211D}}^{r\times 1} $ is the dynamics within the POD basis, $ \Phi (x)=\Theta (x)\Phi =\left({\phi}_1(x),\cdots, {\phi}_r\right)\in {\unicode{x211D}}^{1\times r} $ and $ \overline{u}(x)=\Theta (x)\overline{\mathbf{u}} $ . This model still scales on the $ \tilde{m}\approx N $ and is therefore not yet a reduced-order model.

Then, we have two cases, depending on the reduction technique: either hyper-reduction (Section 2.1.1), which is valid in the general case, or Galerkin projection (Section 2.1.2), which is valid only for polynomial nonlinearities.

2.2. NNLS algorithm and considerations

The ECSW hyper-reduction method relies on solving large-scale NNLS problems to identify sparse quadrature weights. As established in Section 2.1, each hyper-reduction problem has the form:

$$ \underset{\boldsymbol{\omega} \ge 0}{\min }{\left\Vert \mathcal{A}\boldsymbol{\omega } -\mathbf{b}\right\Vert}_2^2, $$

where $ \mathcal{A}\in {\mathrm{\mathbb{R}}}^{rN_{\mathrm{snap}}\times \tilde{m}} $ contains nonlinear evaluations at Gauss points, and $ \boldsymbol{\omega} \in {\mathrm{\mathbb{R}}}^{\tilde{m}} $ is the empirical quadrature weight vector. In our application:

• • The reduced dimension is typically $ r\sim \mathcal{O}(10) $ ;

• • The number of snapshots is $ {N}_{\mathrm{snap}}\sim \mathcal{O}\left({10}^3\right) $ ;

• • The number of Gauss points is $ \tilde{m}\sim \mathcal{O}\left({10}^5\right) $ .

This leads to NNLS problems of size $ \sim {10}^4\times {10}^5 $ , which are challenging to even store in memory—let alone solve using dense methods. Indeed, the classical Lawson and Hanson (Reference Lawson and Hanson1995) active set method for NNLS has worst-case complexity:

$$ \mathcal{O}\left({rN}_{\mathrm{snap}}\hskip0.1em {\tilde{m}}^2+{\tilde{m}}^3\right), $$

which, when inserting our scaling, gives:

$$ \mathcal{O}\left(\left(10\cdot {10}^3\right)\cdot {\left({10}^5\right)}^2+{\left({10}^5\right)}^3\right)=\mathcal{O}\left({10}^{13}+{10}^{15}\right). $$

Such complexity renders direct solvers impractical for our setting.

2.3. Online evaluation after hyper-reduction

Let us now summarize the structure and complexity of the reduced-order model during the online phase after hyper-reduction.

Recall the full-order Navier–Stokes system in semi-discrete perturbed form $ {\mathbf{w}}_t=\overline{\mathbf{w}}+\mathbf{w} $ with

(2.11) $$ \mathbf{M}\hskip0.1em {\partial}_t\mathbf{w}=\mathcal{L}\mathbf{w}+\mathcal{B}\left(\mathbf{w},\mathbf{w}\right), $$

where $ \mathbf{M} $ is the mass matrix, $ \mathcal{L} $ the discrete linear operator, and $ \mathcal{B}\left(\cdot, \cdot \right) $ the quadratic convection operator.

The reduced-order model obtained via projection onto the reduced basis $ \Phi $ reads:

(2.12) $$ {\partial}_t\mathbf{a}={\Phi}^{\top}\mathrm{\mathcal{L}\Phi}\hskip0.1em \mathbf{a}+{\Phi}^{\top}\mathcal{B}\left(\Phi \mathbf{a},\Phi \mathbf{a}\right), $$

with $ \mathbf{a}\in {\mathrm{\mathbb{R}}}^r $ the reduced coordinates. When the nonlinear term is treated explicitly, it becomes a time-dependent source term:

$$ {\mathbf{f}}_{\mathrm{red}}(t):= {\Phi}^{\top}\mathcal{B}\left(\Phi {\mathbf{a}}^{t-\Delta t},\Phi {\mathbf{a}}^{t-\Delta t}\right). $$

Using the matrix formulation derived earlier and applying ECSW, this nonlinear term is computed via:

(2.13) $$ {\displaystyle \begin{array}{c}{\mathbf{f}}_{\mathrm{red}}={\Phi}^{\top }{\Theta}_u^{\top }W\left({\Theta}_u\Phi \hskip0.1em \mathbf{a}\odot {\partial}_x{\Theta}_u\Phi \hskip0.1em \mathbf{a}+{\Theta}_v\Phi \hskip0.1em \mathbf{a}\odot {\partial}_y{\Theta}_u\Phi \hskip0.1em \mathbf{a}\right)\\ {}+{\Phi}^{\top }{\Theta}_v^{\top }W\left({\Theta}_u\Phi \hskip0.1em \mathbf{a}\odot {\partial}_x{\Theta}_v\Phi \hskip0.1em \mathbf{a}+{\Theta}_v\Phi \hskip0.1em \mathbf{a}\odot {\partial}_y{\Theta}_v\Phi \hskip0.1em \mathbf{a}\right).\end{array}} $$

2.3.2. Time integration

Assuming a backward Euler scheme and denoting by $ \mathbf{A}={\Phi}^{\top}\mathrm{\mathcal{L}\Phi } $ the reduced linear operator, the time integration step becomes:

$$ {\mathbf{a}}^{t+\Delta t}={\left(I-\Delta t\mathbf{A}\right)}^{-1}\left({\mathbf{a}}^t+\Delta t\hskip0.1em {\mathbf{f}}_{\mathrm{red}}(t)\right). $$

The cost of solving this reduced system scales as $ \mathcal{O}\left({r}^3\right) $ for a direct solve, but can be reduced to $ \mathcal{O}\left({r}^2\right) $ using an lower-upper (LU) decomposition pre-computed once and reused at each time step.

3.1. Test case selection and setup

To validate the proposed hyper-reduction methodology, we consider the fluidic pinball configuration. This geometry consists of three identical cylinders placed at the vertices of an equilateral triangle, generating a rich variety of wake dynamics as the Reynolds number increases. It provides a compact yet versatile benchmark that captures essential features of vortex shedding, symmetry breaking, and nonlinear coupling between wake modes, while remaining computationally manageable for systematic analysis.

The validation test case corresponds to the periodic regime of the fluidic pinball, obtained at $ \mathit{\operatorname{Re}}=30 $ . Dynamically, this configuration behaves similarly to the classical cylinder flow in the Von Kármán vortex shedding regime. A single trajectory is simulated from the steady base flow toward the periodic limit cycle, and $ 1000 $ snapshots are extracted along this trajectory. These snapshots are used both to construct the reduced POD basis and to assemble the ECSW NNLS hyper-reduction problem. The resulting hyper-reduced model is then evaluated on the same trajectory to verify correct implementation and basic dynamical consistency.

This periodic configuration serves as a controlled validation case. It is simple enough that any properly implemented reduced-order model should reproduce the correct limit cycle dynamics, yet dynamic enough to test stability and temporal accuracy. The regime’s transient growth rate, shedding frequency, and limit cycle amplitude are well documented in the literature and have been used extensively in earlier work to benchmark ROM performance, enabling direct consistency checks with earlier results.

Only the periodic regime is considered here, as it provides a transparent and interpretable sanity check confirming that the ECSW implementation reproduces expected reduced-order behavior. Since the goal of this work is to analyze computational cost rather than model accuracy, extending the study to more complex or chaotic regimes would not provide additional insight into numerical performance. The following sections, therefore, focus exclusively on quantifying the computational efficiency and scaling properties of ECSW hyper-reduction (Figure 1).

Figure 1.

Mesh for the fluidic pinball test case, featuring three cylinders in an equilateral triangle configuration.

Figure 1. Long description

A two-dimensional triangular mesh is displayed on a Cartesian coordinate system with X-coordinate from negative 6 to 20 and Y-coordinate from negative 6 to 6.

* At the center-left of the domain, three white circular voids representing cylinders are arranged in an equilateral triangle. One cylinder is centered at X equals negative 1.5, Y equals 0, while the other two are positioned further right at X equals 0, with one at Y equals 1 and the other at Y equals negative 1.

* The mesh density is highest immediately surrounding these cylinders, where the triangular elements are smallest and appear as a solid blue mass.

* A rectangular region of high-density refinement extends from X equals negative 3 to X equals 10, spanning Y equals negative 2 to 2, capturing the immediate wake area.

* A second, slightly less dense refinement zone continues from X equals 10 to the right boundary at X equals 20, spanning the same vertical range.

* The outermost regions of the domain, particularly toward the top, bottom, and far-left boundaries, feature the coarsest mesh with large, clearly visible triangular elements.

3.2. Validating the implementation on a periodic test case

We evaluate the hyper-reduction method on the periodic regime of the fluidic pinball. This configuration provides a controlled yet representative test case: the dynamics evolve from a steady base flow toward a stable limit cycle, allowing us to assess both transient behavior and long-term accuracy.

The reduced-order model is constructed using a quadratic Galerkin formulation with a basis size chosen as the minimal number of modes required to capture both the transient evolution and the limit cycle dynamics accurately. In this regime, the reference ROM provides a qualitatively correct representation of the full-order dynamics, both in terms of growth toward the attractor and its final periodic behavior.

We then investigate whether the ECSW hyper-reduction is able to recover a similar level of accuracy under different sparsity regimes. Two configurations are considered, corresponding to different tolerances in the NNLS greedy selection of quadrature points.

In the first case, a relative residual threshold of $ \varepsilon =7.5\cdot {10}^{-2} $ is used. This results in retaining 50 Gauss points for the $ u $ -component and 54 for the $ v $ -component in the ECSW method. As shown in Figure 2, the hyper-reduced model at this size does not match the performance of the quadratic ROM or the FOM. It appears to overdamp the system, resulting in a relatively accurate transient but an incorrect limit cycle.

Figure 2.

First three POD coordinates over time for the periodic test case. Comparison between full-order solution (black), quadratic Galerkin ROM (blue), and hyper-reduced model with NNLS tolerance $ \varepsilon =7.5\cdot {10}^{-2} $ (red).

Figure 2. Long description

The figure consists of three vertically aligned panels sharing a common x-axis labeled Time, ranging from 20 to 140. The y-axes represent Mode 1, Mode 2, and Mode 3 respectively, with values ranging from approximately -10 to 10. A legend in the top panel identifies three data series: a solid black line for F O M projected, a solid blue line for Quadratic R O M with 18 modes, and a dashed red line for E C S W R O M with 18 modes and N N L S threshold 0.075.

* Top Panel (Mode 1): All three models show an initial oscillating growth from time 20 to 60. Beyond time 60, the black and blue lines stabilize into a high-amplitude oscillation around a mean of 9. The dashed red line deviates, stabilizing at a lower mean of approximately 7.5 with smaller oscillations.

* Middle Panel (Mode 2): All models show increasing oscillatory amplitude until time 60. From time 60 to 140, the black and blue lines overlap perfectly, oscillating between -11 and 11. The dashed red line shows a phase shift and slightly lower peak amplitudes compared to the other two.

* Bottom Panel (Mode 3): Similar to Mode 2, the black and blue lines show identical high-frequency oscillations between -11 and 11 after time 60. The dashed red line exhibits a noticeable phase lead and reduced amplitude peaks relative to the F O M and Quadratic R O M results.

In the second configuration, the relative tolerance is tightened to $ \varepsilon =5\cdot {10}^{-2} $ , leading to 65 Gauss points for both the $ u $ and $ v $ components. Notably, the two selected sets share only a single common index out of 130 unique points, confirming that the nonlinear terms for each velocity component are supported in largely distinct spatial regions and justifying the use of separate ECSW decompositions. In this case (Figure 3), the hyper-reduced model recovers dynamics that are qualitatively indistinguishable from the quadratic ROM, both during the transient growth and on the limit cycle. This indicates that a moderate increase in quadrature size is sufficient to fully restore the reduced model’s accuracy.

Figure 3.

Same setup as Figure 2, but with NNLS tolerance $ \varepsilon =5\cdot {10}^{-2} $ . The hyper-reduced model and quadratic ROM are nearly indistinguishable.

Figure 3. Long description

A multi-panel line graph with three vertical sections sharing a common x-axis labeled Time, ranging from 20 to 140. The y-axes represent Mode 1, Mode 2, and Mode 3. A legend in the top panel identifies three data series: a solid black line for F O M projected, a solid blue line for Quadratic R O M 18 modes, and a dashed red line for E C S W R O M 18 modes, N N L S threshold equals 0.05.

* The top panel for Mode 1 shows an initial value near 0, followed by a sigmoid-like increase between time 40 and 80, reaching a plateau around 10. The curve exhibits small, high-frequency oscillations throughout the ascent and plateau. The blue and black lines are perfectly overlapped, while the red dashed line shows a slight lag during the steep ascent.

* The middle panel for Mode 2 shows a centered oscillation around y equals 0. The amplitude increases linearly from time 20 to 60, then stabilizes with peaks at positive 11 and troughs at negative 11. All three lines are nearly indistinguishable.

* The bottom panel for Mode 3 shows a similar oscillating pattern to Mode 2, with amplitude growth stabilizing around time 60. The oscillations are centered at 0 with a range between positive 11 and negative 11. The red dashed line shows very minor deviations from the solid lines during the amplitude growth phase.

These results validate ECSW on this classical benchmark and identify the minimal quadrature size required for an accurate representation of this test case with the chosen POD quadratic ROM.

From a computational perspective, we compare the wall-clock time required for full trajectory integration on CPU. The quadratic ROM requires a total simulation time of $ 1.60 $ seconds. The hyper-reduced models significantly reduce this cost, with total simulation times of $ 0.76 $ seconds and $ 0.81 $ seconds for the two configurations, respectively. More precisely, when isolating the right-hand side evaluation, the computational cost drops from $ 5.2\cdot {10}^{-4} $ s to $ 1.45\cdot {10}^{-4} $ s and $ 1.95\cdot {10}^{-4} $ s. This corresponds to an overall speed-up of approximately $ 2\times $ , and up to $ 2.6\times $ when considering only the nonlinear term, while maintaining comparable accuracy in the refined regime.

3.3. Assessment of numerical cost for the considered architecture

One of the main advantages of hyper-reduction is its lower numerical complexity compared to the tensorial quadratic ROM method. In earlier standalone tests on CPU, we showed clear gains, with the hyper-reduced model running about twice as fast in wall-clock time for the chosen number of retained modes and Gauss points. In the present experiments, however, the context is different: the models now run directly on the GPU and are optimized for this hardware. Under these conditions, the results are surprising. Despite using the same reduction sizes as before, the hyper-reduced model no longer shows the expected advantage. To understand this, we measured the runtime of the quadratic ROM alone and the hyper-reduced ROM alone, isolating their computation while keeping the same optimization strategy for the algorithm and hardware. For the full trajectory (18,001 time steps), the quadratic ROM integrates in about 6.1 s, with each nonlinear evaluation costing $ 153\hskip0.22em \mu \mathrm{s} $ on average. By contrast, the hyper-reduced model requires 6.9 s, with nonlinear evaluations taking $ 208\hskip0.22em \mu \mathrm{s} $ each on average. These results run counter to both previous empirical evidence and the theoretical scaling of the method, which should favor hyper-reduction. Section 4 investigates this discrepancy by analyzing cost scaling in more detail and examining how GPU hardware influences the benefits of hyper-reduction.

Although the isolated nonlinear evaluation is about 36% more time-consuming for the hyper-reduced model ( $ 208 $ vs $ 153\hskip0.22em \mu \mathrm{s} $ ), the full trajectory cost increases by only about 13% ( $ 6.9 $ vs $ 6.1\hskip0.1em $ s). This is because the nonlinear right-hand side is only one component of each time step: the remainder is spent on the implicit $ r\times r $ linear solve and on time-stepping bookkeeping, both of which are identical for the two models. Indeed, over the $ 18\hskip0.22em 001 $ steps the right-hand-side evaluation accounts for $ 2.75\hskip0.1em $ s for the quadratic ROM and $ 3.74\hskip0.1em $ s for the hyper-reduced model, leaving $ 3.35\hskip0.1em $ s and $ 3.16\hskip0.1em $ s respectively for the shared remaining workload—consistent up to timing noise, as expected for model-independent operations. The right-hand-side penalty is therefore amortized over these shared costs.

4.1. Theoretical operation count scaling for nonlinearity computation

Before examining wall-clock timings for the tensorized quadratic ROM and the hyper-reduced ROM, we first estimate the exact number of operations required for both approaches. Earlier, we described their asymptotic complexity: $ \mathcal{O}\left({r}^3\right) $ for the tensorial ROM and $ \mathcal{O}(mr) $ for the hyper-reduction, where $ r $ is the size of the retained POD subspace and $ m $ the number of retained Gauss points in ECSW. Here, we go beyond asymptotics and compute the exact non-asymptotic operation counts, which provide a hardware-agnostic baseline for comparing methods and identifying where the trade-off between them shifts.

4.1.1. Numerical cost of the quadratic ROM nonlinearity evaluation

The quadratic ROM computes the quadratic nonlinearity entirely in the reduced space. As expressed in Section 2.1.2, the explicit nonlinear term for one unbatched evaluation is

$$ {f}_k=\sum \limits_{i=1}^r\sum \limits_{j=1}^r{N}_{ijk}\hskip0.1em {a}_i\hskip0.1em {a}_j,\hskip2em N\in {\unicode{x211D}}^{r\times r\times r},\hskip2em a\in {\unicode{x211D}}^r. $$

In practice, this is not implemented with triple loops but with the torch.einsum operator, which supports batching and calls optimized basic linear algebra subprograms (BLAS) kernels. The batched version can be written as

$$ {f}_{lk}=\sum \limits_{i=1}^r\sum \limits_{j=1}^r{N}_{ijk}\hskip0.1em {a}_{li}\hskip0.1em {a}_{lj},\hskip2em a\in {\unicode{x211D}}^{b\times r},\hskip2em \mathbf{f}\in {\unicode{x211D}}^{b\times r}. $$

with $ b $ the batch size.

The contraction is typically performed in two steps. First, compute the intermediate tensor

$$ {T}_{ljk}=\sum \limits_{i=1}^r{N}_{ijk}\hskip0.1em {a}_{li}. $$

This requires $ {br}^3 $ multiplications and $ {br}^3-{br}^2 $ additions.

Then compute

$$ {f}_{lk}=\sum \limits_{j=1}^r{T}_{lj k}\hskip0.1em {a}_{lj}, $$

which requires $ {br}^2 $ multiplications and $ {br}^2- br $ additions.

In total, the quadratic ROM costs

$$ \left({br}^3+{br}^2\right)\;\mathrm{multiplications},\hskip2em \left({br}^3- br\right)\;\mathrm{additions}, $$

for an overall count of $ 2{br}^3+{br}^2- br $ operations (not accounting for MAC optimization).

4.1.2. Numerical cost of hyper-reduction nonlinearity evaluation

For hyper-reduction, the nonlinear evaluation is written as in Section 2.3:

(4.1) $$ {\displaystyle \begin{array}{c}\mathbf{f}={M}_1\left[\left({M}_3\mathbf{a}\right)\odot \left({M}_4\mathbf{a}\right)+\left({M}_5\mathbf{a}\right)\odot \left({M}_6\mathbf{a}\right)\right]\\ {}+{M}_2\left[\left({M}_7\mathbf{a}\right)\odot \left({M}_8\mathbf{a}\right)+\left({M}_9\mathbf{a}\right)\odot \left({M}_{10}\mathbf{a}\right)\right]\end{array}} $$

Where the shorthand matrices are defined from the pre-computed products of the POD basis and FEM sampling operators (see Section 2.3 and Appendix A): $ {M}_1={\Phi}^{\top }{\Theta}_u^{\top }W\in {\unicode{x211D}}^{r\times m} $ , $ {M}_2={\Phi}^{\top }{\Theta}_v^{\top }W\in {\unicode{x211D}}^{r\times m} $ , and the eight projection matrices $ {M}_3,\dots, {M}_{10}\in {\unicode{x211D}}^{m\times r} $ correspond to $ {\Theta}_u\Phi $ , $ {\partial}_x{\Theta}_u\Phi $ , $ {\Theta}_v\Phi $ , $ {\partial}_y{\Theta}_u\Phi $ , $ {\Theta}_u\Phi $ , $ {\partial}_x{\Theta}_v\Phi $ , $ {\Theta}_v\Phi $ , $ {\partial}_y{\Theta}_v\Phi $ , respectively. In batched form:

(4.2) $$ \mathbf{f}=\left[\left(\mathbf{a}{M}_3^{\top}\right)\odot \left(\mathbf{a}{M}_4^{\top}\right)+\left(\mathbf{a}{M}_5^{\top}\right)\odot \left(\mathbf{a}{M}_6^{\top}\right)\right]{M}_1^{\top }+\left[\left(\mathbf{a}{M}_7^{\top}\right)\odot \left(\mathbf{a}{M}_8^{\top}\right)+\left(\mathbf{a}{M}_9^{\top}\right)\odot \left(\mathbf{a}{M}_{10}^{\top}\right)\right]{M}_2^{\top }. $$

Each of the six matrix–vector products costs $ brm $ multiplications and $ b\left(r-1\right)m $ additions. The element-wise products add $ bm $ multiplications, and the element-wise sums add $ bm $ additions. Finally, the two projections back to reduced space each cost $ brm $ multiplications and $ br\left(m-1\right) $ additions.

In total, the hyper-reduction costs

$$ \left(8 brm+4 bm\right)\;\mathrm{multiplications},\hskip2em \left(8 brm-4 bm- br\right)\;\mathrm{additions}, $$

that is $ 16 brm- br $ operations overall.

4.1.3. Trade-off between methods

From these counts, hyper-reduction is cheaper whenever

$$ m\le \frac{r\left(2r+1\right)}{16}. $$

For example, with $ r=20 $ , hyper-reduction is cheaper if $ m<50 $ ; for $ r=40 $ , if $ m<200 $ ; and for $ r=100 $ , if $ m<1256 $ .

These thresholds are broadly consistent with earlier empirical results, although small discrepancies remain. For instance, the 20-mode ROM was observed to be more costly than the 20-mode hyper-reduction with 200 Gauss points, despite the theoretical bound. Such deviations can be explained by hardware effects, parallelization strategies, and optimization details. For a more robust view of performance, Subsection 4.2 reports averaged wall-clock results for isolated nonlinear computations (Table 1).

Table 1.

Exact operation counts for quadratic nonlinearity evaluation using quadratic ROM and ECSW hyper-reduction

Table 1. Long description

The table consists of three columns and four rows including the header.

Column 1: Operation Type

Column 2: Quadratic R O M

Column 3: Hyper-reduction

Row 1 (Multiplications):

- Quadratic R O M: b r cubed plus b r squared

- Hyper-reduction: 8 b r m plus 4 b m

Row 2 (Additions):

- Quadratic R O M: b r cubed minus b r

- Hyper-reduction: 8 b r m minus 4 b m minus b r

Row 3 (Total ops):

- Quadratic R O M: 2 b r cubed plus b r squared minus b r

- Hyper-reduction: 16 b r m minus b r

4.2. Empirical cost scaling and asymptotic behavior on CPU

After establishing the theoretical operation counts, we now evaluate empirical runtimes. The goal is to test how the wall-clock cost of nonlinear evaluations scales with the number of POD modes, and to confirm whether the observed scaling agrees with theoretical predictions: cubic for the quadratic Galerkin ROM and linear for ECSW hyper-reduction.

We run the nonlinear evaluation functions introduced in Section 3 on a single Intel Xeon(R) Gold 5118 core. To obtain a robust estimate of runtime for each configuration, we perform $ 7\times 100 $ function calls. Each block of 100 runs is averaged to reduce measurement noise, and we then take the median of the seven block averages as the final reported value. This procedure combines the precision of averaging with the robustness of a median, limiting the impact of occasional outliers. For consistency with the deep learning code, all computations are batched with batch size $ b=256 $ .

The experiment covers:

• • Quadratic ROMs with POD sizes between 5 and 250 modes (245 data points).

• • Hyper-reduced ROMs over the same POD range, combined with Gauss-point ensemble sizes $ m\in \left\{10,100,200,300,400,500,600,700,800,900,1000\right\} $ , yielding $ 11\times 245=2695 $ configurations.

In total, this corresponds to 171,500 function calls for the quadratic ROM and 1,886,500 for the hyper-reduced models.

Figure 4 plots the results on log–log axes. For the quadratic ROM, the slope of the curve is close to cubic, with an empirical exponent of ≈ 3.25. The slightly super-cubic behavior can be attributed to hierarchical memory effects on the CPU. For hyper-reduction, we illustrate two representative cases ( $ m=100 $ and $ m=1000 $ ). Their asymptotic slopes are ≈ 0.94 and ≈ 1.13, respectively, confirming the predicted linear complexity with respect to the POD size.

Figure 4.

Empirical wall-clock time (log–log scale) for nonlinear evaluation. Comparison between quadratic ROMs and ECSW hyper-reduced models with varying Gauss-point sizes, executed on CPU hardware.

Figure 4. Long description

The x-axis represents the Number of Modes r on a logarithmic scale from 10 super 0 to 10 super 2. The y-axis represents Time per Evaluation in seconds on a logarithmic scale from 10 super minus 4 to 10 super minus 1.

At the top left, a legend lists Quadratic R O M median and E C S W m equals 10 through 1000 median in increments of 100.

Data trends:

* The Quadratic R O M curve is the uppermost line, showing a steep power-law increase. An annotation indicates a slope of Delta y over Delta x equals 3.25.

* Below this, a cluster of E C S W curves shows significantly lower evaluation times. The highest in this group is E C S W m equals 1000, with an annotated slope of Delta y over Delta x equals 1.13.

* As m decreases, the curves shift downward, indicating faster performance.

* The lowest curve, E C S W m equals 10, remains nearly flat across the entire range of r, with an annotated slope of Delta y over Delta x equals 0.94.

* All E C S W models maintain a significant performance advantage over the Quadratic R O M as the number of modes increases beyond 20.

Beyond asymptotics, the point of intersection between the curves is of practical importance. Across all tested Gauss-point sizes, hyper-reduction becomes cheaper than the quadratic ROM at around 50 POD modes. For the case $ m=200 $ , which is most relevant to our applications, hyper-reduction already outperforms the quadratic ROM from about 10 POD modes onward. At 20 modes, the hyper-reduced model requires $ 5\times {10}^{-5} $ s per evaluation, compared to $ 7\times {10}^{-5} $ s for the quadratic ROM. This crossover behavior agrees with our earlier tests.

These experiments confirm that, on CPU hardware, quadratic ROMs scale cubically with POD size, while ECSW hyper-reduction scales linearly. Although hyper-reduction can be more expensive at very small reduced spaces, its advantage grows rapidly with dimension, making it the more efficient choice for realistic POD sizes.

4.3. Empirical cost scaling and asymptotic behavior running on GPU hardware

As a final test, we run the same experiment with both models now running on our NVIDIA Tesla V100S 32GB VRAM GPU chip. This is the hardware we use to run our hybrid model training and inference. The same hardware where we observed slower run times for hyper-reduction compared to the quadratic ROM.

The test follows the same isolated-evaluation setup as on CPU, except that on GPU we report the median over 100 single evaluations together with the inter-quartile range, since the per-call launch jitter is itself of interest and is better shown directly than hidden inside block averages. The Gauss-point ensemble size grid search is also slightly different, with $ m\in \left\{\mathrm{10,20,50,100,500,1000}\right\} $ . This more logarithmic grid-search scale allows us to more clearly see the effects at the extreme ends of the spectrum.

The results for the scaling for the new hardware experiment are very different from the previous ones. Analyzing first the cost of the hyper-reduction models, the cost is nearly constant with respect to both the size of the POD modes and the size of the Gauss-point ensemble, or at least the cost increases only at discrete values.

Concerning the quadratic ROM scaling, initially we see the same near-constant cost scaling with several discrete cost values, but then for models over $ \sim 80 $ POD modes the compute time blows up and follows a new cost scaling that seems exponential or at least asymptotically of high polynomial order of around $ \mathcal{O}\left({r}^{5.5}\right) $ .

These results broadly agree with the timings reported in Section 3.3. In both cases only the nonlinear evaluation itself is timed; what differs is the context surrounding the timed region. During trajectory integration, each right-hand-side evaluation is interleaved with the implicit linear solve and time-stepping bookkeeping, whereas in the present experiment it is called repeatedly in a tight loop. These two settings are not strictly equivalent, so a small difference in absolute cost is to be expected: the $ 18 $ -mode quadratic ROM evaluates in $ 138\hskip0.22em \mu \mathrm{s} $ and the $ 18 $ -mode, $ 200 $ -Gauss-point hyper-reduced model in $ 184\hskip0.32em \mu \mathrm{s} $ , both marginally below the trajectory values of $ 153 $ and $ 208\hskip0.32em \mu \mathrm{s} $ . Crucially, the relative gap between the two methods is preserved: $ 33\% $ here, close to the $ 36\% $ measured during trajectory integration in Section 3.3.

Figure 5 also reports the experimental variability, plotting the median over $ 100 $ evaluations together with the shaded inter-quartile range (IQR), which is visually narrow throughout. The IQR peaks at about $ 28\% $ of the median for the quadratic ROM (at $ r=5 $ ) and at about $ 16\% $ for the hyper-reduced model, confirming that the median is a reliable indicator of practical performance.

Figure 5.

Empirical wall-clock time (log–log scale) showing the computational-cost scaling of quadratic reduction vs hyper-reduction on GPU hardware. Results are reported as medians with inter-quartile ranges.

Figure 5. Long description

A log-log scale line graph titled G P U runtime: quadratic R O M vs E C S W on G P U.

* Axes: The horizontal X-axis represents the Number of modes r, ranging from 10 super 0 to approximately 3 times 10 super 2. The vertical Y-axis represents Time per evaluation in microseconds, ranging from 10 super 2 to 10 super 3.

* Data Series:

- Quadratic R O M (median and I Q R): Represented by a thick blue line. It starts at the lowest point (approx. 40 microseconds) and shows a steep, near-exponential increase as r exceeds 10 super 2, reaching nearly 10 super 3.

- E C S W hyper-reduction (median and I Q R): Multiple lines for m equals 10, 50, 100, 200, 500, and 1000. These lines are clustered horizontally between 200 and 300 microseconds.

* Trends: While the Quadratic R O M is faster for low values of r (below 10 super 2), it scales poorly. The E C S W methods maintain a nearly constant, flat computational cost regardless of the number of modes. The blue Quadratic R O M line intersects and surpasses all E C S W lines at approximately r equals 100, indicating the point where hyper-reduction becomes more efficient on G P U hardware.

To ensure that these observations were not artifacts of a particular implementation, we repeated the profiling experiments with both implementations expressed as torch.nn.Module objects compiled with torch.compile, and also reproduced the benchmark in JAX. In all cases, the relative performance and scaling trends remained unchanged.

Several steps are needed to understand these results: why the hardware scaling appears almost linear at small sizes, and why the cost suddenly explodes at $ r=79 $ for the quadratic ROM. We begin by examining how each method can be parallelized, and then turn to memory scaling and GPU hardware effects.

4.3.1. Parallelization of the nonlinearity evaluation

For the quadratic ROM, the lower compute time can be explained by two factors. First, as shown in the operation counts above, at typical sizes (e.g., $ r=20 $ POD modes and $ m=200 $ Gauss points), the quadratic ROM tensor evaluation requires fewer total operations than hyper-reduction. Second, GPU optimization reduces the effective cost even further.

The computation involves essentially two steps. First,

$$ {T}_{ljk}=\sum \limits_{i=1}^r{N}_{ijk}\hskip0.1em {a}_{li}, $$

followed by

$$ {f}_{lk}=\sum \limits_{j=1}^r{T}_{lj k}\hskip0.1em {a}_{lj}. $$

In both steps, the $ {br}^2 $ (first) and $ br $ (second) repeated operations are non-sequential and can be executed in parallel. In practice, they are likely handled by a batched general matrix–multiplication (GEMM) kernel. Ignoring bandwidth limits, each step can be treated as a near constant-time operation. Thus, the quadratic ROM reduces to two sequential, highly optimized GPU kernels.

The hyper-reduction is more sequential by design. It can be broken down into three main stages: (1) projection into the empirical Gauss-point basis, (2) point-wise products and additions within this space, and (3) projection back into the POD subspace.

For the projection stage, the matrices can be stacked to exploit GPU parallelism. For example, $ {M}_3,{M}_4,{M}_5,{M}_6 $ are grouped into $ {M}_{\mathrm{large}1}\in {\unicode{x211D}}^{4m\times r} $ , and $ {M}_7,{M}_8,{M}_9,{M}_{10} $ into $ {M}_{\mathrm{large}2}\in {\unicode{x211D}}^{4m\times r} $ . In practice, we have 2 Gauss-point ensembles of slightly different sizes, so two separate stacks are required. This allows all projections to be evaluated efficiently using the highly optimized torch.linear routine (perhaps the highest optimized torch GPU function thanks to its role in neural networks), which can again be considered constant-time for small $ m $ and $ r $ .

Once the $ m $ -sized vectors are computed, the element-wise products are trivially parallelizable with compute unified device architecture (CUDA) kernels. The subsequent additions are equally efficient. Together, these operations can be treated as one or perhaps several constant-time extremely fast operations.

Finally, the projection back to POD space is performed by stacking $ {M}_1 $ and $ {M}_2 $ into $ {M}_{\mathrm{large}3}\in {\unicode{x211D}}^{r\times 2m} $ and applying a batched torch.linear. This again corresponds to one constant-time operation.

When memory and bandwidth limits are ignored, the quadratic ROM reduces to two sequential constant-time kernels, while the hyper-reduction requires at least three, and typically more, sequential stages. This explains both the near-constant cost observed for small $ r $ and why hyper-reduction is consistently about 35% slower in this regime: it involves roughly $ 3/2 $ as many sequential steps.

This sequential-stage argument can be complemented by an arithmetic intensity perspective. The quadratic ROM consists entirely of dense GEMM operations; at typical sizes ( $ r=20 $ , $ b=256 $ ), each contraction achieves roughly 9 FLOP/byte, placing it comfortably in the compute-efficient regime of the GPU. The ECSW GEMM stages (forward projections and back-projection) reach a comparable arithmetic intensity of about 8–9 FLOP/byte. However, between these GEMM stages, ECSW requires element-wise Hadamard products over $ b\times m $ vectors, which perform a single floating-point operation per element while reading two operands and writing one result, yielding an arithmetic intensity of approximately 0.08 FLOP/byte—two orders of magnitude lower. Whether these heterogeneous operations are fused by the compiler into a single kernel (as observed with the XLA backend) or dispatched separately, the overall arithmetic intensity of the ECSW evaluation is diluted by the element-wise stages, contributing to its higher wall-clock cost despite comparable GEMM-level performance.

4.3.2. Hardware bandwidth limitations

We now account for hardware limits to explain the near-constant cost at small $ r $ and the sharp increase at $ r\ge 79 $ for the quadratic ROM. The key object is the intermediate tensor $ T\in {\unicode{x211D}}^{b\times r\times r} $ formed in the first contraction. Its size is $ \mid T\mid =4\hskip0.1em b\hskip0.1em {r}^2 $ bytes (float32). While $ N\in {\unicode{x211D}}^{r\times r\times r} $ is streamed slice-by-slice (working set $ \mathcal{O}\left({r}^2\right) $ ) and $ \mathbf{a},\mathbf{f}\in {\unicode{x211D}}^{b\times r} $ are tiny, $ T $ must be repeatedly updated and then consumed. If $ T $ fits entirely in the GPU’s L2 cache, all those updates and reads occur on-chip; dynamic random-access memory (DRAM) traffic is dominated by streaming the tensor $ N $ slice-by-slice from memory, and wall time appears flat. The threshold is set by L2 capacity $ {S}_{L2} $ :

$$ 4\hskip0.1em b\hskip0.1em {r}^2\lesssim {S}_{L2}\hskip1em \Rightarrow \hskip1em r\lesssim \sqrt{\frac{S_{L2}}{4b}}. $$

With $ b=256 $ and $ {S}_{L2}\approx \mathrm{6,291,456} $ B (V100, ~6 MB), this gives $ r\lesssim 78.38 $ , this exactly matches our observed constant cost breaking point. Beyond this point, $ T $ no longer resides in L2. Kernels must tile $ T $ and perform read–modify–write to DRAM for each $ i $ -slice of $ N $ , introducing a dominant bandwidth term $ \propto b\hskip0.1em {r}^3 $ . This is why the cost “cliff” appears precisely around $ r=79 $ : below it, $ T $ reuse is on-chip, and cost is nearly constant; above it, repeated DRAM traffic for $ T $ drives the high-order scaling.

4.4. Guidelines for using hyper-reduction to speed up learning and inference in hybrid deep learning models

The initial motivation for adopting hyper-reduction in hybrid deep learning ROM architectures was its reported computational efficiency. While this remains theoretically true, our results show that scaling arguments alone do not justify its use. In practice, the tight coupling between the ROM and the neural network means both must run on the same hardware. Since deep learning gains enormous speed-ups on GPUs, it is natural to place the ROM there as well. The relevant question is therefore: given a target accuracy, POD size, GPU hardware, and the structure of the nonlinearity, does hyper-reduction actually accelerate the hybrid method?

Before giving practical guidelines, it is useful to discuss how the required quadrature support may evolve for more complex flows. ECSW is an empirical approximation on top of the quadratic POD–Galerkin ROM; in general it cannot be more accurate than the underlying ROM, and coincides with it only when a sufficiently rich quadrature rule is retained. At high Reynolds numbers, two scenarios arise.

In the first, all turbulent scales are resolved through DNS. The Kolmogorov barrier (Ahmed and San, Reference Ahmed and San2020) slows the decay of the POD singular spectrum, forcing the retention of a large number of modes $ r $ . Since Hernández et al. (Reference Hernández, Caicedo and Ferrer2017) show that the number of quadrature points $ m $ scales between $ \mathcal{O}(r) $ and $ \mathcal{O}\left({r}^2\right) $ for greedy non-negative cubature rules (including ECSW), for a given NNLS residual tolerance the number of retained quadrature points grows substantially with turbulence intensity.

In the second, one retains a moderate reduced basis (possibly larger than in the laminar regime studied here, but not drastically so) and closes the unresolved scales through a turbulence model (Lorenzi et al., Reference Lorenzi, Cammi, Luzzi and Rozza2016; Hijazi et al., Reference Hijazi, Stabile, Mola and Rozza2019), a learned closure, or both. Because fewer scales need to be represented, the required number of POD modes $ r $ stays restrained; the same cubature scaling (Hernández et al., Reference Hernández, Caicedo and Ferrer2017) then implies that $ m $ also remains moderate, and ECSW preserves its cost advantage.

Our GPU profiling results support this second scenario. Section 4.3 shows that the hyper-reduced cost remains nearly constant for quadrature sizes up to $ m\approx {10}^3 $ , well above the values needed in our experiments. For moderately richer bases, the dominant operations reduce to batched $ m\times r $ matrix–vector products that better exploit GPU parallelism; the advantage would diminish only when these products become large enough to saturate memory bandwidth or kernel parallelization capacity. These observations remain qualitative projections; dedicated turbulent benchmarks would be needed to confirm them.

From our experiments, we suggest the following:

• • If your nonlinearity is more complex than the quadratic convective term, hyper-reduction is the only practical option.

• • If your application involves a moderately large reduced basis (on the order of a hundred modes, as may arise in industrial turbulent configurations) and your chosen method resolves all retained POD modes through the physics-based operator (Barnett et al., Reference Barnett, Farhat and Maday2023; Menier et al., Reference Menier, Bucci, Yagoubi, Mathelin and Schoenauer2023), then the quadratic ROM computational cost exceeds the efficient threshold. At this scale, the number of quadrature points required remains manageable, making ECSW computable and the preferred option.

• • In contrast, for hybrid ROM–DL approaches where only a small subset of dominant POD modes is evolved through the physics-based operator and the remaining dynamics are captured by a neural correction (Wang et al., Reference Wang, Ripamonti and Hesthaven2020; Menier et al., Reference Menier, Bucci, Yagoubi, Mathelin and Schoenauer2023; Lepage et al., Reference Lepage, Beneddine, Fiorini, Mortazavi, Sipp and Thome2025), the reduced dimension is intentionally kept very small. In this setting, the quadratic nonlinear term can already be evaluated extremely efficiently on GPU. As a result, the hyper-reduced formulation is less efficient than the standard tensorial ROM in this low-dimensional regime.