FEAT Add positivity penalty (#126)

Author: PABannier

skglm/penalties/__init__.py

@@ -1,6 +1,7 @@
 from .base import BasePenalty
 from .separable import (
-    L1_plus_L2, L0_5, L1, L2_3, MCPenalty, SCAD, WeightedL1, IndicatorBox
+    L1_plus_L2, L0_5, L1, L2_3, MCPenalty, SCAD, WeightedL1, IndicatorBox,
+    PositiveConstraint
 )
 from .block_separable import (
     L2_05, L2_1, BlockMCPenalty, BlockSCAD, WeightedGroupL2
@@ -12,5 +13,5 @@
 __all__ = [
     BasePenalty,
     L1_plus_L2, L0_5, L1, L2_3, MCPenalty, SCAD, WeightedL1, IndicatorBox,
-    L2_05, L2_1, BlockMCPenalty, BlockSCAD, WeightedGroupL2, SLOPE
+    PositiveConstraint, L2_05, L2_1, BlockMCPenalty, BlockSCAD, WeightedGroupL2, SLOPE
 ]

skglm/penalties/separable.py

@@ -466,3 +466,48 @@ def is_penalized(self, n_features):
     def generalized_support(self, w):
         """Return a mask with non-zero coefficients."""
         return w != 0
+
+
+class PositiveConstraint(BasePenalty):
+    """Positivity constraint penalty."""
+
+    def __init__(self):
+        pass
+
+    def get_spec(self):
+        return ()
+
+    def params_to_dict(self):
+        return dict()
+
+    def value(self, w):
+        """Compute the value of the PositiveConstraint penalty at w."""
+        return np.inf if (w < 0).any() else 0.
+
+    def prox_1d(self, value, stepsize, j):
+        """Compute the proximal operator of the PositiveConstraint."""
+        return max(0., value)
+
+    def subdiff_distance(self, w, grad, ws):
+        """Compute distance of negative gradient to the subdifferential at w."""
+        subdiff_dist = np.zeros_like(grad)
+        for idx, j in enumerate(ws):
+            if w[j] == 0:
+                # distance of - grad_j to  ]-infty, 0]
+                subdiff_dist[idx] = max(0, -grad[idx])
+            elif w[j] > 0:
+                # distance of - grad_j to 0
+                subdiff_dist[idx] = abs(-grad[idx])
+            else:
+                # subdiff is empty, distance is infinite
+                subdiff_dist[idx] = np.inf
+
+        return subdiff_dist
+
+    def is_penalized(self, n_features):
+        """Return a binary mask with the penalized features."""
+        return np.ones(n_features, bool_)
+
+    def generalized_support(self, w):
+        """Return a mask with non-zero coefficients."""
+        return w != 0

skglm/tests/test_penalties.py

@@ -4,10 +4,12 @@
 from numpy.linalg import norm
 from numpy.testing import assert_array_less
 
+from sklearn.linear_model import LinearRegression
+
 from skglm.datafits import Quadratic, QuadraticMultiTask
 from skglm.penalties import (
     L1, L1_plus_L2, WeightedL1, MCPenalty, SCAD, IndicatorBox, L0_5, L2_3, SLOPE,
-    L2_1, L2_05, BlockMCPenalty, BlockSCAD)
+    PositiveConstraint, L2_1, L2_05, BlockMCPenalty, BlockSCAD)
 from skglm import GeneralizedLinearEstimator, Lasso
 from skglm.solvers import AndersonCD, MultiTaskBCD, FISTA
 from skglm.utils.data import make_correlated_data
@@ -101,5 +103,20 @@ def test_slope():
     np.testing.assert_allclose(ours.coef_, pyslope_out["beta"], rtol=1e-5)
 
 
+@pytest.mark.parametrize("fit_intercept", [True, False])
+def test_nnls(fit_intercept):
+    # compare solutions with sklearn's LinearRegression, note that n_samples >=
+    # n_features for the design matrix to be injective, hence the solution unique
+    clf = GeneralizedLinearEstimator(
+        datafit=Quadratic(),
+        penalty=PositiveConstraint(),
+        solver=AndersonCD(tol=tol, fit_intercept=fit_intercept),
+    ).fit(X, y)
+    reg_nnls = LinearRegression(positive=True, fit_intercept=fit_intercept).fit(X, y)
+
+    np.testing.assert_allclose(clf.coef_, reg_nnls.coef_)
+    np.testing.assert_allclose(clf.intercept_, reg_nnls.intercept_)
+
+
 if __name__ == "__main__":
     pass