splisosm.kernel

splisosm.kernel#

Spatial kernel abstractions.

Classes#

FFTKernel

FFT-based spatial kernel on a periodic 2D raster grid.

IdentityKernel

Identity kernel K = I (placeholder when spatial kernel is skipped).

SpatialCovKernel

Graph-based spatial covariance kernel.

Module Contents#

class splisosm.kernel.FFTKernel(shape, spacing=(1.0, 1.0), rho=0.99, neighbor_degree=1, workers=None)#

Bases: Kernel

FFT-based spatial kernel on a periodic 2D raster grid.

This implementation supports a CAR-style spatial kernel equivalent to a periodic, neighbourhood-graph-based autoregressive model. All operations are \(O(N \log N)\) via the 2-D Fast Fourier Transform.

Parameters:

  • shape (tuple[int, int]) – Grid shape (ny, nx).

  • spacing (tuple[float, float]) – Physical spacing (dy, dx) between neighbouring raster cells.

  • rho (float) – Spatial autocorrelation coefficient in CAR kernel.

  • neighbor_degree (int) – Neighbour ring degree for graph construction. 1 uses nearest neighbours in the periodic metric.

  • workers (int | None) – Number of workers used by scipy.fft.fft2.

  • state. (Initialize kernel-specific)

apply_residual_op(x, epsilon)#

Apply the kernel regression residual operator.

Computed in O(N log N) via FFT as:

R @ v = IFFT2( epsilon / (lambda + epsilon) * FFT2(v) )
Parameters:

  • x (ndarray) – Input with shape (ny, nx) or (ny, nx, m).

  • epsilon (float) – Regularisation / noise level.

Returns:

Residuals of the same shape as x.

Return type:

np.ndarray

eigenvalues(k=None)#

Return kernel eigenvalues in descending order.

Parameters:

k (int | None) – Number of leading eigenvalues to return. None returns all.

Returns:

Eigenvalues in descending order.

Return type:

np.ndarray

power_spectral_density_1d(bins=50)#

Compute the 1D power spectral density (radial profile).

Parameters:

bins (int) – Number of bins for the 1D radial frequency.

Returns:

  • freq_bins (np.ndarray) – The centre frequencies of the valid bins.

  • psd_1d (np.ndarray) – The average power (eigenvalue) in each frequency bin.

Return type:

tuple[ndarray, ndarray]

abstractmethod realization()#

Return the dense kernel matrix (not recommended for large grids).

Return type:

Tensor

square_trace()#

Return tr(K^2).

Return type:

float

trace()#

Return tr(K).

Return type:

float

xtKx(x)#

Compute x^T K x in O(N log N) via FFT.

Parameters:

x (ndarray) – Input with shape (ny, nx) or (ny, nx, m).

Returns:

Scalar for 2D input, or shape (m,) for 3D input.

Return type:

float or np.ndarray

xtKx_approx(x, k=None)#

Same as xtKx() (FFT is already exact for periodic grids).

Parameters:

k (int | None)

n_grid#
neighbor_degree = 1#
rho#
spectrum#
workers = None#
class splisosm.kernel.IdentityKernel(n_spots)#

Bases: Kernel

Identity kernel K = I (placeholder when spatial kernel is skipped).

All spots are treated as independent; no spatial smoothing is applied. This is used when setup_data(skip_spatial_kernel=True) is called and spatial variability testing is not needed.

Parameters:

  • n_spots (int) – Number of spatial locations.

  • state. (Initialize kernel-specific)

eigenvalues(k=None)#

Return k ones (all eigenvalues of I are 1).

Parameters:

k (int | None)

Return type:

Tensor

realization()#

Return the identity matrix of shape (n_spots, n_spots).

Return type:

Tensor

square_trace()#

Return tr(I²) = n_spots.

Return type:

Tensor

trace()#

Return tr(I) = n_spots.

Return type:

Tensor

xtKx(x)#

Compute x^T I x = x^T x.

Parameters:

x (Tensor)

Return type:

Tensor

xtKx_approx(x, k=None)#

Compute x^T I x = x^T x (K = I; no approximation required).

Parameters:
Return type:

Tensor

xtKx_exact(x)#

Compute x^T I x = x^T x (exact, same as xtKx()).

Parameters:

x (Tensor)

Return type:

Tensor

class splisosm.kernel.SpatialCovKernel(coords=None, *, adj_matrix=None, k_neighbors=4, rho=0.99, standardize_cov=True, centering=False)#

Bases: Kernel

Graph-based spatial covariance kernel.

Construct from spot coordinates (from_coordinates() or the coords positional argument) or from a pre-built k-NN adjacency / weight matrix (from_adjacency() or the adj_matrix keyword argument). Exactly one of the two must be provided.

The kernel is stored in one of two modes depending on dataset size:

  • Dense mode (n ≤ DENSE_THRESHOLD): the (n, n) covariance matrix K_sp is materialised, with optional variance standardisation and double-centring applied at construction time.

  • Implicit mode (n > DENSE_THRESHOLD): only the sparse precision inv_cov = K⁻¹ is retained; K x is computed via a cached sparse LU factorisation.

Eigenvalue decomposition and the low-rank Q factor are computed lazily on the first call to eigenvalues().

Construct from spot coordinates or a pre-built adjacency matrix.

Exactly one of coords or adj_matrix must be supplied. Prefer the factory methods from_coordinates() and from_adjacency() for more descriptive call sites.

Parameters:

  • coords (array-like of shape (n_spots, n_dim), optional) – Spot coordinates. A k-NN graph is built from these and used to construct the CAR precision matrix.

  • adj_matrix (sparse or dense matrix of shape (n_spots, n_spots), optional) – Pre-built symmetric k-NN adjacency / weight matrix. The CAR precision is built directly from this matrix.

  • k_neighbors (int) – Number of nearest neighbours when building the graph from coords. Ignored when adj_matrix is provided.

  • rho (float) – CAR spatial autocorrelation coefficient (0 < rho < 1).

  • standardize_cov (bool) – Scale the covariance to unit marginal variance.

  • centering (bool) – Apply double-centring H K H.

eigendecomposition()#

Compute and cache the full eigendecomposition of the dense kernel.

Populates K_eigvals, K_eigvecs, and Q.

Raises:

ValueError – If called on an implicit kernel (no K_sp).

Return type:

None

eigenvalues(k=None)#

Return leading eigenvalues in descending order.

For dense kernels a full eigendecomposition is triggered once and cached. For implicit kernels a partial eigsh via a LinearOperator is used; the low-rank factor Q is cached for subsequent xtKx() calls.

Parameters:

k (int | None) – Number of leading eigenvalues. None returns all cached values.

Returns:

Eigenvalues in descending order.

Return type:

Tensor

classmethod from_adjacency(adj_matrix, rho=0.99, standardize_cov=True, centering=False)#

Build a CAR kernel from a pre-built k-NN adjacency / weight matrix.

Use this when the spatial graph has already been constructed (e.g. from external tools) and you want to skip the coordinate-based k-NN step.

Parameters:

  • adj_matrix (spmatrix | np.ndarray) – Symmetric adjacency or weight matrix of shape (n_spots, n_spots). Can be a dense np.ndarray or any scipy.sparse matrix.

  • rho (float) – CAR spatial autocorrelation coefficient (0 < rho < 1).

  • standardize_cov (bool) – Scale covariance to unit marginal variance.

  • centering (bool) – Apply double-centring H K H.

Return type:

SpatialCovKernel

classmethod from_coordinates(coords, k_neighbors=4, rho=0.99, standardize_cov=True, centering=False)#

Build a CAR spatial kernel from spot coordinates.

Parameters:

  • coords (ndarray | Tensor) – Spot coordinates of shape (n_spots, n_dim).

  • k_neighbors (int) – Number of nearest neighbours for the spatial graph.

  • rho (float) – CAR spatial autocorrelation coefficient (0 < rho < 1).

  • standardize_cov (bool) – Scale covariance to unit marginal variance.

  • centering (bool) – Apply double-centring H K H.

Return type:

SpatialCovKernel

rank()#

Return effective rank of the stored representation.

Return type:

int

realization()#

Return the realised dense kernel matrix.

For implicit kernels this forces dense inversion — prefer calling eigenvalues() first to populate Q, then use Q @ Q.T.

Returns:

Dense matrix of shape (n_spots, n_spots).

Return type:

Tensor

shape()#

Return kernel shape (n_spots, n_spots).

Return type:

tuple of int

square_trace()#

Return tr(K²) [or tr((HKH)²) when centering=True].

Uses ||K||_F² for the dense case. Falls back to a Hutchinson stochastic estimator for implicit kernels.

Returns:

Scalar squared-trace value.

Return type:

Tensor

trace()#

Return tr(K) [or tr(HKH) when centering=True].

Uses exact arithmetic for the dense case (K_sp is already centred). Falls back to a Hutchinson stochastic estimator for implicit kernels.

Returns:

Scalar trace value.

Return type:

Tensor

xtKx(x)#

Compute x^T K x using the cached low-rank factor when available.

When Q has been populated (by a prior eigenvalues() call) the result is x^T Q Q^T x, i.e. consistent with the approximation K ≈ Q Q^T:

  • Dense mode: eigendecomposition() produces a full-rank Q so Q Q^T = K exactly.

  • Implicit mode: eigenvalues() produces a rank-k Q; the result is therefore a rank-k approximation of x^T K x.

Before eigenvalues() is called (Q is None) this falls back to xtKx_exact(). Use xtKx_exact() directly to always obtain the exact quadratic form regardless of whether Q has been computed.

Parameters:

x (Tensor) – Input matrix of shape (n_spots, d).

Returns:

Matrix of shape (d, d).

Return type:

Tensor

See also

xtKx_exact

Always uses the exact path (dense multiply or LU solve).

xtKx_approx(x, k=None)#

Compute x^T K_k x via the top-k low-rank factor.

Calls eigenvalues() to ensure Q is populated for at least k components, then computes (x^T Q_k)(Q_k^T x) where Q_k = Q[:, :k]. The result is consistent with the rank-k eigenvalue approximation used in the Liu null distribution.

Parameters:

  • x (Tensor) – Input matrix of shape (n_spots, d).

  • k (int or None) – Number of leading eigenvalues to include. None uses all available eigenvalues (populates the full Q first).

Returns:

Approximated quadratic form of shape (d, d).

Return type:

Tensor

xtKx_exact(x)#

Compute x^T K x exactly, bypassing any cached low-rank factor.

Always uses the full kernel:

  • Dense mode (K_sp is not None): direct matrix multiply x^T K_{sp} x.

  • Implicit mode (K_sp is None): sparse LU solve M u = x where M = K^{-1}, then returns x^T u = x^T K x.

Unlike xtKx(), this method ignores Q and always returns the exact quadratic form under the full kernel K. This is useful when you need the exact HSIC statistic for reporting while using a low-rank Q for the p-value null distribution.

Parameters:

x (Tensor) – Input matrix of shape (n_spots, d).

Returns:

Matrix of shape (d, d).

Return type:

Tensor

See also

xtKx

Uses cached Q (low-rank approximation) when available.

DENSE_THRESHOLD: int = 5000#

Switch to implicit sparse mode above this number of spots.

K_eigvals: Tensor | None#

Cached eigenvalues in descending order. Populated lazily.

K_eigvecs: Tensor | None#

Cached eigenvectors (dense mode only). Populated lazily.

K_sp: Tensor | None#

Dense centred/standardised kernel (n_spots, n_spots). None in implicit mode.

Q: Tensor | None#

Low-rank factor (n_spots, rank) such that K ≈ Q @ Q.T. Populated lazily after the first eigenvalues() call.

inv_cov: csc_matrix#

Sparse precision matrix M = K⁻¹ of shape (n_spots, n_spots).