quadsv.kernels.base

quadsv.kernels.base#

Abstract bases for the spatial-kernel layer.

This module hosts the two ABCs:

Kernel: the universal interface every kernel implements (Kx / xtKx / xtKy / trace / square_trace / eigenvalues). Subclass this for a custom kernel that can be expressed through a single self._K buffer.
MatrixKernelBase: matrix-form kernels with dense / sparse / sparse-precision auto-switching. Subclass this for a new matrix backend; MatrixKernel is the standard concrete subclass.

Concrete classes live in sibling modules: quadsv.kernels.matrix, quadsv.kernels.fft, and quadsv.kernels.nufft.

Classes#

`Kernel`	Abstract base class shared by all spatial kernels in `quadsv`.
`MatrixKernelBase`	Concrete base for kernels backed by an explicit or implicit `n × n` matrix.

Module Contents#

class quadsv.kernels.base.Kernel(*args, centering=True, **kwargs)[source]#

Bases: abc.ABC

Abstract base class shared by all spatial kernels in quadsv.

Concrete backends:

MatrixKernel — explicit n×n kernel or its sparse precision matrix.
quadsv.fft.FFTKernel — grid kernel via its eigenvalue spectrum.
quadsv.nufft.NUFFTKernel — irregular-point kernel evaluated through a type-1 / type-2 NUFFT round-trip.

Required interface#

Every concrete kernel exposes:

n, method, params, centering — integer / string / dict / bool attributes.
xtKx(), xtKy(), Kx() — quadratic / bilinear / apply primitives.

Note

The empirical data centering (z-scoring) inside quadsv.spatial_q_test() / spatial_r_test() breaks independence across spatial obervations. As a result, the null distribution of the test statistic Q = Zᵀ K Z = Xᵀ (H K H) X / σ² with H = I - 𝟏𝟏ᵀ/n should inspect the spectrum of a centered kernel HKH. Every Kernel carries a centering flag (default True). Set centering=False to recover the raw K moments (useful for diagnostics or theoretical comparison).

abstract Kx(x)[source]#: K x (raw) or HKH x (centered).

abstract eigenvalues(k=None)[source]#

Eigenvalues of K (raw) or HKH (centered).

Parameters:: k (int | None)
Return type:: ndarray

abstract square_trace()[source]#

trace(K²) or trace((HKH)²) = trace(K²) − 2·s₂/n + s₁²/n².

Same contract as trace() re: backend-specific kwargs.

Return type:: float

abstract trace()[source]#

trace(K) or trace(HKH) = trace(K) − s₁/n.

Concrete backends are free to add kwargs for backend-specific options (e.g. MatrixKernelBase exposes n_probes for its precision-stored Hutchinson path).

Return type:: float

abstract xtKx(x)[source]#: xᵀ K x (raw) or xᵀ HKH x (centered).

abstract xtKy(x, y)[source]#: xᵀ K y (raw) or xᵀ HKH y (centered).

centering: bool[source]#

method: str[source]#

n: int[source]#

params: dict[source]#

Parameters:: centering (bool)

class quadsv.kernels.base.MatrixKernelBase(n, method='gaussian', *, centering=True, **kwargs)[source]#

Bases: Kernel

Concrete base for kernels backed by an explicit or implicit n × n matrix.

Subclass this when you want a custom way to construct K (e.g. a bespoke function of coordinates, a learnt kernel, or an ad-hoc adjacency matrix) while inheriting every downstream primitive ready-to-use.

Subclasses must implement _build_kernel() to return either the (n, n) kernel matrix (dense np.ndarray or scipy.sparse) or its precision matrix M = K^{-1}. If a precision matrix is returned, set self.stores_precision = True before calling super().__init__, or flip it inside _build_kernel while the buffer is being built; the MatrixKernel CAR / precomputed-inverse path shows a worked example.

Everything else — cached sparse LU for implicit solves, Hutchinson trace estimation when stores_precision=True, automatic sparse-vs-dense dispatch in xtKx / Kx, standardized-quadratic-form cache (_K_column_sums / xtKx_standardized()), pickle safety for the non-picklable LU factor — is already implemented here.

Handles dense, sparse, and implicit (operator-based) kernels. Switches between an explicit representation (the kernel matrix K is stored and used directly) and an implicit representation (the precision matrix M = K^{-1} is stored and linear systems are solved on demand) based on problem size.

Variables:

n (int) – Number of observations (n).
method (str) – Kernel method name (free-form; used only for diagnostics / provenance).
params (dict) – Resolved kernel parameters — subclasses decide what goes here.
stores_precision (bool) – If True, _K holds the precision matrix M = K^{-1} and linear solves (via a cached LU) are used for xtKx() and trace estimation. If False, _K is the realized kernel matrix (dense or sparse).

Parameters:

n (int)
method (str)
centering (bool)

Notes

The internal buffer _K stores the kernel matrix when stores_precision=False and the precision matrix K^{-1} when stores_precision=True. Public methods (xtKx(), trace(), square_trace(), eigenvalues()) transparently handle both cases; callers should not access _K directly.

Kx(x)[source]#

Apply the kernel operator to x, returning K @ x.

Single public primitive for kernel–vector products. Handles explicit (dense or sparse K) and implicit (precision matrix + cached LU) cases uniformly.

Parameters:: x (np.ndarray or scipy.sparse matrix) – (n,) or (n, M). Sparse inputs are densified internally because scipy.linalg.lu_solve / splu.solve require dense RHS and K @ x typically returns dense anyway.
Returns:: (n,) if x was 1D, else (n, M).
Return type:: np.ndarray

Examples

>>> import numpy as np
>>> from quadsv import MatrixKernel
>>> rng = np.random.default_rng(0)
>>> coords = rng.standard_normal((40, 2))
>>> kernel = MatrixKernel.from_coordinates(coords, method="matern")
>>> kernel.Kx(rng.standard_normal(40)).shape
(40,)

eigenvalues(k=None)[source]#

Eigenvalues of K or HKH (according to self.centering).

Three paths:

Dense ``K`` — form HKH = K − 𝟏 r^T − c 𝟏^T + m · 𝟏𝟏^T in closed form (r / c = row / col means, m = grand mean) and call numpy.linalg.eigvalsh(). Same O(n³) cost as the raw case, no extra memory.
Sparse explicit ``K`` — wrap H K H · v as a scipy.sparse.linalg.LinearOperator and call eigsh(). Preserves K’s sparsity — we never densify.
Implicit precision (sparse ``M = K⁻¹``) — wrap H K H · v = H (M⁻¹ (H v)) where M⁻¹· is the cached sparse-LU solve. Same sparsity benefit as the raw-inverse path.

Results are cached per centering mode; switching self.centering invalidates and recomputes.

Parameters:: k (int, optional) – Number of largest-magnitude eigenvalues to return. If None, returns all (dense) or max(6, n − 2) (sparse/implicit — limited by what eigsh can extract).
Returns:: Eigenvalues sorted in descending order.
Return type:: np.ndarray

realization()[source]#

Return the realized (n, n) kernel matrix.

When self.centering is False, returns K as stored (dense ndarray or sparse matrix). When self.centering is True, returns the centered operator HKH with H = I − 𝟏𝟏ᵀ/n — the operator that xtKx() / Kx() / trace() / eigenvalues() all actually apply. HKH is dense even when K is sparse (each row gets a column-mean subtracted), so the centered path always materializes a dense result.

Notes

If stores_precision is True, this forces expensive dense inversion of the precision matrix. Prefer xtKx() and trace() for implicit kernels.

Return type:: ndarray | spmatrix

square_trace(n_probes=None)[source]#

trace(K²) (raw) or trace((HKH)²) (centered).

Dispatch is automatic based on stores_precision:

Explicit ``K``: Frobenius norm of the stored K — Σᵢⱼ Kᵢⱼ² in O(nnz) for sparse, O(n²) for dense.
Precision-stored: (1/m) Σᵢ ‖K vᵢ‖² with the same cached ±1 probes as trace() (single LU solve shared across both moments).

-2·s₂/n + s₁²/n² is applied on top when centering=True.

Parameters:: n_probes (int, optional) – Probe count for the Hutchinson path (default 15). Ignored on the analytic / explicit-K path. Shared cache with trace().
Return type:: float

trace(n_probes=None)[source]#

trace(K) (raw) or trace(HKH) (centered).

Dispatch is automatic based on stores_precision:

Explicit ``K`` (dense or sparse): diagonal sum of the stored K — exact, O(n) (sparse) or O(n) (dense).
Precision-stored (CAR in the implicit regime, K only available through the LU solve on K⁻¹): Hutchinson estimator (1/m) Σᵢ vᵢᵀ (K vᵢ) with cached ±1 probes solved through the precision.

-s₁/n is applied on top when centering=True.

Parameters:: n_probes (int, optional) – Probe count for the Hutchinson path (default 15). Ignored on the analytic / explicit-K path (no probes involved). On the Hutchinson path the cache is keyed on n_probes; requesting a different count recomputes the LU-solve cache and emits a RuntimeWarning.
Return type:: float

xtKx(x)[source]#

Quadratic form x^T K x (paired diagonal for batched inputs).

Parameters:: x (np.ndarray or scipy.sparse matrix) – (n,) or (n, M). Sparse x is preserved through the final x^T (K x) contraction — only the right side K @ x needs a dense RHS for the solver / BLAS call.
Returns:: Scalar if 1D input, (M,) if batched.
Return type:: float or np.ndarray

xtKx_standardized(x, means, stds)[source]#

Compute z^T K z where z = (x - means) / stds without densifying sparse x.

Sparse-aware expansion using the raw-K primitives:

z^T K z = (x^T K x - 2·μ·(K·𝟏)^T·x + μ²·(𝟏^T K 𝟏)) / σ²

and the cached row sums of K. This is the fast path for standardizing large sparse feature matrices (e.g. scRNA-seq counts) before a Q-test.

Parameters:

x (np.ndarray or scipy.sparse matrix) – (n,) or (n, M). Columns correspond to features.
means (np.ndarray) – (M,) per-feature mean and std (ddof=1 to match quadsv.statistics.spatial_q_test()).
stds (np.ndarray) – (M,) per-feature mean and std (ddof=1 to match quadsv.statistics.spatial_q_test()).

Returns:

(M,) standardized quadratic form values. Columns with std <= 0 are returned as zero.

Return type:

np.ndarray

xtKy(x, y)[source]#

Bilinear form x^T K y (paired diagonal for batched inputs).

For (n, M) batches returns (M,) — the diagonal of X^T K Y in the same column order, matching quadsv.spatial_r_test(). Sparse x is preserved; only K @ y is densified.

Parameters:

x (np.ndarray or scipy.sparse matrix) – (n,) or (n, M) (must share the M).
y (np.ndarray or scipy.sparse matrix) – (n,) or (n, M) (must share the M).

Returns:

Scalar if 1D inputs; (M,) if batched.

Return type:

float or np.ndarray

Examples

>>> import numpy as np
>>> from quadsv import MatrixKernel
>>> rng = np.random.default_rng(0)
>>> coords = rng.standard_normal((40, 2))
>>> kernel = MatrixKernel.from_coordinates(coords, method="matern")
>>> x = rng.standard_normal(40)
>>> y = rng.standard_normal(40)
>>> isinstance(kernel.xtKy(x, y), float)
True

method: str = 'gaussian'[source]#: Kernel method name.

n: int[source]#: Number of observations (samples).

params: dict[source]#: Resolved kernel parameters after defaults are merged with user overrides.

stores_precision: bool = False[source]#: Whether the kernel is stored in precision form (True) or as the realized kernel matrix (False).