splisosm.kernel =============== .. py:module:: splisosm.kernel .. autoapi-nested-parse:: Spatial kernel abstractions. Classes ------- .. autoapisummary:: splisosm.kernel.FFTKernel splisosm.kernel.IdentityKernel splisosm.kernel.SpatialCovKernel Module Contents --------------- .. py:class:: FFTKernel(shape, spacing = (1.0, 1.0), rho = 0.99, neighbor_degree = 1, workers = None, centering = False) Bases: :py:obj:`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 :math:`O(N \log N)` via the 2-D Fast Fourier Transform. :param shape: Grid shape ``(ny, nx)``. :param spacing: Physical spacing ``(dy, dx)`` between neighbouring raster cells. :param rho: Spatial autocorrelation coefficient in CAR kernel. :param neighbor_degree: Neighbour ring degree for graph construction. ``1`` uses nearest neighbours in the periodic metric. :param workers: Number of workers used by ``scipy.fft.fft2``. :param centering: 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; :meth:`xtKx`, :meth:`trace`, :meth:`square_trace`, and :meth:`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. :param Initialize kernel-specific state.: .. py:method:: 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) ) :param x: Input with shape ``(ny, nx)`` or ``(ny, nx, m)``. :param epsilon: Regularisation / noise level. :returns: Residuals of the same shape as ``x``. :rtype: np.ndarray .. py:method:: eigenvalues(k = None) Return kernel eigenvalues in descending order. :param k: Number of leading eigenvalues to return. ``None`` returns all. :returns: Eigenvalues in descending order. :rtype: np.ndarray .. py:method:: power_spectral_density_1d(bins = 50) Compute the 1D power spectral density (radial profile). :param bins: 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. .. py:method:: realization() :abstractmethod: Return the dense kernel matrix (not recommended for large grids). .. py:method:: square_trace() Return ``tr(K^2)``. .. py:method:: trace() Return ``tr(K)``. .. py:method:: xtKx(x) Compute ``x^T K x`` in ``O(N log N)`` via FFT. :param x: Input with shape ``(ny, nx)`` or ``(ny, nx, m)``. :returns: Scalar for 2D input, or shape ``(m,)`` for 3D input. :rtype: float or np.ndarray .. py:method:: xtKx_approx(x, k = None) Same as :meth:`xtKx` (FFT is already exact for periodic grids). .. py:attribute:: n_grid .. py:attribute:: neighbor_degree :value: 1 .. py:attribute:: rho .. py:attribute:: spectrum .. py:attribute:: workers :value: None .. py:class:: IdentityKernel(n_spots, centering = False) Bases: :py:obj:`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. :param n_spots: Number of spatial locations. :type n_spots: int :param centering: 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 :math:`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. :type centering: bool, optional :param Initialize kernel-specific state.: .. py:method:: eigenvalues(k = None) Return the leading eigenvalues of ``K``. For ``K = I`` all :math:`n` eigenvalues are ``1``. For ``H I H = H`` there are :math:`n - 1` unit eigenvalues and a single zero eigenvalue; only the nonzero ones are returned. .. py:method:: realization() Return ``I`` (or ``H = I - (1/n) 1 1^T`` when ``centering=True``). .. py:method:: square_trace() Return tr(K²): ``n`` for ``I`` and ``n - 1`` for ``H`` (idempotent). .. py:method:: trace() Return tr(K): ``n`` for ``I`` and ``n - 1`` for ``H``. .. py:method:: xtKx(x) Compute ``x^T K x`` where ``K = I`` (or ``H I H = H`` when centred). .. py:method:: xtKx_approx(x, k = None) Quadratic form; identical to :meth:`xtKx` (no approximation needed). .. py:method:: xtKx_exact(x) Exact quadratic form; identical to :meth:`xtKx` for this kernel. .. py:class:: SpatialCovKernel(coords = None, *, adj_matrix = None, k_neighbors = 4, rho = 0.99, standardize_cov = True, centering = False) Bases: :py:obj:`Kernel` Graph-based spatial covariance kernel. Construct from spot coordinates (:meth:`from_coordinates` or the ``coords`` positional argument) or from a pre-built k-NN adjacency / weight matrix (:meth:`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 :meth:`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 :meth:`from_coordinates` and :meth:`from_adjacency` for more descriptive call sites. :param coords: Spot coordinates. A k-NN graph is built from these and used to construct the CAR precision matrix. :type coords: array-like of shape ``(n_spots, n_dim)``, optional :param adj_matrix: Pre-built symmetric k-NN adjacency / weight matrix. The CAR precision is built directly from this matrix. :type adj_matrix: sparse or dense matrix of shape ``(n_spots, n_spots)``, optional :param k_neighbors: Number of nearest neighbours when building the graph from ``coords``. Ignored when ``adj_matrix`` is provided. :type k_neighbors: int :param rho: CAR spatial autocorrelation coefficient (0 < rho < 1). :type rho: float :param standardize_cov: Scale the covariance to unit marginal variance. :type standardize_cov: bool :param centering: Apply double-centring H K H. :type centering: bool .. py:method:: eigendecomposition() Compute and cache the full eigendecomposition of the dense kernel. Populates :attr:`K_eigvals`, :attr:`K_eigvecs`, and :attr:`Q`. :raises ValueError: If called on an implicit kernel (no ``K_sp``). .. py:method:: 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 :meth:`xtKx` calls. :param k: Number of leading eigenvalues. ``None`` returns all cached values. :returns: Eigenvalues in descending order. :rtype: torch.Tensor .. py:method:: from_adjacency(adj_matrix, rho = 0.99, standardize_cov = True, centering = False) :classmethod: 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. :param adj_matrix: Symmetric adjacency or weight matrix of shape ``(n_spots, n_spots)``. Can be a dense ``np.ndarray`` or any ``scipy.sparse`` matrix. :param rho: CAR spatial autocorrelation coefficient (0 < rho < 1). :param standardize_cov: Scale covariance to unit marginal variance. :param centering: Apply double-centring H K H. :rtype: SpatialCovKernel .. py:method:: from_coordinates(coords, k_neighbors = 4, rho = 0.99, standardize_cov = True, centering = False) :classmethod: Build a CAR spatial kernel from spot coordinates. :param coords: Spot coordinates of shape ``(n_spots, n_dim)``. :param k_neighbors: Number of nearest neighbours for the spatial graph. :param rho: CAR spatial autocorrelation coefficient (0 < rho < 1). :param standardize_cov: Scale covariance to unit marginal variance. :param centering: Apply double-centring H K H. :rtype: SpatialCovKernel .. py:method:: rank() Return effective rank of the stored representation. :rtype: int .. py:method:: realization() Return the realised dense kernel matrix. For implicit kernels this forces dense inversion — prefer calling :meth:`eigenvalues` first to populate ``Q``, then use ``Q @ Q.T``. :returns: Dense matrix of shape ``(n_spots, n_spots)``. :rtype: torch.Tensor .. py:method:: shape() Return kernel shape ``(n_spots, n_spots)``. :rtype: tuple of int .. py:method:: 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. :rtype: torch.Tensor .. py:method:: 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. :rtype: torch.Tensor .. py:method:: xtKx(x) Compute ``x^T K x`` using the cached low-rank factor when available. When :attr:`Q` has been populated (by a prior :meth:`eigenvalues` call) the result is ``x^T Q Q^T x``, i.e. consistent with the approximation ``K ≈ Q Q^T``: * **Dense mode**: :meth:`eigendecomposition` produces a full-rank Q so ``Q Q^T = K`` exactly. * **Implicit mode**: :meth:`eigenvalues` produces a rank-k Q; the result is therefore a rank-k approximation of ``x^T K x``. Before :meth:`eigenvalues` is called (``Q`` is ``None``) this falls back to :meth:`xtKx_exact`. Use :meth:`xtKx_exact` directly to always obtain the exact quadratic form regardless of whether Q has been computed. :param x: Input matrix of shape ``(n_spots, d)``. :type x: torch.Tensor :returns: Matrix of shape ``(d, d)``. :rtype: torch.Tensor .. seealso:: :py:obj:`xtKx_exact` Always uses the exact path (dense multiply or LU solve). .. py:method:: xtKx_approx(x, k = None) Compute ``x^T K_k x`` via the top-``k`` low-rank factor. Calls :meth:`eigenvalues` to ensure :attr:`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. :param x: Input matrix of shape ``(n_spots, d)``. :type x: torch.Tensor :param k: Number of leading eigenvalues to include. ``None`` uses all available eigenvalues (populates the full Q first). :type k: int or None :returns: Approximated quadratic form of shape ``(d, d)``. :rtype: torch.Tensor .. py:method:: 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 :meth:`xtKx`, this method **ignores** :attr:`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. :param x: Input matrix of shape ``(n_spots, d)``. :type x: torch.Tensor :returns: Matrix of shape ``(d, d)``. :rtype: torch.Tensor .. seealso:: :py:obj:`xtKx` Uses cached Q (low-rank approximation) when available. .. py:attribute:: DENSE_THRESHOLD :type: int :value: 5000 Switch to implicit sparse mode above this number of spots. .. py:attribute:: K_eigvals :type: torch.Tensor | None Cached eigenvalues in descending order. Populated lazily. .. py:attribute:: K_eigvecs :type: torch.Tensor | None Cached eigenvectors (dense mode only). Populated lazily. .. py:attribute:: K_sp :type: torch.Tensor | None Dense centred/standardised kernel ``(n_spots, n_spots)``. ``None`` in implicit mode. .. py:attribute:: Q :type: torch.Tensor | None Low-rank factor ``(n_spots, rank)`` such that K ≈ Q @ Q.T. Populated lazily after the first :meth:`eigenvalues` call. .. py:attribute:: inv_cov :type: scipy.sparse.csc_matrix Sparse precision matrix M = K⁻¹ of shape ``(n_spots, n_spots)``.