Source-linked AI summary

Local and Sparse Kernel Operator Learning for Incompressible Flows: An Operator Lagrange Formulation

Ramansh Sharma, Matthew Lowery, Houman Owhadi, Varun Shankar

arXiv:2610.29009v1cs.LG

TL;DR

Many-query incompressible-flow problems require repeated Navier–Stokes solves, while the earlier property-preserving kernel method relies on global output reconstruction and dense operator regression. This paper combines local divergence-free reconstruction with sparse kernel regression and operator Lagrange weights. Across fluid benchmarks, the selected models are usually sparse and compare favorably with earlier methods.

  • Problem

    Learning incompressible-flow solution operators for repeated queries is costly with the earlier method’s global output reconstruction and dense operator regression.

    Source: §1 · pp. 1–2
  • Method

    The method combines local divergence-free output interpolants with compactly supported sparse kernel regression and operator Lagrange weights for new inputs.

    Source: §1 · p. 2, §6 · p. 28
  • Results

    Across operator benchmarks, test errors decrease roughly like N^-τ with τ between 0.4 and 1.3 except in two degenerate cases; selected kernels are sparse on all but one benchmark and beat reported neural operators by factors from 2 to more than 10^6.

    Source: §6 · p. 28
  • Takeaways & Limitations

    Local reconstruction systems remain fixed-size as output points are added, and storage for most benchmark models grows roughly linearly with training-set size over the tested range.

    Source: §6 · p. 28

    Local interpolants are not coupled across stencils, and the rounding analysis does not cover the affine correction or evaluation before reconstruction.

    Source: §6 · p. 29

Abstract

from arXiv · show

We present a local and sparse kernel method for learning the solution operators of incompressible flows. It builds on our earlier property-preserving kernel method, whose global output reconstruction and dense operator regression scale poorly. Our method replaces the global output reconstruction by local divergence-free interpolants with divergence-free polynomials of any degree, solved in parallel on the GPU. These interpolants are divergence free, periodic on periodic domains, and, with a multiscale kernel, approach Kolmogorov statistics under refinement. A compactly supported kernel, positive semidefinite in any input dimension, makes the operator regression sparse. Its support radius, and hence the sparsity, is fitted to data, and optional affine polynomials in the input capture global trends. Operator Lagrange weights join the two. They cost one sparse solve per new input and act on the training fields before the output reconstruction, so their local coefficients are never formed. We show that, unlike direct evaluation, Lagrange evaluation does not amplify solve errors by the cancellation among spatial coefficients. On ten fluid benchmarks, validation selects sparse matrices nearly everywhere, with densities that mostly decrease as the training set grows, and affine augmentation lowers the error by more than an order of magnitude for forced isotropic turbulence. These sparse models are at least as accurate as our earlier dense method and more accurate than neural operators. Lagrange evaluation is as accurate as direct evaluation with local output reconstructions, and up to nearly three orders of magnitude more accurate with global ones. We also present GPU timings and storage.

1. Introduction

Repeated Navier–Stokes solves make many-query flow problems costly, while the earlier property-preserving kernel method uses global output reconstruction and dense operator regression. This paper joins local divergence-free reconstruction with sparse kernel regression through operator Lagrange weights, and reports strong results across ten fluid benchmarks.

1. Introduction

  • Repeated Navier–Stokes solves dominate workflows requiring velocities for many initial conditions, boundary conditions, or forcing functions, while global reconstruction and dense regression limit scale.The dense operator solve costs O(N^3) operations and O(N^2) storage, and both global stages bind for three-dimensional turbulence with 10,000 or more training functions and thousands of output points.Source: §1 · pp. 1–2
  • Local divergence-free interpolants replace the global output reconstruction, preserving divergence-free predictions while making independent, fixed-size stencil systems suitable for parallel GPU solution.The construction uses divergence-free polynomials of any degree; periodic images give periodic reconstruction, and a multiscale Matérn kernel approaches prescribed Kolmogorov statistics under refinement.Source: §1 · p. 2, §1.1 · p. 4
  • A compactly supported kernel makes operator regression sparse and positive semidefinite in any input dimension; validation fits its support radius, while optional affine polynomials capture input trends.The support radius is selected from validation fields, allowing nearly dense candidates when accuracy requires them; affine augmentation exactly reproduces the included input polynomials.Source: §1 · p. 2
  • Operator Lagrange weights connect regression to reconstruction with one sparse solve per new input and can be applied to training fields before reconstruction, avoiding formation of their local coefficients.The local coefficients would require tens of times the storage of the fields themselves, and Lagrange evaluation avoids the solve-error amplification caused by cancellation among spatial coefficients.Source: §1 · pp. 2–3
  • At earlier-study training sizes, the models are at least as accurate as dense PPKM and the kernel method of Batlle et al. [2] on every data set, and more accurate than the neural operators reported there.Validation selects sparse matrices on nearly all benchmarks; Lagrange evaluation matches direct evaluation with local reconstruction and is up to nearly three orders of magnitude more accurate with global reconstruction.Source: §1 · p. 3

1.1. Problem formulation

  • The learned operator maps sampled initial, boundary, or forcing data to velocity through coefficient prediction and output reconstruction.The output may be spatial or spacetime, and the reconstruction is a collection of local divergence-free interpolants rather than one global interpolant.Source: §1.1 · pp. 3–4

2. Local property-preserving output reconstruction

The output reconstruction replaces one global interpolant with independent local systems combining divergence-free kernels and polynomials. Its stencil construction supports periodic and spacetime outputs, while multiscale Matérn kernels extend the reconstruction to turbulent fields.

2. Local property-preserving output reconstruction

  • Each output point uses a local interpolant combining a divergence-free kernel and divergence-free polynomials, with fixed-size systems independent of output resolution.The local systems are independent and can be assembled, factored, and evaluated in batches on the GPU.Source: §2 · p. 4, §2.1.3 · pp. 6–7

2.1.1. The divergence-free kernel

  • Applying −ΔI+∇∇ᵀ to a scalar kernel yields a matrix-valued kernel with divergence-free columns; the cubic kernel is continuous but not differentiable at its center.The columns have zero divergence distributionally and classical zero divergence away from their center.Source: §2.1.1 · p. 5

2.1.2. Divergence-free polynomials of any degree

  • Differentiated Legendre-product potentials span divergence-free vector polynomials of any degree, and pivoted QR retains an independent subset for each stencil.If stencil geometry forces rank below full dimension, the interpolant reproduces only the retained subspace.Source: §2.1.2 · pp. 5–6

2.1.3. The local system

  • Independent fixed-size stencil systems support parallel GPU solves, but uncoupled local interpolants can jump where stencil assignment changes.Divergence is controlled within each stencil rather than across stencil boundaries.Source: §2.1.3 · pp. 6–7

2.2. Multiscale Matérn output reconstruction for turbulence

  • A weighted mixture of Matérn kernels at dyadic lengths builds prescribed Kolmogorov increment scaling into each stencil’s reconstruction.The finite-stencil reconstruction’s recovery of those statistics also depends on stencil extent and node spacing; the mixture is smooth at its centers.Source: §2.2 · p. 7

2.3. Periodic stencils

  • Periodic stencils use nearest translated images of nodes, shifting all periodic coordinates together so face and corner neighbors are handled consistently.Each image carries its original node’s sampled velocity, and translation preserves divergence.Source: §2.3 · pp. 7–8, Figure 2 · p. 8

2.4. Spacetime outputs

  • For spacetime outputs, snapshots reuse the same spatial stencil system, while temporal Lagrange interpolation adds no further linear system.The spatial basis remains fixed across sampled times; the spacetime reconstruction is divergence free wherever the spatial reconstruction is divergence free, and each added time contributes only another right-hand side.Source: §2.4 · pp. 8–9

3. Sparse operator regression

The section makes operator regression sparse with a compactly supported kernel that is positive semidefinite in any input dimension, then augments it with affine input trends. It combines this regression with local divergence-free output fields and chooses kernel parameters using training data and validation error.

3. Sparse operator regression

  • A compactly supported kernel replaces the dense operator-regression matrix with a sparse one.The earlier global-support kernel required O(N^3) factorization operations and O(N^2) storage.Source: §3 · p. 9

3.1. A compactly supported kernel in any input dimension

  • A low-dimensional projection sets compact support while a Matérn factor handles residual coordinates, yielding a product kernel positive semidefinite in any input dimension.The compactly supported radial factor supplies sparsity without requiring a compactly supported radial kernel to be positive definite in the full input dimension.Source: §3.1 · p. 9

3.2. The operator interpolant

  • Kernel regression predicts all local output coefficients from the same input, and the resulting operator combines fixed divergence-free output fields.The reconstructed output remains divergence free on each stencil regardless of operator-regression accuracy.Source: §3.2 · pp. 10–11

3.3. Affine augmentation

  • Affine functions of up to 15 principal input coordinates are reproduced exactly, while the kernel represents the remainder.This augmentation addresses the poor representation of trends by a compactly supported kernel that returns to zero away from training inputs.Source: §3.3 · p. 11

3.4. Linear systems

  • The systems decouple by stencil, and matrix density is determined by the support radius and data because every nonzero kernel entry is retained.The density is the number of nonzero entries divided by N^2.Source: §3.4 · p. 12

3.5. Sparse Cholesky and whitened QR

  • Sparse Cholesky followed by QR of whitened polynomial columns avoids forming the Schur complement, while nested dissection supports parallel GPU solves.For rank-one inputs, band ordering can reduce memory but leaves prediction solves with less GPU parallelism.Source: §3.5 · p. 12

3.6. Choosing the kernel parameters

  • The support radius trades accuracy against sparsity: fitted shapes and regularization use restricted likelihood, then validation error selects the radius.Candidate radii are first evaluated at densities from 0.10 to 0.90; the restricted-likelihood criterion is used for parameter selection, not uncertainty quantification.Source: §3.6 · pp. 12–13

4. Operator Lagrange prediction

Operator Lagrange prediction replaces output-coefficient solves with one solve for input-dependent weights, applied to the sampled training fields before reconstruction. The section also analyzes rounding-error differences and details the GPU implementation and work-storage trade-offs.

4. Operator Lagrange prediction

  • One solve per new input yields weights that combine the sampled training fields before the linear output reconstruction.This avoids solving separately for every spatial coefficient and forming local coefficients for every training field.Source: §4 · p. 13, §4.1 · p. 14, §4.2 · pp. 14–15

4.1. Operator Lagrange weights

  • The weights reproduce affine functions exactly when affine augmentation is used, and shared linear weights preserve divergence-free constraints on each stencil.The stated reproduction result assumes a positive-definite kernel matrix and full-column-rank input polynomial matrix.Source: §4.1 · p. 14

4.2. Applying the weights before the output reconstruction

  • Applying weights before reconstruction lets the GPU combine sampled fields and apply the reconstruction once per new input, without storing each training field’s local coefficients.The implementation prepares compact evaluation weights on each stencil using adjoint solves with the local matrices.Source: §4.2 · pp. 14–15

4.3. Rounding errors of direct and Lagrange evaluation

  • Lagrange evaluation avoids the spatial-coefficient cancellation factor that amplifies direct solve errors, but the analysis does not establish greater accuracy in general.The analysis excludes affine correction and the order of operations used when weights are applied before output reconstruction.Source: §4.3 · pp. 15–16

4.4. GPU implementation

  • The GPU implementation keeps the factors on the GPU and batches new-input right-hand sides; changing kernel support requires a new ordering and symbolic analysis.Changing kernel values alone requires only numerical refactorization.Source: §4.4 · p. 16

4.5. Work and storage

  • Local reconstruction setup grows linearly with output-point count, while the operator-regression factor’s size depends on kernel fill as well as matrix density.Lagrange weights are dense even when the kernel matrix is sparse; direct evaluation can reuse its solved coefficients, so the cheaper route depends on query count.Source: §4.5 · p. 17, Table 1 · p. 17

5. Numerical experiments

The numerical experiments assess local output reconstruction, compare sparse operator models across ten fluid benchmarks, and examine their storage, cost, and prediction accuracy. They also compare direct and Lagrange evaluation on the same unaugmented models, including with local and global output reconstruction.

5.1.1. Convergence on scattered nodes

  • Scattered-node reconstruction approaches its target convergence orders, while divergence and periodic mismatch remain at rounding level across tested kernels, degrees, and node sets.In 2D, fitted convergence rates are 2.1, 3.9, and 5.8; in 3D they are 1.7–1.9, 3.4–3.7, and 5.0, with sixth-order convergence rising to 5.5 on the three finest sets.Source: §5.1.1 · p. 18, Figure 3 · p. 18, §5.1.2 · p. 18

5.1.3. Kolmogorov statistics

  • The multiscale reconstruction approaches the reference Kolmogorov structure function under refinement, but retains large errors at separations below node spacing.At separation r = 0.8, fitted error-decay exponents are 0.64, 0.73, and 0.67 for degrees 1, 3, and 5; local and global interpolation have nearly the same error, pointing to node resolution rather than stencil locality.Source: §5.1.3 · p. 19, Figure 4 · p. 19

5.1.4. Output reconstruction in the fluid benchmarks

  • Local output reconstruction is not uniformly more accurate than global reconstruction between nodes, but its setup becomes faster at larger node counts and uses less memory at every tested size.Global reconstruction wins on four flows by factors of 1.06–2.7, while local wins on six by factors of 1.04–10; local systems have at most 62 unknowns in 2D and 262 in 3D, versus up to 21000 globally. Setup becomes faster near 3 × 10^3 points in 2D and 5 × 10^3 in 3D.Source: §5.1.4 · p. 19, Table 2 · p. 20, §5.1.5 · p. 20

5.2.2. Sparsity and convergence in 𝑁

  • Validation selects sparse matrices for most benchmarks, their density usually decreases as N grows, and test errors generally decrease with training-set size.At N = 10000, selected densities include 10% for the lid-driven cavity and both cylinder flows; forced isotropic turbulence remains at 90%. Error-decay exceptions include augmented backward-facing step and unaugmented Taylor–Green models.Source: §5.2.2 · p. 22, Figure 6 · p. 22, Figure 7 · p. 23

5.2.3. Kernel shape and affine augmentation

  • Affine augmentation and fitted kernel shape help selectively: augmentation cuts forced-turbulence error by 29-fold, or 39-fold with fitted shape, while fitted shape cuts species-transport error sixfold.For five rank-one flows other than Taylor–Green, model error ratios remain between 0.98 and 1.36; merging-vortices errors are within 12%, so neither option is uniformly better and both are selected by validation.Source: §5.2.3 · p. 23

5.2.4. Comparison with earlier results

  • The sparse local models are at least as accurate as dense PPKM and VKM on all ten benchmarks, and beat the better neural operator by factors from 2 to more than 106.The reported comparisons use the same data sets, but earlier errors average 200 held-out pairs while ours average 100 test pairs; the neural operators were trained in single precision.Source: §5.2.4 · pp. 23, 25, Table 4 · p. 24

5.2.5. Storage

  • Model storage is dominated by training fields over the tested range, while operator Lagrange weights avoid storing local coefficients for each training field.Most model-storage fits have exponents 0.9–1.1, while the Cholesky factor grows as N^1.4 to N^2.4. For degree-four reconstruction, local coefficients would require 31 times training-field storage in 2D and 87 times in 3D.Source: §5.2.5 · p. 25, §5.3.1 · p. 26

5.3. Operator Lagrange weights

  • Lagrange prediction requires one sparse solve per new input and avoids solve-error amplification from cancellation among spatial coefficients.For forced isotropic turbulence at N = 10000, storing training fields takes 1.7 GB versus 147 GB for local coefficients. Global-reconstruction amplification factors exceed 10^7 on all benchmarks except Taylor–Green spacetime, reaching 6 × 10^9; local factors range from 6 to 150 in 2D and are 396 for forced turbulence.Source: §5.3 · p. 26, §5.3.1 · p. 26, §5.3.3 · p. 27

5.3.3. Accuracy of direct and Lagrange evaluation

  • A large spatial-coefficient norm allows but does not guarantee a Lagrange gain, which also requires a non-negligible solve residual.For forced isotropic turbulence, α_spatial ≈ 4×10^9 with global reconstruction, but the gain is at most 1.2.Source: §5.3.3 · p. 27
  • For a global lid-driven-cavity model, the gain grows from 1.3 to 5.1 as N grows, and solving directly for sampled velocities removes it, implicating large spatial coefficients.The comparison follows one global model as N grows; replacing the coefficient solve with a solve for sampled velocities removes the gain.Source: §5.3.3 · p. 27
  • Lagrange evaluation retains cost and storage benefits for local output reconstruction and can protect global or otherwise ill-conditioned reconstruction from amplified solve errors.The authors note this protection also applies to the global output reconstruction of their earlier method.Source: §5.3.3 · p. 27
  • With global output reconstruction, Lagrange evaluation improves accuracy by up to 868× for merging vortices, with median gains above 5 for three benchmarks; with local reconstruction, evaluations differ by at most 3×10^-6.Across all 560 model pairs, the smallest direct-to-Lagrange test-error ratio is 0.99987, so Lagrange evaluation is never less accurate by more than 0.02%.Source: §5.3.3 · p. 27, Table 6 · p. 28

6. Conclusions

The paper combines local divergence-free output reconstruction with sparse operator regression to reduce the scaling costs of property-preserving kernel learning. Across the operator benchmarks, sparse models remain competitive with dense PPKM and outperform the neural-operator results reported in the earlier study.

Output reconstruction

  • The local reconstruction is about as accurate between nodes as the global reconstruction, uses systems about two orders of magnitude smaller, and approaches prescribed Kolmogorov statistics at resolved separations.It is more accurate on six of ten benchmarks and less accurate on the other four; divergence and periodic mismatch remain at rounding level.Source: §6 · p. 28

Benchmark results

  • Test errors decrease roughly like N^-τ, with τ between 0.4 and 1.3 in all but two degenerate cases; the most accurate kernels are sparse on all but one benchmark.Their densities mostly decrease as the training set grows; at the earlier study’s training sizes, sparse models are at least as accurate as dense PPKM and more accurate than the reported neural operators by factors from 2 to more than 10^6.Source: §6 · p. 28

Lagrange evaluation

  • Lagrange evaluation is as accurate as direct evaluation with local reconstructions and up to nearly three orders of magnitude more accurate with global reconstructions.It avoids storing local coefficients of the training fields.Source: §6 · p. 28

Limitations

  • The conclusions are bounded by ten benchmarks; richer input families and comparisons on identical test sets and hardware remain future tests.The local interpolants are divergence free on each stencil but are not coupled across stencils; fixed-density kernel-matrix storage still grows quadratically with training-set size.Source: §6 · pp. 28–29

Credit authorship contribution statement

The appendix records author contributions, vector-coefficient ordering, and implementation choices for constructing local polynomial bases and solving stencil systems. The local basis procedure uses shared preprocessing and batched factorizations, with a fallback for non-finite local solves.

Credit authorship contribution statement

  • Sharma’s listed roles include conceptualization, methodology, software, investigation, validation, visualization, and writing; Lowery’s include methodology and software, while Owhadi’s include conceptualization, methodology, and supervision.Shankar’s listed roles include conceptualization, funding acquisition, methodology, project administration, resources, software, supervision, validation, visualization, and writing.Source: Credit authorship contribution statement · p. 29

A.1. Ordering of vector coefficients

  • The coefficient vectors are ordered by component, while the kernel evaluation groups components by node; a permutation maps between these orderings.The retained polynomial ranks are reported for benchmark stencils at degree four, with full ranks of 20 in 2D and 85 in 3D.Source: §A.1 · pp. 29–30, Table 7 · p. 30

A.2. Polynomial basis and local solves on the GPU

  • A shared SVD removes candidate dependencies before sampling, and batched QR factorizations determine stencil ranks; a pseudoinverse is reserved for non-finite local LU solutions.That fallback preserves divergence freedom but may not interpolate or reproduce the full polynomial space; it was not needed on the ten operator benchmarks or in the manufactured tests.Source: §A.2 · p. 30

A.3. Evaluation away from the nodes

Evaluation away from the data nodes uses a stencil centered at each query point for operator benchmarks, but the manufactured tests use the nearest node’s stencil. Periodic stencils are built from periodic images of nodes.

A.3. Evaluation away from the nodes

  • For operator benchmarks, the GPU implementation centers a newly built stencil basis at each off-node query; manufactured tests instead use the nearest node’s stencil.The Table 2 comparisons use the query-centered rule.Source: §A.3 · pp. 30–31

A.4. Periodic stencils

  • Taylor–Green flows and merging vortices have period 2π in each direction; forced isotropic turbulence has period 2π on the explicitly supplied [0, 2π]^3 box because grid cell centers do not reach its faces.Source: §A.4 · p. 31

A.5. Output reconstruction between the nodes

The between-node comparison holds the operator regression fixed and changes only the output reconstruction. It evaluates reconstruction error and binned structure functions, while the appendix also specifies input normalization and principal-direction selection.

A.5. Output reconstruction between the nodes

  • Validation selects scale factor cσ=16 for both turbulent species transport and forced isotropic turbulence among the tested multiscale-kernel choices.For species transport, error between nodes remains much larger than error at nodes, so nodal accuracy alone does not guarantee accuracy everywhere.Source: §A.5 · p. 31
  • At N=10000, binned structure-function discrepancy is 2.6% for the multiscale kernel versus 1.7% for PHS in species transport, and 0.20% versus 0.40% in forced isotropic turbulence.The multiscale kernel lowers pointwise error between nodes for both flows, but improves binned increments only for forced isotropic turbulence; the references are RANS and coarse LES simulations.Source: §A.5 · p. 31

B.1. Normalization and principal directions

  • At each training size, normalization and principal directions are computed from the first N training inputs; coordinates are divided by their standard deviations, or set to zero when those vanish.The numerical-rank threshold is max{10^-10, 2u_mach max(N, n_da)} relative to the largest singular value.Source: §B.1 · p. 31

B.2. Fitting and validation

The appendix specifies how kernel shape and validation choices are made, then documents computational accounting and prediction checks. It reports bounded-memory assembly but notes that the neighbor search remains quadratic in training-set size.

B.2. Fitting and validation

  • Fitted-shape kernels use 60 Adam steps for each of two polynomial settings, with bounded radius, regularization, anisotropy, and residual-scale parameters.The residual scale is inactive for seven one-dimensional benchmarks and remains one for turbulent species transport.Source: §B.2 · p. 31
  • Validation selects the lowest-error model with density at most 0.90, breaking ties by density, after comparing five density budgets and additional radius searches.The additional search uses four starts and 35 Adam steps on 100 validation fields; fitted-shape models repeat it for both shapes.Source: §B.2 · p. 32

B.3. Work, storage, and timing

  • Sparse matrix assembly bounds memory by processing input pairs in blocks, but its neighbor search still requires quadratic work in the training-set size.Compact support alone does not provide a faster neighbor search.Source: §B.3 · p. 32
  • Reported storage sums stored model and sparse-matrix arrays, estimated solver factors, and solve buffers, but excludes compiled code and temporary workspace.A CSR kernel with N_K nonzeros uses 12N_K + 4(N + 1) bytes.Source: §B.3 · p. 32
  • Warm prediction time is the median of 11 repetitions over 100 test fields after setup and compilation; small between-model timing differences are not significant.Output-reconstruction setup is reported separately from model fitting and validation.Source: §B.3 · p. 32

B.4. Batches and checks

  • All reported models pass residual, affine-reproduction, and batched-versus-compiled prediction checks after factorization.The 100 test inputs are solved in four batches of 32 right-hand sides, with the last batch padded by repeated inputs.Source: §B.4 · p. 32

C. Proofs

The proofs establish completeness of divergence-free polynomial candidates, kernel positive definiteness, and the variational and algebraic properties of operator regression. The appendix also describes ten flow benchmarks spanning varied geometries, inputs, and outputs.

C.1. Completeness of the polynomial candidates

  • In three dimensions, homogeneous divergence-free polynomial components can be represented as curls of polynomial potentials, while rank–nullity gives the space dimension.In two dimensions, rotated gradients of scalar basis functions span divergence-free vector polynomials.Source: §C.1 · pp. 32–33

C.2. Positive definiteness of the kernel

  • The product kernel is positive semidefinite by the Schur product argument, and regularization makes its Gram matrix positive definite.The proof does not require inputs to lie in a low-dimensional affine subspace; regularization also handles standardized inputs that coincide.Source: §C.2 · p. 33

C.3. Variational form of the operator regression

  • The representer argument shows that a regression minimizer lies in the span of data-centered kernel functions.The orthogonal component vanishes at training inputs and only increases the norm.Source: §C.3 · p. 33

C.4. Block elimination and reproduction

  • Block elimination yields a unique Lagrange-weight solution, and direct and Lagrange evaluations agree in exact arithmetic.For zero regularization the weights are cardinal; with positive regularization they are not.Source: §C.4 · pp. 33–34

C.5. Rounding errors of direct and Lagrange evaluation

  • The rounding-error analysis bounds matrix-product errors for both direct and Lagrange evaluation, including products computed in component blocks.The stated bounds use standard matrix-product estimates and Frobenius-norm inequalities.Source: §C.5 · p. 34

D.1. Benchmarks

  • The ten benchmark data sets come from prior simulations and cover flows with varied geometries and boundary conditions, including cylinder flow, Taylor–Green vortices, cavities, and a backward-facing step.The cylinder data vary between vortex-shedding and below-onset regimes; Taylor–Green data include spatial and spacetime mappings.Source: §D.1 · p. 34, §D.1.1 · p. 34, §D.1.2 · p. 34
  • Other benchmarks include merging vortices, turbulent species transport in a Kenics mixer, and forced isotropic turbulence, with inputs ranging from vorticity and inlet velocity to forcing.The forced-turbulence simulations use a Smagorinsky large-eddy model at Reynolds number 6 × 10^6.Source: §D.1.6 · p. 35, §D.1.7 · p. 36, §D.1.8 · pp. 36–38

D.2.1. Domains and fields

The appendix defines divergence-free test fields and node sets for convergence tests, then checks local reconstruction for divergence and periodicity. The reported tests show small divergence and cross-period differences in velocities and Jacobians.

D.2.1. Domains and fields

  • The convergence tests use two divergence-free fields that vanish on the boundary, evaluated on progressively finer node sets in two and three dimensions.The finest sixth-order errors are about 4.5 × 10^-6 in 2D and 1.6 × 10^-3 in 3D for both kernels.Source: §D.2.1 · pp. 38–39

D.2.2. Divergence and periodicity

  • The largest relative divergence across smooth tests is 1.5 × 10^-13 where the output reconstruction is differentiable.Derivatives are taken on each local stencil; the tests exclude cubic-kernel centers.Source: §D.2.2 · p. 39
  • Across periodic copies, the largest differences are 1.4 × 10^-12 for velocities and 1.6 × 10^-12 for Jacobians.Comparisons include box faces, edges, corners, and 128 interior points.Source: §D.2.2 · p. 39

D.2.3. The turbulent ensemble

The turbulent ensemble is constructed to have divergence-free covariance and approximately Kolmogorov structure statistics over a fixed separation range. Degree-five reconstruction errors decrease with refinement at the largest tested separation, while substantial worst-range error remains.

D.2.3. The turbulent ensemble

  • The ensemble covariance is divergence free and its structure function approximately follows the Kolmogorov law for separations 0.03 ≤ r ≤ 0.8.At smaller separations, finite bandwidth makes increments smooth.Source: §D.2.3 · pp. 39–40
  • For degree five, E_S2(0.8) falls from 14% on the coarsest node set to 4.6% on the finest, while the worst error across the range remains 64%.The ensemble has 1514 Fourier modes and its structure functions are computed exactly from the covariance.Source: §D.2.3 · p. 40
  • With the same five-scale Matérn kernel and polynomials on 1770 nodes, degree-five local and global reconstructions have E_S2(0.8) errors of 9.9% and 10.6%, respectively.Their normalized mean square increment errors are 0.168 and 0.166, and their worst-range errors are about 87% for both.Source: §D.2.3 · p. 40

D.3. Field figures

  • The field figures show the first simulation of each data set, with panels of the same quantity sharing a color scale.The transverse direction of the backward-facing step is expanded.Source: §D.3 · p. 40

D.4. Direct and Lagrange evaluation

The comparison tests direct and Lagrange evaluation under local and global output reconstruction. With global reconstruction, large spatial coefficients accompany a growing advantage for Lagrange evaluation; direct evaluation on sampled velocities removes that gap.

D.4. Direct and Lagrange evaluation

  • The evaluation study compares direct and Lagrange predictions for local and global output reconstruction on 100 test inputs, using the same normalization, factorization, support, regularization, and spatial basis.For local reconstruction, direct evaluation solves for stencil coefficients in blocks, while Lagrange evaluation applies weights to those coefficients.Source: §D.4 · p. 40
  • With global output reconstruction, the direct-to-Lagrange error ratio rises from 1.3 to 5.1 as the spatial-coefficient Frobenius norm grows from 1.6 × 10^9 to 7.1 × 10^9.Direct evaluation on sampled velocities is as accurate as Lagrange evaluation, locating the gain in output reconstruction through large coefficients.Source: §D.4 · p. 40, Table 8 · p. 40

D.5. Software

  • The paper’s implementation and result-generation scripts are included in KernelPack-JAX, and released-code predictions agree with original fits to 1.0 × 10^-10.The scripts also produce the reported results and build the tables and figures.Source: §D.5 · pp. 40–41

D.6. All operator benchmarks

The appendix assembles results across ten operator benchmarks, including errors, storage, prediction times, and error-decay rates. It also reports offline and initialization costs, while showing that density alone does not capture storage requirements.

D.6. All operator benchmarks

  • At 9.1% and 21.7% density, the augmented species models have 26 and 35 million factor entries; the backward-facing-step factor has 38 million entries at 33% density.These examples illustrate why the benchmark results report storage separately from density.Source: §D.6 · p. 41
  • Figures 20–23 complete the benchmark comparisons of test error, selected density, accounted storage, and warm prediction time for the other six benchmarks.Table 9 reports least-squares error-decay exponents across training sizes, while Tables 10 and 11 collect errors, storage, and prediction times at N = 10000.Source: §D.6 · p. 41, Table 9 · p. 41, Figure 20 · p. 42

D.6.1. Offline cost

  • At N = 10000, fitting and validating fitted-shape models takes 11–27 minutes, versus 2–7 minutes for fixed-shape models.Rebuilding the selected factorization takes about one to five seconds for most benchmarks and about 43 seconds for dense forced-turbulence models; output reconstruction setup takes 0.03–2.9 seconds after compilation.Source: §D.6.1 · pp. 41, 44
  • With affine augmentation, Taylor–Green models select about 10% density at the two largest training sizes, versus 75% or 90% without augmentation.The selected density is shown for the other six benchmarks.Source: Figure 21 · p. 42
  • Table 22 reports offline selection and rebuild costs alongside initialization costs at N = 10000; initialization includes two synchronized 100-field prediction passes, compilation, and CUDA graph preparation.The four models share output-reconstruction setup and may reuse compiled code.Source: Table 22 · p. 50
Loading 2610.29009v1…