Samplets: Multiscale Analysis for Scattered Data

• Samplets are localized, multiscale signed measures constructed from Dirac point evaluations on scattered data, organized via hierarchical clusters with vanishing moments.
• They enable fast transforms, data and operator compression, and sparse approximation by translating classical wavelet methods to irregular grids.
• Extensions to Banach spaces, RKHS embeddings, and graph signals demonstrate samplets' versatility for adaptive approximation and scalable numerical analysis.

Samplets are localized, multiscale signed measures adapted to scattered data sites. In the classical setting they are constructed from Dirac point evaluations on an unstructured point cloud, organized through a multiresolution analysis on a hierarchical cluster tree, and endowed with vanishing moments against low-degree polynomials. These properties make samplets the scattered-data analogue of wavelets and permit the transfer of fast transforms, data compression, operator compression, and sparse approximation from regular-grid wavelet analysis to irregular data geometries (Harbrecht et al., 2021, Harbrecht et al., 20 Mar 2025). Later work generalized the construction to Banach-space functionals, graph signals, and infinite-data limits leading to multiwavelets, while preserving the central themes of locality, orthogonality or Parseval stability, and moment cancellation (Balazs et al., 2024, Elefante et al., 25 Jul 2025, Giacchi et al., 2 Apr 2026).

1. Definition and conceptual position

In the finite scattered-data setting, one starts from a point set $$X=\{x_1,\dots,x_N\}\subset \Omega\subset \mathbb R^d$$ and the space

$$\mathcal X=\operatorname{span}\{\delta_{x_1},\dots,\delta_{x_N}\},$$

where $\delta_{x_i}$ denotes point evaluation. A samplet is a wavelet-type element of this space, hence a localized signed linear combination of point evaluations,

$$\sigma_{j,k}=\sum_{i=1}^N \omega_{j,k,i}\,\delta_{x_i},$$

with coefficients chosen so that low-degree polynomial moments vanish: $$(p,\sigma_{j,k})_\Omega=0 \qquad \forall p\in\mathcal P_q(\Omega),$$ where $\mathcal P_q(\Omega)$ is the space of polynomials of degree at most $q$ (Harbrecht et al., 2021, Baroli et al., 2023).

This definition encodes the standard wavelet mechanism in a non-grid setting. Samplets are localized because each basis element is supported on a small cluster of sites, multiscale because they arise from a hierarchy of nested spaces, and moment-canceling because they annihilate smooth polynomial backgrounds. The original motivation was to transfer the Tausch–White construction from function spaces on geometries to raw point clouds and sampled data, thereby obtaining a multilevel representation for compression, singularity detection, adaptivity, and kernel-matrix sparsification on arbitrary scattered sites (Harbrecht et al., 2021, Harbrecht et al., 20 Mar 2025).

A recurrent theme in the literature is that samplets should not be identified with continuous basis functions on a regular domain. In their classical form they are discrete signed measures. Only after an RKHS embedding do they become functions such as linear combinations of kernel translates. This distinction is central to later developments in kernel approximation, Gaussian processes, and generalized Banach-space constructions (Harbrecht et al., 2024, Neugebauer, 2024).

2. Finite-data multiresolution construction

The classical construction is based on a multiresolution hierarchy

$$\mathcal X_0\subset \mathcal X_1\subset \cdots \subset \mathcal X_J=\mathcal X, \qquad \mathcal X_{j+1}=\mathcal X_j\overset{\perp}{\oplus}\mathcal S_j,$$

where $\mathcal X_j$ are scaling spaces and $\mathcal S_j$ are detail spaces. The hierarchy is induced by a cluster tree whose root is $$\mathcal X=\operatorname{span}\{\delta_{x_1},\dots,\delta_{x_N}\},$$0, whose non-leaf clusters are disjoint unions of their children, and which is typically taken to be balanced binary in the main complexity statements. At the leaves, the scaling distributions are simply Dirac masses; upward through the tree, finer scaling distributions are transformed into coarser scaling distributions and detail samplets (Harbrecht et al., 2021, Avesani et al., 2024).

The local two-scale relation has the form

$$\mathcal X=\operatorname{span}\{\delta_{x_1},\dots,\delta_{x_N}\},$$1

where $$\mathcal X=\operatorname{span}\{\delta_{x_1},\dots,\delta_{x_N}\},$$2 and $$\mathcal X=\operatorname{span}\{\delta_{x_1},\dots,\delta_{x_N}\},$$3 denote the scaling distributions and samplets associated with a cluster $$\mathcal X=\operatorname{span}\{\delta_{x_1},\dots,\delta_{x_N}\},$$4. Vanishing moments are imposed through a moment matrix built from monomials up to degree $$\mathcal X=\operatorname{span}\{\delta_{x_1},\dots,\delta_{x_N}\},$$5, followed by a QR factorization of its transpose,

$$\mathcal X=\operatorname{span}\{\delta_{x_1},\dots,\delta_{x_N}\},$$6

Because $$\mathcal X=\operatorname{span}\{\delta_{x_1},\dots,\delta_{x_N}\},$$7 is lower triangular, one designates the first $$\mathcal X=\operatorname{span}\{\delta_{x_1},\dots,\delta_{x_N}\},$$8 transformed distributions as scaling distributions and the remainder as samplets, thereby obtaining at least $$\mathcal X=\operatorname{span}\{\delta_{x_1},\dots,\delta_{x_N}\},$$9 vanishing moments (Harbrecht et al., 2021, Harbrecht et al., 20 Mar 2025).

Globally, the construction yields an orthogonal transform matrix $\delta_{x_i}$0 satisfying

$\delta_{x_i}$1

This fast samplet transform maps point-value coordinates to samplet coordinates in linear cost on balanced trees, and its inverse is obtained by transposing the local transforms. The basis construction itself is likewise linear in $\delta_{x_i}$2 under the standard balanced-tree assumptions (Harbrecht et al., 2021, Avesani et al., 2024).

3. Structural properties and coefficient decay

The fundamental basis theorem states that the samplet system is an orthonormal basis of the discrete measure space, with scale-wise dimensions growing dyadically. In representative formulations one has

$\delta_{x_i}$3

together with support localization to single clusters and vanishing moments of order $\delta_{x_i}$4 (Baroli et al., 2023, Avesani et al., 2024).

Localization is quantitative under quasi-uniformity. For quasi-uniform point sets, samplet supports satisfy bounds such as

$\delta_{x_i}$5

and coefficient vectors satisfy $\delta_{x_i}$6-bounds of the form

$\delta_{x_i}$7

These estimates are the discrete analogues of wavelet support and norm control, and they underwrite both data compression and operator compression (Avesani et al., 2024, Neugebauer, 2024).

For smooth data, samplet coefficients decay because vanishing moments remove the low-order Taylor part. A representative estimate is

$\delta_{x_i}$8

for $\delta_{x_i}$9 on a neighborhood of the support. This is the principal reason that smooth regions yield negligible fine-scale coefficients, while nonsmooth regions retain large coefficients and become detectable through the multiscale expansion (Avesani et al., 2024, Baroli et al., 2023).

Later adaptive work makes this interpretation explicit: if $$\sigma_{j,k}=\sum_{i=1}^N \omega_{j,k,i}\,\delta_{x_i},$$0 and the moment order is high enough, then coefficients obey a local estimate of the form

$$\sigma_{j,k}=\sum_{i=1}^N \omega_{j,k,i}\,\delta_{x_i},$$1

so fine-scale coefficients decay rapidly in smooth regions and persist near singularities. This suggests a natural use of samplets as local smoothness indicators for adaptive tree refinement and adaptive center selection (Avesani et al., 2 Apr 2026).

4. RKHS embedding, kernel compression, and multiscale interpolation

A major development is the embedding of samplets into reproducing kernel Hilbert spaces via the Riesz map. If $$\sigma_{j,k}=\sum_{i=1}^N \omega_{j,k,i}\,\delta_{x_i},$$2 is a reproducing kernel and $$\sigma_{j,k}=\sum_{i=1}^N \omega_{j,k,i}\,\delta_{x_i},$$3, then a samplet

$$\sigma_{j,k}=\sum_{i=1}^N \omega_{j,k,i}\,\delta_{x_i},$$4

is sent to the embedded samplet

$$\sigma_{j,k}=\sum_{i=1}^N \omega_{j,k,i}\,\delta_{x_i},$$5

This produces a multiresolution basis of the same finite-dimensional kernel space $$\sigma_{j,k}=\sum_{i=1}^N \omega_{j,k,i}\,\delta_{x_i},$$6, rather than a different approximation space. In this sense samplets are an orthogonal change of coordinates on the discrete data functionals, and after embedding they become a multiscale basis of kernel translates (Baroli et al., 2023, Harbrecht et al., 2024).

In samplet coordinates, the kernel matrix transforms as

$$\sigma_{j,k}=\sum_{i=1}^N \omega_{j,k,i}\,\delta_{x_i},$$7

For asymptotically smooth kernels, off-diagonal entries in $$\sigma_{j,k}=\sum_{i=1}^N \omega_{j,k,i}\,\delta_{x_i},$$8 are small when the supporting clusters are well separated, because vanishing moments cancel low-order polynomial approximations of the kernel. Under the admissibility criterion

$$\sigma_{j,k}=\sum_{i=1}^N \omega_{j,k,i}\,\delta_{x_i},$$9

the transformed matrix can be compressed with controlled Frobenius error, and the compressed matrix has

$$(p,\sigma_{j,k})_\Omega=0 \qquad \forall p\in\mathcal P_q(\Omega),$$0

nonzero entries for quasi-uniform data (Avesani et al., 2024, Harbrecht et al., 2021).

This compression mechanism has several algorithmic consequences. The original samplet paper showed that dense kernel matrices arising in kernel learning, scattered-data approximation, and Gaussian process regression become quasi-sparse in samplet coordinates and can be thresholded to about $$(p,\sigma_{j,k})_\Omega=0 \qquad \forall p\in\mathcal P_q(\Omega),$$1 relevant entries (Harbrecht et al., 2021). In the dual-pair formulation of finite-dimensional RKHS bases, samplets yield a multiresolution orthogonal coordinate system in which both the kernel matrix and its inverse are compressible, and orthogonalization by $$(p,\sigma_{j,k})_\Omega=0 \qquad \forall p\in\mathcal P_q(\Omega),$$2 remains practically local in samplet coordinates (Harbrecht et al., 2024).

For multiscale scattered-data interpolation with Matérn kernels, samplets play a complementary role to levelwise kernel scaling. The approximation is built by residual correction in spaces

$$(p,\sigma_{j,k})_\Omega=0 \qquad \forall p\in\mathcal P_q(\Omega),$$3

where $$(p,\sigma_{j,k})_\Omega=0 \qquad \forall p\in\mathcal P_q(\Omega),$$4 is the fill distance. In that framework, the diagonal blocks $$(p,\sigma_{j,k})_\Omega=0 \qquad \forall p\in\mathcal P_q(\Omega),$$5 remain uniformly well conditioned across levels, while samplet transforms $$(p,\sigma_{j,k})_\Omega=0 \qquad \forall p\in\mathcal P_q(\Omega),$$6 make the multiscale block system compressible. The paper is explicit that the samplet transformation preserves condition numbers,

$$(p,\sigma_{j,k})_\Omega=0 \qquad \forall p\in\mathcal P_q(\Omega),$$7

so samplets do not alleviate ill-conditioning by themselves; favorable conditioning comes from kernel scaling, whereas sparsity and quasi-linear cost come from the samplet basis change (Avesani et al., 2024).

A related misconception concerns interpolation. Samplet-based approximations are generally not Lagrange bases, and the samplet literature explicitly notes that samplet-based approximations “will most likely not satisfy any Lagrange condition.” Their role is algebraic and multiresolution, not nodal-cardinal (Harbrecht et al., 2024).

5. Generalizations and theoretical extensions

The finite-dimensional point-evaluation construction has been generalized in several directions. In Banach spaces, the basic atoms are no longer Dirac masses but arbitrary continuous functionals $$(p,\sigma_{j,k})_\Omega=0 \qquad \forall p\in\mathcal P_q(\Omega),$$8, assumed to form either an $$(p,\sigma_{j,k})_\Omega=0 \qquad \forall p\in\mathcal P_q(\Omega),$$9-frame or a Riesz basis with square-summable coefficients. The construction proceeds on the coefficient side through an isometry on the range of the analysis operator, while localization is induced by supports in $\mathcal P_q(\Omega)$0 and a hierarchy generated by spectral clustering of a similarity graph. Vanishing moments are imposed not only against polynomials but against a user-chosen primitive space $\mathcal P_q(\Omega)$1, and the resulting generalized samplets satisfy an abstract localization estimate

$\mathcal P_q(\Omega)$2

This reformulates samplets as stable multiscale transformations of functionals rather than of point evaluations alone (Balazs et al., 2024).

A second extension concerns the infinite-data limit. In a probabilistic framework with i.i.d. sample sites and polynomial primitives, finite samplets converge weakly to signed measures with broken polynomial densities. These limiting objects form orthonormal multiwavelets associated with a hierarchical partition of the domain, and for congruent partitions one recovers classical multiwavelets with scale- and partition-independent filter coefficients. The construction also extends from total-degree spaces to general downward closed multi-index sets $\mathcal P_q(\Omega)$3, permitting anisotropic moment prescriptions beyond tensor-product settings (Giacchi et al., 2 Apr 2026).

A third extension treats graph signals. There the graph is partitioned into patches, each patch is embedded into a Euclidean space, samplets are constructed in the embedded coordinates, and the basis is pulled back to the graph. This produces a samplet forest: a block-diagonal orthogonal transform with vanishing moments against patchwise polynomial spaces induced from local coordinates. The resulting graph samplets generalize graph Haar constructions from piecewise constants to higher-order local polynomial annihilation, and the paper identifies as compressible those graph signals that are locally approximable by generalized polynomials in the patch coordinates (Elefante et al., 25 Jul 2025).

6. Applications, empirical behavior, and limitations

The earliest applications emphasized direct data compression and feature detection. On one-dimensional signals and images, coefficient thresholding in the orthonormal samplet basis produced strong compression rates: the 2021 paper reports, for example, $\mathcal P_q(\Omega)$4, $\mathcal P_q(\Omega)$5, and $\mathcal P_q(\Omega)$6 compression on a $\mathcal P_q(\Omega)$7 grayscale image, while dominant coefficients localize image features and edges (Harbrecht et al., 2021). The same paper uses samplet-compressed kernel matrices for sparse Cholesky factorizations and Gaussian random field simulation, thereby linking multiresolution data analysis to scalable kernel linear algebra (Harbrecht et al., 2021).

In scattered-data approximation, samplet coordinates support both interpolation and sparsity-constrained approximation. Samplet basis pursuit formulates an $\mathcal P_q(\Omega)$8-regularized least-squares problem in samplet coordinates rather than in the basis of kernel translates, motivated by the claim that the class of signals sparse in the embedded samplet basis is considerably larger than the class sparse in kernel translates. Large-scale experiments on $\mathcal P_q(\Omega)$9 random two-dimensional points showed that a multiresolution semi-smooth Newton solver produced much sparser solutions than ridge regression or FISTA-based baselines, and applications included surface reconstruction and multi-kernel reconstruction of scattered temperature data (Baroli et al., 2023).

Later work extends this logic to adaptive center selection and multikernel sparse regression. In tree-adaptive multiscale kernel lasso, samplet coefficients are used both to compress transformed kernel matrices and to localize important regions of the data through clusterwise energy propagation, after which one representative data site is kept per significant leaf cluster. The reported reductions are substantial: in one $q$0 two-dimensional benchmark, adaptive samplet selection reduced the center set to $q$1, and in a three-dimensional Stanford Bunny reflectance problem to $q$2, after which $q$3-regularized multikernel models were solved by a trust-region semismooth Newton method with online low-rank SVD stabilization (Avesani et al., 2 Apr 2026).

Samplets have also been used as a numerical acceleration layer for Gaussian processes. In that setting the GP prior and posterior formulas remain unchanged, but the dense system

$q$4

is replaced by a transformed sparse approximation

$q$5

The thesis on SampletsGP argues that this can reduce the cubic computational bottleneck to log-linear scale in the best case, while still retaining all observations rather than replacing them by inducing points. The strongest claims are explicitly restricted to low-dimensional problems and to kernels, such as half-integer Matérn kernels, that satisfy the required asymptotic smoothness assumptions (Neugebauer, 2024).

The principal limitations are consistent across the literature. First, the polynomial moment space has dimension

$q$6

so the method is not dimension-robust and is mainly suited to moderate or low dimension (Harbrecht et al., 2021, Neugebauer, 2024). Second, the strongest localization and compression guarantees are typically proved for quasi-uniform point sets, balanced binary trees, and bounded or convex domains, while several numerical studies deliberately go beyond that regime (Avesani et al., 2024, Neugebauer, 2024). Third, matrix sparsity is only quasi-sparsity before truncation, and the practical efficiency of sparse factorizations depends on fill-in and ordering (Harbrecht et al., 2021, Neugebauer, 2024). Finally, the construction and its associated compressed-matrix machinery are considerably more involved than single-scale kernel or localized Lagrange approaches, even though the asymptotic compression can be superior on large problems (Harbrecht et al., 2024, Harbrecht et al., 20 Mar 2025).

Taken together, the literature presents samplets as a general multiresolution technology for irregular data. Their defining structure is stable hierarchical localization combined with prescribed moment cancellation; their principal consequences are sparse representations of data and operators; and their later generalizations show that the construction is not restricted to point clouds, but extends to Banach-space functionals, graph signals, operator-adapted settings, and continuous multiwavelet limits (Balazs et al., 2024, Elefante et al., 25 Jul 2025, Giacchi et al., 2 Apr 2026).