Skip to content

Commit

Permalink
Did all the things from the first JOSS review. Fixed the math notatio…
Browse files Browse the repository at this point in the history 
…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.
  • Loading branch information
Alan Kaptanoglu authored and Alan Kaptanoglu committed Jan 12, 2022
1 parent ff7f436 commit ec6a480
Show file tree
Hide file tree
Showing 11 changed files with 486 additions and 109 deletions.
136 changes: 91 additions & 45 deletions examples/1_feature_overview.ipynb
View file

Large diffs are not rendered by default.

98 changes: 76 additions & 22 deletions pysindy/optimizers/constrained_sr3.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
----------
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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,
Expand All @@ -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
Expand Down Expand Up @@ -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])
Expand All @@ -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):
"""
Expand Down Expand Up @@ -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(
Expand All @@ -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
Expand Down
48 changes: 42 additions & 6 deletions pysindy/optimizers/frols.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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,
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
24 changes: 21 additions & 3 deletions pysindy/optimizers/sindy_pi.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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)
Expand All @@ -101,6 +106,7 @@ def __init__(
thresholds=None,
model_subset=None,
normalize_columns=False,
verbose_cvxpy=False,
):
super(SINDyPI, self).__init__(
threshold=threshold,
Expand All @@ -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.")
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 "
Expand All @@ -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 "
Expand Down
Loading

0 comments on commit ec6a480

Please sign in to comment.