Add seed and sigma parameters to Solver.eigs() #112
Description
Feature Description
Currently, Solver.eigs() calls scipy.sparse.linalg.eigsh() without setting its v0 parameter, denoting the initialisation vector. v0=None by default, leading scipy to create for its ARPACK routines a random intialisation vector of values sampled uniformly from [-1, 1] (see line 360 here). This sampling can be determined by another rng parameter, which accepts a numpy.random.Generator instance. However, since Solver.eigs() again uses the default rng=None, computed eigenvalues and eigenvectors are not exactly the same across multiple calls. Adding some seed: int | np.ndarray | None = None parameter would maintain current behaviour by default, while allowing users to regenerate identical arrays of eigenvectors and eigenvalues:
# Setup intitialization vector
if seed is None or isinstance(seed, int):
rng = np.random.default_rng(seed)
v0 = None
else:
v0 = np.asarray_chkfinite(seed)
if v0.shape != (self.mass.shape[0],):
raise ValueError(f"seed must be None, an integer, or an array of shape ({self.mass.shape[0]},), matching the number of mesh vertices.")
rng = None
eigenvalues, eigenvectors = eigsh(
self.stiffness, k, self.mass, sigma=sigma, OPinv=op_inv, v0=v0, rng=rng
)
The else branch is less important but allows the user to directly set the initialisation vector, for whatever reason. If this is not desired, the setup can just be rng = np.random.default_rng(seed).
Solver.eigs() also sets sigma = -0.01 internally, such that the first eigenvalues and eigenvectors are always returned and are in order. Adding a sigma: float = -0.01 parameter would allow users to compute higher eigenvalues/eigenvectors by increasing sigma, without having to compute the preceding components of the basis. Since this no longer guarantees that outputs are sorted by eigenvalue, sorting afterwards may be desirable:
if sigma > -0.01:
sort_inds = np.argsort(eigenvalues)
eigenvalues = eigenvalues[sort_inds]
eigenvectors = eigenvectors[:, sort_inds]
Motivation
Adding seed would support reproducible analyses and consistent unit tests involving Solver.eigs(). Adding sigma would reduce compute/time when only a few higher eigenvalues/eigenvectors are desired.
Additional Context
If you are interested in either of these parameters, I am happy to write up a pull request and also add tests to ensure that the reproducibility and/or sorting behave as intended.