diff --git a/doc/modules/api.rst b/doc/modules/api.rst index ac72edb05..8a11dcc57 100644 --- a/doc/modules/api.rst +++ b/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 diff --git a/setup.py b/setup.py index fdcf83104..d496c3531 100644 --- a/setup.py +++ b/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/", diff --git a/tick/optim/prox/__init__.py b/tick/optim/prox/__init__.py index 765a52256..097628326 100644 --- a/tick/optim/prox/__init__.py +++ b/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"] \ No newline at end of file + "ProxBinarsity"] diff --git a/tick/optim/prox/prox_oscar.py b/tick/optim/prox/prox_oscar.py index 0ccb3ffea..4e71f2a0a 100644 --- a/tick/optim/prox/prox_oscar.py +++ b/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 diff --git a/tick/optim/prox/prox_slope.py b/tick/optim/prox/prox_slope.py index 8d88e9c2e..f5111f5f5 100644 --- a/tick/optim/prox/prox_slope.py +++ b/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 diff --git a/tick/optim/prox/src/CMakeLists.txt b/tick/optim/prox/src/CMakeLists.txt index 26a5f219a..3ae769d75 100644 --- a/tick/optim/prox/src/CMakeLists.txt +++ b/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 diff --git a/tick/optim/prox/src/prox_oscar.cpp b/tick/optim/prox/src/prox_oscar.cpp new file mode 100644 index 000000000..5d8b76f3d --- /dev/null +++ b/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; + } + } +} diff --git a/tick/optim/prox/src/prox_oscar.h b/tick/optim/prox/src/prox_oscar.h new file mode 100644 index 000000000..cd11d0152 --- /dev/null +++ b/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_ diff --git a/tick/optim/prox/src/prox_slope.cpp b/tick/optim/prox/src/prox_slope.cpp index 257945ae5..7d146250e 100644 --- a/tick/optim/prox/src/prox_slope.cpp +++ b/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; } diff --git a/tick/optim/prox/src/prox_slope.h b/tick/optim/prox/src/prox_slope.h index f71ec8087..7297df328 100644 --- a/tick/optim/prox/src/prox_slope.h +++ b/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; diff --git a/tick/optim/prox/src/prox_sorted_l1.cpp b/tick/optim/prox/src/prox_sorted_l1.cpp index 4bcc12749..dcc457832 100644 --- a/tick/optim/prox/src/prox_sorted_l1.cpp +++ b/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; } diff --git a/tick/optim/prox/src/prox_sorted_l1.h b/tick/optim/prox/src/prox_sorted_l1.h index e4f970f44..26e6c0f69 100644 --- a/tick/optim/prox/src/prox_sorted_l1.h +++ b/tick/optim/prox/src/prox_sorted_l1.h @@ -5,14 +5,8 @@ #include "prox.h" -enum class WeightsType { - bh = 0, - oscar -}; - class ProxSortedL1 : public Prox { protected: - WeightsType weights_type; ArrayDouble weights; bool weights_ready; @@ -22,11 +16,9 @@ class ProxSortedL1 : public Prox { ArrayDouble &x) const; public: - ProxSortedL1(double strength, WeightsType weights_type, - bool positive); + ProxSortedL1(double strength, bool positive); - ProxSortedL1(double strength, WeightsType weights_type, ulong start, - ulong end, bool positive); + ProxSortedL1(double strength, ulong start, ulong end, bool positive); const std::string get_class_name() const override; @@ -35,16 +27,7 @@ class ProxSortedL1 : public Prox { void call(const ArrayDouble &coeffs, double t, ArrayDouble &out, ulong start, ulong end) override; - inline WeightsType get_weights_type() const { - return weights_type; - } - - inline void set_weights_type(WeightsType weights_type) { - this->weights_type = weights_type; - weights_ready = false; - } - - inline double get_weight_i(ulong i) { + inline double get_weight_i(ulong i) const { return weights[i]; } diff --git a/tick/optim/prox/swig/prox_module.i b/tick/optim/prox/swig/prox_module.i index e925307d4..7466bd43a 100644 --- a/tick/optim/prox/swig/prox_module.i +++ b/tick/optim/prox/swig/prox_module.i @@ -15,6 +15,7 @@ %shared_ptr(ProxTV); %shared_ptr(ProxElasticNet); %shared_ptr(ProxSlope); +%shared_ptr(ProxOscar); %shared_ptr(ProxMulti); %shared_ptr(ProxEquality); %shared_ptr(ProxBinarsity); @@ -45,8 +46,10 @@ %include prox_slope.i +%include prox_oscar.i + %include prox_multi.i %include prox_equality.i -%include prox_binarsity.i \ No newline at end of file +%include prox_binarsity.i diff --git a/tick/optim/prox/swig/prox_oscar.i b/tick/optim/prox/swig/prox_oscar.i new file mode 100644 index 000000000..3c0310583 --- /dev/null +++ b/tick/optim/prox/swig/prox_oscar.i @@ -0,0 +1,18 @@ +// License: BSD 3 clause + +%{ +#include "prox_oscar.h" +%} + +class ProxOscar : public Prox { + public: + ProxOscar(double strength, double ratio, bool positive); + + ProxOscar(double strength, double ratio, unsigned long start, unsigned long end, bool positive); + + virtual double get_ratio() const; + + virtual void set_ratio(double ratio); + + inline double get_weight_i(unsigned long i) const; +}; diff --git a/tick/optim/prox/swig/prox_slope.i b/tick/optim/prox/swig/prox_slope.i index 2e1caff3a..9eefd0a2e 100644 --- a/tick/optim/prox/swig/prox_slope.i +++ b/tick/optim/prox/swig/prox_slope.i @@ -14,5 +14,5 @@ class ProxSlope : public Prox { inline void set_false_discovery_rate(double fdr); - inline double get_weight_i(unsigned long i); + inline double get_weight_i(unsigned long i) const; }; diff --git a/tick/optim/prox/tests/prox_oscar_test.py b/tick/optim/prox/tests/prox_oscar_test.py new file mode 100644 index 000000000..972f1e12a --- /dev/null +++ b/tick/optim/prox/tests/prox_oscar_test.py @@ -0,0 +1,168 @@ +# License: BSD 3 clause + +import unittest + +import numpy as np +from scipy.stats import norm as normal + +from tick.optim.prox import ProxSlope +from tick.optim.prox.tests.prox import TestProx + + +class Test(TestProx): + + def get_coeffs_weak(self, n_coeffs): + return (-1) ** np.arange(n_coeffs) * \ + np.sqrt(2 * np.log(np.linspace(1, 10, n_coeffs) * n_coeffs)) + + def get_coeffs_strong(self, n_coeffs): + return 5 * self.get_coeffs_weak(n_coeffs) + + def get_weights_bh(self, strength, fdr, size): + tmp = fdr / (2 * size) + return strength * normal.ppf(1 - tmp * np.arange(1, size + 1)) + + def test_ProxSlope_call(self): + """...Test_of ProxSlope.call + """ + n_coeffs = 30 + np.random.seed(seed=123) + coeffs = np.zeros(n_coeffs) + coeffs[:10] = self.get_coeffs_strong(10) + y = coeffs + 3 * np.random.randn(n_coeffs) + fdr = 0.6 + strength = 4. + step = 1. + prox = ProxSlope(strength=strength, fdr=fdr) + + x_sl1_truth \ + = np.array([2.13923654, -3.53844034, 7.31022147, -9.87610726, + 6.03085298, -3.53844034, 2.13923654, -9.08606547, + 9.87610726, -9.87610726, -0., -0., 0.08661054, -0., + -0., -0., 1.78604443, 1.78604443, 0., 0., 0., + 0.08661054, -0., 0., -0., -0., 0., -0.08661054, -0., + -0.]) + np.testing.assert_almost_equal(prox.call(y, step=step), x_sl1_truth) + + strength = 2. + step = 2. + prox.strength = strength + np.testing.assert_almost_equal(prox.call(y, step=step), x_sl1_truth) + + prox.range = (10, 20) + strength = 4. + step = 1. + prox.strength = strength + + x_sl1_truth \ + = np.array([7.47293832, -9.24669781, 13.88963598, -18.0998993, + 12.24994736, -9.35363322, 7.29476114, -16.08880976, + 18.79749156, -17.7744925, -0., -0., 0., -0., -0., -0., + 0., 0., 0., 0., 2.21210573, 4.47219608, -2.80750161, + 3.52748713, -3.761642, -1.91325451, 2.72131559, + -4.2860421, -0.42020616, -2.58526469]) + np.testing.assert_almost_equal(prox.call(y, step=step), x_sl1_truth) + + strength = 2. + step = 2. + prox.strength = strength + np.testing.assert_almost_equal(prox.call(y, step=step), x_sl1_truth) + + def test_ProxSlope_weights_bh(self): + """...Test of ProxSlope weights computation + """ + n_samples = 5000 + n_outliers = int(0.1 * n_samples) + interc0 = np.zeros(n_samples) + interc0[:n_outliers] = self.get_coeffs_strong(n_outliers) + y = interc0 + 3 * np.random.randn(interc0.shape[0]) + + strength = 2.5 + prox = ProxSlope(strength=strength) + prox.value(y) + prox.strength = 20 + size = y.shape[0] + tmp = prox.fdr / (2 * size) + weights = strength * normal.ppf(1 - tmp * np.arange(1, size + 1)) + + prox_weights = np.array([prox._prox.get_weight_i(i) + for i in range(size)]) + np.testing.assert_almost_equal(prox_weights, weights, decimal=7) + + strength = 2.5 + prox = ProxSlope(strength=strength, range=(300, 3000)) + prox.value(y) + prox.strength = 20 + size = prox.range[1] - prox.range[0] + tmp = prox.fdr / (2 * size) + weights = strength * normal.ppf(1 - tmp * np.arange(1, size + 1)) + + prox_weights = np.array([prox._prox.get_weight_i(i) + for i in range(size)]) + np.testing.assert_almost_equal(prox_weights, weights, decimal=7) + + def test_ProxSlope_value(self): + """...Test of ProxSlope.value + """ + n_samples = 5000 + n_outliers = int(0.1 * n_samples) + interc0 = np.zeros(n_samples) + interc0[:n_outliers] = self.get_coeffs_strong(n_outliers) + y = interc0 + 3 * np.random.randn(interc0.shape[0]) + + fdr = 0.6 + strength = 2.5 + prox = ProxSlope(strength=strength, fdr=fdr) + prox_value = prox.value(y) + + weights = self.get_weights_bh(strength, fdr, size=y.shape[0]) + y_abs = np.abs(y) + sub_y = y_abs + value = sub_y[np.argsort(-sub_y)].dot(weights) + places = 4 + self.assertAlmostEqual(prox_value, value, places=places) + + strength = 7.4 + fdr = 0.2 + prox.strength = strength + prox.fdr = fdr + prox_value = prox.value(y) + + weights = self.get_weights_bh(strength, fdr, size=y.shape[0]) + y_abs = np.abs(y) + sub_y = y_abs + value = sub_y[np.argsort(-sub_y)].dot(weights) + self.assertAlmostEqual(prox_value, value, places=places) + + prox.range = (300, 3000) + fdr = 0.6 + strength = 2.5 + prox.strength = strength + prox.fdr = fdr + prox_value = prox.value(y) + + a, b = prox.range + size = b - a + weights = self.get_weights_bh(strength, fdr, size=size) + y_abs = np.abs(y[a:b]) + sub_y = y_abs + value = sub_y[np.argsort(-sub_y)].dot(weights) + self.assertAlmostEqual(prox_value, value, places=places) + + prox.range = (300, 3000) + strength = 7.4 + fdr = 0.2 + prox.strength = strength + prox.fdr = fdr + prox_value = prox.value(y) + + a, b = prox.range + size = b - a + weights = self.get_weights_bh(strength, fdr, size=size) + y_abs = np.abs(y[a:b]) + sub_y = y_abs + value = sub_y[np.argsort(-sub_y)].dot(weights) + self.assertAlmostEqual(prox_value, value, places=places) + +if __name__ == '__main__': + unittest.main()