Skip to content

Promote weak form SINDy from a feature library to a sibling SINDy class.#678

Merged
Jacob-Stevens-Haas merged 47 commits intomainfrom
weak
Mar 24, 2026
Merged

Promote weak form SINDy from a feature library to a sibling SINDy class.#678
Jacob-Stevens-Haas merged 47 commits intomainfrom
weak

Conversation

@Jacob-Stevens-Haas
Copy link
Copy Markdown
Member

@Jacob-Stevens-Haas Jacob-Stevens-Haas commented Feb 20, 2026

This PR, long in the making, refactors the weak form. Previously, weak form was implemented as a feature library, and the LHS calculation was guarded by an "if weak: do this, else: calculate derivatives". However, the weak form is a lot more than just a feature library, and the implementation forced a lot of "spooky action at a distance" in other classes and functions, which meant that touching those classes could cause weak sindy to break (e.g. #395). See #351 for more overall motivation.

In order for PDELibrary to support this change, it was easier to remove SINDy-PI from the PDELibrary and make it its own class, as discussed in in #192. This PR implements the recommended fix.

What do we get from refactoring WeakSINDy?

  • For one, we instantly gain the ability to fit multiple trajectories (with potentially different spatiotemporal grids), since st_grid in WeakSINDy is analagous to t in SINDy.
  • We also gain the ability to use predict() in weak problems, because they use normal feature libraries. This includes simulate and score for ODEs.
  • We can use smoothing from different differentiation methods on the data before forming weak features.
  • Special guards in GeneralizedLibrary for tensoring the weak form with a regular library are no longer needed. I believe there were also cases where they may also have been wrong in cases.
  • Along the way, I've isolated a lot of the functionality to support unit tests. These give confidence that I'm keeping things mathematically correct for all of the quadrature/integration by parts. That rigor also allows us to accept a user-supplied test function or test PRs of alternate test functions.
  • Clarifies how to use control features and steady-state solutions (e.g. $$f(x) = Du$$)
  • Choosing subdomains is factored into a helper function that returns a SubdomainSpecs dataclass. People can implement alternate strategies (as mentioned in the original Messenger and Bortz paper) for choosing subdomains over the legacy random choice.
  • The use of a separate SINDy class allows one to implement the noise calculations of Messenger and Bortz to pass to a BINDy like algorithm (See Feature: Adding BINDy (based on Greedy search that maximise evidence) to PySINDy #674)

What does it cost?

WeakPDELibrary is maintained but deprecated because the new WeakSINDy class is not compatible with

  • SINDy-PI. To do this correctly requires more explicitly symbolic input feature names, i.e. recognizing that u_1 is a derivative of the first spatial dimension, whereas u1 is the first system coordinate. This PR is a step in that direction, but ultimately we will need to use sympy.
  • GeneralizedLibrary. However, this PR is is correct with odd tensored/concat constructions of libraries, due to the flattening/distributive operation, where GeneralizedLibrary was not. But because of the distributive operation, it's unclear how inputs_per_library would work. Perhaps if inputs_per_library were passed in transform() (after all, input_features come in fit, not __init__) it would be easier.

Developer notes

The _weak.py module begins with a docstring describing the relevant helper objects.

  • Added performance benchmarks in CI for SINDy and WeakSINDy runtime, peakmemory, and fit/predict performance.
  • Weights are no longer stored in a 2D array, indexed by derivative operation and subdomain. Instead, they are stored in a dictionary mapping from derivative op to a 1D list of weights on each subdomain. Removing an indexing convention makes the lookups clearer.
  • Throughout the repo, t and st_grid now follow the same convention as most other arrays. Everything gets promoted into a minimum 2D AxesArray, with the second-to-last axis is time and the last axis is coordinate.
  • All tensored/concatenated libraries are distributed so as to make integration by parts and term ordering clearer.
  • SemiTerm is an internal object to keep track of an integral term symbolically. They are Semi because sometimes multiple need to be summed up to get a single feature - this was always the case, but now it's explicit what operation each term is doing and what order to assemble the regression matrix.
  • In the future, the regular PDE library may become a superclass of this one to enable regular PDE-FIND of multiple trajectories and remove some of the hacks here.

Future steps

  1. Get JointSINDy merged
  2. Find an alternate way to pass _inputs_per_library, enabling GeneralizedLibrary
  3. Implement PDEFIND in a parent class between SINDy and this class, making PDELibrary merely a dataclass
  4. more symbolically-aware function libraries
  5. some private methods shared by SINDy and WeakSINDy that SINDy-PI can invoke in order to correctly construct the feature library (rather than using feature_library.transform(), as now).

Remaining to do:

  • Rebuild SINDy-PI as a separate class that uses pysindy objects, rather than fitting awkwardly into other classes.
  • Find the missing -1 factor
  • Add docstrings
  • Implement partial simulate()/score() for ODEs
  • Add tests for u_dot_wk
  • Add test for eval_semiterm()
  • Add test for local einsum fix
  • Check control
  • Add runnable example to the docstring.
  • Linting (Syntax for tests requires upgrading to 3.12, which requires upgrading pre-commits to conflicting versions. This is a TODO: For afterwards).
  • Add LHS to fit() to allow steady-state solutions (e.g. $$f(x) = Du$$)
  • Remove or modify reshape_samples_to_spatial_grid now that we don't use it for weak (may still be needed for PDELibrary) Fixed PDE prediction shapes across SINDy and WeakSINDy. Still probably need to remove this method, though.
  • Make SemiTerm a class, not just a tuple type alias
  • Verify max derivative order for test function is correctly implemented
  • Check that _linear_weights works for arbitrary test function

Jacob-Stevens-Haas and others added 30 commits February 24, 2026 18:26
Long term, some revision of TensoredLibrary is in order, as the
inputs_per_library is GeneralizedLibrary bleeding in.  This is a consequence of
the fitted/unfitted mutability of scikitlearn estimators.

Along the way, fix a potential bug in GeneralizedLibrary and extract PDE feature
naming function.
Enforces input shape commonality across all models
This ensures that coefficients and feature names will line up.
It does require reassigning the sorted library to self in order
to make predict work.
In addition to having a differentiation method, this change adds simulate()
and score() to WeakSINDy.
Fixed in previous commit, but einsum was seeing the ellipses on
diferent arguments as different (and not shared) axes.

Two cases:

- Ellipsis in two subscripts handled as shared axes
- Conflicting ellipsis axis names across operands raise a ValueError.

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Add two tests exercising _eval_semiterm directly:

- test_eval_semiterm_output_shape: verifies the result has shape
  (n_subdomains, n_coord) for a trivial no-derivative, no-feature term.
- test_eval_semiterm_phi_prime_vanishes: verifies that integrating a
  constant against phi' gives zero, since the bump test function
  vanishes at both boundaries (an analytically exact result).

Also keep more types as AxesArray and handle spatial derivatives inside
semiterms more correctly based on user arguments.

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
This allows discovery of steady-state differential equations, e.g.
Poisson's equation.  This requires some functions to allow eliminating the
time derivative, and perhaps that's an uncomfortable amount of spooky action.
Would love a more elegant way, but the things I do for backwards compatibility.

For such systems, convert_u_dot_integral does *not* move derivative
onto test function, assuming they are handled by the caller
This allows discovery of steady-state differential equations, e.g.
Poisson's equation.  This requires some functions to allow eliminating the
time derivative, and perhaps that's an uncomfortable amount of spooky action.
Would love a more elegant way, but the things I do for backwards compatibility.

For such systems, convert_u_dot_integral does *not* move derivative
onto test function, assuming they are handled by the caller
Add `ParallelImplicitSINDy` in `pysindy/_sindypi.py`, to replace the
removed functionality from SINDyPI optimizer and implicit=True in PDELibrary.
This does not work with WeakSINDy.

In the future, weak SINDy PI can be restored with (a) a new
PDEFIND class that handles spatial_grid the same way as WeakSINDy, then (b)
more symbolic handling of features.  That would allow WeakSINDy and PDEFIND
the opportunity to understand derivative features in the RHS better and provide
and API for ParallelImplicitSINDy to call SINDy subclasses to transform the
derivative and non-derivative features.

Deprecate the `SINDyPI` optimizer (which was a CVXPY-based SR3 variant)
with a `DeprecationWarning` pointing users to the new class.

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Add docstrings to WeakSINDy.fit, WeakSINDy.simulate, and _make_domains.

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
@Jacob-Stevens-Haas Jacob-Stevens-Haas force-pushed the weak branch 3 times, most recently from f6e0948 to 0d9700f Compare February 25, 2026 23:11
@Jacob-Stevens-Haas
Copy link
Copy Markdown
Member Author

Jacob-Stevens-Haas commented Feb 25, 2026

| Change   | Before [b7558be6] <base>   | After [16abb9ee] <pr>   |   Ratio | Benchmark (Parameter)                                   |
|----------|----------------------------|-------------------------|---------|---------------------------------------------------------|
| +        | 1.27±0.01ms                | 1.47±0.01ms             |    1.16 | finite_difference_1D.FiniteDifferenceBM.time_derivative |
| -        | 63.7±20ms                  | 23.1±0.2ms              |    0.36 | special.PDEFindSteady.time_pdefind_steady               |

SOME BENCHMARKS HAVE CHANGED SIGNIFICANTLY.
PERFORMANCE DECREASED.

So it looks like some early work I did to improve the clarity of how different axes are differentiated in FiniteDifference caused a 15% slowdown. The work was necessary because finite difference assumed certain things about axis that may not always be true, and reverting just the finite difference code breaks.

I think it's an acceptable slowdown, given that this is truly the beginning of our benchmarking experience in pysindy and previously these regressions would've gone unnoticed.

Also, the PR speeds up PDEFIND on a steady-state solution by ~64%! It was pretty cool to see how much the overhead of allowing SINDy-PI functionality to occur in PDELibrary was costing.

@Jacob-Stevens-Haas
Copy link
Copy Markdown
Member Author

Jacob-Stevens-Haas commented Feb 25, 2026

@llfung FYSA this touches BINDy (in pretty trivial ways):

  • _check_multiple_trajectories now also looks at t
  • formatting change from pre-commit version
  • replaces non-public call diff_method._differentiate(x, t) with the public diff_method(x, t)

More interestingly, it enables the kind of numerical-accuracy benchmarks you originally wrote for BINDy as a lorenz test file. If you're interested adding something like it back in, see this file for the skeleton - you'd just need to make a benchmark subclass that instantiates the BINDy() model . Example CI output:

[51.04%] ··· defaults.SINDyLorenz.peakmem_fit                              277M
[52.08%] ··· defaults.SINDyLorenz.peakmem_fit_predict                      277M
[53.12%] ··· defaults.SINDyLorenz.time_fit                          10.8±0.07ms
[54.17%] ··· defaults.SINDyLorenz.time_fit_predict                   13.5±0.6ms
[55.21%] ··· defaults.SINDyLorenz.track_score_fit             1.744562254953333
[56.25%] ··· defaults.SINDyLorenz.track_score_fit_predict   0.03302684681199324
[57.29%] ··· defaults.WeakLorenz.peakmem_fit                               279M
[58.33%] ··· defaults.WeakLorenz.peakmem_fit_predict                       279M
[59.38%] ··· defaults.WeakLorenz.time_fit                               149±4ms
[60.42%] ··· defaults.WeakLorenz.time_fit_predict                       151±3ms

@Jacob-Stevens-Haas

This comment was marked as outdated.

This comment was marked as outdated.

Copy link
Copy Markdown

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Pull request overview

Refactors weak-form SINDy into a first-class WeakSINDy model (rather than a special-cased feature library), and splits SINDy-PI into a dedicated ParallelImplicitSINDy class, while also standardizing array/axes conventions across core APIs, differentiation, and tests.

Changes:

  • Added WeakSINDy (pysindy/_weak.py) with supporting weak-form integration utilities and a comprehensive new test suite.
  • Added ParallelImplicitSINDy (pysindy/_sindypi.py) and deprecated the legacy SINDyPI optimizer path.
  • Updated core input validation + differentiation call patterns to prefer differentiators as callables (vs ._differentiate) and aligned CI/tooling to newer Python/pre-commit versions.

Reviewed changes

Copilot reviewed 38 out of 40 changed files in this pull request and generated 5 comments.

Show a summary per file
File Description
test/utils/test_axes.py Adds einsum ellipsis tests for AxesArray behavior.
test/test_weak.py New unit tests covering weak-form utilities and WeakSINDy.
test/test_sindyc.py Updates expected exceptions/messages for new input validation rules.
test/test_pysindy.py Updates tests for revised trajectory validation + adds ParallelImplicitSINDy tests.
test/test_optimizers/test_optimizers.py Adjusts optimizer coef shape expectations to standardized conventions.
test/test_feature_library.py Updates library tests to match PDELibrary API changes and removes some fixtures.
test/test_differentiation.py Migrates tests from ._differentiate to differentiator call syntax and updates shapes.
test/test_bindy.py Aligns BINDy tests with new input validation exception types.
test/conftest.py Reshapes fixtures to match updated array conventions and removes data_1d_bad_shape.
pysindy/utils/bindy.py Switches to callable differentiator interface.
pysindy/utils/base.py Introduces validate_no_reshape for stricter validation + updates control validation.
pysindy/utils/_axes.py Tightens einsum ellipsis handling and erroring on conflicting axis names.
pysindy/optimizers/sindy_pi.py Deprecates SINDyPI optimizer in favor of ParallelImplicitSINDy.
pysindy/optimizers/evidence_greedy.py Minor cleanup/formatting.
pysindy/optimizers/base.py Makes optimizer interface more explicit (abstract fit).
pysindy/feature_library/weak_pde_library.py Deprecation cleanup and internal differentiation call updates.
pysindy/feature_library/pde_library.py Refactors PDELibrary toward derivative-only behavior + new feature name helper.
pysindy/feature_library/generalized_library.py Adjusts get_feature_names input feature handling.
pysindy/feature_library/base.py Removes legacy correct_shape; adds EmptyLibrary; tweaks feature-name handling.
pysindy/differentiation/sindy_derivative.py Adjusts differentiation to new time/axes conventions.
pysindy/differentiation/finite_difference.py Refactors coefficient building around AxesArray time representation.
pysindy/differentiation/base.py Adds BaseDifferentiation.__call__ normalization + AxesArray wrapping.
pysindy/_weak.py New weak-form implementation + WeakSINDy model.
pysindy/_typing.py Updates type aliases and introduces TrajectoryType helper.
pysindy/_sindypi.py New ParallelImplicitSINDy implementation.
pysindy/_core.py Reworks multiple-trajectory validation/adaptation and shape standardization.
pysindy/__init__.py Exposes WeakSINDy and ParallelImplicitSINDy at package top-level.
pyproject.toml Bumps minimum Python to 3.11 and adds dev dependency sindy_exp.
examples/tutorial_1/example.py Minor formatting tweak.
examples/1_feature_overview/example.py Updates examples to use callable differentiators and revised PDE usage.
docs/conf.py Enables sphinx doctest extension.
asv_bench/benchmarks/special.py Adds a PDE steady-state benchmark (currently contains runtime issues).
asv_bench/benchmarks/defaults.py Adds benchmark scaffolding/mixins for SINDy/WeakSINDy.
asv_bench/asv.conf.json Updates ASV install command to include extras.
.pre-commit-config.yaml Updates black/flake8 versions.
.github/workflows/release.yml Bumps release workflow Python to 3.11.
.github/workflows/notebooks.yml Bumps notebooks workflow Python to 3.11.
.github/workflows/main.yml Updates CI Python versions + switches to pre-commit GitHub Action + adds ASV publishing steps.
.flake8 Comment clarifications for ignored rules.

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

Jacob-Stevens-Haas and others added 3 commits February 27, 2026 14:37
Also fix bug in EmptyLibrary.transform, set_params (Weak) return
self, and remove unused (and broken) fixture

Co-authored-by: copilot-swe-agent[bot] <198982749+Copilot@users.noreply.github.com>
Copy link
Copy Markdown
Collaborator

@yb6599 yb6599 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hey Jake, the documentation and type annotations are thorough. Aside from some minor comments, LGTM.

@Jacob-Stevens-Haas Jacob-Stevens-Haas merged commit c4421fc into main Mar 24, 2026
7 of 10 checks passed
@Jacob-Stevens-Haas Jacob-Stevens-Haas deleted the weak branch March 24, 2026 21:53
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants