From b3bc69babd52e7bf62bec46bb0af2ecbda327ea9 Mon Sep 17 00:00:00 2001 From: Antoine Baker Date: Thu, 27 Mar 2025 15:23:25 +0100 Subject: [PATCH 1/9] use full SVD --- sklearn/linear_model/_bayes.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/sklearn/linear_model/_bayes.py b/sklearn/linear_model/_bayes.py index 27ce01d0e75d5..18d7e3215f317 100644 --- a/sklearn/linear_model/_bayes.py +++ b/sklearn/linear_model/_bayes.py @@ -293,8 +293,17 @@ def fit(self, X, y, sample_weight=None): coef_old_ = None XT_y = np.dot(X.T, y) - U, S, Vh = linalg.svd(X, full_matrices=False) + # The full SVD is needed for the covariance matrix + # U_full: (M, M), S: K, Vh_full: (N, N) + U_full, S, Vh_full = linalg.svd(X, full_matrices=True) + K = len(S) eigen_vals_ = S**2 + eigen_vals_full = np.zeros(n_features, dtype=dtype) + eigen_vals_full[0:K] = eigen_vals_ + # The reduced SVD is enough for the updates + # U: (M, K), S: K, Vh: (K, N) + U = U_full[:, 0:K] + Vh = Vh_full[0:K, :] # Convergence loop of the bayesian ridge regression for iter_ in range(self.max_iter): @@ -353,11 +362,10 @@ def fit(self, X, y, sample_weight=None): self.scores_.append(s) self.scores_ = np.array(self.scores_) - # posterior covariance is given by 1/alpha_ * scaled_sigma_ - scaled_sigma_ = np.dot( - Vh.T, Vh / (eigen_vals_ + lambda_ / alpha_)[:, np.newaxis] + # posterior covariance + self.sigma_ = np.dot( + Vh_full.T, Vh_full / (alpha_ * eigen_vals_full + lambda_)[:, np.newaxis] ) - self.sigma_ = (1.0 / alpha_) * scaled_sigma_ self._set_intercept(X_offset_, y_offset_, X_scale_) From 730f0a316679d8d2ceddaf41265619ddffdd9a18 Mon Sep 17 00:00:00 2001 From: Antoine Baker Date: Thu, 27 Mar 2025 15:51:56 +0100 Subject: [PATCH 2/9] add test --- sklearn/linear_model/tests/test_bayes.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/sklearn/linear_model/tests/test_bayes.py b/sklearn/linear_model/tests/test_bayes.py index 6fae1536582c8..af0db731ff72f 100644 --- a/sklearn/linear_model/tests/test_bayes.py +++ b/sklearn/linear_model/tests/test_bayes.py @@ -94,6 +94,17 @@ def test_bayesian_ridge_parameter(): assert_almost_equal(rr_model.intercept_, br_model.intercept_) +def test_bayesian_covariance_matrix(): + # Check the posterior covariance matrix sigma_ + X, y = diabetes.data, diabetes.target + reg = BayesianRidge().fit(X, y) + n_samples, n_features = X.shape + covariance_matrix = np.linalg.inv( + reg.lambda_ * np.identity(n_features) + reg.alpha_ * np.dot(X.T, X) + ) + assert_array_almost_equal(reg.sigma_, covariance_matrix) + + def test_bayesian_sample_weights(): # Test correctness of the sample_weights method X = np.array([[1, 1], [3, 4], [5, 7], [4, 1], [2, 6], [3, 10], [3, 2]]) From 308f10ba31f1da5900568e3c68b98ff6a550da1d Mon Sep 17 00:00:00 2001 From: Antoine Baker Date: Thu, 27 Mar 2025 16:58:49 +0100 Subject: [PATCH 3/9] update test --- sklearn/linear_model/tests/test_bayes.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sklearn/linear_model/tests/test_bayes.py b/sklearn/linear_model/tests/test_bayes.py index af0db731ff72f..5131c7e89c28d 100644 --- a/sklearn/linear_model/tests/test_bayes.py +++ b/sklearn/linear_model/tests/test_bayes.py @@ -96,9 +96,9 @@ def test_bayesian_ridge_parameter(): def test_bayesian_covariance_matrix(): # Check the posterior covariance matrix sigma_ - X, y = diabetes.data, diabetes.target - reg = BayesianRidge().fit(X, y) - n_samples, n_features = X.shape + X, y = datasets.make_regression(n_samples=10, n_features=20) + n_features = X.shape[1] + reg = BayesianRidge(fit_intercept=False).fit(X, y) covariance_matrix = np.linalg.inv( reg.lambda_ * np.identity(n_features) + reg.alpha_ * np.dot(X.T, X) ) From 718751153f0edb9ac73551198dc9ce6aee352510 Mon Sep 17 00:00:00 2001 From: Antoine Baker Date: Thu, 27 Mar 2025 17:24:40 +0100 Subject: [PATCH 4/9] changelog --- .../upcoming_changes/sklearn.linear_model/31094.fix.rst | 4 ++++ 1 file changed, 4 insertions(+) create mode 100644 doc/whats_new/upcoming_changes/sklearn.linear_model/31094.fix.rst diff --git a/doc/whats_new/upcoming_changes/sklearn.linear_model/31094.fix.rst b/doc/whats_new/upcoming_changes/sklearn.linear_model/31094.fix.rst new file mode 100644 index 0000000000000..ea20a96f68599 --- /dev/null +++ b/doc/whats_new/upcoming_changes/sklearn.linear_model/31094.fix.rst @@ -0,0 +1,4 @@ +- :class:`linear_model.BayesianRidge` now uses the full SVD to + correctly estimate the posterior covariance matrix `sigma_` + when `n_samples < n_features`. + By :user:`Antoine Baker ` From e073189f318bb1854489e716270021028e1f689a Mon Sep 17 00:00:00 2001 From: Antoine Baker Date: Thu, 27 Mar 2025 18:17:15 +0100 Subject: [PATCH 5/9] full SVD only for wide dataset --- sklearn/linear_model/_bayes.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/sklearn/linear_model/_bayes.py b/sklearn/linear_model/_bayes.py index 18d7e3215f317..adf515d44d1d9 100644 --- a/sklearn/linear_model/_bayes.py +++ b/sklearn/linear_model/_bayes.py @@ -293,16 +293,18 @@ def fit(self, X, y, sample_weight=None): coef_old_ = None XT_y = np.dot(X.T, y) - # The full SVD is needed for the covariance matrix - # U_full: (M, M), S: K, Vh_full: (N, N) - U_full, S, Vh_full = linalg.svd(X, full_matrices=True) + # Let M, N = n_samples, n_features and K = min(M, N). + # The posterior covariance matrix needs Vh_full: (N, N). + # The full SVD is only required when n_samples < n_features. + # When n_samples < n_features, K=M and full_matrices=True + # U: (M, M), S: M, Vh_full: (N, N), Vh: (M, N) + # When n_samples > n_features, K=N and full_matrices=False + # U: (M, N), S: N, Vh_full: (N, N), Vh: (N, N) + U, S, Vh_full = linalg.svd(X, full_matrices=(n_samples < n_features)) K = len(S) eigen_vals_ = S**2 eigen_vals_full = np.zeros(n_features, dtype=dtype) eigen_vals_full[0:K] = eigen_vals_ - # The reduced SVD is enough for the updates - # U: (M, K), S: K, Vh: (K, N) - U = U_full[:, 0:K] Vh = Vh_full[0:K, :] # Convergence loop of the bayesian ridge regression From 9c4acbe8a3bdb056c8a85ae6ab0b34bb8debe7a8 Mon Sep 17 00:00:00 2001 From: antoinebaker Date: Fri, 28 Mar 2025 15:33:19 +0100 Subject: [PATCH 6/9] Update doc/whats_new/upcoming_changes/sklearn.linear_model/31094.fix.rst Co-authored-by: Olivier Grisel --- .../upcoming_changes/sklearn.linear_model/31094.fix.rst | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/doc/whats_new/upcoming_changes/sklearn.linear_model/31094.fix.rst b/doc/whats_new/upcoming_changes/sklearn.linear_model/31094.fix.rst index ea20a96f68599..b65d96bccd7d2 100644 --- a/doc/whats_new/upcoming_changes/sklearn.linear_model/31094.fix.rst +++ b/doc/whats_new/upcoming_changes/sklearn.linear_model/31094.fix.rst @@ -1,4 +1,3 @@ -- :class:`linear_model.BayesianRidge` now uses the full SVD to - correctly estimate the posterior covariance matrix `sigma_` - when `n_samples < n_features`. +- :class:`linear_model.BayesianRidge` now uses the full SVD to correctly estimate + the posterior covariance matrix `sigma_` when `n_samples < n_features`. By :user:`Antoine Baker ` From afe9c1148120b64a14dad8e5d076ee77b5d9aa73 Mon Sep 17 00:00:00 2001 From: Antoine Baker Date: Fri, 28 Mar 2025 15:34:29 +0100 Subject: [PATCH 7/9] test wide and long data --- sklearn/linear_model/tests/test_bayes.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sklearn/linear_model/tests/test_bayes.py b/sklearn/linear_model/tests/test_bayes.py index 5131c7e89c28d..4b5083c0fcf09 100644 --- a/sklearn/linear_model/tests/test_bayes.py +++ b/sklearn/linear_model/tests/test_bayes.py @@ -94,10 +94,10 @@ def test_bayesian_ridge_parameter(): assert_almost_equal(rr_model.intercept_, br_model.intercept_) -def test_bayesian_covariance_matrix(): +@pytest.mark.parametrize("n_samples, n_features", [(10, 20), (20, 10)]) +def test_bayesian_covariance_matrix(n_samples, n_features): # Check the posterior covariance matrix sigma_ - X, y = datasets.make_regression(n_samples=10, n_features=20) - n_features = X.shape[1] + X, y = datasets.make_regression(n_samples, n_features) reg = BayesianRidge(fit_intercept=False).fit(X, y) covariance_matrix = np.linalg.inv( reg.lambda_ * np.identity(n_features) + reg.alpha_ * np.dot(X.T, X) From 7f64378e997e97d33e4b73e717d774626bd60f4d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9mie=20du=20Boisberranger?= Date: Sun, 13 Apr 2025 15:48:53 +0200 Subject: [PATCH 8/9] use assert_close + ref issue --- sklearn/linear_model/tests/test_bayes.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/sklearn/linear_model/tests/test_bayes.py b/sklearn/linear_model/tests/test_bayes.py index 4b5083c0fcf09..ef32dbec60a6a 100644 --- a/sklearn/linear_model/tests/test_bayes.py +++ b/sklearn/linear_model/tests/test_bayes.py @@ -11,6 +11,7 @@ from sklearn.utils import check_random_state from sklearn.utils._testing import ( _convert_container, + assert_allclose, assert_almost_equal, assert_array_almost_equal, assert_array_less, @@ -96,13 +97,16 @@ def test_bayesian_ridge_parameter(): @pytest.mark.parametrize("n_samples, n_features", [(10, 20), (20, 10)]) def test_bayesian_covariance_matrix(n_samples, n_features): - # Check the posterior covariance matrix sigma_ + """Check the posterior covariance matrix sigma_ + + Non-regression test for https://github.com/scikit-learn/scikit-learn/issues/31093 + """ X, y = datasets.make_regression(n_samples, n_features) reg = BayesianRidge(fit_intercept=False).fit(X, y) covariance_matrix = np.linalg.inv( reg.lambda_ * np.identity(n_features) + reg.alpha_ * np.dot(X.T, X) ) - assert_array_almost_equal(reg.sigma_, covariance_matrix) + assert_allclose(reg.sigma_, covariance_matrix) def test_bayesian_sample_weights(): From cdf92225679b30e9af5676a805a07c57e3786786 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=A9mie=20du=20Boisberranger?= Date: Sun, 13 Apr 2025 16:27:59 +0200 Subject: [PATCH 9/9] add global_random_seed and increase tol [all random seeds] test_bayesian_covariance_matrix --- sklearn/linear_model/tests/test_bayes.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/sklearn/linear_model/tests/test_bayes.py b/sklearn/linear_model/tests/test_bayes.py index ef32dbec60a6a..9f7fabb749f52 100644 --- a/sklearn/linear_model/tests/test_bayes.py +++ b/sklearn/linear_model/tests/test_bayes.py @@ -96,17 +96,19 @@ def test_bayesian_ridge_parameter(): @pytest.mark.parametrize("n_samples, n_features", [(10, 20), (20, 10)]) -def test_bayesian_covariance_matrix(n_samples, n_features): +def test_bayesian_covariance_matrix(n_samples, n_features, global_random_seed): """Check the posterior covariance matrix sigma_ Non-regression test for https://github.com/scikit-learn/scikit-learn/issues/31093 """ - X, y = datasets.make_regression(n_samples, n_features) + X, y = datasets.make_regression( + n_samples, n_features, random_state=global_random_seed + ) reg = BayesianRidge(fit_intercept=False).fit(X, y) covariance_matrix = np.linalg.inv( reg.lambda_ * np.identity(n_features) + reg.alpha_ * np.dot(X.T, X) ) - assert_allclose(reg.sigma_, covariance_matrix) + assert_allclose(reg.sigma_, covariance_matrix, rtol=1e-6) def test_bayesian_sample_weights():