Merge pull request scikit-learn#2882 from larsmans/expit

Speed up RBM training with scipy.special.expit

Author: larsmans

sklearn/neural_network/rbm.py

@@ -4,6 +4,7 @@
 # Authors: Yann N. Dauphin <[email protected]>
 #          Vlad Niculae
 #          Gabriel Synnaeve
+#          Lars Buitinck
 # License: BSD 3 clause
 
 import time
@@ -19,7 +20,8 @@
 from ..utils import gen_even_slices
 from ..utils import issparse
 from ..utils.extmath import safe_sparse_dot
-from ..utils.extmath import logistic_sigmoid
+from ..utils.extmath import log_logistic
+from ..utils.fixes import expit             # logistic function
 
 
 class BernoulliRBM(BaseEstimator, TransformerMixin):
@@ -130,8 +132,9 @@ def _mean_hiddens(self, v):
         h : array-like, shape (n_samples, n_components)
             Corresponding mean field values for the hidden layer.
         """
-        return logistic_sigmoid(safe_sparse_dot(v, self.components_.T)
-                                + self.intercept_hidden_)
+        p = safe_sparse_dot(v, self.components_.T)
+        p += self.intercept_hidden_
+        return expit(p, out=p)
 
     def _sample_hiddens(self, v, rng):
         """Sample from the distribution P(h|v).
@@ -169,8 +172,9 @@ def _sample_visibles(self, h, rng):
         v : array-like, shape (n_samples, n_features)
             Values of the visible layer.
         """
-        p = logistic_sigmoid(np.dot(h, self.components_)
-                             + self.intercept_visible_)
+        p = np.dot(h, self.components_)
+        p += self.intercept_visible_
+        expit(p, out=p)
         p[rng.uniform(size=p.shape) < p] = 1.
         return np.floor(p, p)
 
@@ -274,7 +278,7 @@ def score_samples(self, X):
 
         fe = self._free_energy(v)
         fe_ = self._free_energy(v_)
-        return v.shape[1] * logistic_sigmoid(fe_ - fe, log=True)
+        return v.shape[1] * log_logistic(fe_ - fe)
 
     def fit(self, X, y=None):
         """Fit the model to the data X.

sklearn/utils/extmath.py

@@ -14,7 +14,7 @@
 from scipy.sparse import issparse
 from distutils.version import LooseVersion
 
-from . import check_random_state
+from . import check_random_state, deprecated
 from .fixes import qr_economic
 from ._logistic_sigmoid import _log_logistic_sigmoid
 from ..externals.six.moves import xrange
@@ -567,36 +567,38 @@ def svd_flip(u, v):
     return u, v
 
 
+@deprecated('to be removed in 0.17; use scipy.special.expit or log_logistic')
 def logistic_sigmoid(X, log=False, out=None):
-    """
-    Implements the logistic function, ``1 / (1 + e ** -x)`` and its log.
+    """Logistic function, ``1 / (1 + e ** (-x))``, or its log."""
+    from .fixes import expit
+    fn = log_logistic if log else expit
+    return fn(X, out)
+
 
-    This implementation is more stable by splitting on positive and negative
-    values and computing::
 
-        1 / (1 + exp(-x_i)) if x_i > 0
-        exp(x_i) / (1 + exp(x_i)) if x_i <= 0
+def log_logistic(X, out=None):
+    """Compute the log of the logistic function, ``log(1 / (1 + e ** -x))``.
 
-    The log is computed using::
+    This implementation is numerically stable because it splits positive and
+    negative values::
 
-        -log(1 + exp(-x_i)) if x_i > 0
+        -log(1 + exp(-x_i))     if x_i > 0
         x_i - log(1 + exp(x_i)) if x_i <= 0
 
+    For the ordinary logistic function, use ``sklearn.utils.fixes.expit``.
+
     Parameters
     ----------
     X: array-like, shape (M, N)
         Argument to the logistic function
 
-    log: boolean, default: False
-        Whether to compute the logarithm of the logistic function.
-
     out: array-like, shape: (M, N), optional:
         Preallocated output array.
 
     Returns
     -------
     out: array, shape (M, N)
-        Value of the logistic function evaluated at every point in x
+        Log of the logistic function evaluated at every point in x
 
     Notes
     -----
@@ -611,15 +613,7 @@ def logistic_sigmoid(X, log=False, out=None):
     if out is None:
         out = np.empty_like(X)
 
-    if log:
-        _log_logistic_sigmoid(n_samples, n_features, X, out)
-    else:
-        # logistic(x) = (1 + tanh(x / 2)) / 2
-        out[:] = X
-        out *= .5
-        np.tanh(out, out)
-        out += 1
-        out *= .5
+    _log_logistic_sigmoid(n_samples, n_features, X, out)
 
     if is_1d:
         return np.squeeze(out)

sklearn/utils/fixes.py

@@ -6,6 +6,7 @@
 # Authors: Emmanuelle Gouillart <[email protected]>
 #          Gael Varoquaux <[email protected]>
 #          Fabian Pedregosa <[email protected]>
+#          Lars Buitinck
 #
 # License: BSD 3 clause
 
@@ -109,6 +110,27 @@ def _logaddexp(x1, x2, out=None):
     logaddexp = np.logaddexp
 
 
+try:
+    from scipy.special import expit     # SciPy >= 0.10
+except ImportError:
+    def expit(x, out=None):
+        """Logistic sigmoid function, ``1 / (1 + exp(-x))``.
+
+        See sklearn.utils.extmath.log_logistic for the log of this function.
+        """
+        if out is None:
+            out = np.copy(x)
+
+        # 1 / (1 + exp(-x)) = (1 + tanh(x / 2)) / 2
+        # This way of computing the logistic is both fast and stable.
+        out *= .5
+        np.tanh(out, out)
+        out += 1
+        out *= .5
+
+        return out
+
+
 def _bincount(X, weights=None, minlength=None):
     """Replacing np.bincount in numpy < 1.6 to provide minlength."""
     result = np.bincount(X, weights)

sklearn/utils/tests/test_extmath.py

@@ -3,6 +3,8 @@
 #          Denis Engemann <[email protected]>
 #
 # License: BSD 3 clause
+import warnings
+
 import numpy as np
 from scipy import sparse
 from scipy import linalg
@@ -22,7 +24,7 @@
 from sklearn.utils.extmath import row_norms
 from sklearn.utils.extmath import weighted_mode
 from sklearn.utils.extmath import cartesian
-from sklearn.utils.extmath import logistic_sigmoid
+from sklearn.utils.extmath import log_logistic, logistic_sigmoid
 from sklearn.utils.extmath import fast_dot, _fast_dot
 from sklearn.datasets.samples_generator import make_low_rank_matrix
 
@@ -273,32 +275,16 @@ def test_cartesian():
 
 def test_logistic_sigmoid():
     """Check correctness and robustness of logistic sigmoid implementation"""
-    naive_logsig = lambda x: 1 / (1 + np.exp(-x))
-    naive_log_logsig = lambda x: np.log(naive_logsig(x))
-
-    # Simulate the previous Cython implementations of logistic_sigmoid based on
-    #http://fa.bianp.net/blog/2013/numerical-optimizers-for-logistic-regression
-    def stable_logsig(x):
-        out = np.zeros_like(x)
-        positive = x > 0
-        negative = x <= 0
-        out[positive] = 1. / (1 + np.exp(-x[positive]))
-        out[negative] = np.exp(x[negative]) / (1. + np.exp(x[negative]))
-        return out
+    naive_logistic = lambda x: 1 / (1 + np.exp(-x))
+    naive_log_logistic = lambda x: np.log(naive_logistic(x))
 
     x = np.linspace(-2, 2, 50)
-    assert_array_almost_equal(logistic_sigmoid(x), naive_logsig(x))
-    assert_array_almost_equal(logistic_sigmoid(x, log=True),
-                              naive_log_logsig(x))
-    assert_array_almost_equal(logistic_sigmoid(x), stable_logsig(x),
-                              decimal=16)
-
-    extreme_x = np.array([-100, 100], dtype=np.float)
-    assert_array_almost_equal(logistic_sigmoid(extreme_x), [0, 1])
-    assert_array_almost_equal(logistic_sigmoid(extreme_x, log=True), [-100, 0])
-    assert_array_almost_equal(logistic_sigmoid(extreme_x),
-                              stable_logsig(extreme_x),
-                              decimal=16)
+    with warnings.catch_warnings(record=True):
+        assert_array_almost_equal(logistic_sigmoid(x), naive_logistic(x))
+    assert_array_almost_equal(log_logistic(x), naive_log_logistic(x))
+
+    extreme_x = np.array([-100., 100.])
+    assert_array_almost_equal(log_logistic(extreme_x), [-100, 0])
 
 
 def test_fast_dot():

sklearn/utils/tests/test_fixes.py

@@ -5,9 +5,9 @@
 import numpy as np
 
 from nose.tools import assert_equal
-from numpy.testing import assert_array_equal
+from numpy.testing import assert_almost_equal, assert_array_equal
 
-from ..fixes import _in1d, _copysign, divide
+from ..fixes import _in1d, _copysign, divide, expit
 
 
 def test_in1d():
@@ -16,6 +16,16 @@ def test_in1d():
     assert_equal(_in1d(a, b).sum(), 5)
 
 
+def test_expit():
+    """Check numerical stability of expit (logistic function)."""
+
+    # Simulate our previous Cython implementation, based on
+    #http://fa.bianp.net/blog/2013/numerical-optimizers-for-logistic-regression
+    assert_almost_equal(expit(100.), 1. / (1. + np.exp(-100.)), decimal=16)
+    assert_almost_equal(expit(-100.), np.exp(-100.) / (1. + np.exp(-100.)),
+                        decimal=16)
+
+
 def test_divide():
     assert_equal(divide(.6, 1), .600000000000)