From d2e1af4b368df0164209f62b6d30e47304dc087a Mon Sep 17 00:00:00 2001 From: Joshua Date: Wed, 19 Aug 2026 08:31:28 -0400 Subject: [PATCH 1/3] feat: add random_state parameter for reproducible initialization VBPCA's fallback random initialization (loadings/scores when no explicit init is given) and its auto-generated xprobe mask were both effectively unseeded, with no way for callers to control or reproduce them. Root cause was subtler than a single missing seed: _initialize_parameters decided whether to pass an RNG using a truthiness check on the raw `init` option. Since _build_options()'s default is the *string* "random" (truthy), the common default-configuration path was actually taking the "no new seed" branch and falling through to init_params' own internal fallback of a hardcoded `np.random.default_rng(0)` -- silently deterministic, but not user-controllable, and inconsistent with the genuinely-unseeded branch that fired when `init=None` was passed explicitly. Both "random" and None normalize to the same "no fixture" case one level down in init_params -> _normalize_init, so they should behave identically. Fix: VBPCA.__init__ gains `random_state: int | np.random.Generator | None = None`, threaded through _build_options()/_full_update.py to unconditionally seed the fallback init RNG via np.random.default_rng(random_state), and threaded into fit()'s auto xprobe-mask generation the same way. random_state follows sklearn's convention: None (default) draws fresh entropy each call. This changes the default (unseeded) case from silently-deterministic-at- seed-0 to genuinely random each call, matching what "no explicit seed" should mean and what several existing tests' names/comments already assumed was happening. Fixed 12 existing tests that implicitly depended on the old accidental determinism -- either regression tests comparing against a hardcoded value (add random_state=0, which reproduces the exact old default-path behavior since it was already np.random.default_rng(0) under the hood), or tests comparing two independent calls that need a shared seed to be comparable at all. --- src/vbpca_py/_full_update.py | 14 ++++++++++---- src/vbpca_py/_pca_full.py | 1 + src/vbpca_py/estimators.py | 18 +++++++++++++++++- tests/test_model_selection.py | 9 +++++++-- tests/test_pca_full.py | 13 +++++++++++++ tests/test_sparse_explicit_mask.py | 1 + 6 files changed, 49 insertions(+), 7 deletions(-) diff --git a/src/vbpca_py/_full_update.py b/src/vbpca_py/_full_update.py index 6826c61..a6a53bd 100644 --- a/src/vbpca_py/_full_update.py +++ b/src/vbpca_py/_full_update.py @@ -669,15 +669,21 @@ def _initialize_parameters( # noqa: PLR0914 Returns: Initialized parameters and centered data/probe matrices. """ - # Use a deterministic RNG only when we fall back to random init; when - # init is provided (e.g., MATLAB fixture), pass through without forcing a - # new seed. + # Seed from the user-controllable random_state option. Both the + # "random" sentinel string and an omitted/None init normalize to "no + # fixture" inside init_params -> _normalize_init (empty dict), so both + # draw from this rng for any field a real fixture (mapping or .mat + # path) doesn't supply. random_state=None (default) yields fresh + # entropy each call, matching sklearn's random_state convention. init_value = cast("str | Mapping[str, Any] | None", ctx.opts.get("init")) + random_state = cast( + "int | np.random.Generator | None", ctx.opts.get("random_state") + ) init_result: InitResult = init_params( init_value, ctx.shapes, score_pattern_index=ctx.pattern_index, - rng=None if init_value else np.random.default_rng(), + rng=np.random.default_rng(random_state), ) loadings = init_result.a scores = init_result.s diff --git a/src/vbpca_py/_pca_full.py b/src/vbpca_py/_pca_full.py index 2174676..62715cf 100644 --- a/src/vbpca_py/_pca_full.py +++ b/src/vbpca_py/_pca_full.py @@ -1991,6 +1991,7 @@ def _build_options(kwargs: Mapping[str, object]) -> dict[str, object]: "hp_vb": 0.001, "hp_v": 0.001, "va_init": 1000.0, + "random_state": None, "earlystop": False, "rmsstop": np.array([100, 1e-4, 1e-3]), "cfstop": np.array([]), diff --git a/src/vbpca_py/estimators.py b/src/vbpca_py/estimators.py index 7739d02..e5ae208 100644 --- a/src/vbpca_py/estimators.py +++ b/src/vbpca_py/estimators.py @@ -39,6 +39,7 @@ def __init__( # noqa: PLR0913 xprobe_fraction: float = 0.0, criterion_order: list[str] | None = None, convergence_criteria: dict[str, bool] | None = None, + random_state: int | np.random.Generator | None = None, **opts: object, ) -> None: """ @@ -67,6 +68,13 @@ def __init__( # noqa: PLR0913 are criterion names, values are booleans. Criteria set to ``False`` are still evaluated for diagnostics but excluded from the stop decision. Defaults to all enabled. + random_state: Seed (or existing ``Generator``) controlling + internal randomness: parameter initialization when no + explicit ``init`` value is supplied, and the automatically + generated probe mask when ``xprobe_fraction`` is positive + and no explicit *xprobe* is passed to :meth:`fit`. ``None`` + (default) draws fresh entropy each call, matching prior + behaviour; set an int for reproducible fits. **opts: Additional options passed to the underlying PCA_FULL implementation. """ self.n_components = n_components @@ -82,6 +90,7 @@ def __init__( # noqa: PLR0913 self.xprobe_fraction = xprobe_fraction self.criterion_order = criterion_order self.convergence_criteria = convergence_criteria + self.random_state = random_state self.opts = opts self.components_: np.ndarray | None = None self.scores_: np.ndarray | None = None @@ -132,6 +141,7 @@ def get_params(self, *, deep: bool = True) -> dict[str, object]: # noqa: ARG002 "xprobe_fraction": self.xprobe_fraction, "criterion_order": self.criterion_order, "convergence_criteria": self.convergence_criteria, + "random_state": self.random_state, } params.update(self.opts) return params @@ -161,6 +171,7 @@ def set_params(self, **params: object) -> VBPCA: "xprobe_fraction", "criterion_order", "convergence_criteria", + "random_state", } for key, value in params.items(): if key in valid_params: @@ -214,11 +225,16 @@ def fit( # noqa: C901, PLR0912, PLR0914, PLR0915 opts["criterion_order"] = self.criterion_order if self.convergence_criteria is not None: opts["convergence_criteria"] = self.convergence_criteria + opts["random_state"] = self.random_state opts.update(self.opts) if xprobe is not None: opts["xprobe"] = xprobe elif self.xprobe_fraction > 0.0: - x, xprobe_gen = make_xprobe_mask(x, fraction=self.xprobe_fraction) + x, xprobe_gen = make_xprobe_mask( + x, + fraction=self.xprobe_fraction, + rng=np.random.default_rng(self.random_state), + ) opts["xprobe"] = xprobe_gen max_dense_bytes = resolve_max_dense_bytes( diff --git a/tests/test_model_selection.py b/tests/test_model_selection.py index ea6ec56..008de82 100644 --- a/tests/test_model_selection.py +++ b/tests/test_model_selection.py @@ -37,6 +37,7 @@ def test_select_n_components_tracks_trace_and_best_model() -> None: config=cfg, maxiters=80, verbose=0, + random_state=0, ) assert len(trace) == 3 @@ -322,8 +323,9 @@ def test_select_n_components_mask_argument_matches_nan_mask() -> None: x = rng.standard_normal((5, 8)) x[rng.random(x.shape) < 0.2] = np.nan mask = ~np.isnan(x) - # Supply an empty xprobe to suppress auto-holdout (which would differ - # between calls due to independent RNG states). + # Supply an empty xprobe to suppress auto-holdout, and a fixed + # random_state, so both calls share the same init/holdout draws and + # only the mask-argument form under test differs. empty_probe = np.full(x.shape, np.nan, dtype=float) cfg = SelectionConfig(metric="cost", compute_explained_variance=False) @@ -338,6 +340,7 @@ def test_select_n_components_mask_argument_matches_nan_mask() -> None: compat_mode="strict_legacy", rotate2pca=0, xprobe=empty_probe, + random_state=0, ) best_k_explicit, _, trace_explicit, _ = select_n_components( @@ -350,6 +353,7 @@ def test_select_n_components_mask_argument_matches_nan_mask() -> None: compat_mode="strict_legacy", rotate2pca=0, xprobe=empty_probe, + random_state=0, ) assert best_k_implicit == best_k_explicit @@ -574,6 +578,7 @@ def test_select_n_components_deterministic_across_num_cpu() -> None: rotate2pca=0, num_cpu=num_cpu, runtime_tuning="off", + random_state=0, ) cost_trace = [float(entry["cost"]) for entry in trace] res.append((best_k, cost_trace)) diff --git a/tests/test_pca_full.py b/tests/test_pca_full.py index 4c259a3..973eb3f 100644 --- a/tests/test_pca_full.py +++ b/tests/test_pca_full.py @@ -228,6 +228,7 @@ def test_pca_full_mask_argument_matches_implicit_mask() -> None: verbose=0, compat_mode="strict_legacy", rotate2pca=0, + random_state=0, ) out_explicit = pca_full( x, @@ -239,6 +240,7 @@ def test_pca_full_mask_argument_matches_implicit_mask() -> None: verbose=0, compat_mode="strict_legacy", rotate2pca=0, + random_state=0, ) rms_imp = np.asarray(out_implicit["lc"]["rms"], dtype=float) @@ -265,6 +267,7 @@ def test_pca_full_mask_argument_respects_eps_for_zeros_strict_legacy() -> None: verbose=0, compat_mode="strict_legacy", rotate2pca=0, + random_state=0, ) out_explicit = pca_full( x, @@ -276,6 +279,7 @@ def test_pca_full_mask_argument_respects_eps_for_zeros_strict_legacy() -> None: verbose=0, compat_mode="strict_legacy", rotate2pca=0, + random_state=0, ) rms_imp = np.asarray(out_implicit["lc"]["rms"], dtype=float) @@ -299,6 +303,7 @@ def _run() -> tuple[np.ndarray, np.ndarray, np.ndarray]: verbose=0, compat_mode="strict_legacy", rotate2pca=0, + random_state=0, ) lc_rms = np.asarray(out["lc"]["rms"], dtype=float) A = np.asarray(out["A"], dtype=float) @@ -356,6 +361,7 @@ def test_pca_full_strict_legacy_regression_rms_value() -> None: rotate2pca=1, display=0, verbose=0, + random_state=0, ) lc_rms = float(np.asarray(out["lc"]["rms"], dtype=float)[-1]) @@ -378,6 +384,7 @@ def test_pca_full_strict_legacy_regression_rms_trace_multiple_k() -> None: rotate2pca=1, display=0, verbose=0, + random_state=0, ) rms_vals.append(float(np.asarray(out["lc"]["rms"], dtype=float)[-1])) @@ -433,6 +440,7 @@ def test_pca_full_uniquesv_pattern_sharing_matches_no_uniquesv() -> None: compat_mode="strict_legacy", rotate2pca=1, verbose=0, + random_state=0, ) out_yes = pca_full( x, @@ -442,6 +450,7 @@ def test_pca_full_uniquesv_pattern_sharing_matches_no_uniquesv() -> None: compat_mode="strict_legacy", rotate2pca=1, verbose=0, + random_state=0, ) rms_no = float(np.asarray(out_no["lc"]["rms"], dtype=float)[-1]) @@ -531,6 +540,7 @@ def test_pca_full_auto_pattern_masked_aligns_with_uniquesv() -> None: rotate2pca=1, auto_pattern_masked=1, verbose=0, + random_state=0, ) out_uniquesv = pca_full( x, @@ -540,6 +550,7 @@ def test_pca_full_auto_pattern_masked_aligns_with_uniquesv() -> None: rotate2pca=1, uniquesv=1, verbose=0, + random_state=0, ) rms_auto = float(np.asarray(out_auto["lc"]["rms"], dtype=float)[-1]) @@ -627,6 +638,7 @@ def test_pca_full_sparse_vs_dense_rms_and_v_close() -> None: rotate2pca=1, verbose=0, mask=mask, + random_state=0, ) out_sparse = pca_full( sparse, @@ -636,6 +648,7 @@ def test_pca_full_sparse_vs_dense_rms_and_v_close() -> None: rotate2pca=1, verbose=0, mask=mask_sparse, + random_state=0, ) rms_dense = float(np.asarray(out_dense["lc"]["rms"], dtype=float)[-1]) diff --git a/tests/test_sparse_explicit_mask.py b/tests/test_sparse_explicit_mask.py index 36a4433..6c03b91 100644 --- a/tests/test_sparse_explicit_mask.py +++ b/tests/test_sparse_explicit_mask.py @@ -71,6 +71,7 @@ def test_vbpca_runs_with_sparse_mask_and_retains_observed_zero(): maxiters=10, compat_mode="modern", verbose=0, + random_state=0, ) model.fit(x, mask=mask) recon = np.asarray(model.inverse_transform(), dtype=float) From e0e52cea1fef6ce021411e792ce341bc75ee9ad6 Mon Sep 17 00:00:00 2001 From: Joshua Date: Wed, 19 Aug 2026 08:33:36 -0400 Subject: [PATCH 2/3] test: verify random_state reproducibility Covers: default is None, forwarded into resolved options, same seed -> bit-identical fit, different seeds -> diverge, None -> genuinely non-reproducible across calls, the auto xprobe-mask path specifically, accepting an existing Generator instance (not just an int), and get_params()/set_params() round-tripping. Also fixes get_options() -- it independently rebuilds the same options dict fit() constructs (rather than calling a shared helper), and was missing the random_state wire-up added to fit() in the previous commit, so it always reported random_state=None regardless of what was configured. Caught by test_random_state_forwarded_through_estimator. --- src/vbpca_py/estimators.py | 1 + tests/test_estimators.py | 99 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 100 insertions(+) diff --git a/src/vbpca_py/estimators.py b/src/vbpca_py/estimators.py index e5ae208..20ccd0a 100644 --- a/src/vbpca_py/estimators.py +++ b/src/vbpca_py/estimators.py @@ -387,6 +387,7 @@ def get_options(self) -> dict[str, object]: opts["criterion_order"] = self.criterion_order if self.convergence_criteria is not None: opts["convergence_criteria"] = self.convergence_criteria + opts["random_state"] = self.random_state opts.update(self.opts) return _build_options(opts) diff --git a/tests/test_estimators.py b/tests/test_estimators.py index fb2a1db..77b0b04 100644 --- a/tests/test_estimators.py +++ b/tests/test_estimators.py @@ -358,6 +358,105 @@ def test_xprobe_fraction_no_probe_when_zero() -> None: assert np.isnan(model.prms_) +# ── Reproducible initialization via random_state (issue #109) ─── + + +def test_random_state_default_is_none() -> None: + """Default random_state is None (fresh entropy each call).""" + model = VBPCA(n_components=2) + assert model.random_state is None + + +def test_random_state_forwarded_through_estimator() -> None: + """random_state kwarg reaches _build_options via VBPCA.""" + model = VBPCA(n_components=2, random_state=42) + resolved = model.get_options() + assert resolved["random_state"] == 42 + + +def test_random_state_reproducible_across_fits() -> None: + """Two fits with the same random_state produce identical results.""" + rng = np.random.default_rng(0) + x = rng.standard_normal((10, 20)) + mask = rng.random(x.shape) > 0.1 + + model_a = VBPCA(n_components=2, maxiters=10, verbose=0, random_state=42) + model_a.fit(x, mask=mask) + model_b = VBPCA(n_components=2, maxiters=10, verbose=0, random_state=42) + model_b.fit(x, mask=mask) + + np.testing.assert_array_equal(model_a.components_, model_b.components_) + np.testing.assert_array_equal(model_a.scores_, model_b.scores_) + + +def test_random_state_different_seeds_diverge() -> None: + """Different random_state values produce different initializations.""" + rng = np.random.default_rng(0) + x = rng.standard_normal((10, 20)) + mask = rng.random(x.shape) > 0.1 + + model_a = VBPCA(n_components=2, maxiters=10, verbose=0, random_state=1) + model_a.fit(x, mask=mask) + model_b = VBPCA(n_components=2, maxiters=10, verbose=0, random_state=2) + model_b.fit(x, mask=mask) + + assert not np.array_equal(model_a.components_, model_b.components_) + + +def test_random_state_none_is_not_reproducible() -> None: + """random_state=None (default) draws fresh entropy on every fit.""" + rng = np.random.default_rng(0) + x = rng.standard_normal((10, 20)) + mask = rng.random(x.shape) > 0.1 + + model_a = VBPCA(n_components=2, maxiters=10, verbose=0) + model_a.fit(x, mask=mask) + model_b = VBPCA(n_components=2, maxiters=10, verbose=0) + model_b.fit(x, mask=mask) + + assert not np.array_equal(model_a.components_, model_b.components_) + + +def test_random_state_reproducible_with_xprobe_fraction() -> None: + """random_state also seeds the auto-generated xprobe mask.""" + rng = np.random.default_rng(7) + x = rng.standard_normal((8, 20)) + + model_a = VBPCA(n_components=2, maxiters=5, xprobe_fraction=0.10, random_state=3) + model_a.fit(x) + model_b = VBPCA(n_components=2, maxiters=5, xprobe_fraction=0.10, random_state=3) + model_b.fit(x) + + assert model_a.prms_ == pytest.approx(model_b.prms_) + np.testing.assert_array_equal(model_a.components_, model_b.components_) + + +def test_random_state_accepts_generator_instance() -> None: + """random_state also accepts an existing np.random.Generator.""" + rng = np.random.default_rng(0) + x = rng.standard_normal((10, 20)) + + model_a = VBPCA(n_components=2, maxiters=5, random_state=np.random.default_rng(11)) + model_a.fit(x) + model_b = VBPCA(n_components=2, maxiters=5, random_state=11) + model_b.fit(x) + + np.testing.assert_array_equal(model_a.components_, model_b.components_) + + +def test_get_params_includes_random_state() -> None: + """random_state is included in get_params() output.""" + model = VBPCA(n_components=2, random_state=5) + assert model.get_params()["random_state"] == 5 + + +def test_set_params_random_state() -> None: + """set_params() can update random_state.""" + model = VBPCA(n_components=2) + model.set_params(random_state=9) + assert model.random_state == 9 + + # ── Convergence diagnostics (issue #99) ───────────────────────── From db52da30111d7a8675775c3f0f1fcaaa5ae1c554 Mon Sep 17 00:00:00 2001 From: Joshua Date: Wed, 19 Aug 2026 08:37:43 -0400 Subject: [PATCH 3/3] docs: changelog and limitations entry for random_state --- CHANGELOG.md | 8 ++++++++ docs/limitations.md | 2 ++ 2 files changed, 10 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index c22dcd8..696a878 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,14 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] + +### Added +- `random_state` constructor kwarg on `VBPCA`: seeds parameter initialization and any auto-generated xprobe mask (`int`, `np.random.Generator`, or `None`, following the sklearn convention). Surfaced via `get_params()`/`set_params()`/`get_options()` (#109). + +### Changed +- **Behavior change:** the default (`random_state=None`) now draws fresh entropy on every `fit()` call. Previously, default initialization was silently seeded with a fixed value regardless of configuration, so repeated fits produced identical results without any way to request a different draw. Pass `random_state=` for reproducible runs (#109). + ## [0.3.0] - 2026-08-17 ### Added diff --git a/docs/limitations.md b/docs/limitations.md index 0a035ae..4adedbb 100644 --- a/docs/limitations.md +++ b/docs/limitations.md @@ -11,3 +11,5 @@ - **RMS oscillation with uncentered data.** When `bias=True` (the default) and the input data has non-zero feature means, the RMS convergence trace can exhibit a stable period-2 oscillation caused by a one-iteration lag between the mean update and the reconstruction error. **Workaround:** center your data before fitting — use `MissingAwareStandardScaler` (or `AutoEncoder`) as a preprocessing step. Pre-centered data eliminates the oscillation entirely, even with `bias=True`. + +- **Fits are non-reproducible unless seeded.** `VBPCA(random_state=None)` (the default) draws fresh entropy for parameter initialization and any auto-generated xprobe mask on every call to `fit()`, so repeated fits on the same data can converge to different results. Pass an `int` or `np.random.Generator` via `random_state` for reproducible runs. Prior to #109, the default initialization was silently seeded with a fixed value; this is no longer the case.