[BUG] TrappingSR3 finds unstable model

[Jacob-Stevens-Haas]

Hey, @akaptano , another question for you. Fitting the data in #397 with a TrappingSR3 optimizer finds a model that grows unbounded.

Reproducing code example:

Using the data from #397, fitting a TrappingSR3 Optimizer gives a linear, growing model. Thus, simulation grows exponentially and kills integration:

model = ps.SINDy(
    differentiation_method = ps.SmoothedFiniteDifference(),
    feature_library = ps.PolynomialLibrary(degree=2, include_bias=False),
    optimizer = ps.TrappingSR3(),

)
model.fit(arr, t=5e-4)
model.print()
model.simulate(arr[0], t=np.arange(100)

Error message:

It prints the discovered model ((x0)' = 11.570 x0) and then the integration hangs, finally giving the LSODA issue:

ValueError: Input X contains infinity or a value too large for dtype('float64').

Discussion

It is my understanding that although the data does not contain the guarantees of trapping...

The trapping algorithm only applies to fluid and plasma flows with energy-preserving, quadratic nonlinear structure, so we need to explicitly constrain the coefficients to conform to this structure.

... that nevertheless, trapping would constrain the coefficients to result in a model that is globally stable. Am I incorrect?

Also, and I'm not sure if this is relevant, in the larger problem, setting include_bias=True gives the TrappingSR3 warning that "The feature library is the wrong shape or not quadratic, so please correct this if you are attempting to use the trapping algorithm with the stability term included. Setting PL and PQ tensors to zeros for now." However, it appears to discover a globally stable model so long as this constant term is allowed. This constant appears (implicitly, via default arguments) in the TrappingSR3 docstring example (which also raises the "wrong shape" warning). Is this constant significant?

pysindy version

master

[akaptano]

Sorry I have been so busy this summer and haven't got a chance to look at most of the open issues you want my input in. For this one, the answer is rather simple -- the guarantee only holds if you are using the energy-preserving constraints that are used in the example notebook 8. Moreover, your identified model is purely linear, and the trapping theorem anyways only applies for quadratically nonlinear models. There are a number of ways we should forbid/clarify certain uses of the trappingSR3 optimization to clarify the use cases, but a master's student @MPeng5 will open a pull request for an extension of the trappingSR3 optimizer soon. It is a better version, and we can make some of the fixes (like forbidding certain nonsensical uses of this optimizer) when that pull request opens up. Thanks for your work!

[Jacob-Stevens-Haas]

Hey no worries, I've been able to make it through reasonably well and making pretty good progress. Up until the most recent few issues, the issues/PRs were all things I thought I fully understood but wanted to make sure to give you a chance to voice a concern/veto because they were (small) breaking changes or related to your code.

The linearity in the above model is a consequence of trying to find a smaller example that could be copy-pasted, but the optimizer calculated a zero coefficient for x0^2. In the full system I'm working with, TrappingSR3 calculates nonzero coefficients for the quadratic terms which still blow up in simulation. Moreover the full data appears to have stable, decaying oscillations, which is why Trapping & a quadratic library seemed to make sense:

With the energy-preserving constraints - I thought the optimizer enforced them, regardless of the data, and if the data didn't actually obey those constraints, the simulated model would perform poorly. Is that not the case?

[Jacob-Stevens-Haas]

I'm looking forward to seeing what @MPeng5 is working on! Is it modifying the existing TrappingSR3 class, or adding a new one?

If it's possible, it might be an opportunity to refactor a lot of the common parts of SR3 and subordinate classes so as to reduce the amount of repeated code. There's a lot of copied logic in how trimming is applied, how the cvxpy optimization is formed, how constraints, verbose iteration, &c are applied.

[MPeng5]

Basically we just modified the trappingSR3 and added a new loss term there. This new one will not use constraints though. Hopefully we can get the jupyter notebook done soon.