Did all the things from the first JOSS review. Fixed the math notatio…

…n in the optimizer descriptions. Added a verbose option for all the optimizers, tested it, added unit tests, and added simple example in notebook 1. Standardized a bunch of variables between SR3 variants. Fixed the error in test_complexity.

pysindy/optimizers/constrained_sr3.py

@@ -18,22 +18,22 @@ class ConstrainedSR3(SR3):
 
     .. math::
 
-        0.5\\|y-Xw\\|^2_2 + \\lambda \\times R(v)
-        + (0.5 / nu)\\|w-v\\|^2_2
+        0.5\\|y-Xw\\|^2_2 + \\lambda R(u)
+        + (0.5 / \\nu)\\|w-u\\|^2_2
 
         \\text{subject to}
 
     .. math::
 
-        Cw = d
+        CW = d
 
-    over v and w where :math:`R(v)` is a regularization function, C is a
+    over U and W, where :math:`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."
-        arXiv preprint arXiv:1906.10612 (2019).
+        IEEE Access 8 (2020): 169259-169271.
 
     Parameters
     ----------
@@ -115,6 +115,15 @@ class ConstrainedSR3(SR3):
     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)
@@ -150,6 +159,8 @@ def __init__(
         initial_guess=None,
         thresholds=None,
         inequality_constraints=False,
+        verbose=False,
+        verbose_cvxpy=False,
     ):
         super(ConstrainedSR3, self).__init__(
             threshold=threshold,
@@ -164,8 +175,10 @@ def __init__(
             fit_intercept=fit_intercept,
             copy_X=copy_X,
             normalize_columns=normalize_columns,
+            verbose=verbose,
         )
 
+        self.verbose_cvxpy = verbose_cvxpy
         self.reg = get_regularization(thresholder)
         self.use_constraints = (constraint_lhs is not None) and (
             constraint_rhs is not None
@@ -243,13 +256,18 @@ def _update_coef_cvxpy(self, x, y, coef_sparse):
 
         # default solver is OSQP here but switches to ECOS for L2
         try:
-            prob.solve(max_iter=self.max_iter, eps_abs=self.tol, eps_rel=self.tol)
+            prob.solve(
+                max_iter=self.max_iter,
+                eps_abs=self.tol,
+                eps_rel=self.tol,
+                verbose=self.verbose_cvxpy,
+            )
         # Annoying error coming from L2 norm switching to use the ECOS
         # solver, which uses "max_iters" instead of "max_iter", and
         # similar semantic changes for the other variables.
         except TypeError:
             try:
-                prob.solve(abstol=self.tol, reltol=self.tol)
+                prob.solve(abstol=self.tol, reltol=self.tol, verbose=self.verbose_cvxpy)
             except cp.error.SolverError:
                 print("Solver failed, setting coefs to zeros")
                 xi.value = np.zeros(coef_sparse.shape[0] * coef_sparse.shape[1])
@@ -276,25 +294,48 @@ def _update_sparse_coef(self, coef_full):
         self.history_.append(coef_sparse.T)
         return coef_sparse
 
-    def _objective(self, x, y, coef_full, coef_sparse, trimming_array=None):
+    def _objective(self, x, y, q, coef_full, coef_sparse, trimming_array=None):
         """Objective function"""
+        if q != 0:
+            print_ind = q % (self.max_iter // 10.0)
+        else:
+            print_ind = q
         R2 = (y - np.dot(x, coef_full)) ** 2
         D2 = (coef_full - coef_sparse) ** 2
         if self.use_trimming:
             assert trimming_array is not None
             R2 *= trimming_array.reshape(x.shape[0], 1)
+
         if self.thresholds is None:
-            return (
-                0.5 * np.sum(R2)
-                + self.reg(coef_full, 0.5 * self.threshold ** 2 / self.nu)
-                + 0.5 * np.sum(D2) / self.nu
-            )
+            regularization = self.reg(coef_full, self.threshold ** 2 / self.nu)
+            if print_ind == 0 and self.verbose:
+                row = [
+                    q,
+                    np.sum(R2),
+                    np.sum(D2) / self.nu,
+                    regularization,
+                    np.sum(R2) + np.sum(D2) + regularization,
+                ]
+                print(
+                    "{0:10d} ... {1:10.4e} ... {2:10.4e} ... {3:10.4e}"
+                    " ... {4:10.4e}".format(*row)
+                )
+            return 0.5 * np.sum(R2) + 0.5 * regularization + 0.5 * np.sum(D2) / self.nu
         else:
-            return (
-                0.5 * np.sum(R2)
-                + self.reg(coef_full, 0.5 * self.thresholds ** 2 / self.nu)
-                + 0.5 * np.sum(D2) / self.nu
-            )
+            regularization = self.reg(coef_full, self.thresholds ** 2 / self.nu)
+            if print_ind == 0 and self.verbose:
+                row = [
+                    q,
+                    np.sum(R2),
+                    np.sum(D2) / self.nu,
+                    regularization,
+                    np.sum(R2) + np.sum(D2) + regularization,
+                ]
+                print(
+                    "{0:10d} ... {1:10.4e} ... {2:10.4e} ... {3:10.4e}"
+                    " ... {4:10.4e}".format(*row)
+                )
+            return 0.5 * np.sum(R2) + 0.5 * regularization + 0.5 * np.sum(D2) / self.nu
 
     def _reduce(self, x, y):
         """
@@ -333,12 +374,25 @@ def _reduce(self, x, y):
                 x_expanded, (n_samples * n_targets, n_targets * n_features)
             )
 
+        # Print initial values for each term in the optimization
+        if self.verbose:
+            row = [
+                "Iteration",
+                "|y - Xw|^2",
+                "|w-u|^2/v",
+                "R(u)",
+                "Total Error: |y - Xw|^2 + |w - u|^2 / v + R(u)",
+            ]
+            print(
+                "{: >10} ... {: >10} ... {: >10} ... {: >10} ... {: >10}".format(*row)
+            )
+
         objective_history = []
         if self.inequality_constraints:
             coef_sparse = self._update_coef_cvxpy(x_expanded, y, coef_sparse)
-            objective_history.append(self._objective(x, y, coef_full, coef_sparse))
+            objective_history.append(self._objective(x, y, 0, coef_full, coef_sparse))
         else:
-            for _ in range(self.max_iter):
+            for k in range(self.max_iter):
                 if self.use_trimming:
                     x_weighted = x * trimming_array.reshape(n_samples, 1)
                     H = np.dot(x_weighted.T, x) + np.diag(
@@ -363,11 +417,11 @@ def _reduce(self, x, y):
                     )
 
                     objective_history.append(
-                        self._objective(x, y, coef_full, coef_sparse, trimming_array)
+                        self._objective(x, y, k, coef_full, coef_sparse, trimming_array)
                     )
                 else:
                     objective_history.append(
-                        self._objective(x, y, coef_full, coef_sparse)
+                        self._objective(x, y, k, coef_full, coef_sparse)
                     )
                 if self._convergence_criterion() < self.tol:
                     # TODO: Update this for trimming/constraints

pysindy/optimizers/frols.py

@@ -46,6 +46,10 @@ class FROLS(BaseOptimizer):
     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)
@@ -84,6 +88,7 @@ def __init__(
         max_iter=10,
         alpha=0.05,
         ridge_kw=None,
+        verbose=False,
     ):
         super(FROLS, self).__init__(
             fit_intercept=fit_intercept,
@@ -96,6 +101,7 @@ def __init__(
         self.kappa = kappa
         if self.max_iter <= 0:
             raise ValueError("Max iteration must be > 0")
+        self.verbose = verbose
 
     def _normed_cov(self, a, b):
         return np.vdot(a, b) / np.vdot(a, a)
@@ -135,6 +141,31 @@ def _reduce(self, x, y):
         """
         n_samples, n_features = x.shape
         n_targets = y.shape[1]
+        cond_num = np.linalg.cond(x)
+        if self.kappa is not None:
+            l0_penalty = self.kappa * cond_num
+        else:
+            l0_penalty = 0.0
+
+        # Print initial values for each term in the optimization
+        if self.verbose:
+            row = [
+                "Iteration",
+                "Index",
+                "|y - Xw|^2",
+                "a * |w|_2",
+                "|w|_0",
+                "b * |w|_0",
+                "Total: |y-Xw|^2+a*|w|_2+b*|w|_0",
+            ]
+            print(
+                "{: >10} ... {: >5} ... {: >10} ... {: >10} ... {: >5}"
+                " ... {: >10} ... {: >10}".format(*row)
+            )
+            print(
+                " Note that these error metrics are related but not the same"
+                " as the loss function being minimized by FROLS!"
+            )
 
         # History of selected functions: [iteration x output x coefficients]
         self.history_ = np.zeros((n_features, n_targets, n_features), dtype=x.dtype)
@@ -191,14 +222,19 @@ def _reduce(self, x, y):
 
                 if i >= self.max_iter:
                     break
-
-        if self.kappa is not None:
-            l0_penalty = self.kappa * np.linalg.cond(x)
-        else:
-            l0_penalty = 0.0
+                if self.verbose:
+                    coef = self.history_[i, k, :]
+                    R2 = np.sum((y[:, k] - np.dot(x, coef).T) ** 2)
+                    L2 = self.alpha * np.sum(coef ** 2)
+                    L0 = np.count_nonzero(coef)
+                    row = [i, k, R2, L2, L0, l0_penalty * L0, R2 + L2 + l0_penalty * L0]
+                    print(
+                        "{0:10d} ... {1:5d} ... {2:10.4e} ... {3:10.4e}"
+                        " ... {4:5d} ... {5:10.4e} ... {6:10.4e}".format(*row)
+                    )
 
         # Function selection: L2 error for output k at iteration i is given by
-        #    sum(ERR_global[k, :i]), and the number of nonzero coefficients is (i+1)
+        # sum(ERR_global[k, :i]), and the number of nonzero coefficients is (i+1)
         l2_err = np.cumsum(self.ERR_global, axis=1)
         l0_norm = np.arange(1, n_features + 1)[:, None]
         self.loss_ = l2_err + l0_penalty * l0_norm  # Save for debugging

pysindy/optimizers/sindy_pi.py

@@ -15,7 +15,7 @@ class SINDyPI(SR3):
 
     .. math::
 
-        0.5\\|X-Xw\\|^2_2 + \\lambda \\times R(v)
+        0.5\\|X-Xw\\|^2_2 + \\lambda R(w)
 
     over w where :math:`R(v)` is a regularization function. See the following
     reference for more details:
@@ -76,6 +76,11 @@ class SINDyPI(SR3):
         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)
@@ -101,6 +106,7 @@ def __init__(
         thresholds=None,
         model_subset=None,
         normalize_columns=False,
+        verbose_cvxpy=False,
     ):
         super(SINDyPI, self).__init__(
             threshold=threshold,
@@ -125,6 +131,7 @@ def __init__(
             )
 
         self.unbias = False
+        self.verbose_cvxpy = verbose_cvxpy
         if model_subset is not None:
             if not isinstance(model_subset, list):
                 raise ValueError("Model subset must be in list format.")
@@ -155,6 +162,7 @@ def _update_parallel_coef_constraints(self, x):
                 "of features in the candidate library"
             )
         for i in self.model_subset:
+            print("Model ", i)
             xi = cp.Variable(n_features)
             # Note that norm choice below must be convex,
             # so thresholder must be L1 or L2
@@ -183,7 +191,12 @@ def _update_parallel_coef_constraints(self, x):
                 [xi[i] == 0.0],
             )
             try:
-                prob.solve(max_iter=self.max_iter, eps_abs=self.tol, eps_rel=self.tol)
+                prob.solve(
+                    max_iter=self.max_iter,
+                    eps_abs=self.tol,
+                    eps_rel=self.tol,
+                    verbose=self.verbose_cvxpy,
+                )
                 if xi.value is None:
                     warnings.warn(
                         "Infeasible solve on iteration " + str(i) + ", try "
@@ -195,7 +208,12 @@ def _update_parallel_coef_constraints(self, x):
             # solver, which uses "max_iters" instead of "max_iter", and
             # similar semantic changes for the other variables.
             except TypeError:
-                prob.solve(max_iters=self.max_iter, abstol=self.tol, reltol=self.tol)
+                prob.solve(
+                    max_iters=self.max_iter,
+                    abstol=self.tol,
+                    reltol=self.tol,
+                    verbose=self.verbose_cvxpy,
+                )
                 if xi.value is None:
                     warnings.warn(
                         "Infeasible solve on iteration " + str(i) + ", try "