quadsv.kernels.nufft

Non-uniform FFT (NUFFT) spectra, kernel and spatial tests for irregular data

When data sit on a regular grid (e.g., a rasterized Visium slide), quadsv.power_spectrum_2d() computes \(|\hat{x}(k)|^2\) with a plain 2D FFT. For data whose spatial coordinates are irregular — e.g., imaging-based in situ platforms, Slide-seq, or a Visium slide read straight from adata.obsm['spatial'] without rasterization — power_spectrum_2d_nufft() evaluates the type-1 NUFFT

\[\hat c(k_y, k_x) = \sum_{j=1}^{n} c_j \, \exp\!\bigl[-i(k_y\,y_j + k_x\,x_j)\bigr]\]

on the same uniform (ny, nx) k-space grid that power_spectrum_2d() would produce for a rasterized input of the same physical extent, and returns \(|\hat c|^2\) in the scipy FFT layout (DC at [0, 0]). Anything downstream — quadsv.comparators.multisample.radial_bin_spectrum(), quadsv.ComparatorIrregular — works identically.

Notation (shared across this module)

Dimensions:

• n: number of spots (on the irregular grid).

• (ny, nx): internal uniform k-grid dimensions; n' = ny · nx.

• (dy, dx): physical spacing per k-grid cell, same unit as the spatial coordinates after multiplying unit_scale.

• unit_scale: multiplier that converts the input coordinates S to the same unit as (dy, dx) (e.g., 0.35 if S are in pixels at 0.35 μm/pixel). Samples from different slides and platforms may ship coordinates in different units; this parameter harmonizes them onto the same physical unit for the internal k-grid and all downstream spectra and tests.

Vectors and matrices:

• S: the n × 2 spatial coordinate matrix of the irregular points, ordered as (y, x).

• K: the n × n translation-invariant kernel at the irregular points.

• K': the n' × n' grid kernel with FFT eigenvalues λ(k) = F(K')(k).

• U: the n × n' type-2 NUFFT evaluation matrix; the band-limited approximation is K ≈ (1/n') · U · diag(λ) · Uᴴ.

• x̂ = Uᴴ x: type-1 NUFFT of a length-N signal onto the k-grid (ny, nx) (vectorized).

Classes

NUFFTKernel

Spatial kernel over irregular 2D coordinates evaluated via NUFFTs.

Functions

power_spectrum_2d_nufft(coords, values, grid_shape, ...)

Compute the 2D power spectrum via type-1 NUFFT

Module Contents

class quadsv.kernels.nufft.NUFFTKernel(coords, grid_shape=None, spacing=None, method='matern', unit_scale=1.0, oversample=2.0, eps=1e-06, workers=None, *, centering=True, **kwargs)

Bases: quadsv.kernels.base.Kernel

Spatial kernel over irregular 2D coordinates evaluated via NUFFTs.

Parallels quadsv.fft.FFTKernel (which requires a regular grid) and implements the Kernel interface so it plugs into quadsv.statistics.spatial_q_test() / quadsv.statistics.spatial_r_test() the same way.

The band-limited approximation of the n × n irregular-point operator is K ≈ (1/n') · U · diag(λ) · Uᴴ, where U is the n × n' type-2 NUFFT matrix and λ = F(K') is the grid kernel’s spectrum. Under this approximation, Parseval’s identity gives the fast quadratic form xᵀ K x = (1/n') Σ_k λ(k) |x̂(k)|² with x̂ = Uᴴ x (a single type-1 NUFFT). The matrix-vector primitive Kx() uses the companion two-shot NUFFT K z = (1/n') · U · (λ ⊙ Uᴴ z) and backs the Hutchinson-based cumulant estimator (quadsv.statistics._hutchinson_cumulants()) and the bipartite R-test cross matrix.

quadsv.spatial_q_test() always uses the k-space Parseval path (xtKx()); quadsv.spatial_r_test() dispatches on shape — paired diagonal for M_x == M_y via xtKy(), full (M_x, M_y) cross matrix via Kx() otherwise. The matmul counterparts (xtKx_matmul() / xtKy_matmul()) are exposed for callers that want to compute the same form directly at the n irregular points — they agree with the spectral path to NUFFT precision (eps) on both regular and irregular coordinates.

Parameters:

• coords (np.ndarray) – Spot coordinates of shape (n, 2) in order (y, x).

• grid_shape (tuple[int, int], optional) – (ny, nx) of the internal uniform k-grid. If None (default), auto-inferred from coords: the grid is sized to cover the bounding box and to resolve the sampling Nyquist set by the median nearest-neighbor distance (fully coordinate-driven, kernel-agnostic). Override only when you know you need a finer or coarser grid.

• spacing (tuple[float, float], optional) – (dy, dx) physical spacing per k-grid cell (same unit as coords after unit_scale). If None (default), auto-inferred alongside grid_shape. When both are supplied, users are responsible for ensuring ny · dy, nx · dx covers the coordinate extent.

• method (str, default 'matern') – Kernel method forwarded to FFTKernel. One of 'gaussian', 'matern', 'moran', 'graph_laplacian', 'car'.

• unit_scale (float, default 1.0) – Multiplier applied to coords so they share the same unit as spacing (e.g. 0.35 if coords are in pixels at 0.35 μm/pixel).

• oversample (float, default 2.0) – Auto-grid oversampling factor above the sampling Nyquist. Used only when grid_shape / spacing are auto-derived. Larger values give a finer k-grid (more accurate, slower); 2.0 is safe for all tested kernels.

• eps (float, default 1e-6) – NUFFT tolerance forwarded to finufft.

• workers (int, optional) – Forwarded to scipy.fft (used by Kx_grid()) and reserved for future finufft parallelism. None uses the SciPy default.

• **kwargs – Method-specific kernel hyperparameters forwarded to the internal FFTKernel. bandwidth / nu for gaussian / matern; rho for car; plus a coord-aware k_neighbors for the graph methods (moran / graph_laplacian / car). Note the k_neighbors-to-neighbor_degree mapping differs from that of FFTKernel due to an oversampled internal grid. neighbor_degree is chosen so the internal-grid-ring cutoff matches the median k-th nearest-neighbor distance among the coords — matching MatrixKernel’s mutual-k-NN graph semantic up to the band-limit of the internal grid. See _resolve_k_neighbors_on_coords() for the mapping detail. Pass neighbor_degree directly to bypass this and use the raw grid-ring semantic as in FFTKernel.

• centering (bool)

Variables:

• coords (np.ndarray) – Original (n, 2) coordinates.

• n (int) – Number of spots n.

• grid_shape (tuple[int, int]) – Internal k-grid shape (ny, nx).

• spacing (tuple[float, float]) – Physical spacing per k-grid cell (dy, dx).

• method (str) – Kernel method name.

• params (dict) – Resolved kernel hyperparameters (snapshot of the internal FFT kernel).

• workers (int or None) – scipy.fft worker count used by Kx_grid().

• stores_precision (bool) – Always False — NUFFTKernel never holds an n × n matrix.

Kx(z)

Matrix–vector product K z at the n irregular coordinates.

Implements the band-limited apply

\[K z \;\approx\; \tfrac{1}{n'} \, U \bigl(\lambda \odot U^{\mathsf H} z\bigr),\]

evaluated as type-1 NUFFT → elementwise multiply by λ(k) / n' → type-2 NUFFT. Output length n, same shape as z. Base primitive for xtKx_matmul(), xtKy_matmul(), the Hutchinson cumulant estimator used by Liu’s null approximation, and the bipartite R-test in DetectorIrregular.

Parameters:

z (np.ndarray) – (n,) or (n, M).

Returns:

Same shape as z.

Return type:

np.ndarray

Kx_grid(x)

Grid-domain companion of Kx() — (ny, nx) spatial output.

Whereas Kx() returns the length-n apply at the irregular coordinates, Kx_grid() returns the apply evaluated on the internal uniform grid. Pipeline: type-1 NUFFT → undo the coordinate-centering phase (needed here because we keep complex coefficients; the square-magnitude and adjoint-round-trip paths of xtKx() / Kx() absorb it automatically) → multiply by λ(k) → ifftshift → ifft2 → real.

Parameters:

x (np.ndarray) – (n,) or (n, M).

Returns:

Real (ny, nx) or (ny, nx, M) in the scipy FFT layout (DC at [0, 0]).

Return type:

np.ndarray

eigenvalues(k=None, return_full_layout=False)

Eigenvalues of the n × n irregular-point operator.

Returns the spectrum of the realized operator K_n (raw) or H K_n H (centered) at the irregular coordinates — not the internal FFT-grid spectrum. Purely matrix-free: no dense n × n construction.

Two paths:

• ``k`` given — Lanczos top-k. Wrap Kx() as a LinearOperator and call eigsh(which='LM')() for the top k eigenvalues by magnitude. Cost: O(k × nufft). Not cached.

• ``k=None`` — Toeplitz-M (analytic). See _eigvals_toeplitz_M(). Applies when Λ is PSD and the support size r = #{k : |λ(k)| > 10⁻⁵ · max|λ|} is below _TOEPLITZ_R_THRESHOLD (default 2000). Exact to NUFFT eps; cost dominated by O(r³) eigvalsh on a dense r × r reduced matrix. Cached per centering mode.

Full-spectrum requests that fall outside Toeplitz-M’s reach (indefinite Λ like Moran, or broad-support kernels like CAR at strong coupling) raise NotImplementedError rather than falling back to an approximate density reconstruction. For Liu’s method, use compute_null_params(..., method='liu', liu_n_probes=...)() to get cumulant-based Liu directly from \(2m\) matvecs.

Parameters:

• k (int, optional) – Top-k eigenvalues by magnitude (Lanczos). None returns the full spectrum via Toeplitz-M.

• return_full_layout (bool, default False) – Kept for API compatibility with FFTKernel.eigenvalues(); ignored.

Returns:

Eigenvalues in descending order. Length k when k is given, n for the full-spectrum path.

Return type:

np.ndarray

square_trace()

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

Closed-form (n²/n'²) · λᵀ Ψ λ with Toeplitz Ψ_{k,k'} = |φ(k'-k)|², evaluated as a linear (non-circular) 2D convolution of |φ|² with λ via a doubled-grid FFT in O(n' log n'). φ(j) = (1/n) Σ_i exp(-ij·y_i) is evaluated on a (2·ny, 2·nx) mode grid by a separate type-1 NUFFT of the ones vector (see _coord_power_spectrum_doubled()), and λ is zero-padded to the same doubled layout so the FFT convolution does not wrap values of |φ|² across the n'-period — a silent bias of up to ~45% on broad-spectrum kernels (CAR, graph_laplacian) on the typical oversampled NUFFT grid. On a regular grid where coords coincide with the k-grid, |φ|² = δ and the formula collapses to (n/n')² · Σ_k λ(k)². Adjusts by -2·s₂/n + s₁²/n² when centering=True.

Observed band-limit residuals vs. explicit Kx(I) truth: ≲ 1e-7 on Gaussian / Matern, ~0.1 % on CAR, ~1 % on graph_laplacian, and ~0.05–1.2 % on Moran (indefinite Λ) — accurate across regular, irregular, and clustered coord layouts.

Return type:

float

trace()

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

Closed-form (n/n') · Σ_k λ(k) — exact because the diagonal of G = UᴴU is n regardless of coord arrangement, so trace(K_n) = (n/n') · trace(K_grid) is independent of the coord layout. Adjusts by -s₁/n when centering=True (s₁ = 𝟏ᵀ K 𝟏 via a single K·𝟏 apply in Kernel._ones_stats()).

Return type:

float

xtKx(x)

Quadratic form xᵀ K x via k-space Parseval.

Implements the default path:

\[x^T K x \;=\; \frac{1}{n'} \sum_k \lambda(k) \, |\hat x(k)|^{2}\]

using one type-1 NUFFT of x and an elementwise Parseval sum. Only the real power spectrum |x̂|² of shape (ny, nx) is materialized — no ifft2, no spatial-grid copy of x.

Parameters:

x (np.ndarray) – (n,) for one feature or (n, M) for M features.

Returns:

Scalar if x is 1-D, shape (M,) otherwise.

Return type:

float or np.ndarray

See also

xtKx_matmul

Compute xᵀ · Kx via the length-n matrix product.

xtKx_matmul(x)

Quadratic form xᵀ K x via direct matmul.

Computes Q_B = xᵀ · self.Kx(x) end-to-end at the n irregular points. Sparse-aware on x (x.multiply(Kx).sum). ~2× the NUFFT work of xtKx() per feature; agrees with it to NUFFT precision on regular grids and to the torus-BC band (~1–2 %) on irregular ones.

Parameters:

x (np.ndarray or scipy.sparse matrix) – (n,) or (n, M).

Returns:

Scalar for 1-D input; (M,) for batched.

Return type:

float or np.ndarray

xtKy(x, y)

Bilinear form xᵀ K y via cross Parseval.

Implements the default path:

\[x^T K y \;=\; \frac{1}{n'} \sum_k \lambda(k) \, \overline{\hat x(k)}\, \hat y(k).\]

Paired same-M convention — returns the diagonal of Xᵀ K Y (shape (M,)) for batched inputs, scalar for 1-D inputs. For the bipartite (M_x, M_y) cross matrix use xtKy_matmul() (or build it explicitly via X.T @ self.Kx(Y)).

Parameters:

• x (np.ndarray) – (n,) or (n, M) — must share shape.

• y (np.ndarray) – (n,) or (n, M) — must share shape.

Returns:

Scalar for 1-D inputs; (M,) for batched.

Return type:

float or np.ndarray

xtKy_matmul(x, y)

Bilinear form xᵀ K y via direct matmul.

Returns the paired (M,) diagonal of Xᵀ K Y (sparse-aware on x). For the full (M_x, M_y) bipartite cross matrix build it explicitly as X.T @ self.Kx(Y) — that’s what DetectorIrregular does for compute_rstat when the two feature blocks have different widths.

Parameters:

• x (np.ndarray or scipy.sparse matrix) – (n,) or (n, M).

• y (np.ndarray or scipy.sparse matrix) – (n,) or (n, M).

Returns:

Scalar for 1-D inputs; (M,) for batched.

Return type:

float or np.ndarray

coords: ndarray

grid_shape: tuple[int, int]

method: str = 'matern'

n: int

params: dict

spacing: tuple[float, float]

stores_precision: bool = False

workers: int | None = None

quadsv.kernels.nufft.power_spectrum_2d_nufft(coords, values, grid_shape, spacing, unit_scale=1.0, eps=1e-06, center_coords=True)

Compute the 2D power spectrum via type-1 NUFFT

This function computes the power spectrum \(P(k) = |\hat{c}(k)|^2\) of one or more non-uniform spatial signals via a type-1 NUFFT. The output has the same (ny, nx) layout as quadsv.kernels.fft.power_spectrum_2d() with fft_solver='fft2': DC at [0, 0], Nyquist at [ny/2, nx/2] (when dimensions are even).

Parameters:

• coords (np.ndarray) – Non-uniform spatial coordinates, shape (n, 2) in the order (y, x). Values outside the physical domain implied by grid_shape and spacing are folded into [-π, π) by finufft.

• values (np.ndarray) – Signal strengths at each coordinate. Shape (n,) for a single feature, or (n, M) for M stacked features (e.g., genes) on the same coordinates. Real-valued; promoted to complex internally.

• grid_shape (tuple[int, int]) – (ny, nx) of the target uniform k-space grid. Match whatever grid you use for rasterized samples so the two paths produce comparable spectra.

• spacing (tuple[float, float]) – (dy, dx) physical spacing per cell of the target grid, in the same unit as unit_scale * coords. Together with grid_shape this defines the physical domain extent (ny · dy, nx · dx).

• unit_scale (float, default 1.0) – Multiplier applied to coords before scaling into [-π, π). Use this to convert per-sample coordinate units into the common unit of spacing (e.g., 0.35 if coords are in Visium full-res pixels at 0.35 μm/pixel and spacing is in μm).

• eps (float, default 1e-6) – NUFFT tolerance forwarded to finufft.

• center_coords (bool, default True) – If True, subtract the mean of coords before scaling — avoids wrapping artefacts when coordinates are stored with an arbitrary origin offset (e.g., Visium pixel coordinates start at a few thousand). Power spectra are translation-invariant so recentering does not change the result.

Returns:

Power spectrum. Shape (ny, nx) for 1D values or (ny, nx, M) for 2D values, with DC at index [0, 0].

Return type:

np.ndarray

Raises:

• ImportError – If finufft is not installed.

• ValueError – If input shapes are inconsistent.

Examples

>>> import numpy as np
>>> coords = np.random.default_rng(0).uniform(0, 100, size=(500, 2))
>>> vals = np.random.default_rng(1).standard_normal(500)
>>> P = power_spectrum_2d_nufft(coords, vals, grid_shape=(32, 32), spacing=(4.0, 4.0))
>>> P.shape
(32, 32)