quadsv.fft#

FFT-based spatial testing for regular grid data.

class quadsv.fft.FFTKernel(shape: Tuple[int, int], spacing: Tuple[float, float] = (1.0, 1.0), topology: str = 'square', method: str = 'matern', workers: int | None = None, fft_solver: str = 'fft2', **kwargs)[source]#

Bases: object

FFT-accelerated spatial kernel for dense grid data.

Operates on evenly-spaced grid data (Raster data) with spectral decomposition via FFT with periodic (Torus) boundary conditions.

ny, nx

Grid dimensions.

Type:

int

topology#

Grid topology.

Type:

{‘square’, ‘hex’}

method#

Kernel method.

Type:

str

params#

Kernel parameters (bandwidth, nu, neighbor_degree, rho).

Type:

dict

spectrum#

Precomputed eigenvalues of the kernel matrix. Shape (n_rfft,) where n_rfft = ny * (nx//2 + 1).

Type:

np.ndarray

__init__(shape: Tuple[int, int], spacing: Tuple[float, float] = (1.0, 1.0), topology: str = 'square', method: str = 'matern', workers: int | None = None, fft_solver: str = 'fft2', **kwargs) → None[source]#

Initialize FFT-accelerated spatial kernel for grid data.

Parameters:

  • shape (tuple of int) – Grid dimensions (ny, nx).

  • spacing (tuple of float, default (1.0, 1.0)) – Physical distance between pixels (dy, dx).

  • topology ({'square', 'hex'}, default 'square') – Grid topology. ‘hex’ is for Visium-like hexagonal layouts.

  • method (str, default 'matern') – Kernel method: ‘gaussian’, ‘matern’, ‘moran’, ‘graph_laplacian’, ‘car’.

  • workers (Optional[int], default None) – Number of parallel workers for fft computations.

  • fft_solver ({'fft2', 'rfft2'}, default 'fft2') – FFT solver to use. ‘fft2’ (full FFT) or ‘rfft2’ (real FFT, ~50% memory). Default is ‘fft2’ for better compatibility and robustness on most architectures.

  • **kwargs (dict) – Kernel parameters (bandwidth, nu, neighbor_degree, rho).

Examples

>>> kernel = FFTKernel((64, 64), method='gaussian', bandwidth=2.0)
>>> kernel = FFTKernel((64, 64), topology='hex', method='matern')
__repr__() → str[source]#

Return a detailed, machine-readable representation of the FFTKernel.

Returns:

String representation in angle-bracket format.

Return type:

str

__str__() → str[source]#

Return a human-friendly, multi-line representation of the FFTKernel.

Returns:

Multi-line string summary.

Return type:

str

xtKx(x: ndarray) → float | ndarray[source]#

Compute the quadratic form x^T K x efficiently using FFT.

Uses Parseval’s theorem to compute the result in frequency domain for O(N log N) complexity instead of O(N²).

Parameters:

x (np.ndarray) – Input data tensor. Shape (ny, nx) for single feature or (ny, nx, M) for M features.

Returns:

Quadratic form value(s). Scalar if input was 2D, shape (M,) if input was 3D.

Return type:

float or np.ndarray

eigenvalues(k: int | None = None, return_full: bool = False) → ndarray[source]#

Get the eigenvalues of the kernel matrix.

Parameters:

  • k (int, optional) – Number of largest eigenvalues to return. If None, returns all.

  • return_full (bool, default False.) – Only for fft_solver=’rfft2’. If True, returns eigenvalues in full FFT layout (ny, nx) flattened.

Returns:

Eigenvalues. If return_full=True, shape is (ny * nx,).

If k specified, returns top-k in descending order.

Return type:

np.ndarray

Notes

If fft_solver=’rfft2’, spectrum is stored in rfft2 format, not full FFT format. To convert to full FFT format, use return_full=True.

trace() → float[source]#

Compute the trace of the kernel matrix.

Returns:

Trace of K (sum of eigenvalues).

Return type:

float

square_trace() → float[source]#

Compute the trace of the squared kernel matrix.

Returns:

Trace of K² (sum of squared eigenvalues).

Return type:

float

quadsv.fft.spatial_q_test_fft(Xn: ndarray, kernel: FFTKernel, return_pval: bool = True, is_standardized: bool = False) → float | ndarray | Tuple[float, float] | Tuple[ndarray, ndarray][source]#

FFT-accelerated spatial Q-test for grid data.

Tests whether a spatial variable exhibits significant clustering or dispersion using FFT-based spectral decomposition. This function provides fast approximation via Parseval’s theorem compared to dense kernel methods.

Parameters:

  • Xn (np.ndarray) – Input data tensor. Shape (ny, nx) for single feature or (ny, nx, M) for M features. Order follows kernel dimensions. Will be automatically reshaped to 3D if 2D.

  • kernel (FFTKernel) – Pre-constructed FFT kernel object for grid data.

  • return_pval (bool, default True) – If True, returns (Q, pval) tuple; if False, returns Q only.

  • is_standardized (bool, default False) – If True, skips Z-score standardization internally. Otherwise standardizes per-feature (mean 0, std 1) across spatial dimensions.

Returns:

  • Q (float or np.ndarray) – Test statistic. Scalar if input was 2D; array of shape (M,) if 3D.

  • pval (float or np.ndarray, optional) – Tail probability under null hypothesis. Only returned if return_pval=True. Uses Liu’s method for most kernels; Normal approximation for Moran’s I.

Raises:

AssertionError – If Xn spatial dimensions don’t match kernel shape (ny, nx).

Notes

Under H₀: data is spatially independent. Under H₁: mean-shift present.

Computationally: Q = z^T K z where z is standardized data. Uses FFT via Parseval’s theorem to compute $Q = sum_{i,j} lambda_{i,j} Z^2_{i,j}$ in O(N log N) time instead of O(N³) dense methods.

For Moran’s I kernel (which has negative eigenvalues), uses Normal approximation based on asymptotic theory. For other kernels, uses Liu’s chi-squared mixture approximation.

Examples

>>> ny, nx = 32, 32
>>> kernel = FFTKernel((ny, nx), method='gaussian', bandwidth=1.0)
>>> data = np.random.randn(ny, nx)
>>> Q, pval = spatial_q_test_fft(data, kernel)
quadsv.fft.spatial_r_test_fft(Xn: ndarray, Yn: ndarray, kernel: FFTKernel, return_pval: bool = True, is_standardized: bool = False) → float | ndarray | Tuple[float, float] | Tuple[ndarray, ndarray][source]#

FFT-accelerated spatial R-test (bivariate) for grid data.

Tests for spatial co-variation between two variables using the specified kernel. Computes the cross-variance statistic R = x^T K y via FFT-based spectral methods.

Parameters:

  • Xn (np.ndarray) – First input tensor. Shape (ny, nx) for single feature or (ny, nx, M) for M features. Will be automatically reshaped to 3D if 2D.

  • Yn (np.ndarray) – Second input tensor. Must have the same shape as Xn.

  • kernel (FFTKernel) – Pre-constructed FFT kernel object for grid data.

  • return_pval (bool, default True) – If True, returns (R, pval) tuple; if False, returns R only.

  • is_standardized (bool, default False) – If True, skips standardization. Otherwise standardizes each variable independently (mean 0, std 1) across spatial dimensions.

Returns:

  • R (float or np.ndarray) – Test statistic (cross-variance). Scalar if input was 2D; array of shape (M,) if 3D.

  • pval (float or np.ndarray, optional) – Two-tailed p-value under null hypothesis (no spatial co-variation). Based on Normal approximation: $z = R / sqrt{text{Trace}(K^2)}$. Only returned if return_pval=True.

Raises:

AssertionError – If Xn and Yn shapes don’t match, or spatial dimensions don’t match kernel.

Notes

Under H₀: x and y are spatially independent. Under H₁: spatial co-clustering or co-dispersion present.

Computationally: R = z_x^T K z_y where z_x, z_y are standardized data. Uses FFT via Parseval’s theorem: $R = frac{1}{N} sum_{i,j} overline{Z_x}_{i,j} lambda_{i,j} Z_y_{i,j}$.

P-value calculation assumes asymptotic Normality with variance estimated from kernel trace: $text{Var}(R) approx text{Trace}(K^2) / N^2$. Returns two-tailed probability: $p = 2 P(|Z| > |z\text{-score}|)$.

Examples

>>> ny, nx = 32, 32
>>> kernel = FFTKernel((ny, nx), method='gaussian', bandwidth=1.0)
>>> x_data = np.random.randn(ny, nx)
>>> y_data = np.random.randn(ny, nx)
>>> R, pval = spatial_r_test_fft(x_data, y_data, kernel)

Usage example#

import numpy as np
from quadsv.fft import FFTKernel
from quadsv.statistics import spatial_q_test

# Create FFT kernel for 1000x1000 grid
kernel_fft = FFTKernel(
    shape=(1000, 1000),
    method='car',
    rho=0.9,
    topology='square'
)

# Generate synthetic grid data
data_grid = np.random.randn(1000, 1000)

# Compute Q-test efficiently (O(N log N))
Q, pval = spatial_q_test(data_grid, kernel_fft, return_pval=True)
print(f"Q = {Q:.4f}, p-value = {pval:.4e}")

Supported Topologies#

  • square: 4-neighbor connectivity (default)

  • hexagonal: 6-neighbor connectivity

Performance#

For N = 1,000,000 spots (1000×1000 grid): - Dense eigendecomposition: ~10 hours - FFT eigendecomposition: ~1 minute - Q-test: O(N log N) ≈ 20 seconds