From 0c63de2bda96a73bec565d084edff81ae0306d11 Mon Sep 17 00:00:00 2001
From: Guilherme Henriques <guilherme.henriques@tecnico.ulisboa.pt>
Date: Fri, 28 Mar 2025 10:46:14 +0000
Subject: [PATCH] Fix #23179: added pseudo-likelihood normalization option in
 RBM

---
 sklearn/neural_network/_rbm.py           | 22 ++++++++++++++++---
 sklearn/neural_network/tests/test_rbm.py | 27 ++++++++++++++++++++++++
 2 files changed, 46 insertions(+), 3 deletions(-)

diff --git a/sklearn/neural_network/_rbm.py b/sklearn/neural_network/_rbm.py
index 1e1d3c2e11b7c..a92c9f28ef526 100644
--- a/sklearn/neural_network/_rbm.py
+++ b/sklearn/neural_network/_rbm.py
@@ -68,6 +68,11 @@ class BernoulliRBM(ClassNamePrefixFeaturesOutMixin, TransformerMixin, BaseEstima
         Pass an int for reproducible results across multiple function calls.
         See :term:`Glossary <random_state>`.
 
+    normalize : bool, default=False
+        If True, returns the pseudo-likelihood normalized by the number
+        of visible units (features). This is useful for comparing models
+        trained on inputs of different dimensionalities.
+
     Attributes
     ----------
     intercept_hidden_ : array-like of shape (n_components,)
@@ -135,6 +140,7 @@ class BernoulliRBM(ClassNamePrefixFeaturesOutMixin, TransformerMixin, BaseEstima
         "n_iter": [Interval(Integral, 0, None, closed="left")],
         "verbose": ["verbose"],
         "random_state": ["random_state"],
+        "normalize": ["boolean"],
     }
 
     def __init__(
@@ -146,6 +152,7 @@ def __init__(
         n_iter=10,
         verbose=0,
         random_state=None,
+        normalize=False,
     ):
         self.n_components = n_components
         self.learning_rate = learning_rate
@@ -153,6 +160,7 @@ def __init__(
         self.n_iter = n_iter
         self.verbose = verbose
         self.random_state = random_state
+        self.normalize = normalize
 
     def transform(self, X):
         """Compute the hidden layer activation probabilities, P(h=1|v=X).
@@ -341,7 +349,7 @@ def _fit(self, v_pos, rng):
         h_neg[rng.uniform(size=h_neg.shape) < h_neg] = 1.0  # sample binomial
         self.h_samples_ = np.floor(h_neg, h_neg)
 
-    def score_samples(self, X):
+    def score_samples(self, X, normalize=False):
         """Compute the pseudo-likelihood of X.
 
         Parameters
@@ -349,6 +357,11 @@ def score_samples(self, X):
         X : {array-like, sparse matrix} of shape (n_samples, n_features)
             Values of the visible layer. Must be all-boolean (not checked).
 
+        normalize : bool, default=False
+            If True, returns the pseudo-likelihood normalized by the number
+            of visible units (features). This is useful for comparing models
+            trained on inputs of different dimensionalities.
+
         Returns
         -------
         pseudo_likelihood : ndarray of shape (n_samples,)
@@ -380,7 +393,10 @@ def score_samples(self, X):
         fe = self._free_energy(v)
         fe_ = self._free_energy(v_)
         # log(expit(x)) = log(1 / (1 + exp(-x)) = -np.logaddexp(0, -x)
-        return -v.shape[1] * np.logaddexp(0, -(fe_ - fe))
+        if normalize:
+            return -np.logaddexp(0, -(fe_ - fe))
+        else:
+            return -v.shape[1] * np.logaddexp(0, -(fe_ - fe))
 
     @_fit_context(prefer_skip_nested_validation=True)
     def fit(self, X, y=None):
@@ -430,7 +446,7 @@ def fit(self, X, y=None):
                     % (
                         type(self).__name__,
                         iteration,
-                        self.score_samples(X).mean(),
+                        self.score_samples(X, normalize=self.normalize).mean(),
                         end - begin,
                     )
                 )
diff --git a/sklearn/neural_network/tests/test_rbm.py b/sklearn/neural_network/tests/test_rbm.py
index 8211c9735923d..00e8ef27bed2d 100644
--- a/sklearn/neural_network/tests/test_rbm.py
+++ b/sklearn/neural_network/tests/test_rbm.py
@@ -177,6 +177,33 @@ def test_rbm_verbose():
         sys.stdout = old_stdout
 
 
+def test_rbm_pseudolikelihood():
+    X1 = np.ones([11, 100])
+    rbm1 = BernoulliRBM(n_iter=5, random_state=42, normalize=False)
+    rbm1.fit(X1)
+    score1 = rbm1.score_samples(X1, normalize=False).mean()
+
+    X2 = np.ones([11, 1000])
+    rbm2 = BernoulliRBM(n_iter=5, random_state=42, normalize=False)
+    rbm2.fit(X2)
+    score2 = rbm2.score_samples(X2, normalize=False).mean()
+
+    X3 = np.ones([11, 100])
+    rbm3 = BernoulliRBM(n_iter=5, random_state=42, normalize=True)
+    rbm3.fit(X3)
+    score3 = rbm3.score_samples(X3, normalize=True).mean()
+
+    X4 = np.ones([11, 1000])
+    rbm4 = BernoulliRBM(n_iter=5, random_state=42, normalize=True)
+    rbm4.fit(X4)
+    score4 = rbm4.score_samples(X4, normalize=True).mean()
+
+    assert score1 < score3
+    assert score2 < score4
+    assert abs(score1) > abs(score3) * 10
+    assert abs(score2) > abs(score4) * 10
+
+
 @pytest.mark.parametrize("csc_container", CSC_CONTAINERS)
 def test_sparse_and_verbose(csc_container):
     # Make sure RBM works with sparse input when verbose=True