FEAT Add Gamma regressor (#113)

Author: PABannier

doc/api.rst

@@ -54,9 +54,11 @@ Datafits
 .. autosummary::
    :toctree: generated/
 
+   Gamma
    Huber
    Logistic
    LogisticGroup
+   Poisson
    Quadratic
    QuadraticGroup
    QuadraticSVC
@@ -87,4 +89,4 @@ Experimental
 .. autosummary::
    :toctree: generated/
 
-   SqrtLasso
+   SqrtLasso

skglm/datafits/__init__.py

@@ -1,12 +1,12 @@
 from .base import BaseDatafit, BaseMultitaskDatafit
-from .single_task import Quadratic, QuadraticSVC, Logistic, Huber, Poisson
+from .single_task import Quadratic, QuadraticSVC, Logistic, Huber, Poisson, Gamma
 from .multi_task import QuadraticMultiTask
 from .group import QuadraticGroup, LogisticGroup
 
 
 __all__ = [
     BaseDatafit, BaseMultitaskDatafit,
-    Quadratic, QuadraticSVC, Logistic, Huber, Poisson,
+    Quadratic, QuadraticSVC, Logistic, Huber, Poisson, Gamma,
     QuadraticMultiTask,
     QuadraticGroup, LogisticGroup
 ]

skglm/datafits/single_task.py

@@ -429,10 +429,16 @@ def params_to_dict(self):
         return dict()
 
     def initialize(self, X, y):
-        pass
+        if np.any(y <= 0):
+            raise ValueError(
+                "Target vector `y` should only take positive values " +
+                "when fitting a Poisson model.")
 
     def initialize_sparse(self, X_data, X_indptr, X_indices, y):
-        pass
+        if np.any(y <= 0):
+            raise ValueError(
+                "Target vector `y` should only take positive values " +
+                "when fitting a Poisson model.")
 
     def raw_grad(self, y, Xw):
         """Compute gradient of datafit w.r.t ``Xw``."""
@@ -467,3 +473,58 @@ def gradient_scalar_sparse(self, X_data, X_indptr, X_indices, y, Xw, j):
 
     def intercept_update_self(self, y, Xw):
         pass
+
+
+class Gamma(BaseDatafit):
+    r"""Gamma datafit.
+
+    The datafit reads::
+
+    (1 / n_samples) * \sum_i (Xw_i + y_i * exp(-Xw_i) - 1 - log(y_i))
+
+    Note:
+    ----
+    The class is jit compiled at fit time using Numba compiler.
+    This allows for faster computations.
+    """
+
+    def __init__(self):
+        pass
+
+    def get_spec(self):
+        pass
+
+    def params_to_dict(self):
+        return dict()
+
+    def initialize(self, X, y):
+        if np.any(y <= 0):
+            raise ValueError(
+                "Target vector `y` should only take positive values "
+                "when fitting a Gamma model.")
+
+    def initialize_sparse(self, X_data, X_indptr, X_indices, y):
+        if np.any(y <= 0):
+            raise ValueError(
+                "Target vector `y` should only take positive values "
+                "when fitting a Gamma model.")
+
+    def raw_grad(self, y, Xw):
+        """Compute gradient of datafit w.r.t. ``Xw``."""
+        return (1 - y * np.exp(-Xw)) / len(y)
+
+    def raw_hessian(self, y, Xw):
+        """Compute Hessian of datafit w.r.t. ``Xw``."""
+        return (y * np.exp(-Xw)) / len(y)
+
+    def value(self, y, w, Xw):
+        return (np.sum(Xw + y * np.exp(-Xw) - np.log(y)) - 1) / len(y)
+
+    def gradient_scalar(self, X, y, w, Xw, j):
+        return X[:, j] @ (1 - y * np.exp(-Xw)) / len(y)
+
+    def gradient_scalar_sparse(self, X_data, X_indptr, X_indices, y, Xw, j):
+        pass
+
+    def intercept_update_self(self, y, Xw):
+        pass

skglm/tests/test_datafits.py

@@ -4,7 +4,7 @@
 from sklearn.linear_model import HuberRegressor
 from numpy.testing import assert_allclose, assert_array_less
 
-from skglm.datafits import Huber, Logistic, Poisson
+from skglm.datafits import Huber, Logistic, Poisson, Gamma
 from skglm.penalties import L1, WeightedL1
 from skglm.solvers import AndersonCD, ProxNewton
 from skglm import GeneralizedLinearEstimator
@@ -84,5 +84,35 @@ def test_poisson():
     assert_allclose(model.coef_, w_statsmodels, rtol=1e-4)
 
 
+def test_gamma():
+    try:
+        import statsmodels.api as sm
+    except ImportError:
+        pytest.xfail("`statsmodels` not found. `Gamma` datafit can't be tested.")
+
+    # When n_samples < n_features, the unregularized Gamma objective does not have a
+    # unique minimizer.
+    rho = 1e-2
+    n_samples, n_features = 100, 10
+    X, y, _ = make_correlated_data(n_samples, n_features, random_state=0)
+    y[y <= 0] = 0.1
+    tol = 1e-14
+
+    alpha_max = np.linalg.norm(X.T @ (1 - y), ord=np.inf) / n_samples
+    alpha = rho * alpha_max
+
+    gamma_model = sm.GLM(y, X, family=sm.families.Gamma(sm.families.links.Log()))
+    gamma_results = gamma_model.fit_regularized(
+        method="elastic_net", L1_wt=1, cnvrg_tol=tol, alpha=alpha)
+
+    clf = GeneralizedLinearEstimator(
+        datafit=Gamma(),
+        penalty=L1(alpha),
+        solver=ProxNewton(fit_intercept=False, tol=tol)
+    ).fit(X, y)
+
+    np.testing.assert_allclose(clf.coef_, gamma_results.params, rtol=1e-6)
+
+
 if __name__ == '__main__':
     pass