quadsv.kernels.base =================== .. py:module:: quadsv.kernels.base .. autoapi-nested-parse:: Abstract bases for the spatial-kernel layer. This module hosts the two ABCs: - :class:`Kernel`: the universal interface every kernel implements (``Kx`` / ``xtKx`` / ``xtKy`` / ``trace`` / ``square_trace`` / ``eigenvalues``). Subclass this for a custom kernel that can be expressed through a single ``self._K`` buffer. - :class:`MatrixKernelBase`: matrix-form kernels with dense / sparse / sparse-precision auto-switching. Subclass this for a new matrix backend; :class:`~quadsv.MatrixKernel` is the standard concrete subclass. Concrete classes live in sibling modules: :mod:`quadsv.kernels.matrix`, :mod:`quadsv.kernels.fft`, and :mod:`quadsv.kernels.nufft`. Classes ------- .. autoapisummary:: quadsv.kernels.base.Kernel quadsv.kernels.base.MatrixKernelBase Module Contents --------------- .. py:class:: Kernel(*args, centering = True, **kwargs) Bases: :py:obj:`abc.ABC` Abstract base class shared by all spatial kernels in :mod:`quadsv`. Concrete backends: - :class:`MatrixKernel` — explicit n×n kernel or its sparse precision matrix. - :class:`quadsv.fft.FFTKernel` — grid kernel via its eigenvalue spectrum. - :class:`quadsv.nufft.NUFFTKernel` — irregular-point kernel evaluated through a type-1 / type-2 NUFFT round-trip. Required interface ------------------ Every concrete kernel exposes: - ``n``, ``method``, ``params``, ``centering`` — integer / string / dict / bool attributes. - :meth:`xtKx`, :meth:`xtKy`, :meth:`Kx` — quadratic / bilinear / apply primitives. .. note:: The empirical data centering (z-scoring) inside :func:`quadsv.spatial_q_test` / :func:`~quadsv.spatial_r_test` breaks independence across spatial obervations. As a result, the null distribution of the test statistic ``Q = Zᵀ K Z = Xᵀ (H K H) X / σ²`` with ``H = I - 𝟏𝟏ᵀ/n`` should inspect the spectrum of a centered kernel ``HKH``. Every :class:`Kernel` carries a ``centering`` flag (default ``True``). Set ``centering=False`` to recover the raw ``K`` moments (useful for diagnostics or theoretical comparison). .. py:method:: Kx(x) :abstractmethod: ``K x`` (raw) or ``HKH x`` (centered). .. py:method:: eigenvalues(k = None) :abstractmethod: Eigenvalues of ``K`` (raw) or ``HKH`` (centered). .. py:method:: square_trace() :abstractmethod: ``trace(K²)`` or ``trace((HKH)²) = trace(K²) − 2·s₂/n + s₁²/n²``. Same contract as :meth:`trace` re: backend-specific kwargs. .. py:method:: trace() :abstractmethod: ``trace(K)`` or ``trace(HKH) = trace(K) − s₁/n``. Concrete backends are free to add kwargs for backend-specific options (e.g. :class:`MatrixKernelBase` exposes ``n_probes`` for its precision-stored Hutchinson path). .. py:method:: xtKx(x) :abstractmethod: ``xᵀ K x`` (raw) or ``xᵀ HKH x`` (centered). .. py:method:: xtKy(x, y) :abstractmethod: ``xᵀ K y`` (raw) or ``xᵀ HKH y`` (centered). .. py:attribute:: centering :type: bool .. py:attribute:: method :type: str .. py:attribute:: n :type: int .. py:attribute:: params :type: dict .. py:class:: MatrixKernelBase(n, method = 'gaussian', *, centering = True, **kwargs) Bases: :py:obj:`Kernel` Concrete base for kernels backed by an explicit or implicit ``n × n`` matrix. Subclass this when you want a custom way to construct ``K`` (e.g. a bespoke function of coordinates, a learnt kernel, or an ad-hoc adjacency matrix) while inheriting every downstream primitive ready-to-use. Subclasses must implement :meth:`_build_kernel` to return either the ``(n, n)`` kernel matrix (dense ``np.ndarray`` or ``scipy.sparse``) or its precision matrix ``M = K^{-1}``. If a precision matrix is returned, set ``self.stores_precision = True`` before calling ``super().__init__``, or flip it inside ``_build_kernel`` while the buffer is being built; the :class:`MatrixKernel` CAR / precomputed-inverse path shows a worked example. Everything else — cached sparse LU for implicit solves, Hutchinson trace estimation when ``stores_precision=True``, automatic sparse-vs-dense dispatch in ``xtKx`` / ``Kx``, standardized-quadratic-form cache (``_K_column_sums`` / :meth:`xtKx_standardized`), pickle safety for the non-picklable LU factor — is already implemented here. Handles dense, sparse, and implicit (operator-based) kernels. Switches between an explicit representation (the kernel matrix ``K`` is stored and used directly) and an implicit representation (the precision matrix ``M = K^{-1}`` is stored and linear systems are solved on demand) based on problem size. :ivar n: Number of observations (``n``). :vartype n: int :ivar method: Kernel method name (free-form; used only for diagnostics / provenance). :vartype method: str :ivar params: Resolved kernel parameters — subclasses decide what goes here. :vartype params: dict :ivar stores_precision: If ``True``, ``_K`` holds the precision matrix ``M = K^{-1}`` and linear solves (via a cached LU) are used for :meth:`xtKx` and trace estimation. If ``False``, ``_K`` is the realized kernel matrix (dense or sparse). :vartype stores_precision: bool .. rubric:: Notes The internal buffer ``_K`` stores the kernel matrix when ``stores_precision=False`` and the precision matrix ``K^{-1}`` when ``stores_precision=True``. Public methods (:meth:`xtKx`, :meth:`trace`, :meth:`square_trace`, :meth:`eigenvalues`) transparently handle both cases; callers should not access ``_K`` directly. .. py:method:: Kx(x) Apply the kernel operator to ``x``, returning ``K @ x``. Single public primitive for kernel–vector products. Handles explicit (dense or sparse ``K``) and implicit (precision matrix + cached LU) cases uniformly. :param x: ``(n,)`` or ``(n, M)``. Sparse inputs are densified internally because ``scipy.linalg.lu_solve`` / ``splu.solve`` require dense RHS and ``K @ x`` typically returns dense anyway. :type x: np.ndarray or scipy.sparse matrix :returns: ``(n,)`` if ``x`` was 1D, else ``(n, M)``. :rtype: np.ndarray .. rubric:: Examples >>> import numpy as np >>> from quadsv import MatrixKernel >>> rng = np.random.default_rng(0) >>> coords = rng.standard_normal((40, 2)) >>> kernel = MatrixKernel.from_coordinates(coords, method="matern") >>> kernel.Kx(rng.standard_normal(40)).shape (40,) .. py:method:: eigenvalues(k = None) Eigenvalues of ``K`` or ``HKH`` (according to ``self.centering``). Three paths: - **Dense ``K``** — form ``HKH = K − 𝟏 r^T − c 𝟏^T + m · 𝟏𝟏^T`` in closed form (``r`` / ``c`` = row / col means, ``m`` = grand mean) and call :func:`numpy.linalg.eigvalsh`. Same O(n³) cost as the raw case, no extra memory. - **Sparse explicit ``K``** — wrap ``H K H · v`` as a :class:`scipy.sparse.linalg.LinearOperator` and call :func:`eigsh`. Preserves ``K``'s sparsity — we never densify. - **Implicit precision (sparse ``M = K⁻¹``)** — wrap ``H K H · v = H (M⁻¹ (H v))`` where ``M⁻¹·`` is the cached sparse-LU solve. Same sparsity benefit as the raw-inverse path. Results are cached per centering mode; switching ``self.centering`` invalidates and recomputes. :param k: Number of largest-magnitude eigenvalues to return. If None, returns all (dense) or ``max(6, n − 2)`` (sparse/implicit — limited by what ``eigsh`` can extract). :type k: int, optional :returns: Eigenvalues sorted in descending order. :rtype: np.ndarray .. py:method:: realization() Return the realized ``(n, n)`` kernel matrix. When ``self.centering`` is False, returns ``K`` as stored (dense ndarray or sparse matrix). When ``self.centering`` is True, returns the centered operator ``HKH`` with ``H = I − 𝟏𝟏ᵀ/n`` — the operator that :meth:`xtKx` / :meth:`Kx` / :meth:`trace` / :meth:`eigenvalues` all actually apply. ``HKH`` is dense even when ``K`` is sparse (each row gets a column-mean subtracted), so the centered path always materializes a dense result. .. rubric:: Notes If ``stores_precision`` is True, this forces expensive dense inversion of the precision matrix. Prefer :meth:`xtKx` and :meth:`trace` for implicit kernels. .. py:method:: square_trace(n_probes = None) ``trace(K²)`` (raw) or ``trace((HKH)²)`` (centered). Dispatch is automatic based on :attr:`stores_precision`: - **Explicit ``K``**: Frobenius norm of the stored ``K`` — ``Σᵢⱼ Kᵢⱼ²`` in ``O(nnz)`` for sparse, ``O(n²)`` for dense. - **Precision-stored**: ``(1/m) Σᵢ ‖K vᵢ‖²`` with the same cached ``±1`` probes as :meth:`trace` (single LU solve shared across both moments). ``-2·s₂/n + s₁²/n²`` is applied on top when ``centering=True``. :param n_probes: Probe count for the Hutchinson path (default 15). **Ignored** on the analytic / explicit-``K`` path. Shared cache with :meth:`trace`. :type n_probes: int, optional :rtype: float .. py:method:: trace(n_probes = None) ``trace(K)`` (raw) or ``trace(HKH)`` (centered). Dispatch is automatic based on :attr:`stores_precision`: - **Explicit ``K``** (dense or sparse): diagonal sum of the stored ``K`` — exact, ``O(n)`` (sparse) or ``O(n)`` (dense). - **Precision-stored** (CAR in the implicit regime, ``K`` only available through the LU solve on ``K⁻¹``): Hutchinson estimator ``(1/m) Σᵢ vᵢᵀ (K vᵢ)`` with cached ``±1`` probes solved through the precision. ``-s₁/n`` is applied on top when ``centering=True``. :param n_probes: Probe count for the Hutchinson path (default 15). **Ignored** on the analytic / explicit-``K`` path (no probes involved). On the Hutchinson path the cache is keyed on ``n_probes``; requesting a different count recomputes the LU-solve cache and emits a :class:`RuntimeWarning`. :type n_probes: int, optional :rtype: float .. py:method:: xtKx(x) Quadratic form ``x^T K x`` (paired diagonal for batched inputs). :param x: ``(n,)`` or ``(n, M)``. Sparse ``x`` is preserved through the final ``x^T (K x)`` contraction — only the right side ``K @ x`` needs a dense RHS for the solver / BLAS call. :type x: np.ndarray or scipy.sparse matrix :returns: Scalar if 1D input, ``(M,)`` if batched. :rtype: float or np.ndarray .. py:method:: xtKx_standardized(x, means, stds) Compute ``z^T K z`` where ``z = (x - means) / stds`` *without* densifying sparse ``x``. Sparse-aware expansion using the raw-``K`` primitives:: z^T K z = (x^T K x - 2·μ·(K·𝟏)^T·x + μ²·(𝟏^T K 𝟏)) / σ² and the cached row sums of ``K``. This is the fast path for standardizing large sparse feature matrices (e.g. scRNA-seq counts) before a Q-test. :param x: ``(n,)`` or ``(n, M)``. Columns correspond to features. :type x: np.ndarray or scipy.sparse matrix :param means: ``(M,)`` per-feature mean and std (``ddof=1`` to match :func:`quadsv.statistics.spatial_q_test`). :type means: np.ndarray :param stds: ``(M,)`` per-feature mean and std (``ddof=1`` to match :func:`quadsv.statistics.spatial_q_test`). :type stds: np.ndarray :returns: ``(M,)`` standardized quadratic form values. Columns with ``std <= 0`` are returned as zero. :rtype: np.ndarray .. py:method:: xtKy(x, y) Bilinear form ``x^T K y`` (paired diagonal for batched inputs). For ``(n, M)`` batches returns ``(M,)`` — the diagonal of ``X^T K Y`` in the same column order, matching :func:`quadsv.spatial_r_test`. Sparse ``x`` is preserved; only ``K @ y`` is densified. :param x: ``(n,)`` or ``(n, M)`` (must share the M). :type x: np.ndarray or scipy.sparse matrix :param y: ``(n,)`` or ``(n, M)`` (must share the M). :type y: np.ndarray or scipy.sparse matrix :returns: Scalar if 1D inputs; ``(M,)`` if batched. :rtype: float or np.ndarray .. rubric:: Examples >>> import numpy as np >>> from quadsv import MatrixKernel >>> rng = np.random.default_rng(0) >>> coords = rng.standard_normal((40, 2)) >>> kernel = MatrixKernel.from_coordinates(coords, method="matern") >>> x = rng.standard_normal(40) >>> y = rng.standard_normal(40) >>> isinstance(kernel.xtKy(x, y), float) True .. py:attribute:: method :type: str :value: 'gaussian' Kernel method name. .. py:attribute:: n :type: int Number of observations (samples). .. py:attribute:: params :type: dict Resolved kernel parameters after defaults are merged with user overrides. .. py:attribute:: stores_precision :type: bool :value: False Whether the kernel is stored in precision form (``True``) or as the realized kernel matrix (``False``).