TICK-367 Add Oscar penalization

doc/modules/api.rst

@@ -146,6 +146,7 @@ This module contains all the proximal operators available in tick.
    optim.prox.ProxPositive
    optim.prox.ProxEquality
    optim.prox.ProxSlope
+   optim.prox.ProxOscar
    optim.prox.ProxTV
    optim.prox.ProxBinarsity
 

setup.py

@@ -535,6 +535,7 @@ def add_dir_name(dir_name, filenames):
                   "prox_multi.cpp",
                   "prox_equality.cpp",
                   "prox_slope.cpp",
+                  "prox_oscar.cpp",
                   "prox_binarsity.cpp"],
     "h_files": ["prox.h",
                 "prox_separable.h",
@@ -549,6 +550,7 @@ def add_dir_name(dir_name, filenames):
                 "prox_multi.h",
                 "prox_equality.h",
                 "prox_slope.h",
+                "prox_oscar.h",
                 "prox_binarsity.h"],
     "swig_files": ["prox_module.i"],
     "module_dir": "./tick/optim/prox/",

tick/optim/prox/__init__.py

@@ -10,6 +10,7 @@
 from .prox_tv import ProxTV
 from .prox_nuclear import ProxNuclear
 from .prox_slope import ProxSlope
+from .prox_oscar import ProxOscar
 from .prox_elasticnet import ProxElasticNet
 from .prox_multi import ProxMulti
 from .prox_equality import ProxEquality
@@ -24,7 +25,8 @@
            "ProxTV",
            "ProxNuclear",
            "ProxSlope",
+           "ProxOscar",
            "ProxElasticNet",
            "ProxMulti",
            "ProxEquality",
-           "ProxBinarsity"]
+           "ProxBinarsity"]

tick/optim/prox/prox_oscar.py

@@ -1,52 +1,43 @@
 # License: BSD 3 clause
 
-# -*- coding: utf8 -*-
-
 from tick.optim.prox.base import Prox
 import numpy as np
 
-from tick.optim.prox.build.prox import ProxSortedL1 as _ProxSortedL1
-from tick.optim.prox.build.prox import WeightsType_bh, \
-    WeightsType_oscar
-
-
-# TODO: put also the OSCAR weights
-# TODO: we should be able to put any weights we want...
+from tick.optim.prox.build.prox import ProxOscar as _ProxOscar
 
 
-class ProxSortedL1(Prox):
-    """Proximal operator of sorted L1 penalization
+class ProxOscar(Prox):
+    """Proximal operator of the OSCAR penalization.
+    This penalization combines L1 penalization with a clustering penalization,
+    that induces exact equality of weights corresponding to correlated features,
+    so that clusters can be represented by a single coefficient.
+    This penalization is therefore particularly relevant in high-dimensional
+    problems with strong features correlation.
 
     Parameters
     ----------
     strength : `float`
         Level of penalization
 
-    fdr : `float`, default=0.6
-        Desired False Discovery Rate for detection of non-zeros in
-        the coefficients
-
-    weights_type : "bh" | "oscar", default="bh"
-
-        * If "bh", we use Benjamini-Hochberg weights, under a Gaussian
-          error assumption, and expect a FDR control
-        * If "oscar", there is no FDR control, and we get the OSCAR
-          penalization, see notes below for references
+    ratio : `float`
+        The Oscar ratio parameter, with ratio >= 0. For ratio = 0 this is L1
+        regularization, while a large ratio provides only the clustering effect.
 
     range : `tuple` of two `int`, default=`None`
         Range on which the prox is applied. If `None` then the prox is
         applied on the whole vector
 
-    Attributes
-    ----------
-    weights : `np.array`, shape=(n_coeffs,)
-        The weights used in the penalization. They are automatically
-        setted, depending on the ``weights_type`` and ``fdr``
-        parameters.
+    positive : `bool`, default=`False`
+        If True, apply an extra projection onto the set of vectors with
+        non-negative entries
 
     Notes
     -----
-    Uses the stack-based algorithm for FastProxL1 from
+    This penalization was introduced in
+    * Simultaneous regression shrinkage, variable selection and clustering of
+      predictors with OSCAR, by Bondell H.D. and Reich B.J., Biometrics. 2008
+
+    It uses the stack-based algorithm for FastProxL1 from
 
     * SLOPE--Adaptive Variable Selection via Convex Optimization, by
       Bogdan, M. and Berg, E. van den and Sabatti, C. and Su, W. and Candes, E. J.
@@ -58,59 +49,27 @@ class ProxSortedL1(Prox):
             "writable": True,
             "cpp_setter": "set_strength"
         },
-        "fdr": {
+        "ratio": {
             "writable": True,
-            "cpp_setter": "set_fdr"
-        },
-        "_weights_type": {
-            "writable": False,
-            "cpp_setter": "set_weights_type"
+            "cpp_setter": "set_ratio"
         },
         "positive": {
             "writable": True,
             "cpp_setter": "set_positive"
-        },
-        "weights": {
-            "writable": False,
         }
     }
 
-    def __init__(self, strength: float, fdr: float=0.6,
-                 weights_type: str="bh", range: tuple=None,
+    def __init__(self, strength: float, ratio: float, range: tuple=None,
                  positive: bool=False):
         Prox.__init__(self, range)
         self.strength = strength
-        self.fdr = fdr
-        self.weights_type = weights_type
+        self.ratio = ratio
         self.positive = positive
-        self.weights = None
         if range is None:
-            self._prox = _ProxSortedL1(self.strength, self.fdr,
-                                       self._weights_type,
-                                       self.positive)
+            self._prox = _ProxOscar(self.strength, self.ratio, self.positive)
         else:
-            self._prox = _ProxSortedL1(self.strength, self.fdr,
-                                       self._weights_type,
-                                       self.range[0], self.range[1],
-                                       self.positive)
-
-    @property
-    def weights_type(self):
-        if self._weights_type == WeightsType_bh:
-            return "bh"
-        elif self._weights_type == WeightsType_oscar:
-            return "oscar"
-
-    @weights_type.setter
-    def weights_type(self, val):
-        if val == "bh":
-            self._set("_weights_type", WeightsType_bh)
-        elif val == "oscar":
-            self._set("_weights_type", WeightsType_oscar)
-            raise NotImplementedError("``oscar`` weights.")
-        else:
-            raise ValueError("``weights_type`` must be either 'bh' "
-                             "or 'oscar'")
+            self._prox = _ProxOscar(self.strength, self.ratio,
+                                    self.range[0], self.range[1], self.positive)
 
     def _call(self, coeffs: np.ndarray, t: float, out: np.ndarray):
         self._prox.call(coeffs, t, out)
@@ -130,8 +89,3 @@ def value(self, coeffs: np.ndarray) -> float:
             Value of the penalization at ``coeffs``
         """
         return self._prox.value(coeffs)
-
-    def _as_dict(self):
-        dd = Prox._as_dict(self)
-        del dd["weights"]
-        return dd

tick/optim/prox/prox_slope.py

@@ -25,12 +25,9 @@ class ProxSlope(Prox):
         Range on which the prox is applied. If `None` then the prox is
         applied on the whole vector
 
-    Attributes
-    ----------
-    weights : `np.array`, shape=(n_coeffs,)
-        The weights used in the penalization. They are automatically
-        setted, depending on the ``weights_type`` and ``fdr``
-        parameters.
+    positive : `bool`, default=`False`
+        If True, apply an extra projection onto the set of vectors with
+        non-negative entries
 
     Notes
     -----
@@ -53,9 +50,6 @@ class ProxSlope(Prox):
         "positive": {
             "writable": True,
             "cpp_setter": "set_positive"
-        },
-        "weights": {
-            "writable": False,
         }
     }
 
@@ -65,7 +59,6 @@ def __init__(self, strength: float, fdr: float=0.6, range: tuple=None,
         self.strength = strength
         self.fdr = fdr
         self.positive = positive
-        self.weights = None
         if range is None:
             self._prox = _ProxSlope(self.strength, self.fdr, self.positive)
         else:
@@ -91,8 +84,3 @@ def value(self, coeffs: np.ndarray) -> float:
             Value of the penalization at ``coeffs``
         """
         return self._prox.value(coeffs)
-
-    def _as_dict(self):
-        dd = Prox._as_dict(self)
-        del dd["weights"]
-        return dd

tick/optim/prox/src/CMakeLists.txt

@@ -8,6 +8,7 @@ add_library(tick_prox EXCLUDE_FROM_ALL
         prox_zero.cpp prox_zero.h
         prox_sorted_l1.cpp prox_sorted_l1.h
         prox_slope.cpp prox_slope.h
+        prox_oscar.cpp prox_oscar.h
         prox_tv.cpp prox_tv.h
         prox_l1w.cpp prox_l1w.h
         prox_elasticnet.cpp prox_elasticnet.h

tick/optim/prox/src/prox_oscar.cpp

@@ -0,0 +1,49 @@
+// License: BSD 3 clause
+
+#include "prox_oscar.h"
+
+ProxOscar::ProxOscar(double strength,
+                     double ratio,
+                     bool positive)
+  : ProxSortedL1(strength, positive) {
+  set_ratio(ratio);
+}
+
+ProxOscar::ProxOscar(double strength,
+                     double ratio,
+                     ulong start,
+                     ulong end,
+                     bool positive)
+  : ProxSortedL1(strength, start, end, positive) {
+  set_ratio(ratio);
+}
+
+const std::string ProxOscar::get_class_name() const {
+  return "ProxOscar";
+}
+
+void ProxOscar::compute_weights(void) {
+  if (!weights_ready) {
+    ulong size = end - start;
+    weights = ArrayDouble(size);
+    for (ulong i = 0; i < size; i++) {
+      weights[i] = strength * (ratio * (size - i - 1)  + 1);
+    }
+    weights_ready = true;
+  }
+}
+
+double ProxOscar::get_ratio() const {
+  return ratio;
+}
+
+void ProxOscar::set_ratio(double ratio) {
+  if (ratio < 0) {
+    TICK_ERROR("Ratio should be non-negative");
+  } else {
+    if (ratio != this->ratio) {
+      weights_ready = false;
+      this->ratio = ratio;
+    }
+  }
+}

tick/optim/prox/src/prox_oscar.h

@@ -0,0 +1,29 @@
+#ifndef TICK_OPTIM_PROX_SRC_PROX_OSCAR_H_
+#define TICK_OPTIM_PROX_SRC_PROX_OSCAR_H_
+
+// License: BSD 3 clause
+
+#include "prox_sorted_l1.h"
+
+class ProxOscar : public ProxSortedL1 {
+ protected:
+  void compute_weights(void) override;
+  double ratio;
+
+ public:
+  ProxOscar(double strength, double ratio, bool positive);
+
+  ProxOscar(double strength,
+            double ratio,
+            ulong start,
+            ulong end,
+            bool positive);
+
+  const std::string get_class_name() const override;
+
+  virtual double get_ratio() const;
+
+  virtual void set_ratio(double ratio);
+};
+
+#endif  // TICK_OPTIM_PROX_SRC_PROX_OSCAR_H_

tick/optim/prox/src/prox_slope.cpp

@@ -5,7 +5,7 @@
 ProxSlope::ProxSlope(double strength,
                      double false_discovery_rate,
                      bool positive)
-  : ProxSortedL1(strength, WeightsType::bh, positive) {
+    : ProxSortedL1(strength, positive) {
   this->false_discovery_rate = false_discovery_rate;
 }
 
@@ -14,7 +14,7 @@ ProxSlope::ProxSlope(double strength,
                      ulong start,
                      ulong end,
                      bool positive)
-  : ProxSortedL1(strength, WeightsType::bh, start, end, positive) {
+    : ProxSortedL1(strength, start, end, positive) {
   this->false_discovery_rate = false_discovery_rate;
 }
 

tick/optim/prox/src/prox_slope.h

@@ -29,7 +29,7 @@ class ProxSlope : public ProxSortedL1 {
   inline void set_false_discovery_rate(double false_discovery_rate) {
     if (false_discovery_rate <= 0 || false_discovery_rate >= 1) {
       TICK_ERROR("False discovery rate must be in (0, 1) but received "
-                   << false_discovery_rate)
+                     << false_discovery_rate)
     }
     if (false_discovery_rate != this->false_discovery_rate) {
       weights_ready = false;

tick/optim/prox/src/prox_sorted_l1.cpp

@@ -3,20 +3,16 @@
 #include "prox_sorted_l1.h"
 
 ProxSortedL1::ProxSortedL1(double strength,
-                           WeightsType weights_type,
                            bool positive)
-  : Prox(strength, positive) {
-  this->weights_type = weights_type;
+    : Prox(strength, positive) {
   weights_ready = false;
 }
 
 ProxSortedL1::ProxSortedL1(double strength,
-                           WeightsType weights_type,
                            ulong start,
                            ulong end,
                            bool positive)
-  : Prox(strength, start, end, positive) {
-  this->weights_type = weights_type;
+    : Prox(strength, start, end, positive) {
   weights_ready = false;
 }