Discovering a simple PDE with noisy data -- need suggestion for specifying differentiation_method

[Lawrence5]

Starting with a simple 2nd order PDE:

I solved analytically and generated training data with 100 time intervals and 100 spacial intervals. PYSINDY produces excellent results.

However, when I add a small amount of Gaussian noise to the training data, I need to specify smoothing in the differentiation method. Every method I tried gives me a similar error. Can anyone give me a suggestion? Thank you very much.

model = ps.SINDy(feature_library=pde_lib, optimizer=optimizer,differentiation_method=trend_filtered) model.fit(u_flat, t=dt) model.print()

`

ValueError Traceback (most recent call last)
Cell In[4], line 84
81 model.print()
83 model = ps.SINDy(feature_library=pde_lib, optimizer=optimizer,differentiation_method=trend_filtered)
---> 84 model.fit(u_flat, t=dt)
85 model.print()

File ~/.local/lib/python3.10/site-packages/pysindy/pysindy.py:343, in SINDy.fit(self, x, t, x_dot, u, multiple_trajectories, unbias, quiet, ensemble, library_ensemble, replace, n_candidates_to_drop, n_subset, n_models, ensemble_aggregator)
337 u = validate_control_variables(
338 x,
339 u,
340 trim_last_point=(self.discrete_time and x_dot is None),
341 )
342 self.n_control_features_ = u[0].shape[u[0].ax_coord]
--> 343 x, x_dot = self._process_multiple_trajectories(x, t, x_dot)
345 # Set ensemble variables
346 self.ensemble = ensemble

File ~/.local/lib/python3.10/site-packages/pysindy/pysindy.py:665, in SINDy._process_multiple_trajectories(self, x, t, x_dot)
663 x = [xi[:-1] for xi in x]
664 else:
--> 665 x_dot = [
666 self.feature_library.calc_trajectory(
667 self.differentiation_method, xi, ti
668 )
669 for xi, ti in _zip_like_sequence(x, t)
670 ]
671 return x, x_dot

File ~/.local/lib/python3.10/site-packages/pysindy/pysindy.py:666, in (.0)
663 x = [xi[:-1] for xi in x]
664 else:
665 x_dot = [
--> 666 self.feature_library.calc_trajectory(
667 self.differentiation_method, xi, ti
668 )
669 for xi, ti in _zip_like_sequence(x, t)
670 ]
671 return x, x_dot

File ~/.local/lib/python3.10/site-packages/pysindy/feature_library/base.py:87, in BaseFeatureLibrary.calc_trajectory(self, diff_method, x, t)
85 def calc_trajectory(self, diff_method, x, t):
86 axes = x.dict
---> 87 x_dot = diff_method(x, t=t)
88 return AxesArray(x_dot, axes)

File ~/.local/lib/python3.10/site-packages/pysindy/differentiation/base.py:49, in BaseDifferentiation.call(self, x, t)
48 def call(self, x, t=1):
---> 49 return self._differentiate(x, t)

File ~/.local/lib/python3.10/site-packages/pysindy/differentiation/sindy_derivative.py:73, in SINDyDerivative._differentiate(self, x, t)
70 raise ValueError("t must be a positive constant or an array")
71 t = arange(x.shape[0]) * t
---> 73 return dxdt(x, t, axis=0, **self.kwargs)

File ~/.local/lib/python3.10/site-packages/derivative/differentiation.py:55, in dxdt(x, t, kind, axis, **kwargs)
53 else:
54 method = _gen_method(x, t, kind, axis, **kwargs)
---> 55 return method.d(x, t, axis=axis)

File ~/.local/lib/python3.10/site-packages/derivative/differentiation.py:194, in Derivative.d(self, X, t, axis)
176 def d(self, X, t, axis=1):
177 """
178 Compute the derivative of measurements X taken at times t.
179
(...)
192 ValueError: Requires that X.shape[axis] equals len(t). If X is flat, requires that len(X) equals len(t).
193 """
--> 194 X, flat = _align_axes(X, t, axis)
196 if X.shape[1] == 1:
197 dX = X

File ~/.local/lib/python3.10/site-packages/derivative/differentiation.py:245, in _align_axes(X, t, axis)
243 raise ValueError("Invalid axis.")
244 else:
--> 245 raise ValueError("Invalid shape of X.")
247 if X.shape[1] != len(t):
248 raise ValueError("Desired X axis size does not match t size.")

ValueError: Invalid shape of X.
`

[akaptano]

Given the error, what is the shape of u_flat?

[akaptano]

Also I recommend forget differentiation and use the WeakPDELibrary (see Example 12 jupyter notebook)

[Lawrence5]

Alan, Thanks for the reply. The shape of u_flat is (100,100,1). Below is the entire program (less extraneous stuff). As you can see, even without the noise in the training data, I get the same error any time I include the differentiation method when I instantiate the model. I seems that my syntax is wrong, but I don't know the correct way to specify it.

Yes, I understand that WeakPDELibrary is probably better, but I would like to know what is wrong with my syntax in any case. Thanks again.

`import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import Lasso
from sklearn.metrics import mean_squared_error
from scipy.integrate import solve_ivp
import pysindy as ps
​
rng = np.random.default_rng()
t0=0.01
tend=0.1
tnum=100
dt=(tend-t0)/tnum
x0=-0.2
xend=0.2
xnum=100
dx=(xend-x0)/xnum
​
t_train=np.arange(t0,tend,dt)
x_train=np.arange(x0,xend,dx)
u_train=np.zeros((xnum,tnum))
​
for i in range(xnum):
for j in range(tnum):
u_train[i][j] = 1/np.sqrt(t_train[j]) * np.exp(-1.0*x_train[i]*x_train[i]/(t_train[j]))

u_dot = ps.FiniteDifference(order=2, d=2, axis=0, is_uniform=True)._differentiate(u_train, dt)
ux = ps.FiniteDifference(d=1, axis=0, drop_endpoints=False)._differentiate(u_train, dx)
uxx = ps.SmoothedFiniteDifference(d=2, axis=0)._differentiate(u_train, dx)
​
u_flat = u_train.reshape(len(t_train), len(x_train), 1)
print(np.shape(t_train), np.shape(x_train), np.shape(u_train),np.shape(u_flat))

​library_functions = [lambda x: x_train]
library_function_names = [lambda x: x]
​
pde_lib = ps.PDELibrary(
library_functions=library_functions,
function_names=library_function_names,
derivative_order=2,
spatial_grid=x_train,
include_bias=False,
is_uniform=True
)

pde_lib.fit(u_flat)
print("2nd order derivative library: ")
print(pde_lib.get_feature_names())
print('STLSQ model: ')

LSoptimizer = ps.STLSQ(threshold=0.001, alpha=1e-5, normalize_columns=True)
SRoptimizer = ps.SR3(threshold=0.05, thresholder='l0', tol=1e-10, normalize_columns=True, max_iter=1000)
smoothed_fd = ps.SmoothedFiniteDifference(axis=1)._differentiate(u_train, t=dt)
spline_derivative = ps.SINDyDerivative(kind="spline", s=1e-2)
sg_filtered = ps.SINDyDerivative(kind="savitzky_golay", left=0.5, right=0.5, order=3)
trend_filtered = ps.SINDyDerivative(kind="trend_filtered", order=0, alpha=1e-2)
​
model = ps.SINDy(feature_library=pde_lib, optimizer=LSoptimizer)
model.fit(u_flat, t=dt)
model.print()

model = ps.SINDy(feature_library=pde_lib, optimizer=LSoptimizer,differentiation_method=smoothed_fd)
model.fit(u_flat,t=dt)
model.print()

​

(100,) (100,) (100, 100) (100, 100, 1)
2nd order derivative library:
['x0', 'x0_1', 'x0_11', 'x0x0_1', 'x0x0_11']
STLSQ model:
(x0)' = -0.106 x0 + 0.250 x0_11

---

TypeError Traceback (most recent call last)
Cell In[3], line 63
60 model.print()
62 model = ps.SINDy(feature_library=pde_lib, optimizer=LSoptimizer,differentiation_method=smoothed_fd)
---> 63 model.fit(u_flat,t=dt)
64 #model.fit(u_flat)
65 model.print()

`

[znicolaou]

Hey @Lawrence5 & @akaptano ,

I tried this out of curiosity, and it almost worked for my environment. I think it'll work for you if you install pysindy from the current github repository, as in

git clone https://github.com/dynamicslab/pysindy.git
cd pysindy
pip install -r requirements-dev.txt

The only change I needed to make to your code was

# smoothed_fd = ps.SmoothedFiniteDifference(axis=1)._differentiate(u_train, t=dt)
smoothed_fd = ps.SmoothedFiniteDifference(axis=1)

since the differentiation_method takes an object of type BaseDifferentiation rather than the actual derivative values. You can also use method.fit(u_flat,t=dt,x_dot=ps.SmoothedFiniteDifference(axis=1)._differentiate(u_flat, t=dt) instead, which does the same thing.

Note that the PDELibrary can also take a differentiation_method if you want the smoothed finite differences to apply to the spatial derivatives, since the SINDy class option only applies to the temporal derivatives. The syntax is a little different for the library case--it takes a class of type BaseDifferentiation rather than an object (sorry about the inconsistency). Use both options if you want both smooth spatial and temporal derivatives. Last point: if you want to use the is_uniform=True option with a smooth spatial derivativeoption, use diff_kwargs={'is_uniform':True} in the library declaration, as in

LSoptimizer = ps.STLSQ(threshold=0.001, alpha=1e-5, normalize_columns=True)
smoothed_spatial_fd = ps.SmoothedFiniteDifference
smoothed_fd = ps.SmoothedFiniteDifference(axis=1)

pde_lib = ps.PDELibrary(
library_functions=library_functions,
function_names=library_function_names,
derivative_order=2,
spatial_grid=x_train,
include_bias=False,
differentiation_method=smoothed_spatial_fd,
diff_kwargs={'is_uniform':True}
)
model = ps.SINDy(feature_library=pde_lib, optimizer=LSoptimizer,differentiation_method=smoothed_fd)
model.fit(u_flat,t=dt)
model.print()

[Lawrence5]

Zachary,
Thank you so much. That helps a lot. As I suspected, the error message was leading me to look in the wrong place, but it makes sense now.
Lawrence

[Lawrence5]

Oh and thanks for the tip on the spatial derivatives. I would never have figured that out on my own. You saved me several days no doubt!