Ensemble Stability#

stability #

CLASS DESCRIPTION
EstimatorResult

Results from a single estimator fit.

SpectralStabilityEnsemble

Ensemble-based stability assessment via perturbation analysis.

Attributes#

Classes#

EstimatorResult #

Bases: NamedTuple

Results from a single estimator fit.

PARAMETER DESCRIPTION
eigenvalues

Eigenvalues from the estimator, shape (out_features,) or None.

TYPE: Tensor | None

eigenvectors

Eigenvectors from the estimator, shape (n_samples, out_features) or None.

TYPE: Tensor | None

kernel

Kernel matrix from the estimator, shape (n_samples, n_samples) or None.

TYPE: Tensor | None

SpectralStabilityEnsemble(estimator: ModelType, perturbation_type: Literal['noise', 'subsample', 'bootstrap'] = 'noise', n_trials: int = 10, noise_scale: float = 0.05, subsample_fraction: float = 0.8, metric: Literal['eigenvector_alignment', 'eigenvalue', 'kernel_frobenius', 'kernel_alignment'] = 'eigenvector_alignment', show_progress: bool = True, fit_manager: FitManager | None = None) #

Bases: Ensemble

Ensemble-based stability assessment via perturbation analysis.

Evaluates robustness of spectral methods by fitting multiple perturbed versions of an estimator and comparing results. Useful for determining optimal number of components and validating embedding quality.

PARAMETER DESCRIPTION
estimator

Estimator to analyze.

TYPE: ModelType

perturbation_type

Type of perturbation for stability testing.

  • "noise": Add Gaussian noise to data
  • "subsample": Random subsampling without replacement
  • "bootstrap": Bootstrap resampling with replacement

TYPE: str, by default "noise" DEFAULT: 'noise'

n_trials

Number of perturbation trials to run.

TYPE: int, by default 10 DEFAULT: 10

noise_scale

Standard deviation of Gaussian noise (for perturbation_type="noise").

TYPE: float, by default 0.05 DEFAULT: 0.05

subsample_fraction

Fraction of data to retain in subsampling (for "subsample" and "bootstrap").

TYPE: float, by default 0.8 DEFAULT: 0.8

metric

Stability metric to compute.

  • "eigenvector_alignment": Dominant eigenspace alignment
  • "eigenvalue": Eigenvalue reconstruction error
  • "kernel_frobenius": Kernel Frobenius error
  • "kernel_alignment": Center kernel alignment (CKA)

TYPE: str, by default "eigenvector_alignment" DEFAULT: 'eigenvector_alignment'

show_progress

Whether to display a progress bar during stability analysis.

TYPE: bool, by default True DEFAULT: True

fit_manager

Manager for unified training of both Estimator and Parametric models. Required for using Parametric (PyTorch Lightning) models in stability analysis. See FitManager for configuration options.

TYPE: FitManager | None DEFAULT: None

ATTRIBUTE DESCRIPTION
scores

Stability scores from all trials, shape (n_trials,).

TYPE: Tensor

results

Results from each perturbed trial.

TYPE: list

reference_result

Result from original unperturbed data.

TYPE: dict

Examples:

Analyze stability of diffusion map under noise:

>>> from spectre.decomposition import DiffusionMap
>>> from spectre.kernel import GaussianKernel
>>> from spectre.pairwise_distance import PairwiseDistanceEuclidean
>>> from spectre.ensemble import SpectralStabilityEnsemble
>>>
>>> kernel_fn = GaussianKernel(bw_method="median")
>>> dm = DiffusionMap(
...     kernel_fn=kernel_fn,
...     distance_fn=PairwiseDistanceEuclidean(),
...     out_features=10,
... )
>>>
>>> stability = SpectralStabilityEnsemble(
...     estimator=dm,
...     perturbation_type="noise",
...     n_trials=20,
...     noise_scale=0.05,
... )

Analyze stability with parametric model:

>>> from spectre.parametric import SpectralMap
>>> sm = SpectralMap(model=..., kernel_fn=..., ...)
>>> stability = SpectralStabilityEnsemble(
...     estimator=sm,
...     perturbation_type="bootstrap",
...     n_trials=50,
... )
>>> stability.fit(X)

Compare stability across methods:

>>> dm_stability = SpectralStabilityEnsemble(dm, n_trials=30)
>>> pca_stability = SpectralStabilityEnsemble(pca, n_trials=30)
>>>
>>> dm_stability.fit(X)
>>> pca_stability.fit(X)
>>>
>>> print(f"DM stability: {dm_stability.score():.3f}")
>>> print(f"PCA stability: {pca_stability.score():.3f}")
METHOD DESCRIPTION
fit

Run stability analysis on input data.

predict

Return stability scores as prediction.

score

Return mean stability score.

compute

Compute stability metric between reference and perturbed results.

Source code in spectre/ensemble/stability.py
def __init__(
    self,
    estimator: ModelType,
    perturbation_type: Literal["noise", "subsample", "bootstrap"] = "noise",
    n_trials: int = 10,
    noise_scale: float = 0.05,
    subsample_fraction: float = 0.8,
    metric: Literal[
        "eigenvector_alignment",
        "eigenvalue",
        "kernel_frobenius",
        "kernel_alignment",
    ] = "eigenvector_alignment",
    show_progress: bool = True,
    fit_manager: FitManager | None = None,
):
    super().__init__(estimators=[estimator], fit_manager=fit_manager)

    if not isinstance(estimator, ModelType):
        raise TypeError(
            f"`estimator` must be Estimator or Parametric, "
            f"got {type(estimator).__name__}."
        )
    self.estimator = estimator

    allowed_perturbations = ["noise", "subsample", "bootstrap"]
    if perturbation_type not in allowed_perturbations:
        raise ValueError(
            f"Unknown perturbation_type '{perturbation_type}'. "
            f"Options: {allowed_perturbations}."
        )
    self.perturbation_type = perturbation_type

    if not isinstance(n_trials, int):
        raise TypeError(
            f"n_trials must be an integer, got {type(n_trials).__name__}."
        )
    check_in_interval(n_trials, "[1, inf)")
    self.n_trials = n_trials

    if not isinstance(noise_scale, (int, float)):
        raise TypeError(
            f"noise_scale must be a number, got {type(noise_scale).__name__}."
        )
    check_in_interval(noise_scale, "(0, inf)")
    self.noise_scale = float(noise_scale)

    if not isinstance(subsample_fraction, (int, float)):
        raise TypeError(
            f"subsample_fraction must be a number, "
            f"got {type(subsample_fraction).__name__}."
        )
    check_in_interval(subsample_fraction, "(0, 1]")
    self.subsample_fraction = float(subsample_fraction)

    if metric == "eigenvalue":
        if not hasattr(estimator, "eigenvalues"):
            raise ValueError(
                f"Estimator must have 'eigenvalues' attribute for {metric} metric."
            )
    elif metric == "eigenvector_alignment":
        if not hasattr(estimator, "eigenvectors"):
            raise ValueError(
                f"Estimator must have 'eigenvectors' attribute for {metric} metric."
            )
    elif metric in ["kernel_frobenius", "kernel_alignment"]:
        if not hasattr(estimator, "kernel"):
            raise ValueError(
                f"Estimator must have 'kernel' attribute for {metric} metric."
            )
    else:
        raise ValueError(f"Unknown metric '{metric}'.")
    self.metric = metric

    if metric in ["kernel_frobenius", "kernel_alignment"] and perturbation_type in [
        "subsample",
        "bootstrap",
    ]:
        raise ValueError(
            f"Kernel metrics ('{metric}') are incompatible with "
            f"`perturbation_type='{perturbation_type}'`. Use 'noise' perturbation "
            f"or switch to 'eigenvalue' or 'eigenvector_alignment' metrics."
        )

    if not isinstance(show_progress, bool):
        raise TypeError(
            f"show_progress must be a boolean, got {type(show_progress).__name__}."
        )
    self.show_progress = show_progress

    # Results storage
    self._scores: torch.Tensor | None = None
    self._estimator_reference: ModelType | None = None
    self._estimators: list[EstimatorResult] = []
Attributes#
scores: torch.Tensor property #

Stability scores from all trials.

RETURNS DESCRIPTION
Tensor

Stability scores of shape (n_trials,).

results: list[EstimatorResult] property #

Results from each trial.

estimator_reference: ModelType property #

Result from original data.

Functions#
fit(X: torch.Tensor, weights: torch.Tensor | None = None, target: torch.Tensor | None = None) -> SpectralStabilityEnsemble #

Run stability analysis on input data.

PARAMETER DESCRIPTION
X

Input data of shape (n_samples, in_features).

TYPE: Tensor

weights

Sample weights of shape (n_samples, ).

TYPE: Tensor | None DEFAULT: None

target

Target values for supervised methods.

TYPE: Tensor | None DEFAULT: None

RETURNS DESCRIPTION
SpectralStabilityEnsemble

Returns self for method chaining.

Source code in spectre/ensemble/stability.py
def fit(
    self,
    X: torch.Tensor,
    weights: torch.Tensor | None = None,
    target: torch.Tensor | None = None,
) -> "SpectralStabilityEnsemble":
    """
    Run stability analysis on input data.

    Parameters
    ----------
    X : torch.Tensor
        Input data of shape `(n_samples, in_features)`.

    weights : torch.Tensor | None, optional
        Sample weights of shape `(n_samples, )`.

    target : torch.Tensor | None, optional
        Target values for supervised methods.

    Returns
    -------
    SpectralStabilityEnsemble
        Returns self for method chaining.
    """
    check_2d(X)

    self.in_features = X.shape[1]

    # Fit reference estimator on original data
    estimator_ref = copy.deepcopy(self.estimator)
    self._estimator_reference = self._fit_estimator(
        estimator_ref, X, weights=weights, target=target, estimator_name="reference"
    )

    self.out_features = estimator_ref.out_features

    scores = []
    self._results = []

    trials_range = range(self.n_trials)
    if self.show_progress:
        desc = f"Spectral stability ({self.perturbation_type})"
        trials_range = tqdm(trials_range, desc=desc, unit="trial")

    for trial_idx in trials_range:
        # Generate perturbed data
        X_p, weights_p = self._perturb(X, weights, X.shape[0])

        # Fit estimator on perturbed data
        estimator_p = copy.deepcopy(self.estimator)
        self._fit_estimator(
            estimator_p,
            X_p,
            weights=weights_p,
            target=target,
            estimator_name=f"trial_{trial_idx}",
        )

        result = self.compute(self._estimator_reference, estimator_p)
        self._results.append(result)

        scores.append(result.score)

        if self.show_progress and isinstance(trials_range, tqdm):
            mean_score = torch.tensor(scores).mean().item()
            trials_range.set_postfix(
                stability=f"{mean_score:.3f}", current=f"{result.score:.3f}"
            )

    self._scores = torch.tensor(scores)
    self.is_fitted = True

    return self
predict(X: torch.Tensor) -> torch.Tensor #

Return stability scores as prediction.

PARAMETER DESCRIPTION
X

Input data (not used, stability scores are independent of new data).

TYPE: Tensor

RETURNS DESCRIPTION
Tensor

Stability scores of shape (n_trials,). Higher values indicate better stability.

Source code in spectre/ensemble/stability.py
def predict(self, X: torch.Tensor) -> torch.Tensor:
    """
    Return stability scores as prediction.

    Parameters
    ----------
    X : torch.Tensor
        Input data (not used, stability scores are independent of new data).

    Returns
    -------
    torch.Tensor
        Stability scores of shape `(n_trials,)`. Higher values indicate
        better stability.
    """
    self.check_fitted()
    return self._scores
score(X: torch.Tensor | None = None, weights: torch.Tensor | None = None, target: torch.Tensor | None = None) -> float #

Return mean stability score.

PARAMETER DESCRIPTION
X

Not used (stability is independent of new data).

TYPE: Tensor | None DEFAULT: None

weights

Not used.

TYPE: Tensor | None DEFAULT: None

target

Not used.

TYPE: Tensor | None DEFAULT: None

RETURNS DESCRIPTION
float

Mean stability score across all trials. Higher is more stable.

Source code in spectre/ensemble/stability.py
def score(
    self,
    X: torch.Tensor | None = None,
    weights: torch.Tensor | None = None,
    target: torch.Tensor | None = None,
) -> float:
    """
    Return mean stability score.

    Parameters
    ----------
    X : torch.Tensor | None, optional
        Not used (stability is independent of new data).

    weights : torch.Tensor | None, optional
        Not used.

    target : torch.Tensor | None, optional
        Not used.

    Returns
    -------
    float
        Mean stability score across all trials. Higher is more stable.
    """
    self.check_fitted()
    return self._scores.mean().item()
compute(estimator_reference: ModelType, estimator_p: ModelType) -> EstimatorResult #

Compute stability metric between reference and perturbed results.

RETURNS DESCRIPTION
EstimatorResult

Results including stability score and estimator outputs if applicable.

Source code in spectre/ensemble/stability.py
def compute(
    self,
    estimator_reference: ModelType,
    estimator_p: ModelType,
) -> EstimatorResult:
    """
    Compute stability metric between reference and perturbed results.

    Returns
    -------
    EstimatorResult
        Results including stability score and estimator outputs if applicable.
    """
    if self.metric == "eigenvalue":
        eigenvalues_r = estimator_reference.eigenvalues
        eigenvalues_p = estimator_p.eigenvalues

        min_len = min(len(eigenvalues_r), len(eigenvalues_p))
        eigenvalues_r = eigenvalues_r[:min_len]
        eigenvalues_p = eigenvalues_p[:min_len]

        # Use correlation as stability measure (higher = more stable)
        error = eigenvalue_reconstruction_error(
            eigenvalues_p, eigenvalues_r, reduce="correlation"
        )
        # Correlation is in [-1, 1], map to [0, 1]
        score = (error + 1.0) / 2.0

    elif self.metric == "eigenvector_alignment":
        eigenvectors_r = estimator_reference.eigenvectors
        eigenvectors_p = estimator_p.eigenvectors

        # Match dimensions
        min_samples = min(eigenvectors_r.shape[0], eigenvectors_p.shape[0])
        min_components = min(eigenvectors_r.shape[1], eigenvectors_p.shape[1])

        eigenvectors_r = eigenvectors_r[:min_samples, :min_components]
        eigenvectors_p = eigenvectors_p[:min_samples, :min_components]

        # Subspace overlap (already in [0, 1])
        score = eigenvector_alignment(
            eigenvectors_r,
            eigenvectors_p,
            n_eigen=min_components,
            reduce="subspace_overlap",
        )

    elif self.metric.startswith("kernel"):
        kernel_r = estimator_reference.kernel
        kernel_p = estimator_p.kernel

        check_same_shape(kernel_r, kernel_p)

        if self.metric == "kernel_frobenius":
            # Frobenius error (lower is better, map to [0, 1])
            error = kernel_frobenius_error(kernel_p, kernel_r, normalize=True)
            score = 1.0 - min(error, 1.0)  # Cap at 1.0

        elif self.metric == "kernel_alignment":
            # Kernel alignment (already in [0, 1])
            score = kernel_alignment(kernel_r, kernel_p)
        else:
            raise ValueError(f"Unknown kernel metric: {self.metric}")

    else:
        raise ValueError(f"Unknown metric: {self.metric}")

    return EstimatorResult(
        eigenvalues=estimator_p.eigenvalues,
        eigenvectors=estimator_p.eigenvectors,
        kernel=estimator_p.kernel,
        score=score,
    )

Functions#