FlashDiffusion：融合分块核的谱分解

# FlashDiffusion: Fused Tiled Kernel Spectral Decomposition
Source: [https://arxiv.org/html/2609.38198](https://arxiv.org/html/2609.38198)
Julio Candanedo

###### Abstract

Diffusion maps, and kernel methods more generally, provide an interpretable nonlinear spectral representation basis for geometric learning\. In the geometric limit, small bandwidth, these matrices tend to be high rank and thus require materializing dense Gaussian kernels requires𝒪⁡\(N2\)\\mathcal\{O\}\(N^\{2\}\)memory\. We introduceFlashDiffusion, a matrix\-free method that evaluates dense Gaussian kernel blocks in fused GPU tiles and couples the eigensolver to an empiricalβ\\beta\-flow that selects the finite\-sample resolution scale\. A continuation over sample size and bandwidth warm\-starts increasingly expensive spectral solves from coarser resolutions\.

## 1Introduction

Kernel nonparametric methods and their associated manifold\-learning methods, e\.g\. Laplacian eigenmaps\([Belkin and Niyogi, 2003](https://arxiv.org/html/2609.38198#bib.bib3)\), offer an interpretable route to nonlinear representation learning: nonlinear correlations𝖪=exp⊙⁡\(−β​ℋ\)\\mathsf\{K\}=\\exp\_\{\\odot\}\(\-\\beta\\mathcal\{H\}\)are obtained by exponentiating a matrix of linear correlationsℋ∈ℝN×N\\mathcal\{H\}\\in\\mathbb\{R\}^\{N\\times N\}, e\.g\. pairwise squared distances\. They face two related difficulties\. First, the pairwise kernel requires𝒪⁡\(N2\)\\mathcal\{O\}\(N^\{2\}\)arithmetic and, if materialized,𝒪⁡\(N2\)\\mathcal\{O\}\(N^\{2\}\)memory\. Second, the bandwidthβ\\betafundamentally changes the spectrum and effective rank of the kernel\. Prior work has addressed the computational problem withkkNN graphs,[Muja and Lowe \(2014\)](https://arxiv.org/html/2609.38198#bib.bib42), low\-rank approximations and fast transforms\.111Random Fourier Features\([Rahimi and Recht, 2007](https://arxiv.org/html/2609.38198#bib.bib28);[Rahimi and Recht, 2008](https://arxiv.org/html/2609.38198#bib.bib29)\), Nyström features\([Williams and Seeger, 2000](https://arxiv.org/html/2609.38198#bib.bib22);[Long and Ferguson, 2019](https://arxiv.org/html/2609.38198#bib.bib23)\), and the fast Gauss transform\([Greengard and Strain, 1991](https://arxiv.org/html/2609.38198#bib.bib18);[Raykar and Duraiswami, 2005](https://arxiv.org/html/2609.38198#bib.bib15);[Hertrich, 2024](https://arxiv.org/html/2609.38198#bib.bib16)\)\.We discusskkNN separately in §[A](https://arxiv.org/html/2609.38198#A1), since it sparsifies spatial interactions and hence fundamentally different than truncating spectral rank\. Low\-rank methods face a distinct limitation; at the bandwidths required to resolve a substantial manifold spectrum, the Gaussian kernel is generally high rank, so aggressive compression removes spectral directions the diffusion representation is intended to recover\. Bandwidth selection has meanwhile been approached through median or quantile rules, locally adaptivekk\-th\-neighbor bandwidths\([Zelnik\-Manor and Perona, 2004](https://arxiv.org/html/2609.38198#bib.bib25)\), and the log–log kernel\-sum criterion\([Coifman et al\., 2008](https://arxiv.org/html/2609.38198#bib.bib24)\); these provide useful scale estimates, but depend on the dataset’s sampling and geometry and do not directly optimize finite\-sample spectral resolution\.

We propose that both problems are coupled, and should be treated together\. With finite data, the diffusion matrix approximates an underlying continuum operator\. Singer and Hein\([Singer, 2006](https://arxiv.org/html/2609.38198#bib.bib4);[Hein et al\., 2007](https://arxiv.org/html/2609.38198#bib.bib5)\)showed that the kernel length scale required for consistency depends on sampling density and hence onNN\. In our parameterization, increasingβ\\betadecreases the spatial kernel length scale, so greater sampling can support a largerβ\\betaand resolve finer structure\. We therefore seek an empirical resolution scaleβ⋆​\(N\)\\beta^\{\\star\}\(N\), selected before the leading spectrum reaches the sampling and numerical resolution floors\. Two limits clarify the picture: asβ→0\\beta\\to 0, the kernel approaches its broad\-kernel, low\-rank linear limit\([Candanedo, 2025b](https://arxiv.org/html/2609.38198#bib.bib6)\); along a continuum\-consistent pathN→∞N\\to\\inftywithβ=β⋆​\(N\)\\beta=\\beta^\{\\star\}\(N\), progressively finer manifold structure becomes resolvable and the kernel develops a correspondingly large effective rank\([Clemente and D’Elia, 2018](https://arxiv.org/html/2609.38198#bib.bib26)\)\.

Recent work connects self\-attention with diffusion maps\([Candanedo, 2026](https://arxiv.org/html/2609.38198#bib.bib10)\), suggesting a computational strategy for the quadratic pairwise operator\. FlashAttention\([Dao et al\., 2022](https://arxiv.org/html/2609.38198#bib.bib8);[Zadouri et al\., 2026](https://arxiv.org/html/2609.38198#bib.bib21)\)addresses the analogous𝒪⁡\(N2\)\\mathcal\{O\}\(N^\{2\}\)memory obstruction in attention by evaluating pairwise interactions in GPU tiles rather than materializing the full matrix\.FlashDiffusionapplies the same IO\-aware principle to the Gaussian kernel: pairwise diffusion interactions are evaluated in parallel GPU tiles and used immediately as matrix\-block products inside a block Lanczos eigensolver\([Lanczos, 1950](https://arxiv.org/html/2609.38198#bib.bib20);[Baglama et al\., 2003](https://arxiv.org/html/2609.38198#bib.bib19)\), avoiding explicit materialization of the𝒪⁡\(N2\)\\mathcal\{O\}\(N^\{2\}\)kernel while retaining the dense interaction at the chosen precision\.

FlashDiffusionremains formally𝒪⁡\(N2\)\\mathcal\{O\}\(N^\{2\}\)in pairwise arithmetic, and each Lanczos iteration requires another quadratic kernel application\. At smallβ\\betathe broad Gaussian kernel has rapidly decaying spectrum, so the leading eigenspace is inexpensive to resolve\. Increasingβ\\betareveals additional spectral directions and makes the eigensolve progressively more demanding\. We therefore bootstrap across sample size and bandwidth: starting from a subsamplen≪Nn\\ll N,FlashDiffusionfindsβ⋆​\(n\)\\beta^\{\\star\}\(n\)and its eigenspace at cost𝒪⁡\(n2\)\\mathcal\{O\}\(n^\{2\}\), extends the eigenspace to a larger sample by the Nyström extension\([Bengio et al\., 2003](https://arxiv.org/html/2609.38198#bib.bib27)\), and re\-optimizesβ\\beta\. Repeating this constructs an\(n,β\)\(n,\\beta\)ladder following the empirical resolution curveβ⋆​\(n\)\\beta^\{\\star\}\(n\)toward the full dataset\. For nested subsamples drawn from the same underlying distribution, increased sampling density progressively resolves finer spectral structure, allowing the full\-NNsolve to begin from a well\-adapted eigenspace\.FlashDiffusiontherefore couples bandwidth selection, dense kernel evaluation, and eigenspace continuation while preserving the high\-rank diffusion operator in the regime where its resolved manifold structure is largest\.

## 2FlashDiffusionAlgorithm

### 2\.1Matrix\-free diffusion operator

GivenNNobservationsℛ∈ℝN×D\\mathcal\{R\}\\in\\mathbb\{R\}^\{N\\times D\}, diffusion maps begin with the Gaussian affinity operator:

𝖪=exp⊙⁡\(−β​ℋ\),\\displaystyle\\mathsf\{K\}=\\exp\_\{\\odot\}\\left\(\-\\beta\\,\\mathcal\{H\}\\right\),\(1\)whereℋ\\mathcal\{H\}contains the pairwise squared Euclidean distances\. Letq=𝖪​𝟏q=\\mathsf\{K\}\\mathbf\{1\}\. The Coifman–Lafon normalization\([Coifman and Lafon, 2006](https://arxiv.org/html/2609.38198#bib.bib1)\)defines𝖪α=diag​\(q\)−α​𝖪​diag​\(q\)−α\\mathsf\{K\}\_\{\\alpha\}=\\mathrm\{diag\}\(q\)^\{\-\\alpha\}\\,\\mathsf\{K\}\\,\\mathrm\{diag\}\(q\)^\{\-\\alpha\}, with degreed=𝖪α​𝟏d=\\mathsf\{K\}\_\{\\alpha\}\\mathbf\{1\}and row\-stochastic diffusion operator𝖯=diag​\(d\)−1​𝖪α\\mathsf\{P\}=\\mathrm\{diag\}\(d\)^\{\-1\}\\mathsf\{K\}\_\{\\alpha\}\. Although𝖯\\mathsf\{P\}is generally asymmetric, it is similar to the symmetric operator \(withs=q−α⊙d−1/2s=q^\{\-\\alpha\}\\odot d^\{\-1/2\}, powers taken element\-wise\):

𝖬\\displaystyle\\mathsf\{M\}=diag\(d\)1/2𝖯diag\(d\)−1/2=diag\(s\)𝖪diag\(s\)\.\\displaystyle=\\mathrm\{diag\}\(d\)^\{1/2\}\\,\\mathsf\{P\}\\,\\mathrm\{diag\}\(d\)^\{\-1/2\}=\\mathrm\{diag\}\(s\)\\,\\mathsf\{K\}\\,\\mathrm\{diag\}\(s\)\.\(2\)Thus𝖯\\mathsf\{P\}and𝖬\\mathsf\{M\}share the same eigenvalues, and the diffusion spectrum can be obtained with a symmetric block Lanczos eigensolver\. Note that each freshβ\\betapays an entry fee before any Lanczos iteration can run: the normalization vectors are themselves kernel applications,q=𝖪​𝟏q=\\mathsf\{K\}\\mathbf\{1\}and, forα≠0\\alpha\\neq 0,d=diag​\(q\)−α​𝖪​diag​\(q\)−αd=\\mathrm\{diag\}\(q\)^\{\-\\alpha\}\\,\\mathsf\{K\}\\,\\mathrm\{diag\}\(q\)^\{\-\\alpha\}, two full𝒪⁡\(N2\)\\mathcal\{O\}\(N^\{2\}\)passes that must complete beforessexists\. This is a structural departure from attention, which normalizes each softmax row on\-the\-fly within a single pass,[Milakov and Gimelshein \(2018\)](https://arxiv.org/html/2609.38198#bib.bib44)\. The resulting operator is symmetrizable, eq\.[2](https://arxiv.org/html/2609.38198#S2.E2), unlike self\-attention and this purchase is what licenses the symmetric eigensolver and its guarantees\. Oncessis computed, each block Lanczos iteration reduces to the matrix\-block operation, for a blockV∈ℝN×bV\\in\\mathbb\{R\}^\{N\\times b\}:

𝖬​V=diag⁡\(s\)​𝖪​\(diag⁡\(s\)​V\)\.\\displaystyle\\mathsf\{M\}V=\\mathrm\{diag\}\(s\)\\,\\mathsf\{K\}\\left\(\\mathrm\{diag\}\(s\)V\\right\)\.\(3\)Normalization passes and block Lanczos iterations therefore share a single computational primitive: the Gaussian applicationW↦𝖪​WW\\mapsto\\mathsf\{K\}W, withW=𝟏W=\\mathbf\{1\},q−αq^\{\-\\alpha\}, ordiag⁡\(s\)​V\\mathrm\{diag\}\(s\)V\. This primitive can be evaluated directly from the observations\. With mean\-centered coordinatesℛ¯=ℛ−𝟏​μ⊤\\bar\{\\mathcal\{R\}\}=\\mathcal\{R\}\-\\mathbf\{1\}\\mu^\{\\top\}, define𝒢=ℛ¯​ℛ¯⊤\\mathcal\{G\}=\\bar\{\\mathcal\{R\}\}\\bar\{\\mathcal\{R\}\}^\{\\top\}andr=diag⁡\(𝒢\)r=\\mathrm\{diag\}\(\\mathcal\{G\}\)\. Translation preserves pairwise distances \(and centering boundsrr, taming cancellation in the distance form\), hence:

𝖪=exp⊙⁡\(−β⁡\(r​𝟏⊤\+𝟏​r⊤−2​𝒢\)\)\.\\displaystyle\\mathsf\{K\}=\\exp\_\{\\odot\}\\left\(\-\\beta\\left\(r\\mathbf\{1\}^\{\\top\}\+\\mathbf\{1\}r^\{\\top\}\-2\\mathcal\{G\}\\right\)\\right\)\.\(4\)The kernel is thus a function of Gram blocks of the observations, and never needs to exist as a stored matrix: its tiled evaluation is the subject of §[2\.2](https://arxiv.org/html/2609.38198#S2.SS2)\.

### 2\.2XGEMM: Tiled Kernel Operator

FlashAttention reduces memory traffic by constructing pairwise interactions in fast on\-chip memory rather than materializing the full quadratic matrix in HBM\([Dao et al\., 2022](https://arxiv.org/html/2609.38198#bib.bib8)\);FlashDiffusionapplies the same IO\-aware principle to the Gaussian kernel application in eq\.[3](https://arxiv.org/html/2609.38198#S2.E3)\. Unlike the attention operator, the diffusion operator is scaled on both sides,diag⁡\(s\)​𝖪​diag​\(s\)\\mathrm\{diag\}\(s\)\\,\\mathsf\{K\}\\,\\mathrm\{diag\}\(s\), andsscouples the Lanczos block to global functions of the kernel, and hence no online construction exists\. The costs are nonetheless contained, theqqandddpasses are themselves kernel applications \(§[2\.1](https://arxiv.org/html/2609.38198#S2.SS1)\) and reuse the identical tile schedule withVVreplaced by𝟏\\mathbf\{1\}andq−αq^\{\-\\alpha\}; and within an operator pass the right scaling is the precomputed blockdiag⁡\(s\)​V\\mathrm\{diag\}\(s\)Vwhile the left scaling is a row rescale fused into the accumulation\.FlashDiffusionthus inherits FlashAttention’s tiling at the expense of one or two normalization passes per newβ\\beta\.

For row and column tilesIIandJJ, the observationsℛ¯I\\bar\{\\mathcal\{R\}\}\_\{I\}andℛ¯J\\bar\{\\mathcal\{R\}\}\_\{J\}, their normsrI,rJr\_\{I\},r\_\{J\}, and the Lanczos blockVJV\_\{J\}are loaded from HBM into on\-chip memory\. The tile is evaluated as, with𝒢I​J=ℛ¯I​ℛ¯J⊤\\mathcal\{G\}\_\{IJ\}=\\bar\{\\mathcal\{R\}\}\_\{I\}\\bar\{\\mathcal\{R\}\}\_\{J\}^\{\\top\},

𝖪I​J\\displaystyle\\mathsf\{K\}\_\{IJ\}=exp⊙⁡\(−β⁡\(rI​𝟏⊤\+𝟏​rJ⊤−2​𝒢I​J\)\),\\displaystyle=\\exp\_\{\\odot\}\\left\(\-\\beta\\left\(r\_\{I\}\\mathbf\{1\}^\{\\top\}\+\\mathbf\{1\}r\_\{J\}^\{\\top\}\-2\\mathcal\{G\}\_\{IJ\}\\right\)\\right\),\(5\)followed immediately by:

YI←YI\+𝖪I​J​VJ\.\\displaystyle Y\_\{I\}\\leftarrow Y\_\{I\}\+\\mathsf\{K\}\_\{IJ\}V\_\{J\}\.\(6\)The intermediate matrices𝒢I​J\\mathcal\{G\}\_\{IJ\}and𝖪I​J\\mathsf\{K\}\_\{IJ\}exist only within the working tile and are discarded after accumulation\. Iterating overJJtherefore evaluates the exact dense kernel application without ever storing the𝒪⁡\(N2\)\\mathcal\{O\}\(N^\{2\}\)affinity matrix\. This computation maps directly onto GPU streaming multiprocessors: independent output tiles are distributed across thread blocks, the local Gram product runs on the tensor cores, and the distance transform, exponential, and accumulation are fused into the same GPU kernel as its epilogue\. The implementation targets NVIDIA Ampere and Blackwell architectures\([NVIDIA, 2020](https://arxiv.org/html/2609.38198#bib.bib12);[NVIDIA, 2024](https://arxiv.org/html/2609.38198#bib.bib11)\)usingCUTLASS/CuTeprimitives\([Thakkar and others, 2023](https://arxiv.org/html/2609.38198#bib.bib17)\)\. YieldingXGEMM, an eXponential\-GEMM222General Matrix Multiply\.primitive\.

### 2\.3β\\beta\-function flow

The matrix\-free operator provides the spectrum at a fixed bandwidthβ\\beta;FlashDiffusionadditionally determines the bandwidth at which the requested spectral basis is best resolved\. The algorithm treatslog⁡β\\log\\betaas a one\-dimensional continuation coordinate and advances through a sequence of spectral probes\. Each probe constructs the normalized diffusion operator at the currentβ\\beta, solves for the Perron mode together with the requested nontrivial eigen\-pairs using block implicitly restarted Lanczos, and returns the quantities required to determine the next step and evaluate the bandwidth\-selection criterion\([Sorensen, 1992](https://arxiv.org/html/2609.38198#bib.bib30);[Baglama et al\., 2003](https://arxiv.org/html/2609.38198#bib.bib19)\)\. For centered observations, the initial bandwidth is chosen directly from the data\. Since𝔼i​j​∥xi−xj∥2=2​𝔼i​∥xi∥2\\mathbb\{E\}\_\{ij\}\\lVert x\_\{i\}\-x\_\{j\}\\rVert^\{2\}=2\\mathbb\{E\}\_\{i\}\\lVert x\_\{i\}\\rVert^\{2\},FlashDiffusioninitializes:

β0=min⁡\(0\.1,12​N−1​∑i=1N∥ℛ¯i∥2\)\.\\displaystyle\\beta\_\{0\}=\\min\\left\(0\.1,\\frac\{1\}\{2N^\{\-1\}\\sum\_\{i=1\}^\{N\}\\lVert\\bar\{\\mathcal\{R\}\}\_\{i\}\\rVert^\{2\}\}\\right\)\.\(7\)Hereℛ¯\\bar\{\\mathcal\{R\}\}is the centered data matrix\. The hardcoded cap at0\.10\.1places the initial probe conservatively on the broad\-kernel side when the variance\-based estimate would otherwise begin at a narrow kernel\. The initial scale is therefore determined in𝒪⁡\(N​D\)\\mathcal\{O\}\(ND\)work before any quadratic kernel application is required\.

The potential\-based continuation advances on a dyadic grid,β↦2​β\\beta\\mapsto 2\\beta\. The dyadic structure is required by the effective\-sample estimator below, which compares the off\-diagonal kernel masses at neighboring probesβ\\betaand2​β2\\beta\. Elementwise, the Gaussian kernel satisfies𝖪⁡\(2​β\)=𝖪​\(β\)⊙2\\mathsf\{K\}\(2\\beta\)=\\mathsf\{K\}\(\\beta\)^\{\\odot 2\}, so the two probes correspond exactly to adjacent scales of the same kernel family\. The quantitySoff​\(2​β\)S\_\{\\mathrm\{off\}\}\(2\\beta\)is accumulated during the next probe’s own kernel evaluation, so evaluating the effective\-sample statistic requires no kernel pass beyond those already present in the natural progression of the sweep\. Neighboring probes also have slowly varying eigen\-spaces, so the eigenvectors from the previous probe provide the warm start for the next using the Nyström extension[Nyström \(1930\)](https://arxiv.org/html/2609.38198#bib.bib45)\. The solver follows a gradually rotating spectral subspace rather than restarting from an uninformed initialization\. A non\-dyadic continuation with adaptive ratiosrfast=1\.5r\_\{\\mathrm\{fast\}\}=1\.5,rslow=1\.2r\_\{\\mathrm\{slow\}\}=1\.2, andrcrawl=1\.05r\_\{\\mathrm\{crawl\}\}=1\.05, selected from Lanczos cycle count and Perron\-gap proximity, is retained for API compatibility but does not use the empirical potential selector described below\. Bandwidth selection uses quantities already produced by normalization and the eigensolver\. The first term measures the spectral heat action of the retained nontrivial modes,

𝒱bias​\(β\)=medianj⁡\(−log⁡\(λjλ0\)\)\.\\displaystyle\\mathcal\{V\}\_\{\\mathrm\{bias\}\}\(\\beta\)=\\operatorname\{median\}\_\{j\}\\left\(\-\\log\\left\(\\frac\{\\lambda\_\{j\}\}\{\\lambda\_\{0\}\}\\right\)\\right\)\.\(8\)Asβ→0\\beta\\rightarrow 0, the Gaussian affinity approaches the rank\-one matrix𝟏𝟏⊤\\mathbf\{1\}\\mathbf\{1\}^\{\\top\}and the nontrivial diffusion spectrum collapses, hence𝒱bias→∞\\mathcal\{V\}\_\{\\mathrm\{bias\}\}\\to\\infty\. In the local heat\-kernel regime,𝒱bias​\(β\)∼β−1\\mathcal\{V\}\_\{\\mathrm\{bias\}\}\(\\beta\)\\sim\\beta^\{\-1\}, recording the departure of the requested eigenspace from the broad\-kernel rank\-one limit\. Finite sampling contributes a second scale through an effective sample size measured along the sweep\. LetSoff​\(β\)=∑i≠je−β​∥xi−xj∥2S\_\{\\mathrm\{off\}\}\(\\beta\)=\\sum\_\{i\\neq j\}e^\{\-\\beta\\lVert x\_\{i\}\-x\_\{j\}\\rVert^\{2\}\}\. For neighboring dyadic probes,FlashDiffusiondefines:

neff​\(β\)=Soff​\(β\)2N​Soff​\(2​β\)\.\\displaystyle n\_\{\\mathrm\{eff\}\}\(\\beta\)=\\frac\{S\_\{\\mathrm\{off\}\}\(\\beta\)^\{2\}\}\{NS\_\{\\mathrm\{off\}\}\(2\\beta\)\}\.\(9\)Forkkequally weighted interactions per row this reduces tokk, while for nonuniform weights it gives the corresponding effective interaction count\. In the heat\-kernel regime it scales asneff∼Nβ−d/2n\_\{\\mathrm\{eff\}\}\\sim N\\beta^\{\-d/2\}\. The corresponding relative statistical uncertainty is:

𝒱stat​\(β\)=1neff​\(β\)​𝒱bias​\(β\)\.\\displaystyle\\mathcal\{V\}\_\{\\mathrm\{stat\}\}\(\\beta\)=\\frac\{1\}\{\\sqrt\{n\_\{\\mathrm\{eff\}\}\(\\beta\)\\mathcal\{V\}\_\{\\mathrm\{bias\}\}\(\\beta\)\}\}\.\(10\)Because bothSoff​\(β\)S\_\{\\mathrm\{off\}\}\(\\beta\)andSoff​\(2​β\)S\_\{\\mathrm\{off\}\}\(2\\beta\)are produced by consecutive probes of the dyadic sweep, the statistic requires no auxiliary kernel evaluation\. At largeβ\\beta, the diffusion operator approaches the identity and the leading spectral gaps eventually become comparable to the arithmetic precision of the kernel evaluation and eigensolver\. Define the normalized Perron gap asgP​\(β\)=\(λ0−λ1\)/\|λ0\|g\_\{P\}\(\\beta\)=\(\\lambda\_\{0\}\-\\lambda\_\{1\}\)/\|\\lambda\_\{0\}\|\.FlashDiffusionincorporates a smooth machine\-precision barrier,

𝒱mach​\(β\)=\(κP​ηgP​\(β\)\)p,\\displaystyle\\mathcal\{V\}\_\{\\mathrm\{mach\}\}\(\\beta\)=\\left\(\\frac\{\\kappa\_\{\\mathrm\{P\}\}\\eta\}\{g\_\{P\}\(\\beta\)\}\\right\)^\{p\},\(11\)with exponentp=6p=6and safety factorκP=10\\kappa\_\{\\mathrm\{P\}\}=10, together with the hard admissibility condition:

gP​\(β\)\>κP​η\.\\displaystyle g\_\{P\}\(\\beta\)\>\\kappa\_\{\\mathrm\{P\}\}\\eta\.\(12\)Hereη\\etais the numerical precision scale of the kernel data\-type333Withη=8×10−3\\eta=8\\times 10^\{\-3\}forbf16,10−310^\{\-3\}fortf32,10−710^\{\-7\}forfp32, and2\.22×10−162\.22\\times 10^\{\-16\}forfp64\. The defaults reported here usefp32throughout, withη=10−7\\eta=10^\{\-7\}\. The precision cascade described below is enabled only when multiple data\-types are requested\.The smooth term shapes the approach to the precision boundary from within the admissible region, while the hard condition prevents selection once the Perron and first nontrivial eigenvalues can no longer be reliably separated\. The sweep therefore follows a measured trajectory between two spectral limits\. Increasingβ\\betaprogressively resolves additional nontrivial modes\.FlashDiffusioncombines the heat\-action, finite\-sample, and machine\-precision terms into a dimensionless empirical potential onlog⁡β\\log\\beta,

𝒱⁡\(β\)=𝒱bias​\(β\)\+𝒱stat​\(β\)\+𝒱mach​\(β\)\.\\displaystyle\\mathcal\{V\}\(\\beta\)=\\mathcal\{V\}\_\{\\mathrm\{bias\}\}\(\\beta\)\+\\mathcal\{V\}\_\{\\mathrm\{stat\}\}\(\\beta\)\+\\mathcal\{V\}\_\{\\mathrm\{mach\}\}\(\\beta\)\.\(13\)The ideal bandwidthβ⋆\\beta^\{\\star\}is the minimizer of this finite\-sample resolution potential within the admissible spectral regime\. The discrete sweep returns an achieved bandwidthβ×≤β⋆\\beta^\{\\times\}\\leq\\beta^\{\\star\}, determined by the probes that can be resolved safely at the available sampling and arithmetic precision\. A minimum is accepted once three consecutive probes bracket a stationary point inlog⁡β\\log\\betafor which𝒱bias\\mathcal\{V\}\_\{\\mathrm\{bias\}\}is falling and either𝒱stat\\mathcal\{V\}\_\{\\mathrm\{stat\}\}or𝒱mach\\mathcal\{V\}\_\{\\mathrm\{mach\}\}is rising\. This direction of the flow identifies the bias\-variance well associated with the resolution boundary\. A quadratic fit inlog⁡β\\log\\betaover the bracketing triplet estimates the vertex, which is then resolved by a final probe\. If the Perron wall is reached before the bracket is completed, the last safe probe definesβ×\\beta^\{\\times\}\.

The continuation provides several geometric diagnostics without an additional kernel pass\. In the local\-manifold regime, the log\-derivative of the off\-diagonal kernel mass across a dyadic pair yields the empirical intrinsic dimension:

d^​\(β\)=−2​Δ​log⁡SoffΔ​log⁡β\.\\displaystyle\\widehat\{d\}\(\\beta\)=\-2\\frac\{\\Delta\\log S\_\{\\mathrm\{off\}\}\}\{\\Delta\\log\\beta\}\.\(14\)The corresponding Singer\-type scaling exponent can then be reported diagnostically\. These quantities characterize the observed flow and permit comparison with continuum asymptotics, while the selected bandwidth is determined from the measured finite\-sample spectrum rather than from a prescribed value ofd^\\widehat\{d\}or a hard\-coded scaling law\.

When multiple kernel data\-types are supplied,FlashDiffusionexecutes the continuation as a precision cascade\. A sweep is first performed with the earliest requested data\-type, and its selected state initializes the next precision level\. The final bandwidth and eigenspace are taken from the last data\-type in the cascade\. This separates the precision used to explore theβ\\betatrajectory from the precision required to resolve its final spectral boundary while preserving the same normalization, matrix\-free kernel application, and block Lanczos machinery at every stage\([Bradbury et al\., 2018](https://arxiv.org/html/2609.38198#bib.bib32);[Paszke et al\., 2019](https://arxiv.org/html/2609.38198#bib.bib33)\)\. The default configuration uses a single data\-type,fp32\. The complete algorithm therefore consists of three nested operations:XGEMMsupplies matrix\-free Gaussian applications, block implicitly restarted Lanczos converts those applications into finite spectral probes, and theβ\\betaflow uses successive probes to approach the finite\-sample resolution optimum\. The output is achieved with bandwidthβ×\\beta^\{\\times\}, and its associated nontrivial diffusion eigenvalues and the corresponding symmetric eigenvectors\. From which the Markov right eigenvectors are obtained by diagonal rescaling\.

## 3Experiment:T6T^\{6\}Manifold

We benchmarkFlashDiffusionon a manifold whose continuum spectrum is known in closed form\. The flat six\-torusT6=\(S1\)6T^\{6\}=\(S^\{1\}\)^\{6\}with unequal radii:

R=\(1\.00,1\.08,1\.17,1\.27,1\.38,1\.50\)\\displaystyle R=\(1\.00,\\,1\.08,\\,1\.17,\\,1\.27,\\,1\.38,\\,1\.50\)\(15\)each circle embedded as\(Racosθa,Rasinθa\)\(R\_\{a\}\\cos\\theta\_\{a\},\\,R\_\{a\}\\sin\\theta\_\{a\}\)with independent uniform anglesθa∼Uniform\[0,2π\)\\theta\_\{a\}\\sim\\mathrm\{Uniform\}\[0,2\\pi\), is a curvature\-free product manifold inℝ12\\mathbb\{R\}^\{12\}, zero\-padded toD=32D=32and sample\-centered\.444Only the first twelve ambient coordinates are non\-zero; the padding fixes the ambient dimension for the kernel without adding structure\. Samples are centered by subtracting the sample mean, which preserves the anisotropic radii and matches the centered input expected byFlashDiffusion\. All angles are drawn from a single seed \(4242\) perNNso that the validation diagnostic is reproducible; neither the seed nor the angles enter the solve, only the post\-hoc score of Figure[1](https://arxiv.org/html/2609.38198#S3.F1)c\.The continuum Laplace–Beltrami eigenvalues are∑ama2/Ra2\\sum\_\{a\}m\_\{a\}^\{2\}/R\_\{a\}^\{2\}with integer mode numbersmam\_\{a\}\. The unequal radii break the inter\-circle degeneracies, and becauseRmax/Rmin=1\.50<2R\_\{\\max\}/R\_\{\\min\}=1\.50<2the entire first\-harmonic ladder lies strictly below the first second\-harmonic level: the largest first harmonic is1/Rmin2=1\.001/R\_\{\\min\}^\{2\}=1\.00, while the smallest second harmonic is4/Rmax2≈1\.7784/R\_\{\\max\}^\{2\}\\approx 1\.778, leaving a gap of≈0\.78\\approx 0\.78between the first\-harmonic band and the rest of the spectrum\. Each circle contributes a degenerate doublet\{cos⁡\(m​θa\),sin⁡\(m​θa\)\}\\\{\\cos\(m\\theta\_\{a\}\),\\sin\(m\\theta\_\{a\}\)\\\}per harmonicmm, so the first\-harmonic band consists of six doublets; the leading doublet, that of the largest\-radius circle \(s=5s=5,Rs=1\.50R\_\{s\}=1\.50\), provides an exact target for the recovered modes\.

We compute the leadingk=3k=3nontrivial modes atα=1/2\\alpha=1/2infp32on NVIDIA RTX PRO 6000 \(Blackwell SM120\) GPUs\. Throughout,β⋆​\(N\)\\beta^\{\\star\}\(N\)denotes the ideal resolution scale of §[2\.3](https://arxiv.org/html/2609.38198#S2.SS3): the largest bandwidth at which the requested leading spectrum remains distinguishable from the sampling and numerical resolution floors\. Section[2\.3](https://arxiv.org/html/2609.38198#S2.SS3)shows thatβ⋆\\beta^\{\\star\}exists but is not computable from finite data; the in\-loop heat\-action selector returns an achievable approximationβ×​\(N\)≤β⋆​\(N\)\\beta^\{\\times\}\(N\)\\leq\\beta^\{\\star\}\(N\), subject to the Perron ceilingβ×\\beta^\{\\times\}that enforces the machine\-precision bandwidth floor\. The selector runs as it would on an unlabeled dataset, and the ground\-truth angles enter only the post\-hoc validation metric described below\. The third computed mode \(k=3k=3\) is retained only as a Lanczos buffer, letting the solver certify that the multiplicity\-two target doublet has closed before a possible fourth mode\.

The sample count varies over the power\-of\-two sequenceN=2pN=2^\{p\}, up\-to a dataset of222=4,194,3042^\{22\}=4\{,\}194\{,\}304observations\. Timing variance is largest for small problems while the per\-solve cost grows quadratically, so the number of repetitions decays geometrically from6464at the smallestNNto11at the largest, followingnrep​\(N\)=clip⁡\(⌈64\(22−log2⁡N\)/12⌉,1,64\)n\_\{\\mathrm\{rep\}\}\(N\)=\\mathrm\{clip\}\\left\(\\lceil 64^\{\(22\-\\log\_\{2\}N\)/12\}\\rceil,\\,1,\\,64\\right\), givingnrep​\(210\)=64n\_\{\\mathrm\{rep\}\}\(2^\{10\}\)=64,nrep​\(217\)=6n\_\{\\mathrm\{rep\}\}\(2^\{17\}\)=6, andnrep​\(222\)=1n\_\{\\mathrm\{rep\}\}\(2^\{22\}\)=1\. We report mean wall\-clock time with one standard deviation over repeats\.

2102^\{10\}2142^\{14\}2182^\{18\}2222^\{22\}10010^\{0\}10110^\{1\}10210^\{2\}10310^\{3\}10410^\{4\}10510^\{5\}NNwall time \(s\)\(a\) time to solution1 GPU4 GPUs∝N2\\propto N^\{2\}2102^\{10\}2142^\{14\}2182^\{18\}2222^\{22\}10−0\.510^\{\-0\.5\}10010^\{0\}NNβ⋆\\beta^\{\\star\}\(b\) calibrated bandwidthβ×​\(N\)\\beta^\{\\times\}\(N\)∝N2/\(d\+6\)=N1/6\\propto N^\{2/\(d\+6\)\}=N^\{1/6\}

2102^\{10\}2142^\{14\}2182^\{18\}2222^\{22\}10−410^\{\-4\}10−310^\{\-3\}10−210^\{\-2\}10−110^\{\-1\}NN1−1\-overlap\(c\) first\-harmonic error

Figure 1:FlashDiffusion on uniform samples of the flat torusT6T^\{6\}\(RTX PRO 6000\)\.\(a\)Wall time versus sample count; the dashed guide is exact𝒪⁡\(N2\)\\mathcal\{O\}\(N^\{2\}\)scaling, and the flat region belowN≈214N\\approx 2^\{14\}is the fixed launch/solver overhead\. Four\-GPU points continue the curve toN=222≈4\.2N=2^\{22\}\\approx 4\.2M\.\(b\)The bandwidth emitted by the\(n,β\)\(n,\\beta\)ladder, against the Singer exponent2/\(d\+6\)=1/62/\(d\+6\)=1/6ford=6d=6\(dashed, anchored at the first point\)\.\(c\)Recovery of the first harmonic band: projection error of the leading nontrivial eigenspace onto the exactT6T^\{6\}harmonics, decaying with sample count toward the resolution floor\. Error bars are one standard deviation over repeated runs where multiple runs were performed\.Figure[1](https://arxiv.org/html/2609.38198#S3.F1)a shows time to solution\. BelowN≈214N\\approx 2^\{14\}the runtime is flat at the fixed launch and solver overhead \(≈6\\approx 6s\); above it the measured cost approaches the𝒪⁡\(N2\)\\mathcal\{O\}\(N^\{2\}\)reference of the dense kernel application, eq\.[3](https://arxiv.org/html/2609.38198#S2.E3), with fitted exponentγ=0\.157\\gamma=0\.157over the quadratic regime\. A single GPU reachesN=220N=2^\{20\}in3737min, and four GPUs reachN=222N=2^\{22\}\(≈4\.2\\approx 4\.2M observations\) in3\.53\.5h\. A densefp32affinity matrix atN=222N=2^\{22\}would occupyN2×4​bytes≈7\.0×1013​bytesN^\{2\}\\times 4\\,\\mathrm\{bytes\}\\approx 7\.0\\times 10^\{13\}\\,\\mathrm\{bytes\}, or≈70\\approx 70TB\. AtN=220N=2^\{20\}the four\-GPU configuration completes in781781s against the single\-GPU22252225s, a2\.85×2\.85\\timesspeedup or71%71\\%parallel efficiency, with multi\-GPU runs sharded over the row axis\.

Figure[1](https://arxiv.org/html/2609.38198#S3.F1)b reports the selector’s calibrated bandwidthβ×​\(N\)\\beta^\{\\times\}\(N\), the achievable approximation to the ideal resolution scaleβ⋆​\(N\)\\beta^\{\\star\}\(N\)of §[2\.3](https://arxiv.org/html/2609.38198#S2.SS3), shown against the Singer reference slopeN2/\(d\+6\)=N1/6N^\{2/\(d\+6\)\}=N^\{1/6\}ford=6d=6anchored at the first point\. The reference exponent follows from the graph\-to\-manifold convergence of[Singer \(2006\)](https://arxiv.org/html/2609.38198#bib.bib4);[Hein et al\. \(2007\)](https://arxiv.org/html/2609.38198#bib.bib5), where the bandwidth must shrink with sample count\. For add\-dimensional manifold sampled with a Gaussian kernel, the resulting bandwidth scales asN−2/\(d\+6\)N^\{\-2/\(d\+6\)\}, thusβ⋆\\beta^\{\\star\}grows as the reciprocal\. Sinceβ×≤β⋆\\beta^\{\\times\}\\leq\\beta^\{\\star\}pointwise, the measured curve is a lower estimate of the ideal curve; the visible gap to the Singer guide can reflect either a genuine sub\-Singer scaling ofβ⋆\\beta^\{\\star\}or aβ⋆−β×\\beta^\{\\star\}\-\\beta^\{\\times\}bias whoseNN\-dependence is not separately known\. The measured slope over the quadratic regime is0\.1570\.157, to be compared with1/6≈0\.16671/6\\approx 0\.1667; the measured curve runs somewhat shallower,≈20%\\approx 20\\%below the Singer guide at the largestNN\. A visible break in the per\-NNerror bars nearN=217N=2^\{17\}coincides with the drop ofnrepn\_\{\\mathrm\{rep\}\}into the single\-digit regime and may also reflect a change in the selector’s effective safety margin as the spectrum hardens; it does not shift the central trend ofβ×​\(N\)\\beta^\{\\times\}\(N\)\.

Figure[1](https://arxiv.org/html/2609.38198#S3.F1)c reports recovery of the known leading eigenspace\. Because the continuum doublet is degenerate and may rotate arbitrarily, we score the two\-dimensional span rather than individual eigenvectors: the whitened cross\-covariance between the leading two nontrivial modes and\{cos⁡θs,sin⁡θs\}\\\{\\cos\\theta\_\{s\},\\sin\\theta\_\{s\}\\\}yields a rotation\-invariant overlap in\[0,1\]\[0,1\], accumulated block\-wise infp64so that the diagnostic remains stable at multi\-millionNN\.555WritingAAfor the matrix whose columns are the leading two nontrivial modes,BBfor the columns\(cos⁡θs,sin⁡θs\)\(\\cos\\theta\_\{s\},\\sin\\theta\_\{s\}\), andGA​AG\_\{AA\},GB​BG\_\{BB\},GA​BG\_\{AB\}for their centered covariance blocks, the score is∑i​jCi​j2/2\\sum\_\{ij\}C\_\{ij\}^\{2\}/2withC=GA​A−1/2GA​BGB​B−1/2C=G\_\{AA\}^\{\-1/2\}G\_\{AB\}G\_\{BB\}^\{\-1/2\}, the squared RV coefficient between the two subspaces\. It equals11iff the subspaces coincide and00iff they are orthogonal, up to floating\-point round\-off\. The accumulation is blocked at10610^\{6\}rows and carried infp64so that theN×2N\\times 2blocks are never materialized at full precision\.The subspace error:\(1−overlap\)\(1\-\\mathrm\{overlap\}\)decays from7×10−27\\times 10^\{\-2\}atN=210N=2^\{10\}to≈10−4\\approx 10^\{\-4\}atN=222N=2^\{22\}, confirming that the in\-loop calibration tracks the resolvable geometry as sampling deepens: the leading harmonic is recovered to three significant digits by the time the sample count reaches10610^\{6\}\.

## 4Experiment: Alanine Dipeptide

We evaluateFlashDiffusionon theMDSharealanine\-dipeptide benchmark\([Nüske et al\., 2017](https://arxiv.org/html/2609.38198#bib.bib38);[Wehmeyer and Noé, 2018](https://arxiv.org/html/2609.38198#bib.bib39)\), consisting of three independent250250ns molecular\-dynamics trajectories andN=750,000N=750\{,\}000configurations\. Each configuration contains the1010heavy atoms and is mapped to an E\(3\)\-invariant representation using distances from every heavy atom to the first88heavy atoms in the dataset ordering as anchors\([Candanedo, 2025a](https://arxiv.org/html/2609.38198#bib.bib34)\), givingD=80D=80features\.FlashDiffusionis run withα=1/2\\alpha=1/2infp32and computes256256nontrivial diffusion modes on eight NVIDIA RTX PRO 6000 GPUs\. The eigensolve required73157315s \(2\.032\.03h\), and explicit storage of the corresponding densefp32kernel would require750,0002×4≈2\.25750\{,\}000^\{2\}\\times 4\\approx 2\.25TB, which is once again avoided by the tiled matrix\-free operator\. This experiment uses the fullNNthroughout and exercises theβ\\betacontinuation directly rather than the sample\-size ladder\.

![Refer to caption](https://arxiv.org/html/2609.38198v1/figures/truth.png)−180∘\-180^\{\\circ\}0∘0^\{\\circ\}180∘180^\{\\circ\}−180∘\-180^\{\\circ\}0∘0^\{\\circ\}180∘180^\{\\circ\}ψ\\psiϕ\\phi\(a\) Reference![Refer to caption](https://arxiv.org/html/2609.38198v1/figures/m8.png)−180∘\-180^\{\\circ\}0∘0^\{\\circ\}180∘180^\{\\circ\}ψ\\psi\(b\) 8 DMAP modes![Refer to caption](https://arxiv.org/html/2609.38198v1/figures/m256.png)−180∘\-180^\{\\circ\}0∘0^\{\\circ\}180∘180^\{\\circ\}ψ\\psi\(c\) 256 DMAP modes![Refer to caption](https://arxiv.org/html/2609.38198v1/figures/color1.png)001122334455667788F/kB​TF/k\_\{B\}TFigure 2:Ramachandran free\-energy distributions for alanine dipeptide\.\(a\)Reference distribution,\(b\)reconstruction using 8 diffusion\-map \(DMAP\) modes, and\(c\)reconstruction using 256 DMAP modes\. The horizontal coordinate is the dihedral angleψ\\psiand the vertical coordinate isϕ\\phi; both span\[−180∘,180∘\]\[\-180^\{\\circ\},180^\{\\circ\}\], with\(ϕ,ψ\)=\(0,0\)\(\\phi,\\psi\)=\(0,0\)at the center of each panel\. The shared color bar reports the free energyF/kB​TF/k\_\{B\}Tfrom 0 to 8\.212^\{1\}222^\{2\}232^\{3\}242^\{4\}252^\{5\}262^\{6\}272^\{7\}282^\{8\}10102020mψ=16m\_\{\\psi\}=16D=80D=80Number of diffusion modesmmHeld\-out angular MAE \(deg\)\(a\) Circular decoding atβ⋆=276\.4\\beta^\{\\star\}=276\.4ϕ\\phi\(held\-out traj\.\)ψ\\psi\(held\-out traj\.\)ϕ\\phi\(random holdout\)ψ\\psi\(random holdout\)![Refer to caption](https://arxiv.org/html/2609.38198v1/figures/AlanineDipeptide_traj.png)−150\-150−100\-100−50\-50005050100100150150−150\-150−100\-100−50\-50005050100100150150ϕ\\phi\(deg\)ψ\\psi\(deg\)\(b\) Held\-out Ramachandran surface005510101515202025253030joint error \(deg\)001122334455667788F/kB​TF/k\_\{B\}TFigure 3:Spectral resolution of alanine\-dipeptide coordinates\.\(a\)Mean absolute wrapped error for the backbone dihedralsϕ\\phiandψ\\psias the number of retained diffusion modesmmincreases\. Solid curves evaluate the completely held\-out trajectory, while dashed curves evaluate5,0005\{,\}000randomly held\-out frames from the decoder\-training trajectories\. The vertical markers indicate the ambient feature dimensionD=80D=80and the detected onsetmψ=16m\_\{\\psi\}=16of the principalψ\\psi\-resolution regime\.\(b\)Held\-out Ramachandran free\-energy surface atβ⋆=276\.4\\beta^\{\\star\}=276\.4, with random\-holdout frames colored by joint angular reconstruction error and a strided segment of the decoded held\-out trajectory overlaid as a curve\.Theβ\\betaflow required eight dyadic probes and selectedβ⋆=276\.4\\beta^\{\\star\}=276\.4from the empirical potential\. The sweep continued beyond the selected point to characterize the precision boundary:β=552\.8\\beta=552\.8remained admissible, whileβ=1105\.7\\beta=1105\.7crossed thefp32Perron wall withgP=3\.51×10−9g\_\{P\}=3\.51\\times 10^\{\-9\}\. Atβ⋆\\beta^\{\\star\}, the effective interaction count isneff≈4\.31×104n\_\{\\mathrm\{eff\}\}\\approx 4\.31\\times 10^\{4\}, approximately5\.7%5\.7\\%of the sample count, and the measured local dimension isd^≈4\.85\\widehat\{d\}\\approx 4\.85\. Hered^\\widehat\{d\}is a kernel\-scale diagnostic of the resolved geometry in the E\(3\)\-invariant representation and should not be identified with the two\-dimensional\(ϕ,ψ\)\(\\phi,\\psi\)reaction coordinate alone\. The selectedβ⋆\\beta^\{\\star\}therefore lies comfortably inside the resolved spectral regime rather than at thefp32precision boundary\. Following the diffusion\-map autoencoder construction of[Candanedo \(2024\)](https://arxiv.org/html/2609.38198#bib.bib2), we fit a single linear decoder from the firstmmdiffusion coordinates to the circular targets\(cos⁡ϕ,sin⁡ϕ,cos⁡ψ,sin⁡ψ\)\(\\cos\\phi,\\sin\\phi,\\cos\\psi,\\sin\\psi\)\. The decoder is fit by least squares with ridge regularization10−8​tr⁡\(𝒢\(m\)\)/m10^\{\-8\}\\operatorname\{tr\}\(\\mathcal\{G\}^\{\(m\)\}\)/m, where𝒢\(m\)\\mathcal\{G\}^\{\(m\)\}is the centered decoder Gram matrix\. Two trajectories supply the decoder\-fit population, from which5,0005,000frames are drawn uniformly without replacement as an additional random holdout, while the third250250ns trajectory is excluded entirely from decoder fitting\. The diffusion eigensolve itself remains unsupervised over all configurations; only the decoder uses the split\. We evaluatem∈\{2,4,8,16,32,64,128,256\}m\\in\\\{2,4,8,16,32,64,128,256\\\}using the mean absolute wrapped angular error\|arg⁡exp⁡\(i⁡\(θ^−θ\)\)\|\\big\|\\arg\\exp\(i\(\\widehat\{\\theta\}\-\\theta\)\)\\big\|\. Figure[3](https://arxiv.org/html/2609.38198#S4.F3)shows distinct resolution profiles for the two backbone coordinates\. The held\-outϕ\\phierror falls from26\.69∘26\.69^\{\\circ\}atm=2m=2to7\.71∘7\.71^\{\\circ\}atm=8m=8, whileψ\\psiremains near15∘15^\{\\circ\}throughm=8m=8and begins its principal decrease atmψ=16m\_\{\\psi\}=16\. Atm=256m=256, the completely held\-out trajectory reaches4\.86∘4\.86^\{\\circ\}MAE forϕ\\phiand3\.95∘3\.95^\{\\circ\}forψ\\psi, compared with4\.90∘4\.90^\{\\circ\}and3\.96∘3\.96^\{\\circ\}on the random holdout\. The differences are below0\.05∘0\.05^\{\\circ\}, so the fully held\-out trajectory shows no observable degradation relative to frames drawn from the decoder\-training trajectories at the reported precision\. Reconstruction also continues beyond the ambient feature dimensionD=80D=80: fromm=128m=128tom=256m=256, the held\-out errors decrease from5\.26∘5\.26^\{\\circ\}to4\.86∘4\.86^\{\\circ\}forϕ\\phiand from4\.42∘4\.42^\{\\circ\}to3\.95∘3\.95^\{\\circ\}forψ\\psi\. Figure[3](https://arxiv.org/html/2609.38198#S4.F3)additionally compares the held\-out Ramachandran surface with decoded ensembles atm=8m=8andm=256m=256; the latter recovers the principal metastable regions and their geometry, while a strided held\-out trajectory segment has mean joint angular error7\.1∘7\.1^\{\\circ\}\. The spectrum confirms that these additional modes remain part of a substantial high\-rank operator\. All256256computed eigenvalues are finite, positive, and correctly ordered, withλmax=0\.9985\\lambda\_\{\\max\}=0\.9985andλ256=0\.2127\\lambda\_\{256\}=0\.2127\. The final computed eigenvalue therefore remains approximately21%21\\%of the leading value and lies more than six orders of magnitude above thefp32working precision scale used by the sweep\. Together with the continued reconstruction improvement form\>Dm\>D, this shows that the useful diffusion representation extends well beyond a low\-rank approximation of the8080\-dimensional input features and provides a molecular example of the high\-rank Gaussian regime targeted byFlashDiffusion\.

## 5Conclusion

FlashDiffusionmakes the pragmatic choice to retain the dense𝒪⁡\(N2\)\\mathcal\{O\}\(N^\{2\}\)Gaussian interaction rather than replace it by a sparse or low\-rank surrogate\. GPU tiling removes the need to materialize the kernel, alternatives that achieve this exist in the literature includingKeOps[Feydy et al\. \(2020\)](https://arxiv.org/html/2609.38198#bib.bib43)and we run a comparison with our primitiveXGEMMin §[B](https://arxiv.org/html/2609.38198#A2)\. However, unlike previous work we consider a synergistic connection betweenβ\\beta\-flow and tiling requirements, via the\(n,β\)\(n,\\beta\)continuation introduced in §[2\.3](https://arxiv.org/html/2609.38198#S2.SS3), to reduce the number of expensive full\-data spectral solves\.

This targets the regime in which the Gaussian operator itself is high rank\. Both fixed\-kkNN graphs and low\-rank approximations remain effective when only a small spectral subspace is required\. For data whose useful spectral rank is unknown in advance,FlashDiffusioninstead treats rank and spectral resolution as quantities to be measured from the finite novel sample\.

Theβ\\betaflow also admits a scale interpretation related to renormalization ideas\([Zinn\-Justin, 2019](https://arxiv.org/html/2609.38198#bib.bib7)\)\. For a Gaussian heat kernel, the diffusion time satisfiest∝β−1t\\propto\\beta^\{\-1\}, hence smallβ\\betacorresponds to a long\-time, coarse description and largeβ\\betato progressively shorter\-time, finer structure\.FlashDiffusionfollows this scale dependence in the coarse\-to\-fine direction as sampling increases\. The ideal resolution scaleβ⋆​\(N\)\\beta^\{\\star\}\(N\)characterizes the finite\-sample optimum, while the algorithm returns an achievableβ×\\beta^\{\\times\}constrained by sampling and arithmetic precision\. The accompanying diagnosticsneffn\_\{\\mathrm\{eff\}\}andd^\\widehat\{d\}describe how the effective neighborhood size and local spectral dimension evolve along the same flow\.

The principal remaining computational limitation is the feature dimension\. Each Gaussian application costs𝒪⁡\(N2​D\)\\mathcal\{O\}\(N^\{2\}D\), and therefore becomes cubic whenD∼ND\\sim N\. The experiments considered here remain well inside the quadratic regime, withD/N≈8×10−6D/N\\approx 8\\times 10^\{\-6\}onT6T^\{6\}andD/N≈10−4D/N\\approx 10^\{\-4\}for alanine dipeptide\. Since the kernel depends on the representation through pairwise distances, feature compression provides a natural front end: PCA, Johnson–Lindenstrauss projections\([Johnson and Lindenstrauss, 1984](https://arxiv.org/html/2609.38198#bib.bib40)\), or learned encoders can reduceDDprovided the induced distance distortion is controlled\. This suggests a useful division of labor: the front end determines the geometry presented to the kernel, while diffusion maps perform harmonic analysis on that geometry\. Jointly calibrating representation distortion and spectral resolution is a natural direction for future work\.

Future work includes theoretical guarantees relating the empirical potential minimum toβ⋆\\beta^\{\\star\}, joint calibration of feature\-space compression and diffusion resolution, extensions to other radial kernel families, and applications where high\-rank spectral structure is expected to be essential, e\.g\. in molecular dynamics,[Gowers et al\. \(2016\)](https://arxiv.org/html/2609.38198#bib.bib14)\.

## References

- Baglamaet al\.\(2003\)J\. Baglama, D\. Calvetti, and L\. ReichelIRBL: an implicitly restarted block\-lanczos method for large\-scale hermitian eigenproblems\.[SIAM J\. Sci\. Comput\.](https://doi.org/10.1137/S1064827501397949)24\(5\),pp\.1650–1677\.Cited by:[§1](https://arxiv.org/html/2609.38198#S1.p3.1),[§2\.3](https://arxiv.org/html/2609.38198#S2.SS3.p1.1)\.
- Belkin and Niyogi \(2003\)M\. Belkin and P\. NiyogiLaplacian eigenmaps for dimensionality reduction and data representation\.[Neural Computation](https://dl.acm.org/doi/10.1162/089976603321780317)15\(6\),pp\.1373–1396\.External Links:[Document](https://dx.doi.org/10.1162/089976603321780317)Cited by:[§1](https://arxiv.org/html/2609.38198#S1.p1.1)\.
- Bengioet al\.\(2003\)Y\. Bengio, J\. Paiement, P\. Vincent, O\. Delalleau, N\. Le Roux, and M\. OuimetOut\-of\-sample extensions for LLE, Isomap, MDS, eigenmaps, and spectral clustering\.In[NeurIPS](https://proceedings.neurips.cc/paper/2003/hash/cf05968255451bdefe3c5bc64d550517-Abstract.html),pp\.177–184\.Cited by:[§1](https://arxiv.org/html/2609.38198#S1.p4.1)\.
- Boutsidiset al\.\(2009\)C\. Boutsidis, M\. W\. Mahoney, and P\. DrineasAn improved approximation algorithm for the column subset selection problem\.In[SODA](https://doi.org/10.1137/1.9781611973068.105),pp\.968–977\.Note:[arXiv:0812\.4293](https://arxiv.org/abs/0812.4293)Cited by:[§A\.5](https://arxiv.org/html/2609.38198#A1.SS5.p1.2)\.
- Bradburyet al\.\(2018\)[JAX: composable transformations of Python\+NumPy programs](https://github.com/jax-ml/jax)Cited by:[§2\.3](https://arxiv.org/html/2609.38198#S2.SS3.p4.1)\.
- Calder and García Trillos \(2022\)J\. Calder and N\. García TrillosImproved spectral convergence rates for graph laplacians onε\\varepsilon\-graphs andkk\-NN graphs\.[ACHA](https://doi.org/10.1016/j.acha.2022.02.004)60,pp\.123–175\.Cited by:[§A\.2](https://arxiv.org/html/2609.38198#A1.SS2.p1.3),[§A\.5](https://arxiv.org/html/2609.38198#A1.SS5.p1.1)\.
- Candanedo \(2024\)J\. CandanedoDiffusion map autoencoder\.[arXiv:2409\.05901](https://arxiv.org/abs/2409.05901)\.Cited by:[§4](https://arxiv.org/html/2609.38198#S4.p2.1)\.
- Candanedo \(2025a\)J\. CandanedoE\(3\) invariant representations of biomolecules\.[bioRxiv:2025\.03\.26\.645520](https://www.biorxiv.org/content/10.1101/2025.03.26.645520v1)\.Cited by:[§4](https://arxiv.org/html/2609.38198#S4.p1.1)\.
- Candanedo \(2025b\)J\. CandanedoLinearized diffusion map\.[arXiv:2507\.14257](https://arxiv.org/abs/2507.14257)\.Cited by:[§1](https://arxiv.org/html/2609.38198#S1.p2.1)\.
- Candanedo \(2026\)J\. CandanedoThe diffusion attention connection\.[arXiv:2604\.09560](https://arxiv.org/pdf/2604.09560)\.Cited by:[§1](https://arxiv.org/html/2609.38198#S1.p3.1)\.
- Charlieret al\.\(2021\)B\. Charlier, J\. Feydy, J\. A\. Glaunès, F\. Collin, and G\. DurifKernel operations on the GPU, with autodiff, without memory overflows\.[J\. Mach\. Learn\. Res\.](https://jmlr.org/papers/v22/20-275.html)22\(74\),pp\.1–6\.Note:[arXiv:2004\.11127](https://arxiv.org/abs/2004.11127);[GitHub:getkeops/keops](https://github.com/getkeops/keops/tree/main)Cited by:[Appendix B](https://arxiv.org/html/2609.38198#A2.p1.1)\.
- Clemente and D’Elia \(2018\)G\. Clemente and M\. D’EliaSpectrum of the laplace\-beltrami operator and the phase structure of causal dynamical triangulation\.[Phys\. Rev\. D](https://doi.org/10.1103/PhysRevD.97.124022)97,pp\.124022\.Note:[arXiv:1804\.02294](https://arxiv.org/abs/1804.02294)Cited by:[§1](https://arxiv.org/html/2609.38198#S1.p2.1)\.
- Coifman and Lafon \(2006\)R\. R\. Coifman and S\. LafonDiffusion maps\.[ACHA](https://doi.org/10.1016/j.acha.2006.04.006)21\(1\),pp\.5–30\.Cited by:[§2\.1](https://arxiv.org/html/2609.38198#S2.SS1.p1.2)\.
- Coifmanet al\.\(2008\)R\. R\. Coifman, Y\. Shkolnisky, F\. J\. Sigworth, and A\. SingerGraph laplacian tomography from unknown random projections\.[IEEE Trans\. Image Process\.](https://doi.org/10.1109/TIP.2008.2002305)17\(10\),pp\.1891–1899\.Cited by:[§1](https://arxiv.org/html/2609.38198#S1.p1.1)\.
- Daoet al\.\(2022\)T\. Dao, D\. Y\. Fu, S\. Ermon, A\. Rudra, and C\. RéFlashAttention: fast and memory\-efficient exact attention with IO\-awareness\.In[NeurIPS](https://proceedings.neurips.cc/paper_files/paper/2022/hash/67d57c32e20fd0a7a302cb81d36e40d5-Abstract-Conference.html),Note:[arXiv:2205\.14135](https://arxiv.org/abs/2205.14135)Cited by:[§1](https://arxiv.org/html/2609.38198#S1.p3.1),[§2\.2](https://arxiv.org/html/2609.38198#S2.SS2.p1.1)\.
- Drineaset al\.\(2008\)P\. Drineas, M\. W\. Mahoney, and S\. MuthukrishnanRelative\-error CUR matrix decompositions\.[SIAM J\. Matrix Anal\. Appl\.](https://doi.org/10.1137/07070471X)30\(2\),pp\.844–881\.Cited by:[§A\.5](https://arxiv.org/html/2609.38198#A1.SS5.p1.2)\.
- Feydyet al\.\(2020\)J\. Feydy, J\. A\. Glaunès, B\. Charlier, and M\. M\. BronsteinFast geometric learning with symbolic matrices\.In[NeurIPS](https://proceedings.neurips.cc/paper/2020/hash/a6292668b36ef412fa3c4102d1311a62-Abstract.html),Vol\.33\.Note:[GitHub:getkeops/keops](https://github.com/getkeops/keops)Cited by:[Appendix B](https://arxiv.org/html/2609.38198#A2.p1.1),[§5](https://arxiv.org/html/2609.38198#S5.p1.1)\.
- Gowerset al\.\(2016\)R\. J\. Gowers, M\. Linke, J\. Barnoud, T\. J\. E\. Reddy, M\. N\. Melo, S\. L\. Seyler, J\. Dotson, J\. Domanski, S\. Buchoux, I\. M\. Kenney, and O\. BecksteinMDAnalysis: a python package for the rapid analysis of molecular dynamics simulations\.In[SciPy](https://doi.org/10.25080/Majora-629e541a-00e),pp\.98–105\.Cited by:[§5](https://arxiv.org/html/2609.38198#S5.p5.1)\.
- Greengard and Strain \(1991\)L\. Greengard and J\. StrainThe fast gauss transform\.[SIAM J\. Sci\. Stat\. Comput\.](https://doi.org/10.1137/0912004)12\(1\),pp\.79–94\.Cited by:[footnote 1](https://arxiv.org/html/2609.38198#footnote1)\.
- Heinet al\.\(2007\)M\. Hein, J\. Audibert, and U\. von LuxburgGraph laplacians and their convergence on random neighborhood graphs\.[J\. Mach\. Learn\. Res\.](https://jmlr.org/papers/v8/hein07a.html)8,pp\.1325–1368\.Note:[arXiv:math/0608522](https://arxiv.org/abs/math/0608522)Cited by:[§A\.5](https://arxiv.org/html/2609.38198#A1.SS5.p1.1),[§1](https://arxiv.org/html/2609.38198#S1.p2.1),[§3](https://arxiv.org/html/2609.38198#S3.p5.1)\.
- Hertrich \(2024\)J\. HertrichFast kernel summation in high dimensions via slicing and fourier transforms\.[SIAM J\. Math\. Data Sci\.](https://doi.org/10.1137/24M1632085)6\(4\),pp\.1109–1137\.Note:[arXiv:2401\.08260](https://arxiv.org/abs/2401.08260)Cited by:[footnote 1](https://arxiv.org/html/2609.38198#footnote1)\.
- Johnson and Lindenstrauss \(1984\)W\. B\. Johnson and J\. LindenstraussExtensions of lipschitz mappings into a hilbert space\.In[Contemp\. Math\.](https://doi.org/10.1090/conm/026/737400),Vol\.26,pp\.189–206\.Cited by:[§5](https://arxiv.org/html/2609.38198#S5.p4.1)\.
- Lanczos \(1950\)C\. LanczosAn iteration method for the solution of the eigenvalue problem of linear differential and integral operators\.[J\. Res\. Natl\. Bur\. Stand\.](https://doi.org/10.6028/jres.045.026)45\(4\),pp\.255–282\.Cited by:[§1](https://arxiv.org/html/2609.38198#S1.p3.1)\.
- Long and Ferguson \(2019\)A\. W\. Long and A\. L\. FergusonLandmark diffusion maps \(L\-dMaps\): accelerated manifold learning out\-of\-sample extension\.[ACHA](https://doi.org/10.1016/j.acha.2017.08.004)47\(1\),pp\.190–211\.Note:[arXiv:1706\.09396](https://arxiv.org/abs/1706.09396)Cited by:[footnote 1](https://arxiv.org/html/2609.38198#footnote1)\.
- Malkov and Yashunin \(2020\)Yu\. A\. Malkov and D\. A\. YashuninEfficient and robust approximate nearest neighbor search using hierarchical navigable small world graphs\.[IEEE Trans\. Pattern Anal\. Mach\. Intell\.](https://doi.org/10.1109/TPAMI.2018.2889473)42\(4\),pp\.824–836\.Note:[arXiv:1603\.09320](https://arxiv.org/abs/1603.09320)Cited by:[§A\.5](https://arxiv.org/html/2609.38198#A1.SS5.p1.2)\.
- McQueenet al\.\(2016\)J\. McQueen, M\. Meila, J\. VanderPlas, and Z\. ZhangMegaman: manifold learning with millions of points\.[arXiv:1603\.02763](https://arxiv.org/abs/1603.02763)\.Cited by:[Appendix A](https://arxiv.org/html/2609.38198#A1.p1.1)\.
- Milakov and Gimelshein \(2018\)M\. Milakov and N\. GimelsheinOnline normalizer calculation for softmax\.[arXiv:1805\.02867](https://arxiv.org/abs/1805.02867)\.Cited by:[§2\.1](https://arxiv.org/html/2609.38198#S2.SS1.p1.3)\.
- Muja and Lowe \(2014\)M\. Muja and D\. G\. LoweScalable nearest neighbor algorithms for high dimensional data\.[IEEE Trans\. Pattern Anal\. Mach\. Intell\.](https://doi.org/10.1109/TPAMI.2014.2321376)36\(11\),pp\.2227–2240\.Cited by:[§1](https://arxiv.org/html/2609.38198#S1.p1.1)\.
- Nüskeet al\.\(2017\)F\. Nüske, H\. Wu, J\. Prinz, C\. Wehmeyer, C\. Clementi, and F\. NoéMarkov state models from short non\-equilibrium simulations—analysis and correction of estimation bias\.[J\. Chem\. Phys\.](https://doi.org/10.1063/1.4976518)146\(9\),pp\.094104\.Note:[arXiv:1701\.01665](https://arxiv.org/abs/1701.01665)External Links:[Document](https://dx.doi.org/10.1063/1.4976518)Cited by:[§4](https://arxiv.org/html/2609.38198#S4.p1.1)\.
- NVIDIA \(2020\)NVIDIA[NVIDIA A100 Tensor Core GPU Architecture](https://www.nvidia.com/content/dam/en-zz/Solutions/Data-Center/nvidia-ampere-architecture-whitepaper.pdf)\.Note:NVIDIA architecture whitepaperCited by:[§2\.2](https://arxiv.org/html/2609.38198#S2.SS2.p2.3)\.
- NVIDIA \(2024\)NVIDIA[NVIDIA Blackwell Architecture Technical Brief](https://resources.nvidia.com/en-us-blackwell-architecture/blackwell-architecture-technical-brief)\.Note:NVIDIA technical briefCited by:[§2\.2](https://arxiv.org/html/2609.38198#S2.SS2.p2.3)\.
- Nyström \(1930\)E\. J\. NyströmÜBer die praktische auflösung von integralgleichungen mit anwendungen auf randwertaufgaben\.[Acta Math\.](https://doi.org/10.1007/BF02547521)54\(1\),pp\.185–204\.Cited by:[§2\.3](https://arxiv.org/html/2609.38198#S2.SS3.p2.1)\.
- Paszkeet al\.\(2019\)A\. Paszke, S\. Gross, F\. Massa, A\. Lerer, J\. Bradbury, G\. Chanan, T\. Killeen, Z\. Lin, N\. Gimelshein, L\. Antiga, A\. Desmaison, A\. Köpf, E\. Yang, Z\. DeVito, M\. Raison, A\. Tejani, S\. Chilamkurthy, B\. Steiner, L\. Fang, J\. Bai, and S\. ChintalaPyTorch: an imperative style, high\-performance deep learning library\.In[NeurIPS](https://proceedings.neurips.cc/paper/2019/hash/bdbca288fee7f92f2bfa9f7012727740-Abstract.html),pp\.8024–8035\.Cited by:[§2\.3](https://arxiv.org/html/2609.38198#S2.SS3.p4.1)\.
- Rahimi and Recht \(2007\)A\. Rahimi and B\. RechtRandom features for large\-scale kernel machines\.In[NeurIPS](https://papers.nips.cc/paper/2007/hash/013a006f03dbc5392effeb8f18fda755-Abstract.html),pp\.1177–1184\.Cited by:[footnote 1](https://arxiv.org/html/2609.38198#footnote1)\.
- Rahimi and Recht \(2008\)A\. Rahimi and B\. RechtWeighted sums of random kitchen sinks: replacing minimization with randomization in learning\.In[NeurIPS](https://papers.nips.cc/paper_files/paper/2008/hash/0efe32849d230d7f53049ddc4a4b0c60-Abstract.html),pp\.1313–1320\.Cited by:[footnote 1](https://arxiv.org/html/2609.38198#footnote1)\.
- Raykar and Duraiswami \(2005\)V\. C\. Raykar and R\. DuraiswamiThe improved fast gauss transform with applications to machine learning\.In[NIPS Workshop on Large Scale Kernel Machines](https://users.umiacs.umd.edu/~ramanid/pubs/Raykar_Duraiswami_IFGT_NIPS.pdf),Cited by:[footnote 1](https://arxiv.org/html/2609.38198#footnote1)\.
- Singer \(2006\)A\. SingerFrom graph to manifold laplacian: the convergence rate\.[ACHA](https://doi.org/10.1016/j.acha.2006.03.004)21\(1\),pp\.128–134\.Cited by:[§A\.2](https://arxiv.org/html/2609.38198#A1.SS2.p1.3),[§1](https://arxiv.org/html/2609.38198#S1.p2.1),[§3](https://arxiv.org/html/2609.38198#S3.p5.1)\.
- Sorensen \(1992\)D\. C\. SorensenImplicit application of polynomial filters in akk\-step arnoldi method\.[SIAM J\. Matrix Anal\. Appl\.](https://doi.org/10.1137/0613025)13\(1\),pp\.357–385\.Cited by:[§2\.3](https://arxiv.org/html/2609.38198#S2.SS3.p1.1)\.
- Thakkaret al\.\(2023\)V\. Thakkaret al\.[CUTLASS 3\.0](https://github.com/NVIDIA/cutlass/tree/v3.0.0)\.Note:CUDA Templates for Linear Algebra SubroutinesCited by:[§2\.2](https://arxiv.org/html/2609.38198#S2.SS2.p2.3)\.
- Tinget al\.\(2010\)D\. Ting, L\. Huang, and M\. I\. JordanAn analysis of the convergence of graph laplacians\.In[ICML](https://icml.cc/Conferences/2010/papers/554.pdf),pp\.1079–1086\.Note:[arXiv:1101\.5435](https://arxiv.org/abs/1101.5435)Cited by:[§A\.5](https://arxiv.org/html/2609.38198#A1.SS5.p1.1)\.
- Wehmeyer and Noé \(2018\)C\. Wehmeyer and F\. NoéTime\-lagged autoencoders: deep learning of slow collective variables for molecular kinetics\.[J\. Chem\. Phys\.](https://doi.org/10.1063/1.5011399)148\(24\),pp\.241703\.Note:[arXiv:1710\.11239](https://arxiv.org/abs/1710.11239)External Links:[Document](https://dx.doi.org/10.1063/1.5011399)Cited by:[§4](https://arxiv.org/html/2609.38198#S4.p1.1)\.
- Williams and Seeger \(2000\)C\. K\. I\. Williams and M\. SeegerUsing the nyström method to speed up kernel machines\.In[NeurIPS](https://papers.nips.cc/paper_files/paper/2000/hash/19de10adbaa1b2ee13f77f679fa1483a-Abstract.html),pp\.682–688\.Cited by:[footnote 1](https://arxiv.org/html/2609.38198#footnote1)\.
- Zadouriet al\.\(2026\)T\. Zadouri, M\. Hoehnerbach, J\. Shah, T\. Liu, V\. Thakkar, and T\. DaoFlashAttention\-4: algorithm and kernel pipelining co\-design for asymmetric hardware scaling\.In[MLSys](https://arxiv.org/abs/2603.05451),Note:[arXiv:2603\.05451](https://arxiv.org/abs/2603.05451);[GitHub:Dao\-AILab/flash\-attention](https://github.com/Dao-AILab/flash-attention)Cited by:[§1](https://arxiv.org/html/2609.38198#S1.p3.1)\.
- Zelnik\-Manor and Perona \(2004\)L\. Zelnik\-Manor and P\. PeronaSelf\-tuning spectral clustering\.In[NeurIPS](https://papers.nips.cc/paper_files/paper/2004/hash/40173ea48d9567f1f393b20c855bb40b-Abstract.html),pp\.1601–1608\.Cited by:[§1](https://arxiv.org/html/2609.38198#S1.p1.1)\.
- Zinn\-Justin \(2019\)J\. Zinn\-Justin[From Random Walks to Random Matrices: Selected Topics in Modern Theoretical Physics](https://global.oup.com/academic/product/from-random-walks-to-random-matrices-9780198787754)\.Oxford Graduate Texts,Oxford University Press\.External Links:ISBN 9780198787754Cited by:[§5](https://arxiv.org/html/2609.38198#S5.p3.1)\.

## Appendix ASpectral Consequences of kNN/ANN Sparsification

A common strategy for scaling diffusion maps iskk\-nearest\-neighbor \(kkNN\) sparsification: retain a fixed number of neighbors per sample, store the resulting weighted graph, and apply a sparse eigensolver\. This strategy is representative of many large\-scale spectral\-embedding pipelines that combine approximate neighborhood graphs with iterative solvers, such asmegaman\([McQueen et al\., 2016](https://arxiv.org/html/2609.38198#bib.bib13)\)\. The comparison becomes more consequential when the objective is to recover a substantial spectral basis instead of only a low\-dimensional visualization\. The applications considered here use on the order ofM∼102M\\sim 10^\{2\}modes\. As discussed in the main text, at the bandwidths relevant to this regime the Gaussian kernel generally has substantial effective rank, which limits the effectiveness of very low\-rank approximations\. The matrix\-free formulation therefore preserves the dense Gaussian operator while avoiding both explicit matrix construction and a fixed sparse neighborhood representation\.

The first few eigenfunctions typically vary on spatial scales comparable to the manifold diameter and can therefore be comparatively insensitive to moderate changes in local graph construction\. Higher eigenfunctions vary on progressively shorter scales\. Weyl counting gives a characteristic wavelength of order\(V/M\)1/d\(V/M\)^\{1/d\}for theMM\-th mode on add\-dimensional manifold of volumeVV\. Consequently, a neighborhood graph that gives a visually faithful two\- or three\-dimensional embedding may preserve substantially less information than is required for a larger spectral basis\. We distinguish two questions below\. First, how many neighbors are required for a truncated graph to reproduce the numerically relevant support of the dense Gaussian operator at a prescribed tolerance? Second, what structural changes does hard neighborhood truncation introduce when the retained neighborhood is much smaller than this support scale? The first question leads to a small\-ball scaling law, while the second is addressed through retained row mass, PSD violation, and finite\-sample spectral considerations\.

### A\.1Arithmetic and neighborhood truncation

Consider the dense Gaussian kernel𝖪i​j​\(β\)=exp⁡\(−β​∥xi−xj∥2\)\\mathsf\{K\}\_\{ij\}\(\\beta\)=\\exp\\left\(\-\\beta\\lVert x\_\{i\}\-x\_\{j\}\\rVert^\{2\}\\right\)\. Although every entry is strictly positive in exact arithmetic, numerical computation introduces a practical relative scale below which individual contributions become increasingly insignificant compared withO⁡\(1\)O\(1\)accumulated quantities\. We therefore introduce a working toleranceτ\\tauand regard interactions satisfying𝖪i​j​\(β\)<τ\\mathsf\{K\}\_\{ij\}\(\\beta\)<\\tauas negligible at the requested arithmetic precision\. We useτ≈10−7\\tau\\approx 10^\{\-7\}as anfp32\-scale working tolerance, comparable to the relative machine\-epsilon scale, andτ≈10−15\\tau\\approx 10^\{\-15\}as anfp64\-scale tolerance\. The resulting estimates depend only logarithmically onτ\\tau\. The corresponding support radius is:

Dτ=log⁡\(1/τ\)β\.\\displaystyle D\_\{\\tau\}=\\sqrt\{\\frac\{\\log\(1/\\tau\)\}\{\\beta\}\}\.\(16\)Interactions outside a ball of radiusDτD\_\{\\tau\}therefore have weight below the prescribed working tolerance\. AkkNN construction imposes a different truncation at thekk\-th\-neighbor distanceRk​\(xi\)R\_\{k\}\(x\_\{i\}\)\. This radius is row dependent and varies with the local sampling densityq⁡\(xi\)q\(x\_\{i\}\)\. For the quantitative estimates below we first consider uniform sampling,q⁡\(x\)≡q=1/Vq\(x\)\\equiv q=1/V, and letωd\\omega\_\{d\}denote the volume of the unitdd\-ball\. Locally,kkneighbors occupy a volume of orderk/\(N​q\)k/\(Nq\)\. For fixedβ\\beta, let𝖪~β,k=𝖪β⊙Ak\\widetilde\{\\mathsf\{K\}\}\_\{\\beta,k\}=\\mathsf\{K\}\_\{\\beta\}\\odot A\_\{k\}denote the symmetrically masked kernel, withAkA\_\{k\}the neighborhood mask\. Reproducing all dense\-kernel interactions above toleranceτ\\taurequires approximatelyRk​\(xi\)≳DτR\_\{k\}\(x\_\{i\}\)\\gtrsim D\_\{\\tau\}for every row\. Under nonuniform sampling, this condition becomes spatially heterogeneous, so a single fixedkkgenerally corresponds to different effective support tolerances across the manifold\. This distinction is also relevant to bandwidth selection\. The neighborhood size required to reproduce the dense operator at toleranceτ\\taudepends onβ\\beta, while the preferred bandwidthβ⋆\\beta^\{\\star\}is determined by the spectral calculation\.

A fixed neighborhood graph therefore selects a truncation scale before the bandwidth of interest is known\.

FlashDiffusion evaluates the Gaussian operator directly throughout theβ\\betacontinuation, so the same pairwise interaction model is used during bandwidth selection and eigenspace computation\.

### A\.2Two scales: spectral resolution and kernel support

Two distinct length scales govern the comparison\. The first is the Gaussian smoothing scaleℓβ∼β−1/2\\ell\_\{\\beta\}\\sim\\beta^\{\-1/2\}\. By Weyl counting, the number of spatial modes supported at this scale behaves asM⁡\(β\)∼V​βd/2M\(\\beta\)\\sim V\\beta^\{d/2\}up to geometric constants\. ResolvingMMmodes therefore corresponds to the scaling:

β≳\(MV\)2/d\.\\displaystyle\\beta\\gtrsim\\left\(\\frac\{M\}\{V\}\\right\)^\{2/d\}\.\(17\)The second is the support scaleDτD\_\{\\tau\}in eq\.[16](https://arxiv.org/html/2609.38198#A1.E16)\. Matching the dense Gaussian kernel above the working toleranceτ\\taurequires approximately:

k∗​\(β,τ\)=N​ωdV​\(log⁡\(1/τ\)β\)d/2\.\\displaystyle k^\{\\ast\}\(\\beta,\\tau\)=\\frac\{N\\omega\_\{d\}\}\{V\}\\left\(\\frac\{\\log\(1/\\tau\)\}\{\\beta\}\\right\)^\{d/2\}\.\(18\)The quantityk∗k^\{\\ast\}is therefore a support\-fidelity scale: it is the expected number of kernel interactions per row whose weights exceed the chosen tolerance\. For the Swiss roll atN=106N=10^\{6\},β⋆=15\\beta^\{\\star\}=15, andd=2d=2, this estimate givesk∗≈8×104k^\{\\ast\}\\approx 8\\times 10^\{4\}forτ=10−7\\tau=10^\{\-7\}andk∗≈1\.7×105k^\{\\ast\}\\approx 1\.7\\times 10^\{5\}forτ=10−15\\tau=10^\{\-15\}\. The smallerfp64\-scale tolerance includes a wider range of weak interactions and therefore increases the support count\. These values are much larger than neighborhood sizes such ask∈32,64,128k\\in\{32,64,128\}commonly used to construct sparse embedding graphs\. The comparison establishes that a small fixed\-kkgraph is a strong truncation of the finite\-sample Gaussian operator atβ⋆\\beta^\{\\star\}\. Preservation of a particular leading eigenspace can still occur under stronger truncation because spectral accuracy also depends on the structure of the omitted interactions and on the relevant eigengaps\. Along a calibrated bandwidth path, two reference exponents are useful\. Balancing bias against pointwise fluctuation gives Singer’s scalingβ⋆∼N2/\(d\+6\)\\beta^\{\\star\}\\sim N^\{2/\(d\+6\)\}\([Singer, 2006](https://arxiv.org/html/2609.38198#bib.bib4)\), while balancing against the spectral fluctuation\(Nεd/2\)−1/2\(N\\varepsilon^\{d/2\}\)^\{\-1/2\}givesβ⋆∼N2/\(d\+4\)\\beta^\{\\star\}\\sim N^\{2/\(d\+4\)\}, the corresponding rate for a spectral selection criterion\([Calder and García Trillos, 2022](https://arxiv.org/html/2609.38198#bib.bib35)\)\. On the Swiss roll, the measured bandwidth selection followsβ⋆∼N1/3\\beta^\{\\star\}\\sim N^\{1/3\}ford=2d=2over the tested range\. Substitution into eq\.[18](https://arxiv.org/html/2609.38198#A1.E18)givesk∗∼N4/\(d\+4\)k^\{\\ast\}\\sim N^\{4/\(d\+4\)\}, ork∗∼N2/3k^\{\\ast\}\\sim N^\{2/3\}ford=2d=2, consistent with the measured support counts\. Combining eq\.[17](https://arxiv.org/html/2609.38198#A1.E17)with eq\.[18](https://arxiv.org/html/2609.38198#A1.E18)gives a useful scaling relation\. At the bandwidth associated withMMresolved modes, preserving the dense\-kernel support above toleranceτ\\taurequires:

k≳N​log⁡\(1/τ\)d/2M\.\\displaystyle k\\gtrsim\\frac\{N\\log\(1/\\tau\)^\{d/2\}\}\{M\}\.\(19\)Order\-one geometric constants are omitted in eq\.[19](https://arxiv.org/html/2609.38198#A1.E19)\. Ford=2d=2andN=106N=10^\{6\}, this support\-matching scale is approximately10510^\{5\}neighbors forM=100M=100and10410^\{4\}neighbors forM=1000M=1000at the working tolerances considered here\. A fixed graph with tens or hundreds of neighbors therefore corresponds to a substantially stronger truncation than retaining the numerically relevant support of the dense kernel\. The same relation can be viewed at fixedkk\. Ford=2d=2, the support count scales asβ−1\\beta^\{\-1\}\. Using the measured Swiss\-roll reference point,k∗≈8×104k^\{\\ast\}\\approx 8\\times 10^\{4\}atβ⋆=15\\beta^\{\\star\}=15, a graph withk=64k=64reaches a comparable support radius around:

βsupport≈15​\(8×10464\)≈2×104\.\\displaystyle\\beta\_\{\\mathrm\{support\}\}\\approx 15\\left\(\\frac\{8\\times 10^\{4\}\}\{64\}\\right\)\\approx 2\\times 10^\{4\}\.\(20\)This value is roughly three orders of magnitude larger than the selectedβ⋆\\beta^\{\\star\}\. Ak=64k=64graph therefore corresponds to a much shorter geometric interaction scale than the dense Gaussian operator selected by the spectral criterion\.

### A\.3Numerically relevant support and retained row mass

The support calculation above follows directly from the small\-distance distribution of points on a smooth manifold\. Defineρτ\(β\)=Pr\[𝖪i​j\(β\)≥τ\]\\rho\_\{\\tau\}\(\\beta\)=\\Pr\[\\mathsf\{K\}\_\{ij\}\(\\beta\)\\geq\\tau\]\. Equivalently,ρτ\(β\)=Pr\[d2≤log\(1/τ\)/β\]\\rho\_\{\\tau\}\(\\beta\)=\\Pr\[d^\{2\}\\leq\\log\(1/\\tau\)/\\beta\]\. IfFFdenotes the CDF of pairwise squared distances, thenρτ​\(β\)=F⁡\(log⁡\(1/τ\)/β\)\\rho\_\{\\tau\}\(\\beta\)=F\(\\log\(1/\\tau\)/\\beta\)\. While the bulk ofFFdepends on global manifold geometry, its small\-distance behavior is local\. For a smoothdd\-manifold of volumeVVunder approximately uniform sampling, the interior small\-ball law givesPr\[d≤R\]=\(ωdRd/V\)\(1\+𝒪\(R2\)\)\\Pr\[d\\leq R\]=\(\\omega\_\{d\}R^\{d\}/V\)\(1\+\\mathcal\{O\}\(R^\{2\}\)\)\. Hence

ρτ​\(β\)=ωdV​\(log⁡\(1/τ\)β\)d/2​\(1\+𝒪⁡\(Dτ2\)\)\.\\displaystyle\\rho\_\{\\tau\}\(\\beta\)=\\frac\{\\omega\_\{d\}\}\{V\}\\left\(\\frac\{\\log\(1/\\tau\)\}\{\\beta\}\\right\)^\{d/2\}\\left\(1\+\\mathcal\{O\}\(D\_\{\\tau\}^\{2\}\)\\right\)\.\(21\)The expected support count per row isnτ​\(β\)=N​ρτ​\(β\)n\_\{\\tau\}\(\\beta\)=N\\rho\_\{\\tau\}\(\\beta\), and eq\.[18](https://arxiv.org/html/2609.38198#A1.E18)follows directly from this relation\. Along the spectral scalingβ⋆∼N2/\(d\+4\)\\beta^\{\\star\}\\sim N^\{2/\(d\+4\)\}, the relative support density decreases with sample size while the absolute number of relevant interactions grows asnτ​\(β⋆\)∼N4/\(d\+4\)n\_\{\\tau\}\(\\beta^\{\\star\}\)\\sim N^\{4/\(d\+4\)\}\. The Gaussian operator therefore becomes increasingly localized in relative density while the number of interactions per row associated with a fixed relative tolerance continues to increase\. We verified thed=2d=2scaling on a unit\-variance Swiss roll usingτ=10−7\\tau=10^\{\-7\}\. The prediction isρτ​\(β\)∝β−1\\rho\_\{\\tau\}\(\\beta\)\\propto\\beta^\{\-1\}\. AcrossN=104N=10^\{4\}–10810^\{8\}, the measured productρτ​\(β⋆\)​β⋆\\rho\_\{\\tau\}\(\\beta^\{\\star\}\)\\beta^\{\\star\}remains approximately constant at1\.191\.19, within1%1\\%ofω2​log⁡\(1/τ\)/V\\omega\_\{2\}\\log\(1/\\tau\)/V\. Over the same range, the selected bandwidth followsβ⋆∼N1/3\\beta^\{\\star\}\\sim N^\{1/3\}and the measured support count increases from approximately3\.7×1033\.7\\times 10^\{3\}to1\.7×1061\.7\\times 10^\{6\}interactions per row, consistent withN2/3N^\{2/3\}scaling, while the relative support density decreases from37%37\\%to1\.7%1\.7\\%\. Curvature corrections enter at relative order𝒪⁡\(Dτ2\)\\mathcal\{O\}\(D\_\{\\tau\}^\{2\}\)in the interior, while boundary effects occur within a strip of width𝒪⁡\(Dτ\)\\mathcal\{O\}\(D\_\{\\tau\}\)\. The observed agreement over the tested range indicates that the leading small\-ball scaling captures the dominant dependence relevant to these support estimates\. Support counts depend on the chosen toleranceτ\\tau\. A complementary quantity with no dependence onτ\\tauis the fraction of Gaussian row mass retained by thekkNN mask\. Because the Gaussian weights are nonnegative, define

fk=∑j∈𝒩k​\(i\)exp⁡\(−β​di​j2\)∑jexp⁡\(−β​di​j2\)\.\\displaystyle f\_\{k\}=\\frac\{\\sum\_\{j\\in\\mathcal\{N\}\_\{k\}\(i\)\}\\exp\(\-\\beta d\_\{ij\}^\{2\}\)\}\{\\sum\_\{j\}\\exp\(\-\\beta d\_\{ij\}^\{2\}\)\}\.\(22\)Under the local uniform\-density approximation,

fk=γ⁡\(d/2,β​Rk2\)Γ⁡\(d/2\)\.\\displaystyle f\_\{k\}=\\frac\{\\gamma\(d/2,\\beta R\_\{k\}^\{2\}\)\}\{\\Gamma\(d/2\)\}\.\(23\)SinceRk∼\(k/\(N​q​ωd\)\)1/dR\_\{k\}\\sim\(k/\(Nq\\omega\_\{d\}\)\)^\{1/d\}, the argument scales asβ​Rk2∝β​\(k/N\)2/d\\beta R\_\{k\}^\{2\}\\propto\\beta\(k/N\)^\{2/d\}\. Alongβ⋆∼N2/\(d\+4\)\\beta^\{\\star\}\\sim N^\{2/\(d\+4\)\}, fixedkktherefore givesfk∼kN−4/\(d\+4\)f\_\{k\}\\sim kN^\{\-4/\(d\+4\)\}in the small\-argument regime\. Ford=2d=2, eq\.[23](https://arxiv.org/html/2609.38198#A1.E23)reduces locally tofk=1−e−c2​kf\_\{k\}=1\-e^\{\-c\_\{2\}k\}, wherec2c\_\{2\}contains the sampling\-density and bandwidth factors\. This relation is geometric and independent of the working toleranceτ\\tau\. Table[1](https://arxiv.org/html/2609.38198#A1.T1)verifies the prediction on the Swiss roll\. A fittedc2≈6\.9×10−4c\_\{2\}\\approx 6\.9\\times 10^\{\-4\}reproduces the measured retained fractions to within approximately0\.3%0\.3\\%over the tested range\.

Table 1:Gaussian kernel row mass retained under fixed\-kktruncation on a unit\-variance Swiss roll withN=2×105N=2\\times 10^\{5\}andβ⋆=8\.75\\beta^\{\\star\}=8\.75, averaged over 400 rows\. The locald=2d=2prediction1−e−c2​k1\-e^\{\-c\_\{2\}k\}usesc2≈6\.9×10−4c\_\{2\}\\approx 6\.9\\times 10^\{\-4\}\.At this problem size,k=64k=64retains about4\.4%4\.4\\%of the Gaussian row mass atβ⋆\\beta^\{\\star\}\. Spectral error depends additionally on the organization of the omitted interactions and on the relevant spectral gaps, so the retained\-mass fraction and eigenspace error need not coincide numerically\. The small retained fraction nevertheless shows that fixed\-kkmasking at this scale produces a substantial perturbation of the dense operator\. This distinction is important for bandwidth selection\. FlashDiffusion seeks the largest bandwidth for which a requested spectral subspace remains distinguishable from sampling and numerical fluctuations\. Replacing𝖪β\\mathsf\{K\}\_\{\\beta\}by𝖪~β,k\\widetilde\{\\mathsf\{K\}\}\_\{\\beta,k\}changes the operator, its degree normalization, and its spectral gaps\. The sparsified family therefore has its own finite\-sample resolution boundary, which can occur at a differentβ\\betafrom the corresponding boundary of the dense Gaussian family\.

### A\.4kkNN truncation and PSD

The dense Gaussian Gram matrix satisfies𝖪β⪰0\\mathsf\{K\}\_\{\\beta\}\\succeq 0\. For a Hadamard\-truncated kernel𝖪~β,k=𝖪β⊙Ak\\widetilde\{\\mathsf\{K\}\}\_\{\\beta,k\}=\\mathsf\{K\}\_\{\\beta\}\\odot A\_\{k\}, the Schur product theorem guarantees positive semidefiniteness whenAk⪰0A\_\{k\}\\succeq 0\. A symmetrizedkkNN adjacency mask generally does not satisfy this condition\. Once indefiniteness is introduced, positive diagonal congruence preserves its inertia\. For any positive diagonal matrixDD, Sylvester’s law of inertia gives

inertia⁡\(D⁡\(𝖪β⊙Ak\)​D\)=inertia⁡\(𝖪β⊙Ak\)\.\\displaystyle\\mathrm\{inertia\}\\left\(D\(\\mathsf\{K\}\_\{\\beta\}\\odot A\_\{k\}\)D\\right\)=\\mathrm\{inertia\}\\left\(\\mathsf\{K\}\_\{\\beta\}\\odot A\_\{k\}\\right\)\.\(24\)The symmetric positive\-diagonal rescalings appearing in diffusion\-map normalization therefore preserve any negative eigenvalues introduced by the mask\. Table[2](https://arxiv.org/html/2609.38198#A1.T2)illustrates this effect on the Swiss roll at a problem size for which the complete spectrum can be computed directly\.

Table 2:Spectrum of a symmetricallykkNN\-truncated Gaussian kernel on the Swiss roll withN=3000N=3000,β⋆=2\.158\\beta^\{\\star\}=2\.158, and neighborhood masks symmetrized by union\. Negative spectral mass is∑λ<0\|λ\|/∑\|λ\|\\sum\_\{\\lambda<0\}\|\\lambda\|/\\sum\|\\lambda\|\. The final column reports the smallest eigenvalue of the mask itself\.The table demonstrates directly that hard symmetric neighborhood masking can alter a structural property that the dense Gaussian kernel possesses exactly\. Its purpose is finite\-sample characterization at a size where the complete spectrum is available\. PSD also has a natural diffusion interpretation\. For a positive eigenvalue one may writeλm=e−t​νm\\lambda\_\{m\}=e^\{\-t\\nu\_\{m\}\}withνm≥0\\nu\_\{m\}\\geq 0and interpretνm\\nu\_\{m\}as a real decay rate\. Negative eigenvalues introduced by masking correspond to complex values ofνm\\nu\_\{m\}under this representation\. Portions of the truncated spectrum can therefore lose the real diffusion\-semigroup interpretation carried by the dense Gaussian kernel\.

### A\.5Continuum geometry and ANN perturbations

Furthermore,kkNN methods\([Hein et al\., 2007](https://arxiv.org/html/2609.38198#bib.bib5);[Ting et al\., 2010](https://arxiv.org/html/2609.38198#bib.bib9);[Calder and García Trillos, 2022](https://arxiv.org/html/2609.38198#bib.bib35)\), using exact distances are still a𝒪⁡\(N2\)\\mathcal\{O\}\(N^\{2\}\)theory\. Approximate nearest\-neighbor \(ANN\) methods introduce an additional perturbation\. IfA^k\\widehat\{A\}\_\{k\}is the returned approximate mask, then the effective kernel is:

𝖪^β,k=𝖪β⊙A^k\.\\displaystyle\\widehat\{\\mathsf\{K\}\}\_\{\\beta,k\}=\\mathsf\{K\}\_\{\\beta\}\\odot\\widehat\{A\}\_\{k\}\.\(25\)Missed or substituted edges modify the weighted degrees, theα\\alpha\-normalization, and ultimately the spectrum\. Standard ANN recall\([Malkov and Yashunin, 2020](https://arxiv.org/html/2609.38198#bib.bib31)\)therefore gives only a partial measure of spectral accuracy because the effect of an incorrect edge depends on its Gaussian weight and on the spectral gaps of the modes being estimated\([Boutsidis et al\., 2009](https://arxiv.org/html/2609.38198#bib.bib36);[Drineas et al\., 2008](https://arxiv.org/html/2609.38198#bib.bib37)\)\. Some ANN implementations can additionally introduce run\-to\-run variability, depending on their construction and execution settings\. The broader issue is the additional finite\-sample perturbation introduced by approximate neighborhood construction, whose accuracy must be chosen before or adapted jointly with the bandwidth search\. A neighborhood budget large enough to reproduce the dense\-kernel support at the as\-yet unknownβ⋆\\beta^\{\\star\}must either be selected conservatively in advance or increased as the bandwidth changes\.

### A\.6Synthesis

The comparison above separates several forms of approximation that affect diffusion\-map calculations differently\. A visually stable low\-dimensional embedding can coexist with substantial changes in a larger spectral basis because the first few eigenvectors encode coarse geometry while higher eigenfunctions probe progressively shorter spatial scales\. Matching the dense Gaussian kernel above a relative toleranceτ\\taurequires a neighborhood size scaling as:

k∗​\(β,τ\)∼N​\(log⁡\(1/τ\)β\)d/2\.\\displaystyle k^\{\\ast\}\(\\beta,\\tau\)\\sim N\\left\(\\frac\{\\log\(1/\\tau\)\}\{\\beta\}\\right\)^\{d/2\}\.\(26\)At a bandwidth associated withMMresolved modes, this becomes:

k∗∼N​log⁡\(1/τ\)d/2M\.\\displaystyle k^\{\\ast\}\\sim\\frac\{N\\log\(1/\\tau\)^\{d/2\}\}\{M\}\.\(27\)These relations characterize support fidelity of the dense finite\-sample Gaussian operator\. Spectral fidelity can persist under stronger truncation when the omitted interactions have limited influence on the target eigenspace, making direct eigenspace comparisons an additional empirical diagnostic\. The row\-mass calculation provides a tolerance\-independent measure of truncation strength\. In the Swiss\-roll experiment,k=64k=64retains4\.4%4\.4\\%of the Gaussian row mass atβ⋆\\beta^\{\\star\}, whilek=1024k=1024retains approximately half\. Hard neighborhood masking can also remove the PSD structure of the Gaussian kernel and alter the deeper portion of its diffusion spectrum\.

These observations complement the main\-text analysis of low\-rank approximations\. At bandwidths relevant to resolving many diffusion modes, the Gaussian matrix can simultaneously have substantial effective rank and substantial spatial support\. Low\-rank compression restricts the number of retained spectral directions, while aggressive fixed\-kksparsification removes a large fraction of spatial interactions\. Both operations therefore modify the finite\-sample Gaussian operator in distinct ways\.

## Appendix BKeOpsComparison

Here we demonstrate a comparison of theKeOpslibrary666Usingpip install pykeops==2\.3\(Released: Apr 17, 2025\),[Feydy et al\. \(2020\)](https://arxiv.org/html/2609.38198#bib.bib43);[Charlier et al\. \(2021\)](https://arxiv.org/html/2609.38198#bib.bib41), versus our matrix\-vector \(matvec\) and matrix\-matrix \(matmat\) primitiveXGEMM\. The results are shown in fig\.[4](https://arxiv.org/html/2609.38198#A2.F4), and both achieve similar timing results, and agree tomachine\-precision\(fp32agreement is within∼10−6\\sim 10^\{\-6\}\)\.

2102^\{10\}2142^\{14\}2182^\{18\}2222^\{22\}2242^\{24\}10−310^\{\-3\}10−110^\{\-1\}10110^\{1\}10310^\{3\}NNwall time \(s\)\(a\) matvec \(b=1b=1\)XGEMMKeOps∝N2\\propto N^\{2\}2102^\{10\}2142^\{14\}2182^\{18\}2222^\{22\}2242^\{24\}NN\(b\) matmat \(b=8b=8\)XGEMMKeOps∝N2\\propto N^\{2\}Figure 4:Raw matrix\-free Gaussian\-kernel scaling on the centered anisotropic flat torusT6T^\{6\}\(D=32D=32,β=1\\beta=1\) using one RTX–PRO–6000 GPU\.\(a\)True matrix–vector applicationK​vKv\(b=1b=1\)\.\(b\)Eight\-right\-hand\-side matrix–matrix applicationK​VKV\(b=8b=8\)\.XGEMM–TILT and KeOps evaluate the same dense Gaussian interaction without materializing theN×NN\\times Nkernel matrix\. Points show mean wall time and error bars show one standard deviation over repeated timed applications; the dashed guides indicate exact𝒪⁡\(N2\)\\mathcal\{O\}\(N^\{2\}\)scaling in the large\-NNregime\.