A question concerning examples/12_weakform_SINDy_examples.ipynb

[georgemilosh]

On command line [3] where Lorenz with 20 percent noise is tried I am confused by the comment:

Performance clearly improves as the number of subdomains increases! We can also see that Lorenz is correctly identified despite ~20% noise levels.

Now, first of all, we cannot see Lorenz equation directly because there was no model output statement on line [3] but I added it

# Generate measurement data
dt = 0.002
t_train = np.arange(0, 10, dt)
t_train_span = (t_train[0], t_train[-1])
u0_train = [-8, 8, 27]
u_train = solve_ivp(
    lorenz, t_train_span, u0_train, t_eval=t_train, **integrator_keywords
).y.T

# Instantiate and fit the SINDy model with u_dot
u_dot = ps.FiniteDifference()._differentiate(u_train, t=dt)
model = ps.SINDy()
model.fit(u_train, x_dot=u_dot, t=dt)
model.print()

# Define weak form ODE library
# defaults to derivative_order = 0 if not specified,
# and if spatial_grid is not specified, defaults to None,
# which allows weak form ODEs.
library_functions = [lambda x: x, lambda x: x * x, lambda x, y: x * y]
library_function_names = [lambda x: x, lambda x: x + x, lambda x, y: x + y]
ode_lib = ps.WeakPDELibrary(
    library_functions=library_functions,
    function_names=library_function_names,
    spatiotemporal_grid=t_train,
    is_uniform=True,
    K=100,
)

# Instantiate and fit the SINDy model with the integral of u_dot
optimizer = ps.SR3(
    threshold=0.05, thresholder="l1", max_iter=1000, normalize_columns=True, tol=1e-1
)
model = ps.SINDy(feature_library=ode_lib, optimizer=optimizer)
model.fit(u_train)
model.print()

and it seemed that the models that were discovered where of polynomial order 2, i.e. Lorenz was not "correctly" identified. I do agree with what follows, e.g. the error that is plotted on command line [4] seems to, indeed, indicate that accuracy improves with the number of sub domains.

[georgemilosh]

Burger equation seems fine, but when I try compressible 2D fluid I get with SINDy (not weak):

Noiseless differential fit:
(u)' = -1.000 uv_2 + 0.998 vp_22 + -1.002 uu_1 + -1.002 p^-1p_1 + 1.000 vp_11
(v)' = -1.002 pv_2 + -1.000 p^-1p_2 + 1.000 u^-1p_22 + -0.997 pu_1 + 0.998 u^-1p_11
(p)' = -0.997 pp_2 + -0.999 v^-1v_2 + -1.000 up_1 + -0.998 v^-1u_1
R2 score:  0.9999402527135505
Noisy differential fit:
(u)' = -1.971 u + -0.899 uv_2 + 0.837 u^-1p_2 + -0.593 uu_1
(v)' = -1.624 v + -4.063 vp_2 + -0.838 pv_2
(p)' = -0.995 pp_2 + -0.997 v^-1v_2 + -0.997 up_1 + -0.987 v^-1u_1
R2 score:  0.8679715521429276
CPU times: user 21.5 s, sys: 3.13 s, total: 24.6 s
Wall time: 24.6 s

The R2 score seems reasonable but the expressions are different. Why could this be?

[znicolaou]

Hi @georgemilosh , not sure what you mean by the first question, I see the output from the model.print() statements when I run the notebook from a clone of the repository and also on the docs page: https://pysindy.readthedocs.io/en/latest/examples/12_weakform_SINDy_examples/example.html.

(x0)' = -9.999 x0 + 9.999 x1
(x1)' = 27.992 x0 + -0.999 x1 + -1.000 x0 x2
(x2)' = -2.666 x2 + 1.000 x0 x1
(x0)' = -10.000 x0 + 10.000 x1
(x1)' = 28.000 x0 + -1.000 x1 + -1.000 x0x2
(x2)' = -2.667 x2 + 1.000 x0x1

Also, the Lorenz equations being discovered are second-order polynomials: https://en.wikipedia.org/wiki/Lorenz_system, and the coefficients discovered were those used to generate the data.

For the second question, an R2 score of 0.867 is not very good. SINDy is attempting to minimize the error in the prediction, and it does an okay job of discovering a model with a predicted time derivative that has some correlation with the actual time derivatives for the training data, but the model is not correct (it is overfitting the noise). If we were to evaluate the R2 score test data that differed from the training data, it would be much worse (that would probably be more illustrative, but we didn't generate a test trajectory in the example for simplicity). The example is illustrating that noise causes the nonweak form of SINDy to produce a bad model, while the weak form of SINDy in the next cell does discover the correct model.

[znicolaou]

Ah, I see what you're saying about cell [3]. The label for "correctly identified" was referring to the later output in cell [5], but it isn't entirely accurate either. The model gets the coefficient of the most important terms close, but includes extra small terms too. The "correct" was referring mainly to the better forecasts in the plots in [5], but I'll rephrase it a little.

I'll check the compressible flow output too. There are some updates the notebook needs, including the deprecated quiet=True option.

[znicolaou]

When I run the compressible flow case, my output is correct:

(u)' = -1.000 vu_2 + 0.998 p^-1u_22 + -1.002 uu_1 + -1.002 p^-1p_1 + 1.000 p^-1u_11
(v)' = -1.002 vv_2 + -1.000 p^-1p_2 + 1.000 p^-1v_22 + -0.997 uv_1 + 0.998 p^-1v_11
(p)' = -0.997 pv_2 + -0.999 vp_2 + -1.000 pu_1 + -0.998 up_1
R2 score:  0.9999402527139254
Noisy differential fit:
(u)' = -1.971 u + -0.899 vu_2 + 0.837 p^-1v_2 + -0.593 uu_1
(v)' = -1.624 v + -4.063 p^-1u_2 + -0.838 vv_2
(p)' = -0.995 pv_2 + -0.997 vp_2 + -0.997 pu_1 + -0.987 up_1
R2 score:  0.8679715520180263

It looks like you have an older version that had the library feature names mislabelled. What is your output for

import pysindy
pysindy.__version__

Mine is '1.7.3.dev71+gf6b970b.d20230320'

[znicolaou]

FYI I reinstalled pysindy to be sure in a new environment, and still get the right output for the compressible flow, now with version '1.7.6.dev132+gcddc1e7.d20231006'

[znicolaou]

Ah, okay, I can reproduce your output now--the release of pysindy available in pypi does indeed have the library terms mislabelled. We need to make a new release. In the meanwhile, to get the correct ordering, clone the repository and install it from source

git clone git@github.com:dynamicslab/pysindy.git
cd pysindy
pip install .[dev,cvxpy]

[georgemilosh]

Yes, deprecated quiet and also unbias. Now we get consistent outputs. Thanks

[georgemilosh]

Also, I guess this will be implemented in 2.0 but just in case to make a point about this, with the current build (I think I am on cddc1e7 ) there is also a problem: TypeError: SINDy.fit() got an unexpected keyword argument 'ensemble'.

[georgemilosh]

Btw, for .zsh mac users the solutions would be:

git clone git@github.com:dynamicslab/pysindy.git
cd pysindy
pip install ".[dev,cvxpy]"

I also had to install cmake

I am guessing the recommended version of python is 3.8? I think for me 3.10 also worked.

[znicolaou]

Ah, perfect, thanks for the heads up. FYI the new ensembling approach is in the EnsembleOptimizer class, and a usage example is in the first example notebook.