feat(solvers): add tick solver

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

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}

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)