[josemal98]

During my visit to Devika Narain's lab I needed to run each fit with several seeds to report how robust the embeddings are.

pyRATS fixes its random seed internally: seven calls to np.random.seed(42), plus random_state=42 in KernelPCA and KMeans. This PR exposes that seed as a parameter, RATS(random_state=42).

• By default nothing changes. The default is the 42 used so far, and every random draw takes the same numbers in the same order, so results are identical bit for bit (checked below against e3283f1).
• A different seed can now be passed from outside. Users can check how much an embedding depends on the seed, or reproduce a fit made with another one. Until now the only way was to patch np.random.seed.
• pyRATS no longer modifies NumPy's global random state. Every fit used to leave np.random freshly seeded with 42, so any unseeded random draw made afterwards by the calling code came out the same in every run. This is the only change in behaviour: code that relied on it should seed NumPy itself.

API

random_state : int, RandomState instance, Generator or None, default=42, with scikit-learn's semantics:

value	effect
42 (default)	the results of earlier versions, bit for bit
another int, 0 to 2**32 − 1	that seed everywhere: the same seed gives the same embedding
None	a new seed from NumPy's global generator at every fit
a RandomState or Generator	one seed drawn from it per fit; two fresh generators with the same seed give the same embedding

The default is 42 rather than scikit-learn's None, so that results, and the regression CI, stay deterministic. Generator is accepted although scikit-learn does not accept it yet, since the value is reduced to one int per fit anyway.

What changes

The seed is resolved once at the start of fit_transform and passed explicitly to every place that draws random numbers:

where	random draw	when it matters
lpca	start vector of the truncated SVD (ARPACK) of every neighborhood	always (default kernel=None)
kpca	KernelPCA(random_state=...)	with a kernel
RATS._fit_nbrhd_graph	KMeans(random_state=...)	n_forced_clusters > 1
spectral_alignment	start vector of eigsh	global_init_algo_name="spectral"
compute_far_off_points	the shuffle that picks the far-off points	n_repel > 0
compute_spectral_color_of_pts_on_tear	start vector of eigsh in every tear component	compute_color_of_pts_on_tear

• One stream per fit. The draws that shared NumPy's global stream across functions (the noise seeds of the intermediate views, the far-off shuffle and the spectral start vector) now share one RandomState per fit. It is reseeded wherever the global generator was, so the shuffle continues the same stream as before.
• Workers. The seed reaches the loky workers of the local PCA as an argument, so the result does not depend on n_jobs.
• Tear coloring. compute_color_of_pts_on_tear uses the seed of the last fit, so the gluing colors match that fit.
• Removed calls. The two np.random.seed(42) calls with no draw after them (compute_final_embedding, spectrum_of_laplacian_from_neighbors) are removed. Param's noise draws, active only when noise_var is set, use a local generator with the same seed.

Tests and verification

• tests/test_random_state.py: 13 tests, about 17 s on one core, with small data generated in the file. They check that:

  • the same seed gives the same embedding, and another seed gives another one;
  • the default is 42;
  • RandomState, Generator and None supply one seed per fit;
  • n_jobs=1 and n_jobs=2 agree bit for bit;
  • the tear coloring uses the seed of the fit;
  • a fit leaves NumPy's global state untouched;
  • invalid values raise ValueError.

• ci.yml: one more step runs them. Drop it if you prefer.

• The default is unchanged. I compared this branch with e3283f1, bit for bit (np.array_equal), on these runs:

  • the 56 configurations of the full regression suite (tests/scripts/generate_snapshots.py);
  • three fits with the default constructor on 5 000 points with several cores, so that the local PCA runs in worker processes;
  • nine small configurations: the default, and others that reach the other random draws. These are kernel PCA, the spectral initialization, n_forced_clusters=2, n_repel > 0 with one and with two clusters, min_cluster_size=1, the fast CI configuration and a torn Klein bottle with its tear coloring.

  In all 68 comparisons, y, the tear colors, Utilde, C, n_Utilde_Utilde, Utildeg and c were identical. Two repeated runs of e3283f1 were also identical to each other, so the comparison is not blurred by run-to-run noise. The only combination that fails, n_forced_clusters=2 with the spectral initialization, fails the same way in both versions (see the side notes below).

  Both versions ran in the same environment: Linux, Python 3.13, NumPy 2.5.4, SciPy 1.18.1 and scikit-learn 1.9.1, with 3 cores for the large fits. In every fit, the old code left NumPy's global generator reseeded and the new code left it untouched.

• Seeds other than 42. With random_state=0 and random_state=1, the branch gives bit for bit what the old code gave when its internal 42 was swapped for those seeds by patching np.random.seed in every process. Fits made that way can be repeated with the new parameter.

Example: how much does an embedding depend on the seed?

A different seed always changes the bits. Most often the two embeddings differ only by a rotation or a reflection, but not always. These numbers come from my machine; the details may differ with another BLAS:

• small_spherewithhole (1 986 points), default parameters: seed 0 gives the embedding of seed 42 up to a rotation or a reflection. Seeds 1 and 2 give another embedding, 14% of the diameter away (median per point, after a rigid alignment).
• The Klein bottle of the README (5 000 points, with tearing): seeds 0, 1 and 2 stay within 0.6% of seed 42.

A single seed would not show the first case. The README now has a short section with this snippet:

import numpy as np
from scipy.linalg import orthogonal_procrustes
from scipy.spatial.distance import pdist
from pyRATS import rats

def shift(a, b):
    """Median distance between corresponding points of two embeddings after a
    rigid alignment (rotation, reflection, translation), in units of a's diameter."""
    a, b = a - a.mean(0), b - b.mean(0)
    R, _ = orthogonal_procrustes(b, a)
    return float(np.median(np.linalg.norm(b @ R - a, axis=1)) / pdist(a).max())

Y = {s: rats.RATS(random_state=s).fit_transform(X) for s in (42, 0, 1, 2)}
print({s: round(shift(Y[42], Y[s]), 3) for s in (0, 1, 2)})

Two side notes, not changed here

• With n_forced_clusters > 1 and global_init_algo_name="spectral", the fit fails with TypeError: 'tuple' object does not support item assignment. spectral_alignment returns a tuple (y, Z), and compute_init_embedding assigns it to y and indexes it for the next cluster. I am happy to send a separate fix.
• CONTRIBUTING.md gives 8e27b77 as the regression baseline, while ci.yml uses e862055.

Prepared with the help of Claude Code; reviewed and tested by me.