FIX numerical stability issue in BernoulliRBM

Fixes scikit-learn#2785 by replacing log1p(exp(x)) with logaddexp(0, x).

Author: larsmans

sklearn/neural_network/rbm.py

@@ -1,10 +1,10 @@
 """Restricted Boltzmann Machine
 """
 
-# Main author: Yann N. Dauphin <[email protected]>
-# Author: Vlad Niculae
-# Author: Gabriel Synnaeve
-# License: BSD Style.
+# Authors: Yann N. Dauphin <[email protected]>
+#          Vlad Niculae
+#          Gabriel Synnaeve
+# License: BSD 3 clause
 
 import time
 
@@ -188,8 +188,8 @@ def _free_energy(self, v):
             The value of the free energy.
         """
         return (- safe_sparse_dot(v, self.intercept_visible_)
-                - np.log1p(np.exp(safe_sparse_dot(v, self.components_.T)
-                                  + self.intercept_hidden_)).sum(axis=1))
+                - np.logaddexp(0, safe_sparse_dot(v, self.components_.T)
+                                  + self.intercept_hidden_).sum(axis=1))
 
     def gibbs(self, v):
         """Perform one Gibbs sampling step.

sklearn/neural_network/tests/test_rbm.py

@@ -127,6 +127,11 @@ def test_score_samples():
     s_score = rbm1.score_samples(lil_matrix(X))
     assert_almost_equal(d_score, s_score)
 
+    # Test numerical stability (#2785): would previously generate infinities
+    # and crash with an exception.
+    with np.errstate(under='ignore'):
+        rbm1.score_samples(np.arange(1000) * 100)
+
 
 def test_rbm_verbose():
     rbm = BernoulliRBM(n_iter=2, verbose=10)