Hyperparameter Optimization in the Estimation of PDE and Delay-PDE models from data

This paper proposes an improved method for estimating partial differential and delay partial differential equations from data by integrating time integration with Bayesian optimization and the Bayesian information criterion to automatically optimize hyperparameters, thereby enhancing robustness and expanding the modeling scope across various physical systems.

The Gist

Imagine you are a detective trying to figure out the secret recipe for a delicious soup, but you can't see the ingredients or the chef. All you have is a video of the soup bubbling in the pot over time. You see the steam rise, the vegetables swirl, and the color change, but you don't know why it's happening.

In the world of science, this "soup" is a physical system (like weather patterns, how cells grow, or how a virus spreads), and the "recipe" is a mathematical equation (a Partial Differential Equation, or PDE) that describes exactly how the system behaves.

For a long time, scientists have tried to reverse-engineer these recipes from data. However, it's like trying to guess the recipe while the chef is changing the heat, adding salt, and stirring at the same time. It's messy, and if you get one detail wrong, the whole soup tastes terrible.

This paper introduces a new, smarter detective method called Hyperparameter Optimization to solve this problem. Here is how it works, broken down into simple concepts:

1. The "Giant Ingredient List" (The Library)

Imagine you have a massive cookbook with every possible ingredient and cooking technique imaginable: "add salt," "stir clockwise," "heat to 100 degrees," "wait 5 minutes," "add a pinch of magic dust."

The scientists create a similar "library" of mathematical building blocks. They don't know which ones are in the real soup, so they write down thousands of possibilities.

• The Problem: If you just try to fit all these ingredients to the video, you end up with a recipe that is way too complicated (e.g., "Add salt, but also add a tiny bit of pepper, but also a drop of water, but also..."). This is called overfitting. It works for this specific video, but if you cook the soup again, it fails.

2. The "Smart Filter" (Sparsity and Thresholds)

To fix the over-fitting, the method uses a "Smart Filter." It asks: "Which ingredients are actually doing the heavy lifting?"

• It tries to throw away the tiny, insignificant ingredients (like a pinch of salt that doesn't change the taste) and keeps only the big, important ones.
• The Catch: How do you decide what counts as "tiny"? You need a Threshold.
  • If the threshold is set too high, you throw away important ingredients (like the salt).
  • If it's set too low, you keep the junk (like the magic dust).
  • Old methods used a "guess-and-check" approach to find the right threshold. It was slow and often got it wrong.

3. The "Taste Test" (Bayesian Optimization)

This is the paper's big innovation. Instead of guessing the threshold, they use a Bayesian Optimization system. Think of this as a super-smart AI taste-tester.

• The AI tries a threshold, cooks a "virtual soup" (simulates the equation), and tastes it.
• If the virtual soup looks nothing like the real video, the AI learns: "Okay, that threshold was too high/low. Let's try a different one."
• It repeats this thousands of times, learning from every mistake, until it finds the perfect threshold that makes the virtual soup match the real video perfectly.

4. The "Time Travel" Twist (Delay Equations)

Some systems have a "lag." For example, if you eat a lot of food, you don't get full immediately; it takes 20 minutes.

• The old methods struggled with this "time delay."
• The new method treats the delay time (e.g., "20 minutes") as just another ingredient to be optimized. The AI taste-tester tries different delay times until it finds the one that explains the lag in the data.

5. Why is this better? (The "Robustness" Factor)

The authors tested their method on several "recipes" (mathematical models) that represent real-world physics:

• The "Messy" Data: Real-world data is often "noisy" (like a shaky camera) or "sparse" (missing frames). The old methods often failed here, giving up or giving wrong answers. The new method, by checking how the whole simulation plays out over time (not just a single snapshot), is much more robust. It can handle missing data and noise better.
• The "Complex" Systems: They tested it on systems where different parts behave very differently (some change fast, some slow). The new method can set different "filters" for different parts of the system, whereas old methods used one size fits all.

The Bottom Line

Think of this paper as upgrading from a guessing game to a self-driving car for discovering scientific laws.

• Old Way: You drive a car blindfolded, guessing the road ahead, hoping you don't crash.
• New Way: You have a GPS (Bayesian Optimization) that constantly checks the map, adjusts the speed, and steers you perfectly to the destination (the correct equation), even if the road is bumpy or foggy.

The result is a tool that can look at messy, real-world data and automatically figure out the hidden mathematical laws governing it, without needing a human to tweak the settings manually. This opens the door to discovering new physics in fields like biology, climate science, and engineering much faster than before.

Technical Summary

1. Problem Statement

The paper addresses the challenge of data-driven discovery of governing physical laws, specifically for Partial Differential Equations (PDEs) and Delay Partial Differential Equations (Delay-PDEs).

• Context: While methods for discovering Ordinary Differential Equations (ODEs) from data (e.g., SINDy) are well-established, extending these to spatially extended systems (PDEs) introduces significant challenges regarding scalability, noise sensitivity, and the handling of complex physical constraints (like mass conservation).
• Limitations of Existing Methods:
  • Many existing methods rely on minimizing the residual of the time derivative ($\partial_t u$) directly. This approach is sensitive to noise and often fails to produce models that are stable or accurate when integrated over time.
  • Hyperparameter Tuning: Critical hyperparameters (such as sparsity thresholds, time delays, or regularization strengths) are typically set manually via trial-and-error or computationally expensive grid searches.
  • Model Complexity: Standard libraries often lack the flexibility to handle specific physical constraints (e.g., conservation laws) or non-standard terms (e.g., time delays) without ad-hoc modifications.
  • Overfitting: Models often fail to replicate initial data or become numerically unstable during simulation because they are optimized solely for instantaneous fit rather than long-term predictive capability.

2. Methodology

The authors propose an improved framework that integrates Bayesian Optimization with Sequential Thresholding Least Squares (STLS) and the Bayesian Information Criterion (BIC).

A. Core Optimization Framework

The method seeks to find a sparse set of coefficients $\hat{\sigma}$ for a library of candidate functions $\Theta(u)$ that best describes the dynamics $\partial_t u = \sigma \cdot \Theta(u)$.

1. Initial Estimation: An initial sparse solution is found using STLS to minimize the residual of the time derivative: $R(\sigma) = ||\partial_t u - \sigma \cdot \Theta(u)||_2$.
2. Time Integration: Instead of stopping at the derivative residual, the method performs numerical time integration of the estimated model to generate a reconstructed time series $\hat{u}$.
3. Objective Function (BIC): The optimization targets the Bayesian Information Criterion (BIC) rather than just the derivative residual. The BIC balances the goodness of fit (minimizing the error between the original data $u$ and the integrated reconstruction $\hat{u}$) against model complexity (penalizing the number of non-zero parameters $s$):
   $$\text{BIC} = s \ln(N_t) - 2 \ln(\tilde{L})$$
   where $\tilde{L}$ is the likelihood derived from the integrated error $||u - \hat{u}||_2$.
4. Hyperparameter Search: The authors utilize Hyperopt with a Tree-structured Parzen Estimator (TPE) to automatically optimize hyperparameters. These include:
   • Sparsity Thresholds ($h_i$): Unlike previous methods using a single global threshold, this method allows for individual thresholds for different field variables or groups of terms (e.g., diffusion terms vs. reaction terms).
   • Structural Parameters: The framework can optimize parameters that alter the library itself, such as time delays ($\tau$) in Delay-PDEs.

B. Implementation Details

• Library Construction: A flexible library of ansatz functions is constructed, including polynomial terms, spatial derivatives (up to order $p$), and user-defined terms (e.g., $u(t-\tau)$).
• Numerical Stability: To handle cases where candidate PDEs are numerically unstable during the optimization loop, the method employs a fallback to a trapezoidal integration scheme if the standard Runge-Kutta solver fails.
• Software: The method is implemented in a Python package called TSME (Time Series Model Estimation).

3. Key Contributions

1. Integration of Time Integration into Model Selection: By optimizing based on the deviation of the integrated time series rather than just the instantaneous derivative, the method significantly increases robustness against noise and under-sampling.
2. Automated Hyperparameter Optimization: The use of Bayesian optimization (TPE) to automatically tune sparsity thresholds and structural parameters (like time delays) removes the need for manual tuning and grid searches.
3. Flexible Thresholding: The ability to assign distinct thresholds to different variables or groups of terms allows the method to handle systems where parameters vary by orders of magnitude (e.g., fast vs. slow variables in reaction-diffusion systems).
4. Extension to Delay-PDEs: The framework successfully incorporates time delays as optimizable hyperparameters, enabling the discovery of Delay-PDEs from data, a capability previously limited in sparse identification frameworks.
5. Handling Conservation Laws: The method successfully recovers the Cahn-Hilliard equation (which enforces mass conservation) without requiring explicit constrained optimization, demonstrating flexibility in handling complex physical constraints.

4. Results and Benchmarks

The method was validated on several synthetic benchmarks representing different classes of physical behavior:

• Complex Ginzburg-Landau (cGL) Equation:
  • Performance: The method achieved results comparable to PySINDy on well-sampled data but demonstrated superior robustness on under-sampled data.
  • Sparsity: While PySINDy (using a fixed threshold) often included incorrect terms or failed to identify the correct structure when data was sparse, the proposed BIC-optimized method consistently identified the correct sparse model.
• Phase-Field Models (Allen-Cahn & Cahn-Hilliard):
  • Allen-Cahn: Successfully recovered the equation with <2% error.
  • Cahn-Hilliard: Successfully identified the equation involving second-order spatial derivatives of nonlinear terms ($\nabla^2(u^3)$), a structure that standard PySINDy libraries often struggle with due to restrictive ansatz definitions.
• FitzHugh-Nagumo (FHN) System:
  • Demonstrated the necessity of multiple thresholds. A single global threshold failed to capture the system because the reaction and diffusion terms had vastly different magnitudes. Using variable thresholds ($h_u, h_v$) allowed for the correct identification of both equations simultaneously.
• Chaotic Regimes (cGL with Cross-Diffusion):
  • In chaotic regimes, the method correctly identified the system. When standard single-threshold methods failed to eliminate spurious diffusion terms in "intermittency" regimes, the use of group-specific thresholds for diffusion terms resolved the issue.
• Fisher-KPP with Time Delay:
  • The method successfully identified the time delay $\tau$ as an optimizable parameter, recovering the correct delay value ($\tau \approx 1.0$) and the governing equation from data.

5. Significance and Conclusion

• Robustness: The primary significance is the shift from optimizing for instantaneous fit to optimizing for predictive accuracy via time integration. This makes the method significantly more reliable for real-world applications where data is noisy or sparse.
• Automation: By automating hyperparameter selection, the method reduces the "black box" nature of model discovery, making it accessible to users without deep expertise in tuning regularization parameters.
• Scope Expansion: The ability to handle Delay-PDEs and mass-conserving systems without custom constraints broadens the applicability of data-driven discovery to biological systems (e.g., population dynamics with maturation delays) and materials science (phase separation).
• Future Outlook: The authors suggest that while computational costs are higher than simple least-squares methods, the trade-off is justified by the quality of the resulting models. Future work should focus on applying these methods to experimental data (e.g., microscopy, tokamak discharges) and incorporating sub-sampling or gradient-based optimizations to further improve efficiency.

In summary, this paper presents a robust, automated framework for discovering PDEs and Delay-PDEs that overcomes the limitations of current state-of-the-art methods by integrating time integration into the model selection process and utilizing Bayesian optimization for hyperparameter tuning.