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, centering=False)#

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.
centering (bool) – If True, apply double-centring H K H with H = I - (1/n) 1 1^T. On a periodic grid the all-ones vector is the DFT DC eigenvector, so double-centring is equivalent to zeroing the (0, 0) entry of the spectrum; xtKx(), trace(), square_trace(), and eigenvalues() all reflect the centred kernel. Default False. HSIC-based SV/DU tests should set this to True for consistency with the double-centred HSIC formulation.
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, centering=False)#

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.
centering (bool, optional) – If True, apply double-centring H K H with H = I - (1/n) 1 1^T. For K = I this reduces to H itself (H is idempotent), whose spectrum has \(n - 1\) unit eigenvalues and one zero eigenvalue. Default False. HSIC-based SV/DU tests should set this to True for consistency with the double-centred HSIC formulation.
state. (Initialize kernel-specific)

eigenvalues(k=None)#

Return the leading eigenvalues of K.

For K = I all \(n\) eigenvalues are 1. For H I H = H there are \(n - 1\) unit eigenvalues and a single zero eigenvalue; only the nonzero ones are returned.

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

realization()#

Return I (or H = I - (1/n) 1 1^T when centering=True).

Return type:: Tensor

square_trace()#

Return tr(K²): n for I and n - 1 for H (idempotent).

Return type:: Tensor

trace()#

Return tr(K): n for I and n - 1 for H.

Return type:: Tensor

xtKx(x)#

Compute x^T K x where K = I (or H I H = H when centred).

Parameters:: x (Tensor)
Return type:: Tensor

xtKx_approx(x, k=None)#

Quadratic form; identical to xtKx() (no approximation needed).

Parameters:

x (Tensor)
k (int | None)

Return type:

Tensor

xtKx_exact(x)#

Exact quadratic form; identical to xtKx() for this kernel.

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