diff --git a/example_config.yml b/example_config.yml index 1f58bd3..99dd378 100644 --- a/example_config.yml +++ b/example_config.yml @@ -7,6 +7,7 @@ dataset: solver: - PGD[prox=prox_isotonic] - sortedl1[cd_type=permuted] + - tick # note: use more repetitions to have quantiles on the plot n-repetitions: 1 diff --git a/solvers/tick.py b/solvers/tick.py new file mode 100644 index 0000000..cf0a703 --- /dev/null +++ b/solvers/tick.py @@ -0,0 +1,111 @@ +from typing import ClassVar + +from benchopt import BaseSolver, safe_import_context + +with safe_import_context() as import_ctx: + import numpy as np + from scipy import sparse + from scipy.optimize import brentq + from scipy.stats import norm + from tick.linear_model import ModelLinReg + from tick.prox import ProxSlope + from tick.solver import AGD + + +def _slope_parameters(weights): + weights = np.asarray(weights, dtype=np.float64) + if weights.ndim != 1 or len(weights) == 0 or not np.all(np.isfinite(weights)): + return None + if np.all(weights == 0): + return 0.1, 0.0 + if weights[-1] <= 0 or weights[0] < weights[-1]: + return None + + p = len(weights) + if p == 1: + q = 0.1 + else: + ratio = weights[0] / weights[-1] + + def ratio_error(candidate): + return norm.isf(candidate / (2 * p)) / norm.isf(candidate / 2) - ratio + + lower, upper = 1e-12, 1 - 1e-12 + if ratio_error(lower) > 0 or ratio_error(upper) < 0: + return None + q = brentq(ratio_error, lower, upper, xtol=1e-14) + + bh_weights = norm.isf(q * np.arange(1, p + 1) / (2 * p)) + strength = weights[0] / bh_weights[0] + if not np.allclose(weights, strength * bh_weights, rtol=1e-6, atol=1e-12): + return None + return q, strength + + +class Solver(BaseSolver): + name = "tick" + sampling_strategy = "iteration" + install_cmd = "conda" + requirements: ClassVar[list[str]] = [ + "numpy", + "scipy", + "numpydoc", + "pip::tick>=0.8.0.2", + ] + references: ClassVar[list[str]] = [ + ( + "Bacry, E., Bompaire, M., Gaiffas, S., & Poulsen, S. V. (2017). " + "tick: a Python library for statistical learning, with a particular " + "emphasis on time-dependent modeling. arXiv:1707.03003." + ) + ] + + def skip(self, X, y, alphas, fit_intercept): + if _slope_parameters(alphas) is None: + return True, "tick requires Benjamini-Hochberg SLOPE weights" + return False, None + + def set_objective(self, X, y, alphas, fit_intercept): + self.n_features = X.shape[1] + self.fit_intercept = fit_intercept + parameters = _slope_parameters(alphas) + if parameters is None: + raise ValueError("tick requires Benjamini-Hochberg SLOPE weights") + q, strength = parameters + + if sparse.issparse(X): + X = sparse.csr_matrix(X, dtype=np.float64) + else: + X = np.ascontiguousarray(X, dtype=np.float64) + y = np.ascontiguousarray(y, dtype=np.float64) + model = ModelLinReg(fit_intercept=fit_intercept).fit(X, y) + penalty_range = (0, self.n_features) if fit_intercept else None + prox = ProxSlope(strength=strength, fdr=q, range=penalty_range) + self.solver = AGD( + # The average sample Lipschitz constant bounds the full loss and + # works for sparse designs, unlike tick's dense-only best bound. + step=1 / model.get_lip_mean(), + max_iter=1, + tol=0.0, + linesearch=False, + verbose=False, + ) + self.solver.set_model(model).set_prox(prox) + self.n_coeffs = model.n_coeffs + + @staticmethod + def get_next(previous): + return max(1, 2 * previous) + + def run(self, n_iter): + self.coef = np.zeros(self.n_coeffs) + if n_iter: + self.solver.max_iter = n_iter + self.coef = self.solver.solve(x0=self.coef) + + def get_result(self): + if self.fit_intercept: + beta = np.concatenate(([self.coef[-1]], self.coef[:-1])) + else: + beta = np.concatenate(([0.0], self.coef)) + return {"beta": beta} diff --git a/tests/test_tick.py b/tests/test_tick.py new file mode 100644 index 0000000..9faccf0 --- /dev/null +++ b/tests/test_tick.py @@ -0,0 +1,81 @@ +import numpy as np +import pytest +from scipy import sparse +from scipy.stats import norm + +from objective import Objective +from solvers.tick import Solver, _slope_parameters + +pytest.importorskip("tick") + + +@pytest.mark.parametrize("fit_intercept", [False, True]) +@pytest.mark.parametrize("sparse_input", [False, True]) +def test_tick_solves_bh_problem(fit_intercept, sparse_input): + X = np.array([[-2.0, 1.0], [-1.0, 0.0], [0.0, 1.0], [1.0, 0.0], [2.0, 1.0]]) + if sparse_input: + X = sparse.csc_array(X) + y = np.arange(5.0) + q = 0.2 + penalty_strength = 0.05 + alphas = penalty_strength * norm.ppf(1 - q * np.arange(1, 3) / 4) + + solver = Solver() + assert solver.skip(X, y, alphas, fit_intercept) == (False, None) + solver.set_objective(X, y, alphas, fit_intercept) + solver.run(0) + np.testing.assert_array_equal(solver.get_result()["beta"], np.zeros(3)) + + solver.run(300) + result = solver.get_result() + objective = Objective(q=q, fit_intercept=fit_intercept) + objective.X, objective.y, objective.alphas = X, y, alphas + assert objective.evaluate_result(**result)["rel_duality_gap"] < 1e-6 + if fit_intercept: + np.testing.assert_allclose(result["beta"][0], 2.0, atol=1e-6) + np.testing.assert_allclose( + solver.solver.prox.value(np.array([1.0, 0.5, 7.0])), + alphas @ np.array([1.0, 0.5]), + ) + else: + assert result["beta"][0] == 0.0 + + +def test_tick_skips_non_bh_weights(): + solver = Solver() + skip, reason = solver.skip(None, None, np.ones(2), False) + assert skip + assert "Benjamini-Hochberg" in reason + + +def test_tick_accepts_generated_objective(): + rng = np.random.default_rng(42) + X = rng.normal(size=(30, 5)) + y = 1.5 + X[:, 0] - X[:, 1] + objective = Objective(q=0.2, reg=0.1, fit_intercept=True) + objective.set_data(X, y) + assert set(objective.get_objective()) == { + "X", + "y", + "alphas", + "fit_intercept", + } + + solver = Solver() + assert solver._set_objective(objective) == (False, None) + solver.run(1) + assert solver.get_result()["beta"].shape == (6,) + + +@pytest.mark.parametrize("q", [0.05, 0.2]) +@pytest.mark.parametrize("n_features", [1, 5, 200]) +def test_tick_recovers_penalty_from_weights(q, n_features): + weights = 0.05 * norm.ppf(1 - q * np.arange(1, n_features + 1) / (2 * n_features)) + recovered_q, strength = _slope_parameters(weights) + reconstructed = strength * norm.ppf( + 1 - recovered_q * np.arange(1, n_features + 1) / (2 * n_features) + ) + np.testing.assert_allclose(reconstructed, weights, rtol=1e-7) + if n_features > 1: + np.testing.assert_allclose(recovered_q, q, rtol=1e-7) + np.testing.assert_allclose(strength, 0.05, rtol=1e-7)