MNT replace log_logistic by logaddexp (scikit-learn#27544)

Co-authored-by: Guillaume Lemaitre <g.lemaitre58@gmail.com>
glemaitre · Oct 31, 2023 · 66872ff · 66872ff
1 parent 19cbd13
commit 66872ff
Show file tree

Hide file tree

Showing 6 changed files with 23 additions and 48 deletions.
diff --git a/doc/whats_new/v1.4.rst b/doc/whats_new/v1.4.rst
@@ -426,6 +426,10 @@ Changelog
   array and does not return the correct result for SciPy sparse matrix.
   :pr:`27336` by :user:`Guillaume Lemaitre <glemaitre>`.
 
+- |API| :func:`sklearn.extmath.log_logistic` is deprecated and will be removed in 1.6.
+  Use `-np.logaddexp(0, -x)` instead.
+  :pr:`27544` by :user:`Christian Lorentzen <lorentzenchr>`.
+
 Code and Documentation Contributors
 -----------------------------------
 

diff --git a/setup.py b/setup.py
@@ -405,7 +405,6 @@ def check_package_status(package, min_version):
             "include_np": True,
         },
         {"sources": ["_random.pyx"], "include_np": True},
-        {"sources": ["_logistic_sigmoid.pyx"], "include_np": True},
         {"sources": ["_typedefs.pyx"]},
         {"sources": ["_heap.pyx"]},
         {"sources": ["_sorting.pyx"]},

diff --git a/sklearn/neural_network/_rbm.py b/sklearn/neural_network/_rbm.py
@@ -22,7 +22,7 @@
 )
 from ..utils import check_random_state, gen_even_slices
 from ..utils._param_validation import Interval
-from ..utils.extmath import log_logistic, safe_sparse_dot
+from ..utils.extmath import safe_sparse_dot
 from ..utils.validation import check_is_fitted
 
 
@@ -380,7 +380,8 @@ def score_samples(self, X):
 
         fe = self._free_energy(v)
         fe_ = self._free_energy(v_)
-        return v.shape[1] * log_logistic(fe_ - fe)
+        # log(expit(x)) = log(1 / (1 + exp(-x)) = -np.logaddexp(0, -x)
+        return -v.shape[1] * np.logaddexp(0, -(fe_ - fe))
 
     @_fit_context(prefer_skip_nested_validation=True)
     def fit(self, X, y=None):

diff --git a/sklearn/utils/_logistic_sigmoid.pyx b/sklearn/utils/_logistic_sigmoid.pyx
diff --git a/sklearn/utils/extmath.py b/sklearn/utils/extmath.py
@@ -19,10 +19,10 @@
 import numpy as np
 from scipy import linalg, sparse
 
+from ..utils import deprecated
 from ..utils._param_validation import Interval, StrOptions, validate_params
 from . import check_random_state
 from ._array_api import _is_numpy_namespace, device, get_namespace
-from ._logistic_sigmoid import _log_logistic_sigmoid
 from .sparsefuncs_fast import csr_row_norms
 from .validation import check_array
 
@@ -885,14 +885,15 @@ def svd_flip(u, v, u_based_decision=True):
     return u, v
 
 
+# TODO(1.6): remove
+@deprecated(  # type: ignore
+    "The function `log_logistic` is deprecated and will be removed in 1.6. "
+    "Use `-np.logaddexp(0, -x)` instead."
+)
 def log_logistic(X, out=None):
     """Compute the log of the logistic function, ``log(1 / (1 + e ** -x))``.
 
-    This implementation is numerically stable because it splits positive and
-    negative values::
-
-        -log(1 + exp(-x_i))     if x_i > 0
-        x_i - log(1 + exp(x_i)) if x_i <= 0
+    This implementation is numerically stable and uses `-np.logaddexp(0, -x)`.
 
     For the ordinary logistic function, use ``scipy.special.expit``.
 
@@ -914,19 +915,13 @@ def log_logistic(X, out=None):
     See the blog post describing this implementation:
     http://fa.bianp.net/blog/2013/numerical-optimizers-for-logistic-regression/
     """
-    is_1d = X.ndim == 1
-    X = np.atleast_2d(X)
-    X = check_array(X, dtype=np.float64)
-
-    n_samples, n_features = X.shape
+    X = check_array(X, dtype=np.float64, ensure_2d=False)
 
     if out is None:
         out = np.empty_like(X)
 
-    _log_logistic_sigmoid(n_samples, n_features, X, out)
-
-    if is_1d:
-        return np.squeeze(out)
+    np.logaddexp(0, -X, out=out)
+    out *= -1
     return out
 
 

diff --git a/sklearn/utils/tests/test_extmath.py b/sklearn/utils/tests/test_extmath.py
@@ -687,16 +687,20 @@ def test_cartesian_mix_types(arrays, output_dtype):
     assert output.dtype == output_dtype
 
 
+# TODO(1.6): remove this test
 def test_logistic_sigmoid():
     # Check correctness and robustness of logistic sigmoid implementation
     def naive_log_logistic(x):
         return np.log(expit(x))
 
     x = np.linspace(-2, 2, 50)
-    assert_array_almost_equal(log_logistic(x), naive_log_logistic(x))
+    warn_msg = "`log_logistic` is deprecated and will be removed"
+    with pytest.warns(FutureWarning, match=warn_msg):
+        assert_array_almost_equal(log_logistic(x), naive_log_logistic(x))
 
     extreme_x = np.array([-100.0, 100.0])
-    assert_array_almost_equal(log_logistic(extreme_x), [-100, 0])
+    with pytest.warns(FutureWarning, match=warn_msg):
+        assert_array_almost_equal(log_logistic(extreme_x), [-100, 0])
 
 
 @pytest.fixture()