Choose a numerical compatibility mode¶

The 0.5 native estimator and public component-selection helpers default to compat_mode="modern". Use that policy for new native workflows, including observed zeros and small finite measurements. The explicit vbpca_py.legacy helpers, low-level solver and a directly constructed PCAOptions() retain strict_legacy as their default.

Mode and layout are separate choices. Native estimators always accept samples × features; the legacy adapter and low-level solver accept features × samples. Setting compat_mode="strict_legacy" on a native estimator does not restore the historical interface or serialized-model layout. The 0.5 migration guide covers those changes.

What changes between modes¶

Solver behavior	`strict_legacy`	`modern`
Sparse empty-axis detection without an explicit solver mask	Signed row/column sums; observed values can cancel to zero and remove an axis.	Stored non-NaN observations determine occupied axes, including stored zeros.
Zero and near-zero normalization	Historical `numpy.isclose(value, 0)` matches are replaced by float64 machine epsilon.	Finite values retain their original magnitude, including exact zero; no epsilon substitution is needed.
Observation coordinates in update paths	Retains historical nonzero-based operations, so zeros may be treated differently across preparation/probe paths.	Uses the observation mask through updates, probe exclusion and NaN cleanup.
`bias_update_order="auto"`	Resolves to `legacy`.	Resolves to `post_factor`, updating the mean/RMS after the factor sweep.
`variance_update_order="auto"`	Resolves to `legacy`.	Resolves to `post_rotation`, updating prior variances after component rotation.

An explicit bias_update_order or variance_update_order overrides that mode's ordering default. The latter controls prior variances, not observation noise variance. The algorithm description explains the relationship to returned components.

Machine epsilon here is a historical numerical replacement, not a missing-data code or a measurement threshold for new data. In modern mode, retain observed sparse zeros explicitly in storage. A sparse constructor applied to a dense array can drop those zeros; eliminate_zeros() would also remove their observation status. Use coordinates as in the executed example below. Without an explicit mask, absent sparse entries are missing in both modes; neither mode can infer that a discarded entry was an observed zero. A modern explicit mask can request an observed zero at an absent coordinate.

Dense NaNs are missing at the native estimator boundary. Native masks are consumed before the solver, and native NaNs remain missing even at positions marked true. The historical explicit dense-mask interface in the legacy adapter/low-level solver is authoritative: combine a supplied mask with ~numpy.isnan(X) to exclude existing NaNs. Changing numerical mode does not replace that interface contract. In either mode, observed probes discarded by empty training-axis removal are refused rather than silently dropped.

Split automatic probes without losing observed zeros¶

make_xprobe_mask(X, ..., compat_mode="modern") is the public helper's default. It removes exactly the sampled sparse coordinates and preserves every other stored zero. Stored NaNs are excluded from modern sampling; an explicit mask can request a zero at an absent coordinate. The mask only restricts probe sampling: it does not filter other entries from the returned training copy. Apply the same observation mask when fitting that copy. Neither input buffers nor the input mask are modified. The returned matrices retain the input orientation.

For an older helper-based reproduction, explicitly pass compat_mode="strict_legacy". That retains historical sparse sampling and subtraction, including the possible loss of unselected stored zeros and sampling of stored NaNs. Dense probe sampling is unchanged. The legacy estimator passes its configured mode to automatic probe generation, so its default strict-legacy path retains that behavior. Modern legacy fits also materialize mask-requested zeros directly without epsilon substitution. The native estimator keeps its existing dedicated sparse split in both modes. Component selection passes its requested mode when creating probes.

The executed example below checks the reported nine-entry sparse case: seed 19 holds out a nonzero observation, leaving eight training entries and the unselected stored zero. Using the strict helper instead leaves seven training entries. Compatibility mode does not qualify convergence or held-out coverage.

Migrate an existing workflow¶

For a frozen reproduction, retain the package version, explicit legacy entry point, mode, initialization/seed and all fitting options. Switching numerical mode can change both the observation set and update order; it is not a promise of identical fitted factors or convergence decisions.

For a new native workflow, arrange the observation table as samples × features and use VBPCA(...) or explicitly set compat_mode="modern". Pass aligned masks/probes in that layout. Keep the masking Pipeline responsible for missingness before feature expansion. Compare original-unit predictions and fitting diagnostics when migrating an established application. The stopping guide is a separate choice from compatibility mode; four-iteration examples below do not qualify convergence or calibration.

The following standalone source uses both entry points explicitly, checks that modern masked/unmasked tiny measurements agree, preserves a stored sparse zero and prints the resolved update-order policies. Its modern examples use the same numerical policy; the historical strict fit is shown without an equivalence assertion.

python scripts/example_compatibility_modes.py

"""Keep layout migration separate from numerical compatibility policy."""

from __future__ import annotations

import numpy as np
from scipy.sparse import csr_matrix

from vbpca_py import VBPCA, PCAOptions, make_xprobe_mask
from vbpca_py.legacy import VBPCA as LegacyVBPCA  # noqa: N811


def run_example():
    """Demonstrate modern observations and explicit feature-first migration.

    Returns:
        Models fitted to the same small synthetic observations.
    """
    x = np.random.default_rng(19).normal(size=(8, 3))
    x[:, 0] *= 1e-10
    x[0, 1], x[1, 2] = 0.0, np.nan
    options = {
        "maxiters": 4,
        "niter_broadprior": 0,
        "random_state": 19,
        "num_cpu": 1,
        "runtime_tuning": "off",
        "runtime_report": True,
        "verbose": 0,
    }
    native = VBPCA(1, **options).fit(x)
    modern_legacy = LegacyVBPCA(1, compat_mode="modern", **options).fit(x.T)
    # Explicit legacy masks are authoritative: exclude existing NaNs too.
    mask = ~np.isnan(x.T)
    modern_masked = LegacyVBPCA(1, compat_mode="modern", **options).fit(x.T, mask=mask)
    historical = LegacyVBPCA(1, **options).fit(x.T)
    # Construct coordinates so the observed numeric zero remains stored.
    rows, columns = np.nonzero(~np.isnan(x))
    sparse = csr_matrix((x[rows, columns], (rows, columns)), shape=x.shape)
    modern_sparse = VBPCA(1, **options).fit(sparse)
    np.testing.assert_allclose(native.components_, modern_legacy.components_.T)
    np.testing.assert_allclose(
        modern_masked.mean_, modern_legacy.mean_, rtol=1e-12, atol=0
    )
    np.testing.assert_allclose(
        native.mean_, modern_sparse.mean_, rtol=1e-10, atol=1e-14
    )
    probe_values = np.arange(9.0).reshape(3, 3)
    probe_rows, probe_cols = np.indices(probe_values.shape)
    probe_input = csr_matrix((
        probe_values.ravel(),
        (probe_rows.ravel(), probe_cols.ravel()),
    ))
    train, probe = make_xprobe_mask(probe_input, 0.1, np.random.default_rng(19))
    strict_train, _ = make_xprobe_mask(
        probe_input, 0.1, np.random.default_rng(19), compat_mode="strict_legacy"
    )
    np.testing.assert_equal((train.nnz, probe.nnz, strict_train.nnz), (8, 1, 7))
    np.testing.assert_equal(train.data[0], 0)
    return {
        "native": native,
        "modern_legacy": modern_legacy,
        "modern_masked": modern_masked,
        "historical": historical,
        "modern_sparse": modern_sparse,
    }


def main() -> None:
    models = run_example()
    print("Native default:", models["native"].compat_mode)
    print("Legacy default:", models["historical"].compat_mode)
    print("Low-level configuration default:", PCAOptions().compat_mode.value)
    for name in ["native", "historical"]:
        report = models[name].runtime_report_
        print(name, report["bias_update_order"], report["variance_update_order"])
    print("Modern masked values and explicit native/legacy layouts agree")


if __name__ == "__main__":
    main()