SINDY-PI LinAlgError: Singular matrix #620

balcieb · 2025-05-08T02:25:34Z

balcieb
May 8, 2025

Hi,
I have 3D simulation results for a flame problem. I am trying to obtain a model for the mixture fraction (Z) field multiplied by the progress variable (Y) field using SINDY-PI. The results include 186 timesteps for a spatial grid of (150, 30, 100) in (z,y,x). I attempted to follow the example in 9_sindypi_with_sympy, but I am encountering a Singular matrix error. I believe this is related to the candidate library function I am creating, and I assume that the matrix formed from this library is ill-conditioned. I tried to modify the library but keep getting the same error. Can you help me solve this problem, please?
The code is below:

import numpy as np
import matplotlib.pyplot as plt
# import pandas as pd
import vtkmodules.all as vtk
from vtkmodules.util.numpy_support import vtk_to_numpy
import matplotlib.cm as cm
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

## OBTAIN Z (MIXTURE FRACTION)


Y_arr_BC_Z = np.load("Y_arr_BC_Z.npy")
print("Y_arr_BC_Z.shape =", Y_arr_BC_Z.shape)

## read time_arr from txt file and assign it to time_arr_v2:
time_arr = np.loadtxt('time_array.txt', delimiter=',')
time_arr = np.array(time_arr)
x = np.loadtxt('x_array.txt', delimiter=',')
x = np.array(x)
## READ y_array from txt file and assign it to y_array:
y = np.loadtxt('y_array.txt', delimiter=',')
y = np.array(y)
## READ z_array from txt file and assign it to z_array:
z = np.loadtxt('z_array.txt', delimiter=',')
z = np.array(z)

total_time = time_arr[-1] - time_arr[0]
print("total time = ", total_time)



Y_arr_BC_Z = np.reshape(Y_arr_BC_Z,(len(z), len(y), len(x), len(time_arr)))
print("Y_arr_BC_Z.shape =", Y_arr_BC_Z.shape)


## OBTAIN PROGRESS VARIABLE (Y)

Y_arr_BC_Y = np.load("Y_arr_BC_Y.npy")
print("Y_arr_BC_Y.shape =", Y_arr_BC_Y.shape)


Y_arr_BC_Y = np.reshape(Y_arr_BC_Y,(len(z), len(y), len(x), len(time_arr)))
print("Y_arr_BC_Y.shape = ", Y_arr_BC_Y.shape)


## SINDY MODEL


import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from sklearn.linear_model import Lasso
from scipy.io import loadmat
from sklearn.metrics import mean_squared_error
from scipy.integrate import solve_ivp

import pysindy as ps

# Ignore matplotlib deprecation warnings
import warnings
warnings.filterwarnings("ignore", category=UserWarning)

# Seed the random number generators for reproducibility
np.random.seed(100)

integrator_keywords = {}
integrator_keywords['rtol'] = 1e-12
integrator_keywords['method'] = 'LSODA'
integrator_keywords['atol'] = 1e-12



dt = time_arr[1] - time_arr[0]
dx = x[1] - x[0]
dy = y[1] - y[0]
dz = z[1] - z[0]
print("dx = ", dx)
print("dy = ", dy)
print("dz = ", dz)
print("dt = ", dt)
# print("time_arr =", time_arr)
time_arr = time_arr - time_arr[0]
print("min time_arr=", np.min(time_arr))
print("max time_arr=", np.max(time_arr))
print("len(time_arr)=", len(time_arr))


max_time_arr = np.max(time_arr)
mean_dt = max_time_arr/len(time_arr)
print("mean_dt = ", mean_dt)


dt = mean_dt


# Compute u_t from generated solution
# u = np.zeros((len(z), len(y[y_start:y_end]), len(x[x_start:x_end]), len(time_arr), 2))
u = np.zeros((len(z), len(y[y_start:y_end]), len(x[x_start:x_end]), len(time_arr), 1))
print("u.shape=", u.shape)


# Assuming Y_arr_BC_Z and Y_arr_BC_Y are already defined and have the same shape
# Compute element-wise product for each (x, y, t)
Y_arr_BC_ZY = Y_arr_BC_Z * Y_arr_BC_Y

# Optional: print shape to confirm
print("Y_arr_BC_ZY.shape =", Y_arr_BC_ZY.shape)
print("min Y_arr_BC_ZY =", np.min(Y_arr_BC_ZY))
print("max Y_arr_BC_ZY =", np.max(Y_arr_BC_ZY))


yslice = 75
y_start = 60
y_end = 90
x_start = 25
x_end = 125

# u[:, :, :, 0] = Y_arr_BC_Z[0:m, :, :]
u[:, :, :, :, 0] = Y_arr_BC_ZY[:, y_start:y_end , x_start:x_end, :]
# u[:, :, :, :, 1] = Y_arr_BC_Y[:, y_start:y_end , x_start:x_end, :]
# u_dot = ps.FiniteDifference(axis=3)._differentiate(u, dt)
u_dot = ps.SmoothedFiniteDifference(axis=3,drop_endpoints=False, smoother_kws={'window_length': 30})._differentiate(u, dt)
print("u.shape=", u.shape)
print("u_dot.shape=", u_dot.shape)

## print an error if u_dot incldes NaN
if np.isnan(u_dot).any():
    print("u_dot contains NaN values")
else:
    print("u_dot does not contain NaN values")
## print an error if u_dot incldes inf
if np.isinf(u_dot).any():
    print("u_dot contains inf values")
else:
    print("u_dot does not contain inf values")



# x = x - np.min(x)
x = x - x[0]
y = y - y[0]
print("x min = ", np.min(x))
print("x max = ", np.max(x))
print("y min = ", np.min(y))
print("y max = ", np.max(y))


print("y[y_start:y_end] min = ", np.min(y[y_start:y_end]))
print("y[y_start:y_end] max = ", np.max(y[y_start:y_end]))
X, Y, Z = np.meshgrid(x[x_start:x_end], y[y_start:y_end], z)
print("X.shape = ", X.shape)
print("Y.shape = ", Y.shape)
print("Z.shape = ", Z.shape)



# Choose 60 % of data for training because data is big... 
# can only randomly subsample if you are passing u_dot to model.fit!!!
training_amount = 0.6 # 0.6
train = int(len(t) * training_amount)
print("train =", train)
test = int(len(t)-train)
print("test=", test)
print("len(t)=", len(t))
print("train+test", train+test)
u_train = u[:, :, :, 0:train, :]
print("u_train.shape=", u_train.shape)

u_test = u[:, :, :, train:t_len, :]
print("u_test.shape=", u_test.shape)

u_dot_train = u_dot[:, :, :, 0:train, :]
print("u_dot_train.shape=", u_dot_train.shape)

u_dot_test = u_dot[:, :, :, train:t_len, :]
print("u_dot_test.shape=", u_dot_test.shape)

t_train = t[0:train]
print("t_train.shape=", t_train.shape)

t_test = t[train:t_len]
print("t_test.shape=", t_test.shape)

spatial_grid = np.asarray([X, Y, Z])
spatial_grid = np.transpose(spatial_grid, [3, 1, 2, 0])
print("spatial_grid.shape =", spatial_grid.shape)



## print the shape of the u_train, u_dot_train, t_train:
print("u_train.shape:", u_train.shape)
print("u_dot_train.shape:", u_dot_train.shape)
print("t_train.shape:", t_train.shape)
t_train = np.sort(t_train)
# print("t_train=", t_train)
t_test = np.sort(t_test)
# print("t_test=", t_test)



X_2, Y_2, Z_2, T = np.meshgrid(
	x[x_start:x_end], y[y_start:y_end], z, t_train, indexing="ij"
)
# print("X_2.shape=", X_2.shape)
# print("Y_2.shape=", Y_2.shape)
# print("Z_2.shape=", Z_2.shape)
# print("T.shape=", T.shape)
X_2 = np.transpose(X_2, [2, 1, 0, 3])
print("X_2.shape=", X_2.shape)
Y_2 = np.transpose(Y_2, [2, 1, 0, 3])
print("Y_2.shape=", Y_2.shape)
Z_2 = np.transpose(Z_2, [2, 1, 0, 3])
print("Z_2.shape=", Z_2.shape)
T = np.transpose(T, [2, 1, 0, 3])
print("T.shape=", T.shape)

spatiotemporal = np.asarray([X_2, Y_2, Z_2, T])
# spatial_grid = np.transpose(spatial_grid, [2, 0, 3, 1])
# ## make spatiotemporal shape 150,30,100,111,4
spatiotemporal = np.transpose(spatiotemporal, [1, 2, 3, 4, 0])
print("spatiotemporal.shape =", spatiotemporal.shape)

## MODEL WITH PDE LIBRARY

library_functions = [lambda x: x, lambda x: x * x, lambda x: x * x * x]
library_function_names = [lambda x: x, lambda x: x + x, lambda x: x + x + x]
np.random.seed(1)

pde_lib = ps.PDELibrary(
    library_functions=library_functions,
    function_names=library_function_names,
    derivative_order=3,
    spatial_grid=spatial_grid,
    temporal_grid=t_train,
    include_bias=False,
    is_uniform=True,
    implicit_terms=True,
)

sindy_opt = ps.SINDyPI(
    threshold=1e-1, max_iter=1000, tol=1e-10,
    thresholder='l1', normalize_columns=False
)
model = ps.SINDy(feature_library=pde_lib, optimizer=sindy_opt)
model.fit(u_train, t=dt)
model.print()

LinAlgError                               Traceback (most recent call last)
Cell In[90], [line 21](vscode-notebook-cell:?execution_count=90&line=21)
     [16](vscode-notebook-cell:?execution_count=90&line=16) sindy_opt = ps.SINDyPI(
     [17](vscode-notebook-cell:?execution_count=90&line=17)     threshold=1e-1, max_iter=1000, tol=1e-10,
     [18](vscode-notebook-cell:?execution_count=90&line=18)     thresholder='l1', normalize_columns=False
     [19](vscode-notebook-cell:?execution_count=90&line=19) )
     [20](vscode-notebook-cell:?execution_count=90&line=20) model = ps.SINDy(feature_library=pde_lib, optimizer=sindy_opt)
---> [21](vscode-notebook-cell:?execution_count=90&line=21) model.fit(u_train, t=dt)
     [22](vscode-notebook-cell:?execution_count=90&line=22) model.print()

File /opt/homebrew/anaconda3/lib/python3.11/site-packages/pysindy/pysindy.py:414, 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)
    [412](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/pysindy/pysindy.py:412)     warnings.filterwarnings(action, category=LinAlgWarning)
    [413](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/pysindy/pysindy.py:413)     warnings.filterwarnings(action, category=UserWarning)
--> [414](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/pysindy/pysindy.py:414)     self.model.fit(x, x_dot)
    [416](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/pysindy/pysindy.py:416) # New version of sklearn changes attribute name
    [417](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/pysindy/pysindy.py:417) if float(__version__[:3]) >= 1.0:

File /opt/homebrew/anaconda3/lib/python3.11/site-packages/sklearn/base.py:1151, in _fit_context.<locals>.decorator.<locals>.wrapper(estimator, *args, **kwargs)
   [1144](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/sklearn/base.py:1144)     estimator._validate_params()
   [1146](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/sklearn/base.py:1146) with config_context(
   [1147](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/sklearn/base.py:1147)     skip_parameter_validation=(
   [1148](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/sklearn/base.py:1148)         prefer_skip_nested_validation or global_skip_validation
   [1149](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/sklearn/base.py:1149)     )
   [1150](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/sklearn/base.py:1150) ):
...
File /opt/homebrew/anaconda3/lib/python3.11/site-packages/numpy/linalg/linalg.py:89, in _raise_linalgerror_singular(err, flag)
     [88](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/numpy/linalg/linalg.py:88) def _raise_linalgerror_singular(err, flag):
---> [89](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/numpy/linalg/linalg.py:89)     raise LinAlgError("Singular matrix")

LinAlgError: Singular matrix



## MODEL WITH WEAKPDE LIBRARY:


sindy_opt = ps.SINDyPI(
    threshold=1e-2, max_iter=1000, tol=1e-14,
    thresholder='l1', normalize_columns=False
)

pde_lib = ps.WeakPDELibrary(library_functions=library_functions, 
                        function_names=library_function_names, 
                        derivative_order=2, spatiotemporal_grid=spatiotemporal,
                        include_bias=False, is_uniform=True,implicit_terms=True)

model = ps.SINDy(feature_library=pde_lib, optimizer=sindy_opt)
model.fit(u_train, t=dt)
model.print()

LinAlgError                               Traceback (most recent call last)
Cell In[88], [line 12](vscode-notebook-cell:?execution_count=88&line=12)
      [6](vscode-notebook-cell:?execution_count=88&line=6) pde_lib = ps.WeakPDELibrary(library_functions=library_functions, 
      [7](vscode-notebook-cell:?execution_count=88&line=7)                         function_names=library_function_names, 
      [8](vscode-notebook-cell:?execution_count=88&line=8)                         derivative_order=2, spatiotemporal_grid=spatiotemporal,
      [9](vscode-notebook-cell:?execution_count=88&line=9)                         include_bias=False, is_uniform=True,implicit_terms=True)
     [11](vscode-notebook-cell:?execution_count=88&line=11) model = ps.SINDy(feature_library=pde_lib, optimizer=sindy_opt)
---> [12](vscode-notebook-cell:?execution_count=88&line=12) model.fit(u_train, t=dt)
     [13](vscode-notebook-cell:?execution_count=88&line=13) model.print()

File /opt/homebrew/anaconda3/lib/python3.11/site-packages/pysindy/pysindy.py:414, 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)
    [412](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/pysindy/pysindy.py:412)     warnings.filterwarnings(action, category=LinAlgWarning)
    [413](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/pysindy/pysindy.py:413)     warnings.filterwarnings(action, category=UserWarning)
--> [414](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/pysindy/pysindy.py:414)     self.model.fit(x, x_dot)
    [416](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/pysindy/pysindy.py:416) # New version of sklearn changes attribute name
    [417](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/pysindy/pysindy.py:417) if float(__version__[:3]) >= 1.0:

File /opt/homebrew/anaconda3/lib/python3.11/site-packages/sklearn/base.py:1151, in _fit_context.<locals>.decorator.<locals>.wrapper(estimator, *args, **kwargs)
   [1144](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/sklearn/base.py:1144)     estimator._validate_params()
   [1146](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/sklearn/base.py:1146) with config_context(
   [1147](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/sklearn/base.py:1147)     skip_parameter_validation=(
   [1148](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/sklearn/base.py:1148)         prefer_skip_nested_validation or global_skip_validation
   [1149](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/sklearn/base.py:1149)     )
   [1150](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/sklearn/base.py:1150) ):
...
File /opt/homebrew/anaconda3/lib/python3.11/site-packages/numpy/linalg/linalg.py:89, in _raise_linalgerror_singular(err, flag)
     [88](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/numpy/linalg/linalg.py:88) def _raise_linalgerror_singular(err, flag):
---> [89](https://file+.vscode-resource.vscode-cdn.net/opt/homebrew/anaconda3/lib/python3.11/site-packages/numpy/linalg/linalg.py:89)     raise LinAlgError("Singular matrix")

LinAlgError: Singular matrix

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

SINDY-PI LinAlgError: Singular matrix #620

Uh oh!

{{title}}

Uh oh!

Replies: 0 comments

Select a reply

Uh oh!

SINDY-PI LinAlgError: Singular matrix #620

Uh oh!

balcieb May 8, 2025

Replies: 0 comments

balcieb
May 8, 2025