PaperlayerView exact revision on arXiv

Source-linked AI summary

Weak SINDy For Partial Differential Equations

Daniel A. Messenger, David M. Bortz

arXiv:2007.02848v3math.NAcs.LGstat.ML

TL;DR

PDE discovery methods commonly rely on pointwise derivative approximations, limiting robustness to noisy data. This paper extends WSINDy to PDEs with a weak, convolutional, FFT-based approach, reporting robust recovery across nonlinear waves and spatiotemporal chaos while identifying weaker performance for limit-cycle data.

  • Problem

    Most PDE discovery methods require pointwise derivative approximations, creating a central challenge for robust identification from data.

  • Method

    WSINDy extends weak-form sparse regression to PDEs using integration by parts, convolutional discretization, FFT acceleration, spectral test-function selection, and learned sparsity thresholds.

  • Results

    WSINDy is very robust for nonlinear waves and spatiotemporal chaos, while identification is less effective for reaction-diffusion and Navier–Stokes data with limit cycles.

  • Takeaways & Limitations

    Weak identification provides a practical route to discovering spatiotemporal PDE dynamics without pointwise derivative approximations across several challenging physical settings.

  • Takeaways & Limitations

    Identification is less robust for data with limit cycles, and the accuracy of recovered coefficients remains incompletely understood for deriving recovery guarantees.

Abstract

from arXiv · show

Sparse Identification of Nonlinear Dynamics (SINDy) is a method of system discovery that has been shown to successfully recover governing dynamical systems from data (Brunton et al., PNAS, '16; Rudy et al., Sci. Adv. '17). Recently, several groups have independently discovered that the weak formulation provides orders of magnitude better robustness to noise. Here we extend our Weak SINDy (WSINDy) framework introduced in (arXiv:2005.04339) to the setting of partial differential equations (PDEs). The elimination of pointwise derivative approximations via the weak form enables effective machine-precision recovery of model coefficients from noise-free data (i.e. below the tolerance of the simulation scheme) as well as robust identification of PDEs in the large noise regime (with signal-to-noise ratio approaching one in many well-known cases). This is accomplished by discretizing a convolutional weak form of the PDE and exploiting separability of test functions for efficient model identification using the Fast Fourier Transform. The resulting WSINDy algorithm for PDEs has a worst-case computational complexity of $\mathcal{O}(N^{D+1}\log(N))$ for datasets with $N$ points in each of $D+1$ dimensions (i.e. $\mathcal{O}(\log(N))$ operations per datapoint). Furthermore, our Fourier-based implementation reveals a connection between robustness to noise and the spectra of test functions, which we utilize in an \textit{a priori} selection algorithm for test functions. Finally, we introduce a learning algorithm for the threshold in sequential-thresholding least-squares (STLS) that enables model identification from large libraries, and we utilize scale-invariance at the continuum level to identify PDEs from poorly-scaled datasets. We demonstrate WSINDy's robustness, speed and accuracy on several challenging PDEs.

1. Introduction

The paper extends weak-form SINDy from ODEs to PDEs, avoiding pointwise derivatives through convolutional weak formulations. WSINDy combines FFT acceleration, spectrally informed test-function selection, automatic threshold learning, and scale-invariant identification for robust PDE discovery.

  • Contribution: WSINDy extends the Weak SINDy framework from recovering ODEs to identifying PDEs from data.It is a Galerkin-based sparse regression method using the weak form of the dynamics.
  • Motivation: Pointwise derivative approximations are a major disadvantage of most PDE discovery methods because they reduce robustness to noise.WSINDy instead integrates by parts in space and time against smooth compactly supported test functions.
  • Results: WSINDy recovers PDEs from highly noisy data, including truly weak solutions, while maintaining high-accuracy recovery in the low-noise regime.Its robustness remains effective as the signal-to-noise ratio approaches one.
  • Computational method: The convolutional weak form enables efficient model identification with the Fast Fourier Transform at worst-case cost O(N^(D+1) log(N)), or O(log(N)) operations per datapoint.Subsampling can further reduce computational cost.
  • Algorithmic extensions: Test-function spectral decay can damp noise-dominated modes, motivating an automatic learning algorithm for test-function hyperparameters.The framework also introduces automatic threshold selection for large libraries and uses PDE scale invariance for poorly scaled data.
  • Experiments: Numerical experiments evaluate WSINDy on nonlinear PDEs spanning waves, spatiotemporal chaos, reaction-diffusion, and Navier-Stokes.The article reports several improvements over existing results in the literature.

2. Problem Statement and Notation

The paper formulates PDE discovery as identifying nonlinear functions and differential operators governing a weak solution from a noisy spatiotemporal dataset. It specifies the computational grid, data model, candidate PDE form, and example systems used for evaluation.

  • Discretization: The spatial grid is a tensor product of equally spaced one-dimensional grids, while time is equally spaced with step Δt.The resulting dataset is a (D + 1)-dimensional array with dimensions N_1 × ··· × N_D+1.
  • Data model: The assumed data model is U = u(X, t) + ε, where u is a weak PDE solution and ε is independent identically distributed noise.The noise array has the same dimensions as U.
  • PDE representation: The candidate PDE form expresses D^α_0u as a sum of differential operators applied to nonlinear functions of u.The multi-index collection enumerates the possible true differential operators.
  • Problem formulation: The objective is to identify functions g_s and differential operators D^α_s governing the evolution of u from dataset U and grid (X, t).The evolution term D^α_0u is treated as the reference side of the PDE.
  • Examples: The method is demonstrated on inviscid Burgers, Korteweg–de Vries, Kuramoto–Sivashinsky, nonlinear Schrödinger, Sine–Gordon, reaction-diffusion, and Navier–Stokes systems.The admissible class also includes PDEs such as Allen–Cahn, Cahn–Hilliard, and Boussinesq equations.

3. Weak Formulation

The PDE is converted into a weak, finite-dimensional coefficient-recovery problem by testing against compactly supported functions and integrating by parts. Discretization yields a convolutional linear system whose dimensions depend on the test-function count and library size, while separable FFT evaluation reduces computational cost.

  • Weak formulation: Recovering the PDE becomes a finite-dimensional problem of finding a sparse coefficient vector w⋆ over candidate differential operators and trial functions.The trial functions are known beforehand, and the governing operators are enumerated by a multi-index set.
  • Weak formulation: The PDE is multiplied by a smooth, compactly supported test function and integrated over spacetime before repeated integration by parts.Compact support makes the integration-by-parts transformation valid.
  • Weak formulation: Discretizing the weak-form integrals with observed data produces the linear system b = Gw⋆, where b and G are the left-hand side and Gram matrix.The trapezoidal rule is used and has nearly negligible quadrature error for the employed test functions.
  • Convolutional weak form and discretization: The discrete weak form is expressed as convolutions between derivative-weighted test functions Ψs and library evaluations fj(U) at query points.Translated test functions generate the convolutional form, and the resulting system can be assembled efficiently with the FFT.
  • Convolutional weak form and discretization: The reference grid is controlled by support integers m, whose values significantly affect the algorithm and can be selected automatically from spectral properties of U.The automatic selection uses the dataset’s spectral information.
  • FFT-based implementation and complexity for separable ψ: For separable test functions, convolutions are performed sequentially across coordinates, avoiding direct construction of multidimensional test-function arrays.The 1D components can be precomputed, and individual convolutions can be parallelized.
  • FFT-based implementation and complexity for separable ψ: For n = N/4, FFT-based separability changes the cost from O(N^(2D+2)) to O(N^(D+1) log(N)), a reduction by N^(D+1)/log(N).This comparison applies to the typical choice n = N/4.
  • FFT-based implementation and complexity for separable ψ: When N = 512 and D + 1 = 3, the separable FFT-based convolution is 104 times faster for 100 ≤ n ≤ 450.Figure 1 compares the ratio of naive to separable FFT-based convolution costs in two- and three-dimensional space-time settings.

4. WSINDy Algorithm for PDEs and Hyperparameter Selection

The PDE WSINDy algorithm combines test-function selection, sparsification, regularization, query-point and library choices, and scaling procedures. Its Fourier interpretation links noise robustness to test-function spectra, while MSTLS learns thresholds for large libraries.

  • Algorithm components: WSINDy requires choices for the reference test function, sparsification, regularization, convolution query points, and model library.The algorithm has many hyperparameters, motivating simplifications and automatic selection procedures.
  • Test-function selection: Fourier properties of the test function determine noise robustness and inform automatic selection of its hyperparameters.The convolutional weak form can act as an approximate low-pass filter without altering the data’s frequency content.
  • Sparsification: MSTLS modifies sequential-thresholding least-squares by automatically selecting the threshold λ, enabling PDE discovery from large libraries.Its loss balances projection distance against the number of nonzero terms, including terms outside the dominant data balance and terms with small coefficients.
  • Scaling and remaining choices: Scale-invariance is used to rescale data and coordinates to regularize model recovery for poorly scaled datasets.Query-point and model-library selection are discussed but left for more thorough future investigation.
  • Test-function selection: Larger real-space support reduces variance but can suppress signal-dominated modes, whereas smaller support may retain noise and provide insufficient variance reduction.Test-function design must balance variance reduction against resolution of the underlying dynamics.
  • Limitations: The paper leaves systematic investigation of different test-function performance and query-point and model-library choices to future research.These choices remain scope boundaries for the presented algorithm.
  • Piecewise-polynomial test functions: Piecewise-polynomial test functions are assembled from tensor products of one-dimensional functions with compact support, controlled smoothness, and selected support lengths and degrees.Their parameters determine decay in real and Fourier space and must provide enough smoothness for the required integrations by parts.
  • Piecewise-polynomial test functions: The implementation sets τ = 10^-10 and can choose support sizes m by estimating the critical wavenumber separating noise- and signal-dominated modes.The same spatial test function and a separate temporal test function are used in the examples.

5. Examples

WSINDy is evaluated on seven PDEs across noise levels, amplitudes, and library sizes, with robust recovery even for discontinuous and weak solutions. The experiments also examine learned thresholding, scaling, model-specific failure modes, and coefficient accuracy.

  • Seven PDEs are tested across a range of noise levels, amplitudes, and model libraries.
  • Orders of magnitude improvements over derivative-based methods yield reliable recovery for four of seven PDEs at 100% noise and all PDEs under 20% noise.
  • WSINDy combines a convolutional weak form, FFT implementation, automatic sparsity-threshold selection, and continuum-level rescaling.
  • Model Identification: The average TPR remains above 0.95 as noise approaches 1 for several PDEs, while nonlinear Schrödinger remains above 0.95 until 50% noise.
  • Chosen Examples: WSINDy recovers nonlinear transport models from shock-forming data and remains robust to noise and scaling without special algorithmic modifications.
  • Model Identification: Automatic threshold learning is crucial in large-noise regimes because the feasible threshold range changes across systems and can become prohibitively narrow for Kuramoto-Sivashinsky.
  • Model Identification: Reaction-diffusion and Navier-Stokes exhibit model misidentification above 22% and 27% noise, respectively, through reduced or averaged dynamics.
  • Coefficient Accuracy: Noise-free Navier-Stokes coefficients are recovered to less than 3 significant digits, while inviscid Burgers achieves E∞ = 4.3 × 10^-5 on non-smooth data.

6. Conclusion

WSINDy extends weak-form system identification to PDEs without pointwise derivative approximations, combining efficient computation with data-driven sparsification and scaling strategies. The conclusion reports strong robustness across several PDE classes while identifying limits in model selection, statistical error characterization, and test-function optimality.

  • PDE extension: WSINDy extends weak identification to PDE spatiotemporal dynamics without relying on pointwise derivative approximations or black-box closure models.The method is presented as a PDE system-identification framework built around a weak formulation.
  • Automatic tuning: Learning algorithm components reduce reliance on manual hyperparameter tuning, including automatic learning of sparsification-related parameters.The conclusion specifically highlights learning many hyperparameters and the threshold b_λ from the dataset.
  • Computational design: The convolutional weak form enables FFT-based implementation, while spectral decay provides a criterion for selecting admissible test functions.The Fourier formulation connects computational efficiency with test-function design.
  • Challenging data and libraries: Scale invariance and modified sparsification support recovery from large candidate libraries and poorly scaled data.These strategies are presented as extensions for difficult identification settings.
  • Limitations and future directions: Open limitations include the lack of optimality guarantees for the test-function class, incomplete understanding of coefficient accuracy, and the need for generalized least squares.The conclusion also calls for more precise statistical information about the linear system and recovery-error dependence.
  • Empirical scope: WSINDy is reported as robust for nonlinear waves and spatiotemporal chaos, but less robust for limit-cycle data such as reaction-diffusion and Navier–Stokes.The stated examples include Burgers, Korteweg–de Vries, nonlinear Schrödinger, Sine–Gordon, and Kuramoto–Sivashinsky.

Appendix A. Learning Test Functions From Data

The test-function learning algorithm detects noise-dominated Fourier modes and selects coordinate test-function hyperparameters from spectral and real-space constraints.

  • Spectral changepoint detection: Critical wavenumbers separate signal-dominated from noise-dominated modes in each coordinate.The algorithm estimates these changepoints from power spectra after averaging over the other coordinate.
  • Coordinate-wise analysis: The procedure applies separately to spatial and temporal Fourier transforms, using Fx and Ft for coordinate-wise analysis and F for the full transform.The method is illustrated on a one-dimensional spatial grid sampled at multiple timepoints.
  • Spectral changepoint detection: The cumulative-sum spectrum Hx is approximately linear over noise-dominated modes because the Fourier coefficients of white noise are identically distributed.A best two-piece approximation identifies the changepoint as the corner between signal- and noise-dominated regions.
  • Implementation: Changepoints are detected algorithmically using piecewise-linear approximation, with MATLAB implementation through findchangepts.The method uses an algorithm described in prior work to fit the two-piece approximation.
  • Hyperparameter selection: Test-function hyperparameters are chosen by enforcing spectral decay beyond the detected critical wavenumbers and decay in real space.The resulting constraints determine the support and degree parameters of the coordinate test functions.

2. Enforcing Decay.

The decay conditions select test functions whose Fourier spectra retain signal-dominated modes while suppressing noise-dominated modes, with parameters determined by coupled spectral and real-space constraints.

  • Decay conditions: The test functions are selected so their spectra decay into the tails at the detected critical wavenumbers and their real-space values decay at support boundaries.These conditions are imposed for both spatial and temporal coordinate test functions.
  • Fourier-space approximation: For small frequencies, a Gaussian spectrum serves as a proxy for the Fourier transform of the chosen test function.The approximation is justified using agreement of low-order moments and is applied for small ξ, small a, and large p.
  • Parameter coupling: The spectral constraint relates the test-function degree px to its discrete support mx through the detected critical wavenumber.The spatial scaling uses a = mx∆x, linking discrete support to the physical domain length and grid size.
  • Parameter coupling: The real-space decay condition supplies a second relation between px and mx, and combining both constraints yields an equation for mx.After determining mx, px can be obtained from either constraint.
  • Existence and uniqueness: Under N1 > 4, 0 < τ < 1, and the stated critical-wavenumber condition, the support equation has a unique root mx ≥ 2.The relevant function is monotonically decreasing and changes sign on the specified interval.
  • KS illustration: The learned test function filters noise-dominated modes even when the noisy squared data have substantial bias relative to clean squared data.The demonstration uses the Burgers-type nonlinearity ∂x(u^2) on the same noisy KS dataset.
  • KS illustration: For the 50%-noise KS example, the detected changepoint is kx* = 24, and the noise-dominated region is approximately linear.The learned test-function spectrum agrees with a Gaussian reference over signal-dominated modes k ≤ 24.

Appendix B. Numerical Simulation Methods

The numerical experiments use simulated noise-free datasets for the PDEs, with Navier–Stokes generated by an immersed boundary projection method and the other computations performed in MATLAB.

  • Simulation setup: Noise-free datasets are simulated for each PDE used in the numerical experiments.Dataset dimensions are reported separately in Table 3.
  • Simulation setup: Navier–Stokes data are generated with an immersed boundary projection method implemented in C++.This is the exception to the MATLAB-based computational workflow.
  • Implementation: The remaining computations are performed in MATLAB 2019b.

B.1. Inviscid Burgers.

The inviscid Burgers experiment uses an exact shock-forming solution sampled on a large-amplitude space–time grid to test WSINDy under challenging data scales.

  • Exact solution: The Burgers experiment uses an exact shock-forming solution that becomes discontinuous at t = α^-1.The shock travels along x = A.
  • Test-function illustration: The figure compares clean and noisy data with their power spectra and the resulting test-function construction along a fixed-time spatial slice.
  • Parameter regime: The experiment sets α = 0.5 and A = 1000 to evaluate performance on large-amplitude data.
  • Discretization: Noise-free data are sampled on a 256 × 256 space–time grid with ∆x = 31.25 and ∆t = 0.0157.The spatial points begin at x = −4000 and are indexed by (xi, tj).

B.2. Korteweg-de Vries.

The supplied passages describe numerical datasets and subsampling choices used for system identification, including periodic one-dimensional simulations, a complex-valued NLS dataset, and a two-dimensional wave simulation.

  • Korteweg-de Vries: Periodic one-dimensional datasets use ETDRK4 timestepping and Fourier-spectral differentiation before system identification subsampling.One dataset uses 400 spatial and 2400 temporal points with 25% of timepoints retained; another uses 256 spatial and 1500 temporal points with 20% retained.
  • Korteweg-de Vries: The first one-dimensional dataset uses a two-soliton initial condition with A = 25 and B = 16.The initial condition combines two squared sech profiles centered at shifted spatial locations.
  • NLS: The NLS dataset contains 512 spatial points and 502 timepoints, with both dimensions subsampled by 50%.The data are represented through real and imaginary components, w = u + iv, for system identification.

B.6. Reaction-Diffusion.

The reaction-diffusion passages describe simulations of the λ-ω system and show that noise changes which model WSINDy selects: the full model below a threshold and a purely oscillatory reduction above it.

  • Simulation setup: The λ-ω reaction-diffusion system is simulated on a doubly-periodic square using Fourier-spectral differentiation and ode45 time integration.The computational domain is [−10, 10] × [−10, 10] with t ∈ [0, 10], using N1 = N2 = 256 and N3 = 201.
  • Spiral-wave data: The unstable spiral initial data breaks apart over time but still settles into a limit cycle.This dataset differs from the traditional stable spiral wave data referenced from prior work.
  • Model selection: For σNR ≤ 0.1, WSINDy returns the full model, whereas for σNR > 0.1 it detects the oscillatory reduced model.The reduced model provides a good approximation away from the origin, and the full and reduced models are compared from the same spiral-wave initial conditions.
  • Limitation: The stable spiral wave is not an ideal candidate for identifying the full model from data.This limitation is stated specifically from the perspective of data-driven model selection.

B.7. Navier-Stokes.

The supplied Navier–Stokes passages describe an immersed-boundary simulation of flow around a cylinder and the resulting velocity–vorticity dataset used for identification.

  • Simulation setup: Flow around a cylinder is simulated with an immersed boundary projection method and third-order Runge–Kutta timestepping.The simulation uses spatial and temporal resolutions Δx = Δt = 0.02 for 2000 timesteps after vortex shedding reaches its limit cycle.
  • Dataset: The dataset contains velocity components and vorticity away from the cylinder and boundaries.Samples are taken in the rectangle (x, y) ∈ [1, 7.5] × [−1.5, 1.5].
  • Dataset: Ten percent of the data are retained in time, producing a final temporal resolution of Δt = 0.2.The spatial resolution remains Δx = 0.02.
  • Comparison: Figure 9 compares the full reaction-diffusion model with an oscillatory reduced model at T = 10 from the same spiral-wave initial conditions.Only the v component is shown in the comparison described by the caption passages.
  • Comparison: The reduced model provides a good approximation away from the cylinder.This comparison is stated for the flow setting represented in the supplied passages.
Loading 2007.02848v3…