From d75efbf8747b51e6f71599435b9caa4103c8403e Mon Sep 17 00:00:00 2001 From: Mark Colley <39366189+M-Colley@users.noreply.github.com> Date: Wed, 7 Oct 2026 22:31:51 +0100 Subject: [PATCH 1/5] fix(denoise): learn the signal amplitude and normalize targets The Gaussian-process kernel used by `denoise=True` was `RBF + WhiteKernel + ConstantKernel`, with no amplitude on the RBF term and no `normalize_y`. The smooth component therefore had a fixed prior variance of 1, so for targets far from unit scale the fitted model attributed most of the signal to noise and over-smoothed it. For `y = A sin(2x) + noise` (noise sd = 0.1 A), the RMSE of the "denoised" targets against the true signal was 0.25x the noise level at A = 1, but 7.3x at A = 100 or 1e4, i.e. far worse than not denoising at all. Use `ConstantKernel() * RBF(...) + WhiteKernel(...)` with `normalize_y=True`, and standardize `X` (applying the same transform to `Xresampled`) so the initial length scales suit any input units. The error is now 0.23-0.25x the noise level at every amplitude tested. Co-Authored-By: Claude Opus 5.5 --- pysr/denoising.py | 24 +++++++++++++++++------- pysr/test/test_main.py | 21 +++++++++++++++++++++ 2 files changed, 38 insertions(+), 7 deletions(-) diff --git a/pysr/denoising.py b/pysr/denoising.py index fbf17468b..4d24b3bca 100644 --- a/pysr/denoising.py +++ b/pysr/denoising.py @@ -18,16 +18,26 @@ def denoise( from sklearn.gaussian_process import GaussianProcessRegressor from sklearn.gaussian_process.kernels import RBF, ConstantKernel, WhiteKernel - gp_kernel = RBF(np.ones(X.shape[1])) + WhiteKernel(1e-1) + ConstantKernel() + # Standardize the inputs so the initial length scales suit any units of + # `X`; `normalize_y` below does the same for the targets. + X_mean = X.mean(axis=0) + X_scale = X.std(axis=0) + X_scale[X_scale == 0] = 1.0 + + # The `ConstantKernel` factor learns the signal variance. Without it the + # smooth component has unit variance, so targets far from unit scale are + # mostly attributed to noise. + gp_kernel = ConstantKernel() * RBF(np.ones(X.shape[1])) + WhiteKernel(1e-1) gpr = GaussianProcessRegressor( - kernel=gp_kernel, n_restarts_optimizer=50, random_state=random_state + kernel=gp_kernel, + normalize_y=True, + n_restarts_optimizer=50, + random_state=random_state, ) - gpr.fit(X, y) + gpr.fit((X - X_mean) / X_scale, y) - if Xresampled is not None: - return Xresampled, cast(ndarray, gpr.predict(Xresampled)) - - return X, cast(ndarray, gpr.predict(X)) + X_out = X if Xresampled is None else Xresampled + return X_out, cast(ndarray, gpr.predict((X_out - X_mean) / X_scale)) def multi_denoise( diff --git a/pysr/test/test_main.py b/pysr/test/test_main.py index 453eaa50a..f2c77ff06 100644 --- a/pysr/test/test_main.py +++ b/pysr/test/test_main.py @@ -37,6 +37,7 @@ jl, load_all_packages, ) +from pysr.denoising import denoise from pysr.export_latex import sympy2latex from pysr.export_sympy import pysr2sympy from pysr.feature_selection import _handle_feature_selection, run_feature_selection @@ -1607,6 +1608,26 @@ def test_all_selection_strategies(self): self.assertEqual(idx_score, 1) +class TestDenoising(unittest.TestCase): + def test_denoise_learns_signal_amplitude(self): + rstate = np.random.RandomState(0) + X = rstate.uniform(-3, 3, size=(200, 1)) + Xresampled = np.linspace(-2.5, 2.5, 50)[:, None] + for amplitude in [1.0, 100.0]: + with self.subTest(amplitude=amplitude): + noise = 0.1 * amplitude + y = amplitude * np.sin(2 * X[:, 0]) + noise * rstate.randn(len(X)) + X_out, y_out = denoise( + X, y, Xresampled=Xresampled, random_state=np.random.RandomState(0) + ) + self.assertIs(X_out, Xresampled) + rmse = np.sqrt( + np.mean((y_out - amplitude * np.sin(2 * Xresampled[:, 0])) ** 2) + ) + # Denoising should beat the raw data's error (`noise`) comfortably: + self.assertLess(rmse, 0.5 * noise) + + class TestFeatureSelection(unittest.TestCase): def setUp(self): self.rstate = np.random.RandomState(0) From a901dbb669c3ca3e5c49b2b2f4cf22026f5890ef Mon Sep 17 00:00:00 2001 From: Mark Colley <39366189+M-Colley@users.noreply.github.com> Date: Thu, 8 Oct 2026 20:05:05 +0100 Subject: [PATCH 2/5] perf(denoise): use 5 optimizer restarts instead of 50 Each restart re-runs the marginal-likelihood optimization from random hyperparameters, at O(n^3) per optimizer step, so denoising took minutes for a few hundred rows (about 55-63 s at n=250 and 10 min at n=500 with three features). With standardized inputs and normalized targets, the fits were identical across 0, 2, 10, and 50 restarts on every case tested (RMSE against the true signal of 0.377/0.397 at n=250 for all four, and 0.327/0.330 at n=500 for 2 and 10), so keep a few restarts as a guard against poor local optima and drop the rest. Co-Authored-By: Claude Opus 5.5 --- pysr/denoising.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pysr/denoising.py b/pysr/denoising.py index 4d24b3bca..ed3873d83 100644 --- a/pysr/denoising.py +++ b/pysr/denoising.py @@ -31,7 +31,10 @@ def denoise( gpr = GaussianProcessRegressor( kernel=gp_kernel, normalize_y=True, - n_restarts_optimizer=50, + # With standardized inputs and targets, extra restarts from random + # hyperparameters did not change the fit, but each one costs a full + # optimization, at O(n^3) per step: + n_restarts_optimizer=5, random_state=random_state, ) gpr.fit((X - X_mean) / X_scale, y) From e880e4e0bafa8ef07ed35e02de4b1bb7aee23d76 Mon Sep 17 00:00:00 2001 From: Mark Colley <39366189+M-Colley@users.noreply.github.com> Date: Thu, 8 Oct 2026 20:05:06 +0100 Subject: [PATCH 3/5] fix(denoise): match `Xresampled` columns to `X` by name `check_array(Xresampled)` dropped the column names, so when `X` and `Xresampled` were DataFrames with columns in different orders, the denoised targets were evaluated at, and the search was run on, rows whose columns were silently swapped (e.g. a temperature column read as a fraction). Select the columns of `X` from `Xresampled` by name, and raise a clear error when some are missing. Also reject `weights` together with `Xresampled` when denoising: the denoised targets have one row per row of `Xresampled`, so the weights no longer line up, which previously failed with a bare `AssertionError`. Co-Authored-By: Claude Opus 5.5 --- pysr/sr.py | 15 +++++++++++++++ pysr/test/test_main.py | 20 ++++++++++++++++++++ 2 files changed, 35 insertions(+) diff --git a/pysr/sr.py b/pysr/sr.py index 961cb5225..755265a1a 100644 --- a/pysr/sr.py +++ b/pysr/sr.py @@ -2102,6 +2102,21 @@ def _validate_and_set_fit_params( y_units, ) + if Xresampled is not None and self.denoise and weights is not None: + raise ValueError( + "`weights` cannot be combined with `Xresampled`, as the denoised " + "targets are evaluated at the rows of `Xresampled` instead of `X`." + ) + if isinstance(X, pd.DataFrame) and isinstance(Xresampled, pd.DataFrame): + # Match the columns of `Xresampled` to those of `X` by name: + missing_columns = X.columns.difference(Xresampled.columns) + if len(missing_columns) > 0: + raise ValueError( + "`Xresampled` is missing columns of `X`: " + f"{list(missing_columns)}." + ) + Xresampled = Xresampled[X.columns] + if isinstance(X, pd.DataFrame): if variable_names: variable_names = None diff --git a/pysr/test/test_main.py b/pysr/test/test_main.py index f2c77ff06..e42d6d626 100644 --- a/pysr/test/test_main.py +++ b/pysr/test/test_main.py @@ -1627,6 +1627,26 @@ def test_denoise_learns_signal_amplitude(self): # Denoising should beat the raw data's error (`noise`) comfortably: self.assertLess(rmse, 0.5 * noise) + def test_xresampled_columns_are_matched_by_name(self): + rstate = np.random.RandomState(0) + X = pd.DataFrame({"T": rstate.uniform(300, 400, 20), "x": rstate.rand(20)}) + y = X["T"] * X["x"] + Xresampled = pd.DataFrame( + {"x": rstate.rand(5), "T": rstate.uniform(300, 400, 5)} + ) + model = PySRRegressor(denoise=True) + Xresampled_validated = model._validate_and_set_fit_params( + X, y, Xresampled, None, None, None, None, None + )[2] + np.testing.assert_array_equal( + Xresampled_validated, Xresampled[["T", "x"]].to_numpy() + ) + + with self.assertRaisesRegex(ValueError, "missing columns"): + model.fit(X, y, Xresampled=Xresampled[["x"]]) + with self.assertRaisesRegex(ValueError, "`weights` cannot be combined"): + model.fit(X, y, Xresampled=Xresampled, weights=np.ones(len(y))) + class TestFeatureSelection(unittest.TestCase): def setUp(self): From 41831f319daf04bc92161ce139fa30f7024871e3 Mon Sep 17 00:00:00 2001 From: MilesCranmerBot Date: Thu, 8 Oct 2026 22:11:29 +0100 Subject: [PATCH 4/5] fix(denoise): trim to amplitude, column-order, and weights fixes Learn the GP signal amplitude with ConstantKernel() * RBF and normalize_y=True, keeping 50 optimizer restarts and unscaled inputs. Match a DataFrame Xresampled to X by column name, and reject weights combined with Xresampled under denoise. Register TestDenoising in runtests() so CI runs it. Co-authored-by: Miles Cranmer --- pysr/denoising.py | 22 ++++++---------------- pysr/sr.py | 5 +++-- pysr/test/test_main.py | 5 ++--- 3 files changed, 11 insertions(+), 21 deletions(-) diff --git a/pysr/denoising.py b/pysr/denoising.py index ed3873d83..bfff0978e 100644 --- a/pysr/denoising.py +++ b/pysr/denoising.py @@ -18,29 +18,19 @@ def denoise( from sklearn.gaussian_process import GaussianProcessRegressor from sklearn.gaussian_process.kernels import RBF, ConstantKernel, WhiteKernel - # Standardize the inputs so the initial length scales suit any units of - # `X`; `normalize_y` below does the same for the targets. - X_mean = X.mean(axis=0) - X_scale = X.std(axis=0) - X_scale[X_scale == 0] = 1.0 - - # The `ConstantKernel` factor learns the signal variance. Without it the - # smooth component has unit variance, so targets far from unit scale are - # mostly attributed to noise. gp_kernel = ConstantKernel() * RBF(np.ones(X.shape[1])) + WhiteKernel(1e-1) gpr = GaussianProcessRegressor( kernel=gp_kernel, normalize_y=True, - # With standardized inputs and targets, extra restarts from random - # hyperparameters did not change the fit, but each one costs a full - # optimization, at O(n^3) per step: - n_restarts_optimizer=5, + n_restarts_optimizer=50, random_state=random_state, ) - gpr.fit((X - X_mean) / X_scale, y) + gpr.fit(X, y) + + if Xresampled is not None: + return Xresampled, cast(ndarray, gpr.predict(Xresampled)) - X_out = X if Xresampled is None else Xresampled - return X_out, cast(ndarray, gpr.predict((X_out - X_mean) / X_scale)) + return X, cast(ndarray, gpr.predict(X)) def multi_denoise( diff --git a/pysr/sr.py b/pysr/sr.py index 755265a1a..02e962edb 100644 --- a/pysr/sr.py +++ b/pysr/sr.py @@ -2108,7 +2108,6 @@ def _validate_and_set_fit_params( "targets are evaluated at the rows of `Xresampled` instead of `X`." ) if isinstance(X, pd.DataFrame) and isinstance(Xresampled, pd.DataFrame): - # Match the columns of `Xresampled` to those of `X` by name: missing_columns = X.columns.difference(Xresampled.columns) if len(missing_columns) > 0: raise ValueError( @@ -2840,12 +2839,14 @@ def fit( Resampled training data, of shape (n_resampled, n_features), to generate a denoised data on. This will be used as the training data, rather than `X`. + If both are DataFrames, columns are matched to `X` by name. weights : ndarray | pandas.DataFrame Weight array of the same shape as `y`. Each element is how to weight the mean-square-error loss for that particular element of `y`. Alternatively, if a custom `loss` was set, it will can be used - in arbitrary ways. + in arbitrary ways. Cannot be combined with `Xresampled` + when `denoise=True`. variable_names : list[str] A list of names for the variables, rather than "x0", "x1", etc. If `X` is a pandas dataframe, the column names will be used diff --git a/pysr/test/test_main.py b/pysr/test/test_main.py index e42d6d626..39a13acbd 100644 --- a/pysr/test/test_main.py +++ b/pysr/test/test_main.py @@ -1617,14 +1617,12 @@ def test_denoise_learns_signal_amplitude(self): with self.subTest(amplitude=amplitude): noise = 0.1 * amplitude y = amplitude * np.sin(2 * X[:, 0]) + noise * rstate.randn(len(X)) - X_out, y_out = denoise( + _, y_out = denoise( X, y, Xresampled=Xresampled, random_state=np.random.RandomState(0) ) - self.assertIs(X_out, Xresampled) rmse = np.sqrt( np.mean((y_out - amplitude * np.sin(2 * Xresampled[:, 0])) ** 2) ) - # Denoising should beat the raw data's error (`noise`) comfortably: self.assertLess(rmse, 0.5 * noise) def test_xresampled_columns_are_matched_by_name(self): @@ -3086,6 +3084,7 @@ def runtests(just_tests=False): test_cases = [ TestPipeline, TestBest, + TestDenoising, TestFeatureSelection, TestMiscellaneous, TestHelpMessages, From 3a63b4c1a9b6a82e23d2daf7d2227c4abad262a0 Mon Sep 17 00:00:00 2001 From: MilesCranmerBot Date: Thu, 8 Oct 2026 22:23:04 +0100 Subject: [PATCH 5/5] fix(denoise): keep the original kernel; normalize_y alone fixes large-scale targets Co-authored-by: Miles Cranmer --- pysr/denoising.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pysr/denoising.py b/pysr/denoising.py index bfff0978e..97847cc94 100644 --- a/pysr/denoising.py +++ b/pysr/denoising.py @@ -18,7 +18,7 @@ def denoise( from sklearn.gaussian_process import GaussianProcessRegressor from sklearn.gaussian_process.kernels import RBF, ConstantKernel, WhiteKernel - gp_kernel = ConstantKernel() * RBF(np.ones(X.shape[1])) + WhiteKernel(1e-1) + gp_kernel = RBF(np.ones(X.shape[1])) + WhiteKernel(1e-1) + ConstantKernel() gpr = GaussianProcessRegressor( kernel=gp_kernel, normalize_y=True,