pysindy.optimizers package¶

Submodules¶

pysindy.optimizers.base module¶

Base class for SINDy optimizers.

class pysindy.optimizers.base.ComplexityMixin[source]¶

Bases: object

property complexity¶

class pysindy.optimizers.base.BaseOptimizer(max_iter=20, normalize_columns=False, fit_intercept=False, initial_guess=None, copy_X=True)[source]¶

Bases: LinearRegression, ComplexityMixin

Base class for SINDy optimizers. Subclasses must implement a _reduce method for carrying out the bulk of the work of fitting a model.

Parameters

fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed.
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
initial_guess (np.ndarray, shape (n_features,) or (n_targets, n_features),) – optional (default None) Initial guess for coefficients coef_. If None, the initial guess is obtained via a least-squares fit.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s).
ind_ (array, shape (n_features,) or (n_targets, n_features)) – Array of 0s and 1s indicating which coefficients of the weight vector have not been masked out.
history_ (list) – History of coef_ over iterations of the optimization algorithm.
Theta_ (np.ndarray, shape (n_samples, n_features)) – The Theta matrix to be used in the optimization. We save it as an attribute because access to the full library of terms is sometimes needed for various applications.

fit(x_, y, sample_weight=None, **reduce_kws)[source]¶

Fit the model.

Parameters

x (array-like, shape (n_samples, n_features)) – Training data
y (array-like, shape (n_samples,) or (n_samples, n_targets)) – Target values
sample_weight (float or numpy array of shape (n_samples,), optional) – Individual weights for each sample
reduce_kws (dict) – Optional keyword arguments to pass to the _reduce method (implemented by subclasses)

Returns

self

Return type

returns an instance of self

class pysindy.optimizers.base.EnsembleOptimizer(opt: BaseOptimizer, bagging: bool = False, library_ensemble: bool = False, n_models: int = 20, n_subset: Optional[int] = None, n_candidates_to_drop: int = 1, replace: bool = True, ensemble_aggregator: Optional[Callable] = None)[source]¶

Bases: BaseOptimizer

Wrapper class for ensembling methods.

Parameters

opt (BaseOptimizer) – The underlying optimizer to run on each ensemble
bagging (boolean, optional (default False)) – This parameter is used to allow for “ensembling”, i.e. the generation of many SINDy models (n_models) by choosing a random temporal subset of the input data (n_subset) for each sparse regression. This often improves robustness because averages (bagging) or medians (bragging) of all the models are usually quite high-performing. The user can also generate “distributions” of many models, and calculate how often certain library terms are included in a model.
library_ensemble (boolean, optional (default False)) – This parameter is used to allow for “library ensembling”, i.e. the generation of many SINDy models (n_models) by choosing a random subset of the candidate library terms to truncate. So, n_models are generated by solving n_models sparse regression problems on these “reduced” libraries. Once again, this often improves robustness because averages (bagging) or medians (bragging) of all the models are usually quite high-performing. The user can also generate “distributions” of many models, and calculate how often certain library terms are included in a model.
n_models (int, optional (default 20)) – Number of models to generate via ensemble
n_subset (int, optional (default len(time base))) – Number of time points to use for ensemble
n_candidates_to_drop (int, optional (default 1)) – Number of candidate terms in the feature library to drop during library ensembling.
replace (boolean, optional (default True)) – If ensemble true, whether or not to time sample with replacement.
ensemble_aggregator (callable, optional (default numpy.median)) – Method to aggregate model coefficients across different samples. This method argument is only used if ensemble or library_ensemble is True. The method should take in a list of 2D arrays and return a 2D array of the same shape as the arrays in the list. Example: lambda x: np.median(x, axis=0)

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Regularized weight vector(s). This is the v in the objective function.
coef_full_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s) that are not subjected to the regularization. This is the w in the objective function.
unbias (boolean) – Whether to perform an extra step of unregularized linear regression to unbias the coefficients for the identified support. unbias is automatically set to False if a constraint is used and is otherwise left uninitialized.

pysindy.optimizers.constrained_sr3 module¶

class pysindy.optimizers.constrained_sr3.ConstrainedSR3(threshold=0.1, nu=1.0, tol=1e-05, thresholder='l0', max_iter=30, trimming_fraction=0.0, trimming_step_size=1.0, constraint_lhs=None, constraint_rhs=None, constraint_order='target', normalize_columns=False, fit_intercept=False, copy_X=True, initial_guess=None, thresholds=None, equality_constraints=False, inequality_constraints=False, constraint_separation_index=0, verbose=False, verbose_cvxpy=False)[source]¶

Bases: SR3

Sparse relaxed regularized regression with linear (in)equality constraints.

Attempts to minimize the objective function

\[0.5\|y-Xw\|^2_2 + \lambda R(u) + (0.5 / \nu)\|w-u\|^2_2\]

\[\text{subject to } Cw = d\]

over u and w, where \(R(u)\) is a regularization function, C is a constraint matrix, and d is a vector of values. See the following reference for more details:

Champion, Kathleen, et al. “A unified sparse optimization framework to learn parsimonious physics-informed models from data.” IEEE Access 8 (2020): 169259-169271.

Parameters

threshold (float, optional (default 0.1)) – Determines the strength of the regularization. When the regularization function R is the l0 norm, the regularization is equivalent to performing hard thresholding, and lambda is chosen to threshold at the value given by this parameter. This is equivalent to choosing lambda = threshold^2 / (2 * nu).
nu (float, optional (default 1)) – Determines the level of relaxation. Decreasing nu encourages w and v to be close, whereas increasing nu allows the regularized coefficients v to be farther from w.
tol (float, optional (default 1e-5)) – Tolerance used for determining convergence of the optimization algorithm.
thresholder (string, optional (default 'l0')) – Regularization function to use. Currently implemented options are ‘l0’ (l0 norm), ‘l1’ (l1 norm), ‘l2’ (l2 norm), ‘cad’ (clipped absolute deviation), ‘weighted_l0’ (weighted l0 norm), ‘weighted_l1’ (weighted l1 norm), and ‘weighted_l2’ (weighted l2 norm).
max_iter (int, optional (default 30)) – Maximum iterations of the optimization algorithm.
fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
constraint_lhs (numpy ndarray, optional (default None)) – Shape should be (n_constraints, n_features * n_targets), The left hand side matrix C of Cw <= d. There should be one row per constraint.
constraint_rhs (numpy ndarray, shape (n_constraints,), optional (default None)) – The right hand side vector d of Cw <= d.
constraint_order (string, optional (default "target")) – The format in which the constraints constraint_lhs were passed. Must be one of “target” or “feature”. “target” indicates that the constraints are grouped by target: i.e. the first n_features columns correspond to constraint coefficients on the library features for the first target (variable), the next n_features columns to the library features for the second target (variable), and so on. “feature” indicates that the constraints are grouped by library feature: the first n_targets columns correspond to the first library feature, the next n_targets columns to the second library feature, and so on.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed. Note that this parameter is incompatible with the constraints!
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
initial_guess (np.ndarray, optional (default None)) – Shape should be (n_features) or (n_targets, n_features). Initial guess for coefficients coef_, (v in the mathematical equations) If None, least-squares is used to obtain an initial guess.
thresholds (np.ndarray, shape (n_targets, n_features), optional (default None)) – Array of thresholds for each library function coefficient. Each row corresponds to a measurement variable and each column to a function from the feature library. Recall that SINDy seeks a matrix \(\Xi\) such that \(\dot{X} \approx \Theta(X)\Xi\). thresholds[i, j] should specify the threshold to be used for the (j + 1, i + 1) entry of \(\Xi\). That is to say it should give the threshold to be used for the (j + 1)st library function in the equation for the (i + 1)st measurement variable.
inequality_constraints (bool, optional (default False)) – If True, CVXPY methods are used to solve the problem.
verbose (bool, optional (default False)) – If True, prints out the different error terms every max_iter / 10 iterations.
verbose_cvxpy (bool, optional (default False)) – Boolean flag which is passed to CVXPY solve function to indicate if output should be verbose or not. Only relevant for optimizers that use the CVXPY package in some capabity.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Regularized weight vector(s). This is the v in the objective function.
coef_full_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s) that are not subjected to the regularization. This is the w in the objective function.
unbias (boolean) – Whether to perform an extra step of unregularized linear regression to unbias the coefficients for the identified support. unbias is automatically set to False if a constraint is used and is otherwise left uninitialized.

pysindy.optimizers.frols module¶

class pysindy.optimizers.frols.FROLS(normalize_columns=False, fit_intercept=False, copy_X=True, kappa=None, max_iter=10, alpha=0.05, ridge_kw=None, verbose=False)[source]¶

Bases: BaseOptimizer

Forward Regression Orthogonal Least-Squares (FROLS) optimizer.

Attempts to minimize the objective function \(\|y - Xw\|^2_2 + \alpha \|w\|^2_2\) by iteractively selecting the most correlated function in the library. This is a greedy algorithm.

See the following reference for more details:

Billings, Stephen A. Nonlinear system identification: NARMAX methods in the time, frequency, and spatio-temporal domains. John Wiley & Sons, 2013.

Parameters

fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed.
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
kappa (float, optional (default None)) – If passed, compute the MSE errors with an extra L0 term with strength equal to kappa times the condition number of Theta.
max_iter (int, optional (default 10)) – Maximum iterations of the optimization algorithm. This determines the number of nonzero terms chosen by the FROLS algorithm.
alpha (float, optional (default 0.05)) – Optional L2 (ridge) regularization on the weight vector.
ridge_kw (dict, optional (default None)) – Optional keyword arguments to pass to the ridge regression.
verbose (bool, optional (default False)) – If True, prints out the different error terms every iteration.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s).
history_ (list) – History of coef_. history_[k] contains the values of coef_ at iteration k of FROLS.

Examples

>>> import numpy as np
>>> from scipy.integrate import odeint
>>> from pysindy import SINDy
>>> from pysindy.optimizers import FROLS
>>> lorenz = lambda z,t : [10 * (z[1] - z[0]),
>>>                        z[0] * (28 - z[2]) - z[1],
>>>                        z[0] * z[1] - 8 / 3 * z[2]]
>>> t = np.arange(0, 2, .002)
>>> x = odeint(lorenz, [-8, 8, 27], t)
>>> opt = FROLS(threshold=.1, alpha=.5)
>>> model = SINDy(optimizer=opt)
>>> model.fit(x, t=t[1] - t[0])
>>> model.print()
x0' = -9.999 1 + 9.999 x0
x1' = 27.984 1 + -0.996 x0 + -1.000 1 x1
x2' = -2.666 x1 + 1.000 1 x0

pysindy.optimizers.miosr module¶

class pysindy.optimizers.miosr.MIOSR(target_sparsity=5, group_sparsity=None, alpha=0.01, regression_timeout=10, fit_intercept=False, constraint_lhs=None, constraint_rhs=None, constraint_order='target', normalize_columns=False, copy_X=True, initial_guess=None, verbose=False)[source]¶

Bases: BaseOptimizer

Mixed-Integer Optimized Sparse Regression.

Solves the sparsity constrained regression problem to provable optimality .. math:

\|y-Xw\|^2_2 + \lambda R(u)

\[\text{subject to } \|w\|_0 \leq k\]

by using type-1 specially ordered sets (SOS1) to encode the support of the coefficients. Can optionally add additional constraints on the coefficients or access the gurobi model directly for advanced usage. See the following reference for additional details:

Bertsimas, D. and Gurnee, W., 2022. Learning Sparse Nonlinear Dynamics via Mixed-Integer Optimization. arXiv preprint arXiv:2206.00176.

Parameters

target_sparsity (int, optional (default 5)) – The maximum number of nonzero coefficients across all dimensions. If set, the model will fit all dimensions jointly, potentially reducing statistical efficiency.
group_sparsity (int tuple, optional (default None)) – Tuple of length n_targets constraining the number of nonzero coefficients for each target dimension.
alpha (float, optional (default 0.01)) – Optional L2 (ridge) regularization on the weight vector.
regression_timeout (int, optional (default 10)) – The timeout (in seconds) of the gurobi optimizer to solve and prove optimality (either per dimension or jointly depending on the above sparsity settings).
fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
constraint_lhs (numpy ndarray, optional (default None)) – Shape should be (n_constraints, n_features * n_targets), The left hand side matrix C of Cw <= d. There should be one row per constraint.
constraint_rhs (numpy ndarray, shape (n_constraints,), optional (default None)) – The right hand side vector d of Cw <= d.
constraint_order (string, optional (default "target")) – The format in which the constraints constraint_lhs were passed. Must be one of “target” or “feature”. “target” indicates that the constraints are grouped by target: i.e. the first n_features columns correspond to constraint coefficients on the library features for the first target (variable), the next n_features columns to the library features for the second target (variable), and so on. “feature” indicates that the constraints are grouped by library feature: the first n_targets columns correspond to the first library feature, the next n_targets columns to the second library feature, and so on.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed. Note that this parameter is incompatible with the constraints!
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
initial_guess (np.ndarray, shape (n_features) or (n_targets, n_features), optional (default None)) – Initial guess for coefficients coef_ to warmstart the optimizer.
verbose (bool, optional (default False)) – If True, prints out the Gurobi solver log.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s).
ind_ (array, shape (n_features,) or (n_targets, n_features)) – Array of 0s and 1s indicating which coefficients of the weight vector have not been masked out, i.e. the support of self.coef_.
model (gurobipy.model) – The raw gurobi model being solved.

property complexity¶

pysindy.optimizers.sindy_optimizer module¶

class pysindy.optimizers.sindy_optimizer.SINDyOptimizer(optimizer, unbias=True)[source]¶

Bases: BaseEstimator

Wrapper class for optimizers/sparse regression methods passed into the SINDy object.

Enables single target regressors (i.e. those whose predictions are 1-dimensional) to perform multi target regression (i.e. predictions are 2-dimensional). Also enhances an _unbias function to reduce bias when regularization is used.

Parameters

optimizer (estimator object) – The optimizer/sparse regressor to be wrapped, implementing fit and predict. optimizer should also have the attributes coef_, fit_intercept, and intercept_. Note that attribute normalize is deprecated as of sklearn versions >= 1.0 and will be removed in future versions.
unbias (boolean, optional (default True)) – Whether to perform an extra step of unregularized linear regression to unbias the coefficients for the identified support. For example, if optimizer=STLSQ(alpha=0.1) is used then the learned coefficients will be biased toward 0 due to the L2 regularization. Setting unbias=True will trigger an additional step wherein the nonzero coefficients learned by the optimizer object will be updated using an unregularized least-squares fit.

fit(x, y)[source]¶

predict(x)[source]¶

property coef_¶

property intercept_¶

property complexity¶

pysindy.optimizers.sindy_pi module¶

class pysindy.optimizers.sindy_pi.SINDyPI(threshold=0.1, tol=1e-05, thresholder='l1', max_iter=10000, fit_intercept=False, copy_X=True, thresholds=None, model_subset=None, normalize_columns=False, verbose_cvxpy=False)[source]¶

Bases: SR3

SINDy-PI optimizer

Attempts to minimize the objective function

\[0.5\|X-Xw\|^2_2 + \lambda R(w)\]

over w where \(R(v)\) is a regularization function. See the following reference for more details:

Kaheman, Kadierdan, J. Nathan Kutz, and Steven L. Brunton. SINDy-PI: a robust algorithm for parallel implicit sparse identification of nonlinear dynamics. Proceedings of the Royal Society A 476.2242 (2020): 20200279.

Parameters

threshold (float, optional (default 0.1)) – Determines the strength of the regularization. When the regularization function R is the l0 norm, the regularization is equivalent to performing hard thresholding, and lambda is chosen to threshold at the value given by this parameter. This is equivalent to choosing lambda = threshold^2 / (2 * nu).
tol (float, optional (default 1e-5)) – Tolerance used for determining convergence of the optimization algorithm.
thresholder (string, optional (default 'l1')) – Regularization function to use. Currently implemented options are ‘l1’ (l1 norm), ‘weighted_l1’ (weighted l1 norm), l2, and ‘weighted_l2’ (weighted l2 norm)
max_iter (int, optional (default 10000)) – Maximum iterations of the optimization algorithm.
fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
normalize_columns (boolean, optional (default False)) – This parameter normalizes the columns of Theta before the optimization is done. This tends to standardize the columns to similar magnitudes, often improving performance.
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
thresholds (np.ndarray, shape (n_targets, n_features), optional (default None)) – Array of thresholds for each library function coefficient. Each row corresponds to a measurement variable and each column to a function from the feature library. Recall that SINDy seeks a matrix \(\Xi\) such that \(\dot{X} \approx \Theta(X)\Xi\). thresholds[i, j] should specify the threshold to be used for the (j + 1, i + 1) entry of \(\Xi\). That is to say it should give the threshold to be used for the (j + 1)st library function in the equation for the (i + 1)st measurement variable.
model_subset (np.ndarray, shape(n_models), optional (default None)) – List of indices to compute models for. If list is not provided, the default is to compute SINDy-PI models for all possible candidate functions. This can take a long time for 4D systems or larger.
verbose_cvxpy (bool, optional (default False)) – Boolean flag which is passed to CVXPY solve function to indicate if output should be verbose or not. Only relevant for optimizers that use the CVXPY package in some capabity.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Regularized weight vector(s). This is the v in the objective function.
unbias (boolean) – Whether to perform an extra step of unregularized linear regression to unbias the coefficients for the identified support. unbias is automatically set to False if a constraint is used and is otherwise left uninitialized.

pysindy.optimizers.sr3 module¶

class pysindy.optimizers.sr3.SR3(threshold=0.1, thresholds=None, nu=1.0, tol=1e-05, thresholder='L0', trimming_fraction=0.0, trimming_step_size=1.0, max_iter=30, fit_intercept=False, copy_X=True, initial_guess=None, normalize_columns=False, verbose=False)[source]¶

Bases: BaseOptimizer

Sparse relaxed regularized regression.

Attempts to minimize the objective function

\[0.5\|y-Xw\|^2_2 + \lambda R(u) + (0.5 / \nu)\|w-u\|^2_2\]

where \(R(u)\) is a regularization function. See the following references for more details:

Zheng, Peng, et al. “A unified framework for sparse relaxed regularized regression: SR3.” IEEE Access 7 (2018): 1404-1423.

Champion, K., Zheng, P., Aravkin, A. Y., Brunton, S. L., & Kutz, J. N. (2020). A unified sparse optimization framework to learn parsimonious physics-informed models from data. IEEE Access, 8, 169259-169271.

Parameters

threshold (float, optional (default 0.1)) – Determines the strength of the regularization. When the regularization function R is the L0 norm, the regularization is equivalent to performing hard thresholding, and lambda is chosen to threshold at the value given by this parameter. This is equivalent to choosing lambda = threshold^2 / (2 * nu).
nu (float, optional (default 1)) – Determines the level of relaxation. Decreasing nu encourages w and v to be close, whereas increasing nu allows the regularized coefficients v to be farther from w.
tol (float, optional (default 1e-5)) – Tolerance used for determining convergence of the optimization algorithm.
thresholder (string, optional (default 'L0')) – Regularization function to use. Currently implemented options are ‘L0’ (L0 norm), ‘L1’ (L1 norm), ‘L2’ (L2 norm) and ‘CAD’ (clipped absolute deviation). Note by ‘L2 norm’ we really mean the squared L2 norm, i.e. ridge regression
trimming_fraction (float, optional (default 0.0)) – Fraction of the data samples to trim during fitting. Should be a float between 0.0 and 1.0. If 0.0, trimming is not performed.
trimming_step_size (float, optional (default 1.0)) – Step size to use in the trimming optimization procedure.
max_iter (int, optional (default 30)) – Maximum iterations of the optimization algorithm.
initial_guess (np.ndarray, shape (n_features) or (n_targets, n_features), optional (default None)) – Initial guess for coefficients coef_. If None, least-squares is used to obtain an initial guess.
fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed.
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
thresholds (np.ndarray, shape (n_targets, n_features), optional (default None)) – Array of thresholds for each library function coefficient. Each row corresponds to a measurement variable and each column to a function from the feature library. Recall that SINDy seeks a matrix \(\Xi\) such that \(\dot{X} \approx \Theta(X)\Xi\). thresholds[i, j] should specify the threshold to be used for the (j + 1, i + 1) entry of \(\Xi\). That is to say it should give the threshold to be used for the (j + 1)st library function in the equation for the (i + 1)st measurement variable.
verbose (bool, optional (default False)) – If True, prints out the different error terms every max_iter / 10 iterations.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Regularized weight vector(s). This is the v in the objective function.
coef_full_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s) that are not subjected to the regularization. This is the w in the objective function.
history_ (list) – History of sparse coefficients. history_[k] contains the sparse coefficients (v in the optimization objective function) at iteration k.

Examples

>>> import numpy as np
>>> from scipy.integrate import odeint
>>> from pysindy import SINDy
>>> from pysindy.optimizers import SR3
>>> lorenz = lambda z,t : [10 * (z[1] - z[0]),
>>>                        z[0] * (28 - z[2]) - z[1],
>>>                        z[0] * z[1] - 8 / 3 * z[2]]
>>> t = np.arange(0, 2, .002)
>>> x = odeint(lorenz, [-8, 8, 27], t)
>>> opt = SR3(threshold=0.1, nu=1)
>>> model = SINDy(optimizer=opt)
>>> model.fit(x, t=t[1] - t[0])
>>> model.print()
x0' = -10.004 1 + 10.004 x0
x1' = 27.994 1 + -0.993 x0 + -1.000 1 x1
x2' = -2.662 x1 + 1.000 1 x0

enable_trimming(trimming_fraction)[source]¶

Enable the trimming of potential outliers.

Parameters: trimming_fraction (float) – The fraction of samples to be trimmed. Must be between 0 and 1.

disable_trimming()[source]¶: Disable trimming of potential outliers.

pysindy.optimizers.ssr module¶

class pysindy.optimizers.ssr.SSR(alpha=0.05, max_iter=20, ridge_kw=None, normalize_columns=False, fit_intercept=False, copy_X=True, criteria='coefficient_value', kappa=None, verbose=False)[source]¶

Bases: BaseOptimizer

Stepwise sparse regression (SSR) greedy algorithm.

Attempts to minimize the objective function \(\|y - Xw\|^2_2 + \alpha \|w\|^2_2\) by iteratively eliminating the smallest coefficient

See the following reference for more details:

Boninsegna, Lorenzo, Feliks Nüske, and Cecilia Clementi. “Sparse learning of stochastic dynamical equations.” The Journal of chemical physics 148.24 (2018): 241723.

Parameters

max_iter (int, optional (default 20)) – Maximum iterations of the optimization algorithm.
fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed.
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
kappa (float, optional (default None)) – If passed, compute the MSE errors with an extra L0 term with strength equal to kappa times the condition number of Theta.
criteria (string, optional (default "coefficient_value")) – The criteria to use for truncating a coefficient each iteration. Must be “coefficient_value” or “model_residual”. “coefficient_value”: zero out the smallest coefficient). “model_residual”: choose the N-1 term model with the smallest residual error.
alpha (float, optional (default 0.05)) – Optional L2 (ridge) regularization on the weight vector.
ridge_kw (dict, optional (default None)) – Optional keyword arguments to pass to the ridge regression.
verbose (bool, optional (default False)) – If True, prints out the different error terms every iteration.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s).
history_ (list) – History of coef_. history_[k] contains the values of coef_ at iteration k of SSR
err_history_ (list) – History of coef_. history_[k] contains the MSE of each coef_ at iteration k of SSR

Examples

>>> import numpy as np
>>> from scipy.integrate import odeint
>>> from pysindy import SINDy
>>> from pysindy.optimizers import SSR
>>> lorenz = lambda z,t : [10 * (z[1] - z[0]),
>>>                        z[0] * (28 - z[2]) - z[1],
>>>                        z[0] * z[1] - 8 / 3 * z[2]]
>>> t = np.arange(0, 2, .002)
>>> x = odeint(lorenz, [-8, 8, 27], t)
>>> opt = SSR(alpha=.5)
>>> model = SINDy(optimizer=opt)
>>> model.fit(x, t=t[1] - t[0])
>>> model.print()
x0' = -9.999 1 + 9.999 x0
x1' = 27.984 1 + -0.996 x0 + -1.000 1 x1
x2' = -2.666 x1 + 1.000 1 x0

pysindy.optimizers.stable_linear_sr3 module¶

class pysindy.optimizers.stable_linear_sr3.StableLinearSR3(threshold=0.1, nu=1.0, tol=1e-05, thresholder='l1', max_iter=30, trimming_fraction=0.0, trimming_step_size=1.0, constraint_lhs=None, constraint_rhs=None, constraint_order='target', normalize_columns=False, fit_intercept=False, copy_X=True, initial_guess=None, thresholds=None, equality_constraints=False, inequality_constraints=False, constraint_separation_index=0, verbose=False, verbose_cvxpy=False, gamma=-1e-08)[source]¶

Bases: ConstrainedSR3

Sparse relaxed regularized regression for building a-priori stable linear models. This requires making a matrix negative definite, which can be challenging. Here we use a similar method to the TrappingOptimizer algorithm. Linear equality and linear inequality constraints are both allowed, as in the ConstrainedSR3 optimizer.

Attempts to minimize the objective function

\[0.5\|y-Xw\|^2_2 + \lambda R(u) + (0.5 / \nu)\|w-u\|^2_2\]

\[\text{subject to } Cu = d, Du = e, w negative definite\]

over u and w, where \(R(u)\) is a regularization function, C and D are constraint matrices, and d and e are vectors of values. NOTE: This optimizer is intended for building purely linear models that are guaranteed to be stable.

Parameters

threshold (float, optional (default 0.1)) – Determines the strength of the regularization. When the regularization function R is the l0 norm, the regularization is equivalent to performing hard thresholding, and lambda is chosen to threshold at the value given by this parameter. This is equivalent to choosing lambda = threshold^2 / (2 * nu).
nu (float, optional (default 1)) – Determines the level of relaxation. Decreasing nu encourages w and v to be close, whereas increasing nu allows the regularized coefficients v to be farther from w.
tol (float, optional (default 1e-5)) – Tolerance used for determining convergence of the optimization algorithm.
thresholder (string, optional (default 'l1')) – Regularization function to use. Currently implemented options are ‘l1’ (l1 norm), ‘l2’ (l2 norm), ‘cad’ (clipped absolute deviation), ‘weighted_l1’ (weighted l1 norm), and ‘weighted_l2’ (weighted l2 norm). Note that the thresholder must be convex here.
max_iter (int, optional (default 30)) – Maximum iterations of the optimization algorithm.
fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
constraint_lhs (numpy ndarray, optional (default None)) – Shape should be (n_constraints, n_features * n_targets), The left hand side matrix C of Cw <= d. There should be one row per constraint.
constraint_rhs (numpy ndarray, shape (n_constraints,), optional (default None)) – The right hand side vector d of Cw <= d.
constraint_order (string, optional (default "target")) – The format in which the constraints constraint_lhs were passed. Must be one of “target” or “feature”. “target” indicates that the constraints are grouped by target: i.e. the first n_features columns correspond to constraint coefficients on the library features for the first target (variable), the next n_features columns to the library features for the second target (variable), and so on. “feature” indicates that the constraints are grouped by library feature: the first n_targets columns correspond to the first library feature, the next n_targets columns to the second library feature, and so on.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed. Note that this parameter is incompatible with the constraints!
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
initial_guess (np.ndarray, optional (default None)) – Shape should be (n_features) or (n_targets, n_features). Initial guess for coefficients coef_, (v in the mathematical equations) If None, least-squares is used to obtain an initial guess.
thresholds (np.ndarray, shape (n_targets, n_features), optional (default None)) – Array of thresholds for each library function coefficient. Each row corresponds to a measurement variable and each column to a function from the feature library. Recall that SINDy seeks a matrix \(\Xi\) such that \(\dot{X} \approx \Theta(X)\Xi\). thresholds[i, j] should specify the threshold to be used for the (j + 1, i + 1) entry of \(\Xi\). That is to say it should give the threshold to be used for the (j + 1)st library function in the equation for the (i + 1)st measurement variable.
inequality_constraints (bool, optional (default False)) – If True, CVXPY methods are used to solve the problem.
verbose (bool, optional (default False)) – If True, prints out the different error terms every max_iter / 10 iterations.
verbose_cvxpy (bool, optional (default False)) – Boolean flag which is passed to CVXPY solve function to indicate if output should be verbose or not. Only relevant for optimizers that use the CVXPY package in some capabity.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Regularized weight vector(s). This is the v in the objective function.
coef_full_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s) that are not subjected to the regularization. This is the w in the objective function.
unbias (boolean) – Whether to perform an extra step of unregularized linear regression to unbias the coefficients for the identified support. unbias is automatically set to False if a constraint is used and is otherwise left uninitialized.

pysindy.optimizers.stlsq module¶

class pysindy.optimizers.stlsq.STLSQ(threshold=0.1, alpha=0.05, max_iter=20, ridge_kw=None, normalize_columns=False, fit_intercept=False, copy_X=True, initial_guess=None, verbose=False)[source]¶

Bases: BaseOptimizer

Sequentially thresholded least squares algorithm. Defaults to doing Sequentially thresholded Ridge regression.

Attempts to minimize the objective function \(\|y - Xw\|^2_2 + \alpha \|w\|^2_2\) by iteratively performing least squares and masking out elements of the weight array w that are below a given threshold.

See the following reference for more details:

Brunton, Steven L., Joshua L. Proctor, and J. Nathan Kutz. “Discovering governing equations from data by sparse identification of nonlinear dynamical systems.” Proceedings of the national academy of sciences 113.15 (2016): 3932-3937.

Parameters

threshold (float, optional (default 0.1)) – Minimum magnitude for a coefficient in the weight vector. Coefficients with magnitude below the threshold are set to zero.
alpha (float, optional (default 0.05)) – Optional L2 (ridge) regularization on the weight vector.
max_iter (int, optional (default 20)) – Maximum iterations of the optimization algorithm.
ridge_kw (dict, optional (default None)) – Optional keyword arguments to pass to the ridge regression.
fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed.
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
initial_guess (np.ndarray, shape (n_features) or (n_targets, n_features),) – optional (default None) Initial guess for coefficients coef_. If None, least-squares is used to obtain an initial guess.
verbose (bool, optional (default False)) – If True, prints out the different error terms every iteration.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s).
ind_ (array, shape (n_features,) or (n_targets, n_features)) – Array of 0s and 1s indicating which coefficients of the weight vector have not been masked out, i.e. the support of self.coef_.
history_ (list) – History of coef_. history_[k] contains the values of coef_ at iteration k of sequentially thresholded least-squares.

Examples

>>> import numpy as np
>>> from scipy.integrate import odeint
>>> from pysindy import SINDy
>>> from pysindy.optimizers import STLSQ
>>> lorenz = lambda z,t : [10*(z[1] - z[0]),
>>>                        z[0]*(28 - z[2]) - z[1],
>>>                        z[0]*z[1] - 8/3*z[2]]
>>> t = np.arange(0,2,.002)
>>> x = odeint(lorenz, [-8,8,27], t)
>>> opt = STLSQ(threshold=.1, alpha=.5)
>>> model = SINDy(optimizer=opt)
>>> model.fit(x, t=t[1]-t[0])
>>> model.print()
x0' = -9.999 1 + 9.999 x0
x1' = 27.984 1 + -0.996 x0 + -1.000 1 x1
x2' = -2.666 x1 + 1.000 1 x0

property complexity¶

pysindy.optimizers.trapping_sr3 module¶

class pysindy.optimizers.trapping_sr3.TrappingSR3(evolve_w=True, threshold=0.1, eps_solver=1e-07, relax_optim=True, inequality_constraints=False, eta=None, alpha_A=None, alpha_m=None, gamma=-0.1, tol=1e-05, tol_m=1e-05, thresholder='l1', thresholds=None, max_iter=30, accel=False, normalize_columns=False, fit_intercept=False, copy_X=True, m0=None, A0=None, objective_history=None, constraint_lhs=None, constraint_rhs=None, constraint_order='target', verbose=False, verbose_cvxpy=False)[source]¶

Bases: SR3

Trapping variant of sparse relaxed regularized regression. This optimizer can be used to identify systems with globally stable (bounded) solutions.

Attempts to minimize one of two related objective functions

\[0.5\|y-Xw\|^2_2 + \lambda R(w) + 0.5\|Pw-A\|^2_2/\eta + \delta_0(Cw-d) + \delta_{\Lambda}(A)\]

or

\[0.5\|y-Xw\|^2_2 + \lambda R(w) + \delta_0(Cw-d) + 0.5 * maximumeigenvalue(A)/\eta\]

where \(R(w)\) is a regularization function, which must be convex, \(\delta_0\) is an indicator function that provides a hard constraint of CW = d, and :math:delta_{Lambda} is a term to project the \(A\) matrix onto the space of negative definite matrices. See the following references for more details:

Kaptanoglu, Alan A., et al. “Promoting global stability in data-driven models of quadratic nonlinear dynamics.” arXiv preprint arXiv:2105.01843 (2021).

Zheng, Peng, et al. “A unified framework for sparse relaxed regularized regression: Sr3.” IEEE Access 7 (2018): 1404-1423.

Champion, Kathleen, et al. “A unified sparse optimization framework to learn parsimonious physics-informed models from data.” IEEE Access 8 (2020): 169259-169271.

Parameters

evolve_w (bool, optional (default True)) – If false, don’t update w and just minimize over (m, A)
threshold (float, optional (default 0.1)) – Determines the strength of the regularization. When the regularization function R is the L0 norm, the regularization is equivalent to performing hard thresholding, and lambda is chosen to threshold at the value given by this parameter. This is equivalent to choosing lambda = threshold^2 / (2 * nu).
eta (float, optional (default 1.0e20)) – Determines the strength of the stability term ||Pw-A||^2 in the optimization. The default value is very large so that the algorithm default is to ignore the stability term. In this limit, this should be approximately equivalent to the ConstrainedSR3 method.
alpha_m (float, optional (default eta * 0.1)) – Determines the step size in the prox-gradient descent over m. For convergence, need alpha_m <= eta / ||w^T * PQ^T * PQ * w||. Typically 0.01 * eta <= alpha_m <= 0.1 * eta.
alpha_A (float, optional (default eta)) – Determines the step size in the prox-gradient descent over A. For convergence, need alpha_A <= eta, so typically alpha_A = eta is used.
gamma (float, optional (default 0.1)) – Determines the negative interval that matrix A is projected onto. For most applications gamma = 0.1 - 1.0 works pretty well.
tol (float, optional (default 1e-5)) – Tolerance used for determining convergence of the optimization algorithm over w.
tol_m (float, optional (default 1e-5)) – Tolerance used for determining convergence of the optimization algorithm over m.
thresholder (string, optional (default 'L1')) – Regularization function to use. For current trapping SINDy, only the L1 and L2 norms are implemented. Note that other convex norms could be straightforwardly implemented, but L0 requires reformulation because of nonconvexity.
thresholds (np.ndarray, shape (n_targets, n_features), optional (default None)) – Array of thresholds for each library function coefficient. Each row corresponds to a measurement variable and each column to a function from the feature library. Recall that SINDy seeks a matrix \(\Xi\) such that \(\dot{X} \approx \Theta(X)\Xi\). thresholds[i, j] should specify the threshold to be used for the (j + 1, i + 1) entry of \(\Xi\). That is to say it should give the threshold to be used for the (j + 1)st library function in the equation for the (i + 1)st measurement variable.
eps_solver (float, optional (default 1.0e-7)) – If threshold != 0, this specifies the error tolerance in the CVXPY (OSQP) solve. Default is 1.0e-3 in OSQP.
relax_optim (bool, optional (default True)) – If relax_optim = True, use the relax-and-split method. If False, try a direct minimization on the largest eigenvalue.
inequality_constraints (bool, optional (default False)) – If True, relax_optim must be false or relax_optim = True AND threshold != 0, so that the CVXPY methods are used.
max_iter (int, optional (default 30)) – Maximum iterations of the optimization algorithm.
accel (bool, optional (default False)) – Whether or not to use accelerated prox-gradient descent for (m, A).
m0 (np.ndarray, shape (n_targets), optional (default None)) – Initial guess for vector m in the optimization. Otherwise each component of m is randomly initialized in [-1, 1].
A0 (np.ndarray, shape (n_targets, n_targets), optional (default None)) – Initial guess for vector A in the optimization. Otherwise A is initialized as A = diag(gamma).
fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed.
verbose (bool, optional (default False)) – If True, prints out the different error terms every max_iter / 10 iterations.
verbose_cvxpy (bool, optional (default False)) – Boolean flag which is passed to CVXPY solve function to indicate if output should be verbose or not. Only relevant for optimizers that use the CVXPY package in some capabity.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Regularized weight vector(s). This is the v in the objective function.
history_ (list) – History of sparse coefficients. history_[k] contains the sparse coefficients (v in the optimization objective function) at iteration k.
objective_history_ (list) – History of the value of the objective at each step. Note that the trapping SINDy problem is nonconvex, meaning that this value may increase and decrease as the algorithm works.
A_history_ (list) – History of the auxiliary variable A that approximates diag(PW).
m_history_ (list) – History of the shift vector m that determines the origin of the trapping region.
PW_history_ (list) – History of PW = A^S, the quantity we are attempting to make negative definite.
PWeigs_history_ (list) – History of diag(PW), a list of the eigenvalues of A^S at each iteration. Tracking this allows us to ascertain if A^S is indeed being pulled towards the space of negative definite matrices.
PL_unsym_ (np.ndarray, shape (n_targets, n_targets, n_targets, n_features)) – Unsymmetrized linear coefficient part of the P matrix in ||Pw - A||^2
PL_ (np.ndarray, shape (n_targets, n_targets, n_targets, n_features)) – Linear coefficient part of the P matrix in ||Pw - A||^2
PQ_ (np.ndarray, shape (n_targets, n_targets,) – n_targets, n_targets, n_features) Quadratic coefficient part of the P matrix in ||Pw - A||^2

Examples

>>> import numpy as np
>>> from scipy.integrate import odeint
>>> from pysindy import SINDy
>>> from pysindy.optimizers import TrappingSR3
>>> lorenz = lambda z,t : [10*(z[1] - z[0]),
>>>                        z[0]*(28 - z[2]) - z[1],
>>>                        z[0]*z[1] - 8/3*z[2]]
>>> t = np.arange(0,2,.002)
>>> x = odeint(lorenz, [-8,8,27], t)
>>> opt = TrappingSR3(threshold=0.1)
>>> model = SINDy(optimizer=opt)
>>> model.fit(x, t=t[1]-t[0])
>>> model.print()
x0' = -10.004 1 + 10.004 x0
x1' = 27.994 1 + -0.993 x0 + -1.000 1 x1
x2' = -2.662 x1 + 1.000 1 x0

Module contents¶

class pysindy.optimizers.BaseOptimizer(max_iter=20, normalize_columns=False, fit_intercept=False, initial_guess=None, copy_X=True)[source]¶

Bases: LinearRegression, ComplexityMixin

Base class for SINDy optimizers. Subclasses must implement a _reduce method for carrying out the bulk of the work of fitting a model.

Parameters

fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed.
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
initial_guess (np.ndarray, shape (n_features,) or (n_targets, n_features),) – optional (default None) Initial guess for coefficients coef_. If None, the initial guess is obtained via a least-squares fit.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s).
ind_ (array, shape (n_features,) or (n_targets, n_features)) – Array of 0s and 1s indicating which coefficients of the weight vector have not been masked out.
history_ (list) – History of coef_ over iterations of the optimization algorithm.
Theta_ (np.ndarray, shape (n_samples, n_features)) – The Theta matrix to be used in the optimization. We save it as an attribute because access to the full library of terms is sometimes needed for various applications.

fit(x_, y, sample_weight=None, **reduce_kws)[source]¶

Fit the model.

Parameters

x (array-like, shape (n_samples, n_features)) – Training data
y (array-like, shape (n_samples,) or (n_samples, n_targets)) – Target values
sample_weight (float or numpy array of shape (n_samples,), optional) – Individual weights for each sample
reduce_kws (dict) – Optional keyword arguments to pass to the _reduce method (implemented by subclasses)

Returns

self

Return type

returns an instance of self

class pysindy.optimizers.EnsembleOptimizer(opt: BaseOptimizer, bagging: bool = False, library_ensemble: bool = False, n_models: int = 20, n_subset: Optional[int] = None, n_candidates_to_drop: int = 1, replace: bool = True, ensemble_aggregator: Optional[Callable] = None)[source]¶

Bases: BaseOptimizer

Wrapper class for ensembling methods.

Parameters

opt (BaseOptimizer) – The underlying optimizer to run on each ensemble
bagging (boolean, optional (default False)) – This parameter is used to allow for “ensembling”, i.e. the generation of many SINDy models (n_models) by choosing a random temporal subset of the input data (n_subset) for each sparse regression. This often improves robustness because averages (bagging) or medians (bragging) of all the models are usually quite high-performing. The user can also generate “distributions” of many models, and calculate how often certain library terms are included in a model.
library_ensemble (boolean, optional (default False)) – This parameter is used to allow for “library ensembling”, i.e. the generation of many SINDy models (n_models) by choosing a random subset of the candidate library terms to truncate. So, n_models are generated by solving n_models sparse regression problems on these “reduced” libraries. Once again, this often improves robustness because averages (bagging) or medians (bragging) of all the models are usually quite high-performing. The user can also generate “distributions” of many models, and calculate how often certain library terms are included in a model.
n_models (int, optional (default 20)) – Number of models to generate via ensemble
n_subset (int, optional (default len(time base))) – Number of time points to use for ensemble
n_candidates_to_drop (int, optional (default 1)) – Number of candidate terms in the feature library to drop during library ensembling.
replace (boolean, optional (default True)) – If ensemble true, whether or not to time sample with replacement.
ensemble_aggregator (callable, optional (default numpy.median)) – Method to aggregate model coefficients across different samples. This method argument is only used if ensemble or library_ensemble is True. The method should take in a list of 2D arrays and return a 2D array of the same shape as the arrays in the list. Example: lambda x: np.median(x, axis=0)

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Regularized weight vector(s). This is the v in the objective function.
coef_full_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s) that are not subjected to the regularization. This is the w in the objective function.
unbias (boolean) – Whether to perform an extra step of unregularized linear regression to unbias the coefficients for the identified support. unbias is automatically set to False if a constraint is used and is otherwise left uninitialized.

class pysindy.optimizers.SINDyOptimizer(optimizer, unbias=True)[source]¶

Bases: BaseEstimator

Wrapper class for optimizers/sparse regression methods passed into the SINDy object.

Enables single target regressors (i.e. those whose predictions are 1-dimensional) to perform multi target regression (i.e. predictions are 2-dimensional). Also enhances an _unbias function to reduce bias when regularization is used.

Parameters

optimizer (estimator object) – The optimizer/sparse regressor to be wrapped, implementing fit and predict. optimizer should also have the attributes coef_, fit_intercept, and intercept_. Note that attribute normalize is deprecated as of sklearn versions >= 1.0 and will be removed in future versions.
unbias (boolean, optional (default True)) – Whether to perform an extra step of unregularized linear regression to unbias the coefficients for the identified support. For example, if optimizer=STLSQ(alpha=0.1) is used then the learned coefficients will be biased toward 0 due to the L2 regularization. Setting unbias=True will trigger an additional step wherein the nonzero coefficients learned by the optimizer object will be updated using an unregularized least-squares fit.

fit(x, y)[source]¶

predict(x)[source]¶

property coef_¶

property intercept_¶

property complexity¶

class pysindy.optimizers.SR3(threshold=0.1, thresholds=None, nu=1.0, tol=1e-05, thresholder='L0', trimming_fraction=0.0, trimming_step_size=1.0, max_iter=30, fit_intercept=False, copy_X=True, initial_guess=None, normalize_columns=False, verbose=False)[source]¶

Bases: BaseOptimizer

Sparse relaxed regularized regression.

Attempts to minimize the objective function

\[0.5\|y-Xw\|^2_2 + \lambda R(u) + (0.5 / \nu)\|w-u\|^2_2\]

where \(R(u)\) is a regularization function. See the following references for more details:

Zheng, Peng, et al. “A unified framework for sparse relaxed regularized regression: SR3.” IEEE Access 7 (2018): 1404-1423.

Champion, K., Zheng, P., Aravkin, A. Y., Brunton, S. L., & Kutz, J. N. (2020). A unified sparse optimization framework to learn parsimonious physics-informed models from data. IEEE Access, 8, 169259-169271.

Parameters

threshold (float, optional (default 0.1)) – Determines the strength of the regularization. When the regularization function R is the L0 norm, the regularization is equivalent to performing hard thresholding, and lambda is chosen to threshold at the value given by this parameter. This is equivalent to choosing lambda = threshold^2 / (2 * nu).
nu (float, optional (default 1)) – Determines the level of relaxation. Decreasing nu encourages w and v to be close, whereas increasing nu allows the regularized coefficients v to be farther from w.
tol (float, optional (default 1e-5)) – Tolerance used for determining convergence of the optimization algorithm.
thresholder (string, optional (default 'L0')) – Regularization function to use. Currently implemented options are ‘L0’ (L0 norm), ‘L1’ (L1 norm), ‘L2’ (L2 norm) and ‘CAD’ (clipped absolute deviation). Note by ‘L2 norm’ we really mean the squared L2 norm, i.e. ridge regression
trimming_fraction (float, optional (default 0.0)) – Fraction of the data samples to trim during fitting. Should be a float between 0.0 and 1.0. If 0.0, trimming is not performed.
trimming_step_size (float, optional (default 1.0)) – Step size to use in the trimming optimization procedure.
max_iter (int, optional (default 30)) – Maximum iterations of the optimization algorithm.
initial_guess (np.ndarray, shape (n_features) or (n_targets, n_features), optional (default None)) – Initial guess for coefficients coef_. If None, least-squares is used to obtain an initial guess.
fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed.
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
thresholds (np.ndarray, shape (n_targets, n_features), optional (default None)) – Array of thresholds for each library function coefficient. Each row corresponds to a measurement variable and each column to a function from the feature library. Recall that SINDy seeks a matrix \(\Xi\) such that \(\dot{X} \approx \Theta(X)\Xi\). thresholds[i, j] should specify the threshold to be used for the (j + 1, i + 1) entry of \(\Xi\). That is to say it should give the threshold to be used for the (j + 1)st library function in the equation for the (i + 1)st measurement variable.
verbose (bool, optional (default False)) – If True, prints out the different error terms every max_iter / 10 iterations.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Regularized weight vector(s). This is the v in the objective function.
coef_full_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s) that are not subjected to the regularization. This is the w in the objective function.
history_ (list) – History of sparse coefficients. history_[k] contains the sparse coefficients (v in the optimization objective function) at iteration k.

Examples

>>> import numpy as np
>>> from scipy.integrate import odeint
>>> from pysindy import SINDy
>>> from pysindy.optimizers import SR3
>>> lorenz = lambda z,t : [10 * (z[1] - z[0]),
>>>                        z[0] * (28 - z[2]) - z[1],
>>>                        z[0] * z[1] - 8 / 3 * z[2]]
>>> t = np.arange(0, 2, .002)
>>> x = odeint(lorenz, [-8, 8, 27], t)
>>> opt = SR3(threshold=0.1, nu=1)
>>> model = SINDy(optimizer=opt)
>>> model.fit(x, t=t[1] - t[0])
>>> model.print()
x0' = -10.004 1 + 10.004 x0
x1' = 27.994 1 + -0.993 x0 + -1.000 1 x1
x2' = -2.662 x1 + 1.000 1 x0

enable_trimming(trimming_fraction)[source]¶

Enable the trimming of potential outliers.

Parameters: trimming_fraction (float) – The fraction of samples to be trimmed. Must be between 0 and 1.

disable_trimming()[source]¶: Disable trimming of potential outliers.

class pysindy.optimizers.STLSQ(threshold=0.1, alpha=0.05, max_iter=20, ridge_kw=None, normalize_columns=False, fit_intercept=False, copy_X=True, initial_guess=None, verbose=False)[source]¶

Bases: BaseOptimizer

Sequentially thresholded least squares algorithm. Defaults to doing Sequentially thresholded Ridge regression.

Attempts to minimize the objective function \(\|y - Xw\|^2_2 + \alpha \|w\|^2_2\) by iteratively performing least squares and masking out elements of the weight array w that are below a given threshold.

See the following reference for more details:

Brunton, Steven L., Joshua L. Proctor, and J. Nathan Kutz. “Discovering governing equations from data by sparse identification of nonlinear dynamical systems.” Proceedings of the national academy of sciences 113.15 (2016): 3932-3937.

Parameters

threshold (float, optional (default 0.1)) – Minimum magnitude for a coefficient in the weight vector. Coefficients with magnitude below the threshold are set to zero.
alpha (float, optional (default 0.05)) – Optional L2 (ridge) regularization on the weight vector.
max_iter (int, optional (default 20)) – Maximum iterations of the optimization algorithm.
ridge_kw (dict, optional (default None)) – Optional keyword arguments to pass to the ridge regression.
fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed.
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
initial_guess (np.ndarray, shape (n_features) or (n_targets, n_features),) – optional (default None) Initial guess for coefficients coef_. If None, least-squares is used to obtain an initial guess.
verbose (bool, optional (default False)) – If True, prints out the different error terms every iteration.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s).
ind_ (array, shape (n_features,) or (n_targets, n_features)) – Array of 0s and 1s indicating which coefficients of the weight vector have not been masked out, i.e. the support of self.coef_.
history_ (list) – History of coef_. history_[k] contains the values of coef_ at iteration k of sequentially thresholded least-squares.

Examples

>>> import numpy as np
>>> from scipy.integrate import odeint
>>> from pysindy import SINDy
>>> from pysindy.optimizers import STLSQ
>>> lorenz = lambda z,t : [10*(z[1] - z[0]),
>>>                        z[0]*(28 - z[2]) - z[1],
>>>                        z[0]*z[1] - 8/3*z[2]]
>>> t = np.arange(0,2,.002)
>>> x = odeint(lorenz, [-8,8,27], t)
>>> opt = STLSQ(threshold=.1, alpha=.5)
>>> model = SINDy(optimizer=opt)
>>> model.fit(x, t=t[1]-t[0])
>>> model.print()
x0' = -9.999 1 + 9.999 x0
x1' = 27.984 1 + -0.996 x0 + -1.000 1 x1
x2' = -2.666 x1 + 1.000 1 x0

property complexity¶

class pysindy.optimizers.ConstrainedSR3(threshold=0.1, nu=1.0, tol=1e-05, thresholder='l0', max_iter=30, trimming_fraction=0.0, trimming_step_size=1.0, constraint_lhs=None, constraint_rhs=None, constraint_order='target', normalize_columns=False, fit_intercept=False, copy_X=True, initial_guess=None, thresholds=None, equality_constraints=False, inequality_constraints=False, constraint_separation_index=0, verbose=False, verbose_cvxpy=False)[source]¶

Bases: SR3

Sparse relaxed regularized regression with linear (in)equality constraints.

Attempts to minimize the objective function

\[0.5\|y-Xw\|^2_2 + \lambda R(u) + (0.5 / \nu)\|w-u\|^2_2\]

\[\text{subject to } Cw = d\]

over u and w, where \(R(u)\) is a regularization function, C is a constraint matrix, and d is a vector of values. See the following reference for more details:

Champion, Kathleen, et al. “A unified sparse optimization framework to learn parsimonious physics-informed models from data.” IEEE Access 8 (2020): 169259-169271.

Parameters

threshold (float, optional (default 0.1)) – Determines the strength of the regularization. When the regularization function R is the l0 norm, the regularization is equivalent to performing hard thresholding, and lambda is chosen to threshold at the value given by this parameter. This is equivalent to choosing lambda = threshold^2 / (2 * nu).
nu (float, optional (default 1)) – Determines the level of relaxation. Decreasing nu encourages w and v to be close, whereas increasing nu allows the regularized coefficients v to be farther from w.
tol (float, optional (default 1e-5)) – Tolerance used for determining convergence of the optimization algorithm.
thresholder (string, optional (default 'l0')) – Regularization function to use. Currently implemented options are ‘l0’ (l0 norm), ‘l1’ (l1 norm), ‘l2’ (l2 norm), ‘cad’ (clipped absolute deviation), ‘weighted_l0’ (weighted l0 norm), ‘weighted_l1’ (weighted l1 norm), and ‘weighted_l2’ (weighted l2 norm).
max_iter (int, optional (default 30)) – Maximum iterations of the optimization algorithm.
fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
constraint_lhs (numpy ndarray, optional (default None)) – Shape should be (n_constraints, n_features * n_targets), The left hand side matrix C of Cw <= d. There should be one row per constraint.
constraint_rhs (numpy ndarray, shape (n_constraints,), optional (default None)) – The right hand side vector d of Cw <= d.
constraint_order (string, optional (default "target")) – The format in which the constraints constraint_lhs were passed. Must be one of “target” or “feature”. “target” indicates that the constraints are grouped by target: i.e. the first n_features columns correspond to constraint coefficients on the library features for the first target (variable), the next n_features columns to the library features for the second target (variable), and so on. “feature” indicates that the constraints are grouped by library feature: the first n_targets columns correspond to the first library feature, the next n_targets columns to the second library feature, and so on.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed. Note that this parameter is incompatible with the constraints!
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
initial_guess (np.ndarray, optional (default None)) – Shape should be (n_features) or (n_targets, n_features). Initial guess for coefficients coef_, (v in the mathematical equations) If None, least-squares is used to obtain an initial guess.
thresholds (np.ndarray, shape (n_targets, n_features), optional (default None)) – Array of thresholds for each library function coefficient. Each row corresponds to a measurement variable and each column to a function from the feature library. Recall that SINDy seeks a matrix \(\Xi\) such that \(\dot{X} \approx \Theta(X)\Xi\). thresholds[i, j] should specify the threshold to be used for the (j + 1, i + 1) entry of \(\Xi\). That is to say it should give the threshold to be used for the (j + 1)st library function in the equation for the (i + 1)st measurement variable.
inequality_constraints (bool, optional (default False)) – If True, CVXPY methods are used to solve the problem.
verbose (bool, optional (default False)) – If True, prints out the different error terms every max_iter / 10 iterations.
verbose_cvxpy (bool, optional (default False)) – Boolean flag which is passed to CVXPY solve function to indicate if output should be verbose or not. Only relevant for optimizers that use the CVXPY package in some capabity.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Regularized weight vector(s). This is the v in the objective function.
coef_full_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s) that are not subjected to the regularization. This is the w in the objective function.
unbias (boolean) – Whether to perform an extra step of unregularized linear regression to unbias the coefficients for the identified support. unbias is automatically set to False if a constraint is used and is otherwise left uninitialized.

class pysindy.optimizers.StableLinearSR3(threshold=0.1, nu=1.0, tol=1e-05, thresholder='l1', max_iter=30, trimming_fraction=0.0, trimming_step_size=1.0, constraint_lhs=None, constraint_rhs=None, constraint_order='target', normalize_columns=False, fit_intercept=False, copy_X=True, initial_guess=None, thresholds=None, equality_constraints=False, inequality_constraints=False, constraint_separation_index=0, verbose=False, verbose_cvxpy=False, gamma=-1e-08)[source]¶

Bases: ConstrainedSR3

Sparse relaxed regularized regression for building a-priori stable linear models. This requires making a matrix negative definite, which can be challenging. Here we use a similar method to the TrappingOptimizer algorithm. Linear equality and linear inequality constraints are both allowed, as in the ConstrainedSR3 optimizer.

Attempts to minimize the objective function

\[0.5\|y-Xw\|^2_2 + \lambda R(u) + (0.5 / \nu)\|w-u\|^2_2\]

\[\text{subject to } Cu = d, Du = e, w negative definite\]

over u and w, where \(R(u)\) is a regularization function, C and D are constraint matrices, and d and e are vectors of values. NOTE: This optimizer is intended for building purely linear models that are guaranteed to be stable.

Parameters

threshold (float, optional (default 0.1)) – Determines the strength of the regularization. When the regularization function R is the l0 norm, the regularization is equivalent to performing hard thresholding, and lambda is chosen to threshold at the value given by this parameter. This is equivalent to choosing lambda = threshold^2 / (2 * nu).
nu (float, optional (default 1)) – Determines the level of relaxation. Decreasing nu encourages w and v to be close, whereas increasing nu allows the regularized coefficients v to be farther from w.
tol (float, optional (default 1e-5)) – Tolerance used for determining convergence of the optimization algorithm.
thresholder (string, optional (default 'l1')) – Regularization function to use. Currently implemented options are ‘l1’ (l1 norm), ‘l2’ (l2 norm), ‘cad’ (clipped absolute deviation), ‘weighted_l1’ (weighted l1 norm), and ‘weighted_l2’ (weighted l2 norm). Note that the thresholder must be convex here.
max_iter (int, optional (default 30)) – Maximum iterations of the optimization algorithm.
fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
constraint_lhs (numpy ndarray, optional (default None)) – Shape should be (n_constraints, n_features * n_targets), The left hand side matrix C of Cw <= d. There should be one row per constraint.
constraint_rhs (numpy ndarray, shape (n_constraints,), optional (default None)) – The right hand side vector d of Cw <= d.
constraint_order (string, optional (default "target")) – The format in which the constraints constraint_lhs were passed. Must be one of “target” or “feature”. “target” indicates that the constraints are grouped by target: i.e. the first n_features columns correspond to constraint coefficients on the library features for the first target (variable), the next n_features columns to the library features for the second target (variable), and so on. “feature” indicates that the constraints are grouped by library feature: the first n_targets columns correspond to the first library feature, the next n_targets columns to the second library feature, and so on.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed. Note that this parameter is incompatible with the constraints!
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
initial_guess (np.ndarray, optional (default None)) – Shape should be (n_features) or (n_targets, n_features). Initial guess for coefficients coef_, (v in the mathematical equations) If None, least-squares is used to obtain an initial guess.
thresholds (np.ndarray, shape (n_targets, n_features), optional (default None)) – Array of thresholds for each library function coefficient. Each row corresponds to a measurement variable and each column to a function from the feature library. Recall that SINDy seeks a matrix \(\Xi\) such that \(\dot{X} \approx \Theta(X)\Xi\). thresholds[i, j] should specify the threshold to be used for the (j + 1, i + 1) entry of \(\Xi\). That is to say it should give the threshold to be used for the (j + 1)st library function in the equation for the (i + 1)st measurement variable.
inequality_constraints (bool, optional (default False)) – If True, CVXPY methods are used to solve the problem.
verbose (bool, optional (default False)) – If True, prints out the different error terms every max_iter / 10 iterations.
verbose_cvxpy (bool, optional (default False)) – Boolean flag which is passed to CVXPY solve function to indicate if output should be verbose or not. Only relevant for optimizers that use the CVXPY package in some capabity.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Regularized weight vector(s). This is the v in the objective function.
coef_full_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s) that are not subjected to the regularization. This is the w in the objective function.
unbias (boolean) – Whether to perform an extra step of unregularized linear regression to unbias the coefficients for the identified support. unbias is automatically set to False if a constraint is used and is otherwise left uninitialized.

class pysindy.optimizers.TrappingSR3(evolve_w=True, threshold=0.1, eps_solver=1e-07, relax_optim=True, inequality_constraints=False, eta=None, alpha_A=None, alpha_m=None, gamma=-0.1, tol=1e-05, tol_m=1e-05, thresholder='l1', thresholds=None, max_iter=30, accel=False, normalize_columns=False, fit_intercept=False, copy_X=True, m0=None, A0=None, objective_history=None, constraint_lhs=None, constraint_rhs=None, constraint_order='target', verbose=False, verbose_cvxpy=False)[source]¶

Bases: SR3

Trapping variant of sparse relaxed regularized regression. This optimizer can be used to identify systems with globally stable (bounded) solutions.

Attempts to minimize one of two related objective functions

\[0.5\|y-Xw\|^2_2 + \lambda R(w) + 0.5\|Pw-A\|^2_2/\eta + \delta_0(Cw-d) + \delta_{\Lambda}(A)\]

or

\[0.5\|y-Xw\|^2_2 + \lambda R(w) + \delta_0(Cw-d) + 0.5 * maximumeigenvalue(A)/\eta\]

where \(R(w)\) is a regularization function, which must be convex, \(\delta_0\) is an indicator function that provides a hard constraint of CW = d, and :math:delta_{Lambda} is a term to project the \(A\) matrix onto the space of negative definite matrices. See the following references for more details:

Kaptanoglu, Alan A., et al. “Promoting global stability in data-driven models of quadratic nonlinear dynamics.” arXiv preprint arXiv:2105.01843 (2021).

Zheng, Peng, et al. “A unified framework for sparse relaxed regularized regression: Sr3.” IEEE Access 7 (2018): 1404-1423.

Champion, Kathleen, et al. “A unified sparse optimization framework to learn parsimonious physics-informed models from data.” IEEE Access 8 (2020): 169259-169271.

Parameters

evolve_w (bool, optional (default True)) – If false, don’t update w and just minimize over (m, A)
threshold (float, optional (default 0.1)) – Determines the strength of the regularization. When the regularization function R is the L0 norm, the regularization is equivalent to performing hard thresholding, and lambda is chosen to threshold at the value given by this parameter. This is equivalent to choosing lambda = threshold^2 / (2 * nu).
eta (float, optional (default 1.0e20)) – Determines the strength of the stability term ||Pw-A||^2 in the optimization. The default value is very large so that the algorithm default is to ignore the stability term. In this limit, this should be approximately equivalent to the ConstrainedSR3 method.
alpha_m (float, optional (default eta * 0.1)) – Determines the step size in the prox-gradient descent over m. For convergence, need alpha_m <= eta / ||w^T * PQ^T * PQ * w||. Typically 0.01 * eta <= alpha_m <= 0.1 * eta.
alpha_A (float, optional (default eta)) – Determines the step size in the prox-gradient descent over A. For convergence, need alpha_A <= eta, so typically alpha_A = eta is used.
gamma (float, optional (default 0.1)) – Determines the negative interval that matrix A is projected onto. For most applications gamma = 0.1 - 1.0 works pretty well.
tol (float, optional (default 1e-5)) – Tolerance used for determining convergence of the optimization algorithm over w.
tol_m (float, optional (default 1e-5)) – Tolerance used for determining convergence of the optimization algorithm over m.
thresholder (string, optional (default 'L1')) – Regularization function to use. For current trapping SINDy, only the L1 and L2 norms are implemented. Note that other convex norms could be straightforwardly implemented, but L0 requires reformulation because of nonconvexity.
thresholds (np.ndarray, shape (n_targets, n_features), optional (default None)) – Array of thresholds for each library function coefficient. Each row corresponds to a measurement variable and each column to a function from the feature library. Recall that SINDy seeks a matrix \(\Xi\) such that \(\dot{X} \approx \Theta(X)\Xi\). thresholds[i, j] should specify the threshold to be used for the (j + 1, i + 1) entry of \(\Xi\). That is to say it should give the threshold to be used for the (j + 1)st library function in the equation for the (i + 1)st measurement variable.
eps_solver (float, optional (default 1.0e-7)) – If threshold != 0, this specifies the error tolerance in the CVXPY (OSQP) solve. Default is 1.0e-3 in OSQP.
relax_optim (bool, optional (default True)) – If relax_optim = True, use the relax-and-split method. If False, try a direct minimization on the largest eigenvalue.
inequality_constraints (bool, optional (default False)) – If True, relax_optim must be false or relax_optim = True AND threshold != 0, so that the CVXPY methods are used.
max_iter (int, optional (default 30)) – Maximum iterations of the optimization algorithm.
accel (bool, optional (default False)) – Whether or not to use accelerated prox-gradient descent for (m, A).
m0 (np.ndarray, shape (n_targets), optional (default None)) – Initial guess for vector m in the optimization. Otherwise each component of m is randomly initialized in [-1, 1].
A0 (np.ndarray, shape (n_targets, n_targets), optional (default None)) – Initial guess for vector A in the optimization. Otherwise A is initialized as A = diag(gamma).
fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed.
verbose (bool, optional (default False)) – If True, prints out the different error terms every max_iter / 10 iterations.
verbose_cvxpy (bool, optional (default False)) – Boolean flag which is passed to CVXPY solve function to indicate if output should be verbose or not. Only relevant for optimizers that use the CVXPY package in some capabity.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Regularized weight vector(s). This is the v in the objective function.
history_ (list) – History of sparse coefficients. history_[k] contains the sparse coefficients (v in the optimization objective function) at iteration k.
objective_history_ (list) – History of the value of the objective at each step. Note that the trapping SINDy problem is nonconvex, meaning that this value may increase and decrease as the algorithm works.
A_history_ (list) – History of the auxiliary variable A that approximates diag(PW).
m_history_ (list) – History of the shift vector m that determines the origin of the trapping region.
PW_history_ (list) – History of PW = A^S, the quantity we are attempting to make negative definite.
PWeigs_history_ (list) – History of diag(PW), a list of the eigenvalues of A^S at each iteration. Tracking this allows us to ascertain if A^S is indeed being pulled towards the space of negative definite matrices.
PL_unsym_ (np.ndarray, shape (n_targets, n_targets, n_targets, n_features)) – Unsymmetrized linear coefficient part of the P matrix in ||Pw - A||^2
PL_ (np.ndarray, shape (n_targets, n_targets, n_targets, n_features)) – Linear coefficient part of the P matrix in ||Pw - A||^2
PQ_ (np.ndarray, shape (n_targets, n_targets,) – n_targets, n_targets, n_features) Quadratic coefficient part of the P matrix in ||Pw - A||^2

Examples

>>> import numpy as np
>>> from scipy.integrate import odeint
>>> from pysindy import SINDy
>>> from pysindy.optimizers import TrappingSR3
>>> lorenz = lambda z,t : [10*(z[1] - z[0]),
>>>                        z[0]*(28 - z[2]) - z[1],
>>>                        z[0]*z[1] - 8/3*z[2]]
>>> t = np.arange(0,2,.002)
>>> x = odeint(lorenz, [-8,8,27], t)
>>> opt = TrappingSR3(threshold=0.1)
>>> model = SINDy(optimizer=opt)
>>> model.fit(x, t=t[1]-t[0])
>>> model.print()
x0' = -10.004 1 + 10.004 x0
x1' = 27.994 1 + -0.993 x0 + -1.000 1 x1
x2' = -2.662 x1 + 1.000 1 x0

class pysindy.optimizers.SSR(alpha=0.05, max_iter=20, ridge_kw=None, normalize_columns=False, fit_intercept=False, copy_X=True, criteria='coefficient_value', kappa=None, verbose=False)[source]¶

Bases: BaseOptimizer

Stepwise sparse regression (SSR) greedy algorithm.

Attempts to minimize the objective function \(\|y - Xw\|^2_2 + \alpha \|w\|^2_2\) by iteratively eliminating the smallest coefficient

See the following reference for more details:

Boninsegna, Lorenzo, Feliks Nüske, and Cecilia Clementi. “Sparse learning of stochastic dynamical equations.” The Journal of chemical physics 148.24 (2018): 241723.

Parameters

max_iter (int, optional (default 20)) – Maximum iterations of the optimization algorithm.
fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed.
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
kappa (float, optional (default None)) – If passed, compute the MSE errors with an extra L0 term with strength equal to kappa times the condition number of Theta.
criteria (string, optional (default "coefficient_value")) – The criteria to use for truncating a coefficient each iteration. Must be “coefficient_value” or “model_residual”. “coefficient_value”: zero out the smallest coefficient). “model_residual”: choose the N-1 term model with the smallest residual error.
alpha (float, optional (default 0.05)) – Optional L2 (ridge) regularization on the weight vector.
ridge_kw (dict, optional (default None)) – Optional keyword arguments to pass to the ridge regression.
verbose (bool, optional (default False)) – If True, prints out the different error terms every iteration.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s).
history_ (list) – History of coef_. history_[k] contains the values of coef_ at iteration k of SSR
err_history_ (list) – History of coef_. history_[k] contains the MSE of each coef_ at iteration k of SSR

Examples

>>> import numpy as np
>>> from scipy.integrate import odeint
>>> from pysindy import SINDy
>>> from pysindy.optimizers import SSR
>>> lorenz = lambda z,t : [10 * (z[1] - z[0]),
>>>                        z[0] * (28 - z[2]) - z[1],
>>>                        z[0] * z[1] - 8 / 3 * z[2]]
>>> t = np.arange(0, 2, .002)
>>> x = odeint(lorenz, [-8, 8, 27], t)
>>> opt = SSR(alpha=.5)
>>> model = SINDy(optimizer=opt)
>>> model.fit(x, t=t[1] - t[0])
>>> model.print()
x0' = -9.999 1 + 9.999 x0
x1' = 27.984 1 + -0.996 x0 + -1.000 1 x1
x2' = -2.666 x1 + 1.000 1 x0

class pysindy.optimizers.FROLS(normalize_columns=False, fit_intercept=False, copy_X=True, kappa=None, max_iter=10, alpha=0.05, ridge_kw=None, verbose=False)[source]¶

Bases: BaseOptimizer

Forward Regression Orthogonal Least-Squares (FROLS) optimizer.

Attempts to minimize the objective function \(\|y - Xw\|^2_2 + \alpha \|w\|^2_2\) by iteractively selecting the most correlated function in the library. This is a greedy algorithm.

See the following reference for more details:

Billings, Stephen A. Nonlinear system identification: NARMAX methods in the time, frequency, and spatio-temporal domains. John Wiley & Sons, 2013.

Parameters

fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed.
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
kappa (float, optional (default None)) – If passed, compute the MSE errors with an extra L0 term with strength equal to kappa times the condition number of Theta.
max_iter (int, optional (default 10)) – Maximum iterations of the optimization algorithm. This determines the number of nonzero terms chosen by the FROLS algorithm.
alpha (float, optional (default 0.05)) – Optional L2 (ridge) regularization on the weight vector.
ridge_kw (dict, optional (default None)) – Optional keyword arguments to pass to the ridge regression.
verbose (bool, optional (default False)) – If True, prints out the different error terms every iteration.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s).
history_ (list) – History of coef_. history_[k] contains the values of coef_ at iteration k of FROLS.

Examples

>>> import numpy as np
>>> from scipy.integrate import odeint
>>> from pysindy import SINDy
>>> from pysindy.optimizers import FROLS
>>> lorenz = lambda z,t : [10 * (z[1] - z[0]),
>>>                        z[0] * (28 - z[2]) - z[1],
>>>                        z[0] * z[1] - 8 / 3 * z[2]]
>>> t = np.arange(0, 2, .002)
>>> x = odeint(lorenz, [-8, 8, 27], t)
>>> opt = FROLS(threshold=.1, alpha=.5)
>>> model = SINDy(optimizer=opt)
>>> model.fit(x, t=t[1] - t[0])
>>> model.print()
x0' = -9.999 1 + 9.999 x0
x1' = 27.984 1 + -0.996 x0 + -1.000 1 x1
x2' = -2.666 x1 + 1.000 1 x0

class pysindy.optimizers.SINDyPI(threshold=0.1, tol=1e-05, thresholder='l1', max_iter=10000, fit_intercept=False, copy_X=True, thresholds=None, model_subset=None, normalize_columns=False, verbose_cvxpy=False)[source]¶

Bases: SR3

SINDy-PI optimizer

Attempts to minimize the objective function

\[0.5\|X-Xw\|^2_2 + \lambda R(w)\]

over w where \(R(v)\) is a regularization function. See the following reference for more details:

Kaheman, Kadierdan, J. Nathan Kutz, and Steven L. Brunton. SINDy-PI: a robust algorithm for parallel implicit sparse identification of nonlinear dynamics. Proceedings of the Royal Society A 476.2242 (2020): 20200279.

Parameters

threshold (float, optional (default 0.1)) – Determines the strength of the regularization. When the regularization function R is the l0 norm, the regularization is equivalent to performing hard thresholding, and lambda is chosen to threshold at the value given by this parameter. This is equivalent to choosing lambda = threshold^2 / (2 * nu).
tol (float, optional (default 1e-5)) – Tolerance used for determining convergence of the optimization algorithm.
thresholder (string, optional (default 'l1')) – Regularization function to use. Currently implemented options are ‘l1’ (l1 norm), ‘weighted_l1’ (weighted l1 norm), l2, and ‘weighted_l2’ (weighted l2 norm)
max_iter (int, optional (default 10000)) – Maximum iterations of the optimization algorithm.
fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
normalize_columns (boolean, optional (default False)) – This parameter normalizes the columns of Theta before the optimization is done. This tends to standardize the columns to similar magnitudes, often improving performance.
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
thresholds (np.ndarray, shape (n_targets, n_features), optional (default None)) – Array of thresholds for each library function coefficient. Each row corresponds to a measurement variable and each column to a function from the feature library. Recall that SINDy seeks a matrix \(\Xi\) such that \(\dot{X} \approx \Theta(X)\Xi\). thresholds[i, j] should specify the threshold to be used for the (j + 1, i + 1) entry of \(\Xi\). That is to say it should give the threshold to be used for the (j + 1)st library function in the equation for the (i + 1)st measurement variable.
model_subset (np.ndarray, shape(n_models), optional (default None)) – List of indices to compute models for. If list is not provided, the default is to compute SINDy-PI models for all possible candidate functions. This can take a long time for 4D systems or larger.
verbose_cvxpy (bool, optional (default False)) – Boolean flag which is passed to CVXPY solve function to indicate if output should be verbose or not. Only relevant for optimizers that use the CVXPY package in some capabity.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Regularized weight vector(s). This is the v in the objective function.
unbias (boolean) – Whether to perform an extra step of unregularized linear regression to unbias the coefficients for the identified support. unbias is automatically set to False if a constraint is used and is otherwise left uninitialized.

class pysindy.optimizers.MIOSR(target_sparsity=5, group_sparsity=None, alpha=0.01, regression_timeout=10, fit_intercept=False, constraint_lhs=None, constraint_rhs=None, constraint_order='target', normalize_columns=False, copy_X=True, initial_guess=None, verbose=False)[source]¶

Bases: BaseOptimizer

Mixed-Integer Optimized Sparse Regression.

Solves the sparsity constrained regression problem to provable optimality .. math:

\|y-Xw\|^2_2 + \lambda R(u)

\[\text{subject to } \|w\|_0 \leq k\]

by using type-1 specially ordered sets (SOS1) to encode the support of the coefficients. Can optionally add additional constraints on the coefficients or access the gurobi model directly for advanced usage. See the following reference for additional details:

Bertsimas, D. and Gurnee, W., 2022. Learning Sparse Nonlinear Dynamics via Mixed-Integer Optimization. arXiv preprint arXiv:2206.00176.

Parameters

target_sparsity (int, optional (default 5)) – The maximum number of nonzero coefficients across all dimensions. If set, the model will fit all dimensions jointly, potentially reducing statistical efficiency.
group_sparsity (int tuple, optional (default None)) – Tuple of length n_targets constraining the number of nonzero coefficients for each target dimension.
alpha (float, optional (default 0.01)) – Optional L2 (ridge) regularization on the weight vector.
regression_timeout (int, optional (default 10)) – The timeout (in seconds) of the gurobi optimizer to solve and prove optimality (either per dimension or jointly depending on the above sparsity settings).
fit_intercept (boolean, optional (default False)) – Whether to calculate the intercept for this model. If set to false, no intercept will be used in calculations.
constraint_lhs (numpy ndarray, optional (default None)) – Shape should be (n_constraints, n_features * n_targets), The left hand side matrix C of Cw <= d. There should be one row per constraint.
constraint_rhs (numpy ndarray, shape (n_constraints,), optional (default None)) – The right hand side vector d of Cw <= d.
constraint_order (string, optional (default "target")) – The format in which the constraints constraint_lhs were passed. Must be one of “target” or “feature”. “target” indicates that the constraints are grouped by target: i.e. the first n_features columns correspond to constraint coefficients on the library features for the first target (variable), the next n_features columns to the library features for the second target (variable), and so on. “feature” indicates that the constraints are grouped by library feature: the first n_targets columns correspond to the first library feature, the next n_targets columns to the second library feature, and so on.
normalize_columns (boolean, optional (default False)) – Normalize the columns of x (the SINDy library terms) before regression by dividing by the L2-norm. Note that the ‘normalize’ option in sklearn is deprecated in sklearn versions >= 1.0 and will be removed. Note that this parameter is incompatible with the constraints!
copy_X (boolean, optional (default True)) – If True, X will be copied; else, it may be overwritten.
initial_guess (np.ndarray, shape (n_features) or (n_targets, n_features), optional (default None)) – Initial guess for coefficients coef_ to warmstart the optimizer.
verbose (bool, optional (default False)) – If True, prints out the Gurobi solver log.

Attributes

coef_ (array, shape (n_features,) or (n_targets, n_features)) – Weight vector(s).
ind_ (array, shape (n_features,) or (n_targets, n_features)) – Array of 0s and 1s indicating which coefficients of the weight vector have not been masked out, i.e. the support of self.coef_.
model (gurobipy.model) – The raw gurobi model being solved.

property complexity¶