[BUG] Finite difference shape mismatch

[NolanBrb]

When trying to run the basic examples found at https://pysindy.readthedocs.io/en/latest/examples/3_original_paper/example.html , I run into multiple errors. It seems that there are some issues with the finite difference code, with numpy functions called with arrays of the wrong shapes.

A first error occurs when np.linalg.solve(matrices, b) is called to compute finite differences coefficients. I fixed it by changing b for its transpose b.T. However, I run into a second error in the np.einsum() function.

Reproducing code example:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.cm import rainbow
import numpy as np
from scipy.integrate import solve_ivp
from scipy.io import loadmat
from pysindy.utils import linear_damped_SHO
from pysindy.utils import cubic_damped_SHO
from pysindy.utils import linear_3D
from pysindy.utils import hopf
from pysindy.utils import lorenz

import pysindy as ps

# ignore user warnings
import warnings
warnings.filterwarnings("ignore", category=UserWarning)

np.random.seed(1000)  # Seed for reproducibility

# Integrator keywords for solve_ivp
integrator_keywords = {}
integrator_keywords['rtol'] = 1e-12
integrator_keywords['method'] = 'LSODA'
integrator_keywords['atol'] = 1e-12

# Generate training data

dt = 0.01
t_train = np.arange(0, 25, dt)
t_train_span = (t_train[0], t_train[-1])
x0_train = [2, 0]
x_train = solve_ivp(linear_damped_SHO, t_train_span,
                    x0_train, t_eval=t_train, **integrator_keywords).y.T

# Fit the model

poly_order = 5
threshold = 0.05

model = ps.SINDy(
    optimizer=ps.STLSQ(threshold=threshold),
    feature_library=ps.PolynomialLibrary(degree=poly_order),
)
model.fit(x_train, t=dt)
model.print()

Error message #1

---

ValueError Traceback (most recent call last)
Cell In[1], line 45
39 threshold = 0.05
41 model = ps.SINDy(
42 optimizer=ps.STLSQ(threshold=threshold),
43 feature_library=ps.PolynomialLibrary(degree=poly_order),
44 )
---> 45 model.fit(x_train, t=dt)
46 model.print()

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\pysindy.py:246, in SINDy.fit(self, x, t, x_dot, u)
240 u = validate_control_variables(
241 x,
242 u,
243 trim_last_point=(self.discrete_time and x_dot is None),
244 )
245 self.n_control_features_ = u[0].shape[u[0].ax_coord]
--> 246 x, x_dot = self._process_trajectories(x, t, x_dot)
248 # Append control variables
249 if u is not None:

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\pysindy.py:491, in SINDy._process_trajectories(self, x, t, x_dot)
488 x = [xi[:-1] for xi in x]
489 else:
490 x, x_dot = zip(
--> 491 *[
492 self.feature_library.calc_trajectory(
493 self.differentiation_method, xi, ti
494 )
495 for xi, ti in _zip_like_sequence(x, t)
496 ]
497 )
498 return x, x_dot

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\pysindy.py:492, in (.0)
488 x = [xi[:-1] for xi in x]
489 else:
490 x, x_dot = zip(
491 *[
--> 492 self.feature_library.calc_trajectory(
493 self.differentiation_method, xi, ti
494 )
495 for xi, ti in _zip_like_sequence(x, t)
496 ]
497 )
498 return x, x_dot

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\feature_library\base.py:66, in BaseFeatureLibrary.calc_trajectory(self, diff_method, x, t)
65 def calc_trajectory(self, diff_method, x, t):
---> 66 x_dot = diff_method(x, t=t)
67 x = AxesArray(diff_method.smoothed_x_, x.axes)
68 return x, AxesArray(x_dot, x.axes)

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\differentiation\base.py:53, in BaseDifferentiation.call(self, x, t)
52 def call(self, x, t=1):
---> 53 return self._differentiate(x, t)

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\differentiation\finite_difference.py:278, in FiniteDifference._differentiate(self, x, t)
275 if not self.drop_endpoints:
276 # Forward differences on boundary
277 if not self.periodic:
--> 278 coeffs = self._coefficients_boundary_forward(t)
279 boundary = self._accumulate(coeffs, x)
281 if self.order % 2 == 0:

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\differentiation\finite_difference.py:148, in FiniteDifference._coefficients_boundary_forward(self, t)
146 b = np.zeros(self.stencil_inds.shape).T
147 b[:, self.d] = factorial(self.d)
--> 148 return np.linalg.solve(matrices, b)

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\numpy\linalg_linalg.py:413, in solve(a, b)
410 signature = 'DD->D' if isComplexType(t) else 'dd->d'
411 with errstate(call=_raise_linalgerror_singular, invalid='call',
412 over='ignore', divide='ignore', under='ignore'):
--> 413 r = gufunc(a, b, signature=signature)
415 return wrap(r.astype(result_t, copy=False))

ValueError: solve: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (m,m),(m,n)->(m,n) (size 2 is different from 3)

Fix for error #1 : replace all lines with np.linalg.solve(matrices, b) with np.linalg.solve(matrices, b.T)

Error message #2

AFTER replacing all np.linalg.solve(matrices, b) with np.linalg.solve(matrices, b.T), I get the following error

---

ValueError Traceback (most recent call last)
Cell In[13], line 45
39 threshold = 0.05
41 model = ps.SINDy(
42 optimizer=ps.STLSQ(threshold=threshold),
43 feature_library=ps.PolynomialLibrary(degree=poly_order),
44 )
---> 45 model.fit(x_train, t=dt)
46 model.print()

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\pysindy.py:246, in SINDy.fit(self, x, t, x_dot, u)
240 u = validate_control_variables(
241 x,
242 u,
243 trim_last_point=(self.discrete_time and x_dot is None),
244 )
245 self.n_control_features_ = u[0].shape[u[0].ax_coord]
--> 246 x, x_dot = self._process_trajectories(x, t, x_dot)
248 # Append control variables
249 if u is not None:

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\pysindy.py:491, in SINDy._process_trajectories(self, x, t, x_dot)
488 x = [xi[:-1] for xi in x]
489 else:
490 x, x_dot = zip(
--> 491 *[
492 self.feature_library.calc_trajectory(
493 self.differentiation_method, xi, ti
494 )
495 for xi, ti in _zip_like_sequence(x, t)
496 ]
497 )
498 return x, x_dot

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\pysindy.py:492, in (.0)
488 x = [xi[:-1] for xi in x]
489 else:
490 x, x_dot = zip(
491 *[
--> 492 self.feature_library.calc_trajectory(
493 self.differentiation_method, xi, ti
494 )
495 for xi, ti in _zip_like_sequence(x, t)
496 ]
497 )
498 return x, x_dot

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\feature_library\base.py:66, in BaseFeatureLibrary.calc_trajectory(self, diff_method, x, t)
65 def calc_trajectory(self, diff_method, x, t):
---> 66 x_dot = diff_method(x, t=t)
67 x = AxesArray(diff_method.smoothed_x_, x.axes)
68 return x, AxesArray(x_dot, x.axes)

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\differentiation\base.py:53, in BaseDifferentiation.call(self, x, t)
52 def call(self, x, t=1):
---> 53 return self._differentiate(x, t)

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\differentiation\finite_difference.py:279, in FiniteDifference._differentiate(self, x, t)
277 if not self.periodic:
278 coeffs = self._coefficients_boundary_forward(t)
--> 279 boundary = self._accumulate(coeffs, x)
281 if self.order % 2 == 0:
282 right_len = (self.n_stencil - 1) // 2

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\pysindy\differentiation\finite_difference.py:227, in FiniteDifference._accumulate(self, coeffs, x)
222 x = AxesArray(x, {"ax_unk": list(range(x.ndim))})
223 x_expanded = AxesArray(
224 np.transpose(x[tuple(s)], axes=trans), x.insert_axis(0, "ax_offset")
225 )
226 return np.transpose(
--> 227 np.einsum(
228 "ij...,ij->j...",
229 x_expanded,
230 np.transpose(coeffs),
231 ),
232 np.roll(np.arange(x.ndim), self.axis),
233 )

File c:\Users\nolan\miniconda3\envs\numerical-physics\Lib\site-packages\numpy_core\einsumfunc.py:1429, in einsum(out, optimize, *operands, **kwargs)
1427 if specified_out:
1428 kwargs['out'] = out
-> 1429 return c_einsum(*operands, **kwargs)
1431 # Check the kwargs to avoid a more cryptic error later, without having to
1432 # repeat default values here
1433 valid_einsum_kwargs = ['dtype', 'order', 'casting']

ValueError: operand has more dimensions than subscripts given in einstein sum, but no '...' ellipsis provided to broadcast the extra dimensions.

PySINDy/Python version information:

1.7.6.dev465+g69e51e8 3.11.2 | packaged by conda-forge | (main, Mar 31 2023, 17:45:57) [MSC v.1934 64 bit (AMD64)]

[chaitanya94]

I had the same problem. I tried both the stable version 1.7.5 available from pip as well as compiling from source. The bug is present in both.

An easy way to reproduce the error is:

import numpy as np
import pysindy as ps

t = np.linspace(0, 1, 100)
x = 3 * np.exp(-2 * t)
y = 0.5 * np.exp(t)
X = np.stack((x, y), axis=-1)

model = ps.SINDy(feature_names=["x", "y"])
model.fit(X, t=t)

Solution
For pysindy version 1.7.5 the issue is resolved when I downgrade my numpy to version 1.26.4.

[Jacob-Stevens-Haas]

Thanks @chaitanya94! @NolanBrb, can you report your numpy version?