Got the inequality constraints working in the cvxpy subroutines of tr…

…apping SR3. Requires constraint_order=feature for now.

pysindy/optimizers/trapping_sr3.py

@@ -93,6 +93,10 @@ class TrappingSR3(SR3):
         If relax_optim = True, use the relax-and-split method. If False,
         try a direct minimization on the largest eigenvalue.
 
+    inequality_constraints : bool, optional
+        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.
 
@@ -185,6 +189,7 @@ def __init__(
         threshold=0.1,
         eps_solver=1.0e-7,
         relax_optim=True,
+        inequality_constraints=False,
         eta=1.0e20,
         alpha_A=1.0e20,
         alpha_m=5e19,
@@ -229,11 +234,20 @@ def __init__(
             raise ValueError("tol and tol_m must be positive")
         if not evolve_w and not relax_optim:
             raise ValueError("If doing direct solve, must evolve w")
+        if inequality_constraints and relax_optim and threshold == 0.0:
+            raise ValueError(
+                "Ineq. constr. -> threshold!=0 + relax_optim=True or relax_optim=False."
+            )
+        if inequality_constraints and not evolve_w:
+            raise ValueError(
+                "Use of inequality constraints requires solving for xi (evolve_w=True)."
+            )
 
         self.evolve_w = evolve_w
         self.threshold = threshold
         self.eps_solver = eps_solver
         self.relax_optim = relax_optim
+        self.inequality_constraints = inequality_constraints
         self.m0 = m0
         self.A0 = A0
         self.alpha_A = alpha_A
@@ -311,7 +325,7 @@ def _check_P_matrix(self, r, n_features, N):
         if (np.any(self.PL != 0.0) or np.any(self.PQ != 0.0)) and n_features != N:
             raise ValueError("Feature library is wrong shape or not quadratic")
 
-    def _update_coef_constraints(self, H, x_transpose_y, P_transpose_A):
+    def _update_coef_constraints(self, H, x_transpose_y, P_transpose_A, coef_sparse):
         g = x_transpose_y + P_transpose_A / self.eta
         inv1 = np.linalg.pinv(H, rcond=1e-15)
         inv2 = np.linalg.pinv(
@@ -370,10 +384,16 @@ def _solve_sparse_relax_and_split(self, r, N, x_expanded, y, Pmatrix, A, coef_pr
         ) + self.threshold * cp.norm1(xi)
         cost = cost + cp.sum_squares(Pmatrix @ xi - A.flatten()) / self.eta
         if self.use_constraints:
-            prob = cp.Problem(
-                cp.Minimize(cost),
-                [self.constraint_lhs @ xi == self.constraint_rhs],
-            )
+            if self.inequality_constraints:
+                prob = cp.Problem(
+                    cp.Minimize(cost),
+                    [self.constraint_lhs @ xi <= self.constraint_rhs],
+                )
+            else:
+                prob = cp.Problem(
+                    cp.Minimize(cost),
+                    [self.constraint_lhs @ xi == self.constraint_rhs],
+                )
         else:
             prob = cp.Problem(cp.Minimize(cost))
 
@@ -416,9 +436,9 @@ def _solve_m_relax_and_split(self, r, N, m_prev, m, A, coef_sparse, tk_previous)
 
     def _solve_nonsparse_relax_and_split(self, H, xTy, P_transpose_A, coef_prev):
         if self.use_constraints:
-            coef_sparse = self._update_coef_constraints(H, xTy, P_transpose_A).reshape(
-                coef_prev.shape
-            )
+            coef_sparse = self._update_coef_constraints(
+                H, xTy, P_transpose_A, coef_prev
+            ).reshape(coef_prev.shape)
         else:
             cho = cho_factor(H)
             coef_sparse = cho_solve(cho, xTy + P_transpose_A / self.eta).reshape(
@@ -433,10 +453,16 @@ def _solve_direct_cvxpy(self, r, N, x_expanded, y, Pmatrix, coef_prev):
         ) + self.threshold * cp.norm1(xi)
         cost = cost + cp.lambda_max(cp.reshape(Pmatrix @ xi, (r, r))) / self.eta
         if self.use_constraints:
-            prob = cp.Problem(
-                cp.Minimize(cost),
-                [self.constraint_lhs @ xi == self.constraint_rhs],
-            )
+            if self.inequality_constraints:
+                prob = cp.Problem(
+                    cp.Minimize(cost),
+                    [self.constraint_lhs @ xi <= self.constraint_rhs],
+                )
+            else:
+                prob = cp.Problem(
+                    cp.Minimize(cost),
+                    [self.constraint_lhs @ xi == self.constraint_rhs],
+                )
         else:
             prob = cp.Problem(cp.Minimize(cost))
 
@@ -448,23 +474,26 @@ def _solve_direct_cvxpy(self, r, N, x_expanded, y, Pmatrix, coef_prev):
             return None, None
         coef_sparse = (xi.value).reshape(coef_prev.shape)
 
-        m_cp = cp.Variable(r)
-        L = np.tensordot(self.PL, coef_sparse, axes=([3, 2], [0, 1]))
-        Q = np.reshape(
-            np.tensordot(self.PQ, coef_sparse, axes=([4, 3], [0, 1])), (r, r * r)
-        )
-        Ls = 0.5 * (L + L.T).flatten()
-        cost_m = cp.lambda_max(cp.reshape(Ls - m_cp @ Q, (r, r)))
-        prob_m = cp.Problem(cp.Minimize(cost_m))
+        if np.all(self.PL == 0.0) and np.all(self.PQ == 0.0):
+            return np.zeros(r), coef_sparse  # no optimization over m
+        else:
+            m_cp = cp.Variable(r)
+            L = np.tensordot(self.PL, coef_sparse, axes=([3, 2], [0, 1]))
+            Q = np.reshape(
+                np.tensordot(self.PQ, coef_sparse, axes=([4, 3], [0, 1])), (r, r * r)
+            )
+            Ls = 0.5 * (L + L.T).flatten()
+            cost_m = cp.lambda_max(cp.reshape(Ls - m_cp @ Q, (r, r)))
+            prob_m = cp.Problem(cp.Minimize(cost_m))
 
-        # default solver is SCS here
-        prob_m.solve(eps=self.eps_solver)
+            # default solver is SCS here
+            prob_m.solve(eps=self.eps_solver)
 
-        m = m_cp.value
-        if m is None:
-            print("Infeasible solve over m, increase/decrease eta")
-            return None, None
-        return m, coef_sparse
+            m = m_cp.value
+            if m is None:
+                print("Infeasible solve over m, increase/decrease eta")
+                return None, coef_sparse
+            return m, coef_sparse
 
     def _reduce(self, x, y):
         """

test/optimizers/test_optimizers.py

@@ -4,13 +4,16 @@
 import numpy as np
 import pytest
 from numpy.linalg import norm
+from scipy.integrate import odeint
 from sklearn.base import BaseEstimator
 from sklearn.exceptions import ConvergenceWarning
 from sklearn.exceptions import NotFittedError
 from sklearn.linear_model import ElasticNet
 from sklearn.linear_model import Lasso
 from sklearn.utils.validation import check_is_fitted
 
+from pysindy import FiniteDifference
+from pysindy import PolynomialLibrary
 from pysindy import SINDy
 from pysindy.feature_library import CustomLibrary
 from pysindy.optimizers import ConstrainedSR3
@@ -21,6 +24,10 @@
 from pysindy.utils import supports_multiple_targets
 
 
+def lorenz(z, t):
+    return 10 * (z[1] - z[0]), z[0] * (28 - z[2]) - z[1], z[0] * z[1] - 8 / 3 * z[2]
+
+
 class DummyLinearModel(BaseEstimator):
     # Does not natively support multiple targets
     def fit(self, x, y):
@@ -483,3 +490,58 @@ def test_target_format_constraints(data_linear_combination, optimizer, target_va
     model.fit(x, x_dot)
     print(model.coef_, target_value)
     np.testing.assert_allclose(model.coef_[:, 1], target_value, atol=1e-8)
+
+
+@pytest.mark.parametrize("thresholds", [0.005, 0.05, 0.5])
+@pytest.mark.parametrize("relax_optim", [False, True])
+@pytest.mark.parametrize("noise_levels", [0.0, 0.5, 2.0])
+def test_trapping_inequality_constraints(thresholds, relax_optim, noise_levels):
+    t = np.arange(0, 40, 0.05)
+    x = odeint(lorenz, [-8, 8, 27], t)
+    x = x + np.random.normal(0.0, noise_levels, x.shape)
+    # if order is "feature"
+    constraint_rhs = np.array([-10.0, -2.0])
+    constraint_matrix = np.zeros((2, 30))
+    constraint_matrix[0, 6] = 1.0
+    constraint_matrix[1, 17] = 1.0
+    feature_names = ["x", "y", "z"]
+    opt = TrappingSR3(
+        threshold=thresholds,
+        constraint_lhs=constraint_matrix,
+        constraint_rhs=constraint_rhs,
+        constraint_order="feature",
+        inequality_constraints=True,
+        relax_optim=relax_optim,
+    )
+    poly_lib = PolynomialLibrary(degree=2)
+    model = SINDy(
+        optimizer=opt,
+        feature_library=poly_lib,
+        differentiation_method=FiniteDifference(drop_endpoints=True),
+        feature_names=feature_names,
+    )
+    model.fit(x, t=t[1] - t[0])
+    print(
+        np.dot(constraint_matrix, (model.coefficients()).flatten("F")), constraint_rhs
+    )
+    assert np.all(
+        np.dot(constraint_matrix, (model.coefficients()).flatten("F")) <= constraint_rhs
+    ) or np.allclose(
+        np.dot(constraint_matrix, (model.coefficients()).flatten("F")), constraint_rhs
+    )
+
+
+def test_inequality_constraints_reqs():
+    constraint_rhs = np.array([-10.0, -2.0])
+    constraint_matrix = np.zeros((2, 30))
+    constraint_matrix[0, 6] = 1.0
+    constraint_matrix[1, 17] = 1.0
+    with pytest.raises(ValueError):
+        TrappingSR3(
+            threshold=0.0,
+            constraint_lhs=constraint_matrix,
+            constraint_rhs=constraint_rhs,
+            constraint_order="feature",
+            inequality_constraints=True,
+            relax_optim=True,
+        )