fix: seed VBPCA initialization via random_state (#109)

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=<int>` for reproducible runs (#109).
+
 ## [0.3.0] - 2026-08-17
 
 ### Added

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.

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

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([]),

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(
@@ -371,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)
 

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) ─────────────────────────
 
 

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))

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])

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)