Operator-Adapted Wavelet Decomposition

• Operator-adapted wavelet decomposition is a multiscale method that customizes wavelet bases to a given operator for improved computational accuracy.
• The approach integrates techniques like wavelet-vaguelette, gamblet transforms, diagonal frame decompositions, and Fourier multiplier designs to achieve sparse, localized representations.
• This tailored framework enhances efficiency in inversion, denoising, eigenanalysis, and finite-element corrections while significantly reducing computational complexity.

Operator-adapted wavelet decomposition denotes a class of multiscale constructions in which the decomposition is tailored to a given operator rather than chosen independently of it. In the literature represented here, this adaptation appears in several mathematically distinct forms: the wavelet-vaguelette decomposition for linear inverse problems, the gamblet transform for symmetric positive operators, diagonal frame decompositions that generalize the SVD to redundant frames, generalized wavelets built from a Fourier multiplier operator $L$, and fast product-convolution-to-wavelet algorithms that construct sparse wavelet matrices of spatially varying operators. A common objective is to obtain sparse, localized, or block-diagonal representations in which approximation, inversion, denoising, eigenpair computation, or multiscale finite-element correction become more efficient than in operator-agnostic coordinates (Johnstone et al., 2013, Yoo et al., 2018, Khalidov et al., 2012, Escande et al., 2020).

1. Principal meanings of operator adaptation

A first meaning of operator adaptation is the wavelet-vaguelette decomposition (WVD). For a linear operator $K$, the decomposition assumes wavelets $(\phi_{j_0k},\psi_{jk})$, vaguelette systems $\mathcal U=\{u_{jk}\}$ and $\mathcal V=\{v_{jk}\}$, and scale-dependent pseudo-singular values $\kappa_j$ such that

$$K\psi_{jk}=\kappa_j v_{jk}, \qquad Ku_{jk}=\kappa_j \psi_{jk},$$

with

$$\langle u_{jk},v_{j'k'}\rangle=\delta_{j-j'}\delta_{k-k'}.$$

This converts the white-noise inverse problem

$Y_\epsilon(dt)=(Kf)(t)\,dt+\epsilon W(dt)$

into the heteroscedastic Gaussian sequence model

$$y_{jk}=\theta_{jk}+\epsilon_j z_{jk}, \qquad \epsilon_j=\epsilon \kappa_j^{-1},$$

or, under the normalization $K$0,

$K$1

Here the decomposition is operator-adapted because the analysis functions are not raw wavelets of the data but vaguelettes tied to $K$2 (Johnstone et al., 2013).

A second meaning is the gamblet construction for symmetric positive operators $K$3. Starting from nested measurement functions $K$4, one defines operator-adapted pre-wavelets

$K$5

and operator-adapted wavelets

$K$6

The corresponding spaces satisfy

$K$7

orthogonally in the energy inner product

$K$8

Here adaptation is with respect to the operator-induced energy geometry rather than the $K$9 inner product (Yoo et al., 2018, Xie et al., 2018).

A third meaning is diagonalization by operator-adapted frames. A diagonal frame decomposition (DFD) for $(\phi_{j_0k},\psi_{jk})$0 is a triple $(\phi_{j_0k},\psi_{jk})$1 such that $(\phi_{j_0k},\psi_{jk})$2 is a frame of $(\phi_{j_0k},\psi_{jk})$3, $(\phi_{j_0k},\psi_{jk})$4 is a frame of $(\phi_{j_0k},\psi_{jk})$5, and

$(\phi_{j_0k},\psi_{jk})$6

This yields the reconstruction identity

$(\phi_{j_0k},\psi_{jk})$7

and a direct thresholding rule in operator-adapted coordinates. The construction generalizes the SVD to redundant systems such as wavelets, curvelets, and shearlets (Frikel et al., 2019).

A fourth meaning is basis construction from a Fourier multiplier operator. For an admissible operator $(\phi_{j_0k},\psi_{jk})$8 with symbol $(\phi_{j_0k},\psi_{jk})$9, the generalized scaling functions are

$\mathcal U=\{u_{jk}\}$0

the $\mathcal U=\{u_{jk}\}$1-interpolating spline is

$\mathcal U=\{u_{jk}\}$2

and the wavelet is defined by

$\mathcal U=\{u_{jk}\}$3

In this setting, the wavelet acts like the operator across scales (Khalidov et al., 2012).

A fifth meaning is representation-theoretic decomposition for wavelet representations associated with a finite-to-one endomorphism $\mathcal U=\{u_{jk}\}$4. After lifting to the solenoid $\mathcal U=\{u_{jk}\}$5, one obtains a covariant representation $\mathcal U=\{u_{jk}\}$6 and a direct-integral decomposition

$\mathcal U=\{u_{jk}\}$7

where $\mathcal U=\{u_{jk}\}$8 is a fundamental domain and the fiber representations are irreducible for $\mathcal U=\{u_{jk}\}$9-almost every $\mathcal V=\{v_{jk}\}$0. This is an operator-adapted decomposition in a representation-theoretic sense rather than a numerical sparse-matrix sense (Dutkay et al., 2011).

2. Fast construction of sparse wavelet matrices from product-convolution structure

A particularly explicit algorithmic realization is given for spatially varying operators through product-convolution expansions. Let $\mathcal V=\{v_{jk}\}$1 be the space of discrete $\mathcal V=\{v_{jk}\}$2-dimensional signals on

$\mathcal V=\{v_{jk}\}$3

with periodic boundary conditions, and let

$\mathcal V=\{v_{jk}\}$4

The target object is the wavelet-domain matrix

$\mathcal V=\{v_{jk}\}$5

for a product-convolution approximation

$\mathcal V=\{v_{jk}\}$6

The naive cost of generic wavelet decomposition is stated as $\mathcal V=\{v_{jk}\}$7 in general. The product-convolution model replaces this by a structurally adapted factorization

$\mathcal V=\{v_{jk}\}$8

with

$\mathcal V=\{v_{jk}\}$9

The adaptation occurs because the local kernels $\kappa_j$0 and coefficient maps $\kappa_j$1 are derived from the operator’s spatially varying impulse responses, by direct interpolation, adaptive sampling, PCA on sampled impulse responses, or randomized SVD of the SVIR (Escande et al., 2020).

The algorithm computes an $\kappa_j$2-accurate sparse approximation $\kappa_j$3 of $\kappa_j$4 such that

$\kappa_j$5

Its pipeline is explicit. One first builds or obtains the product-convolution expansion. For each term $\kappa_j$6, one sets

$\kappa_j$7

The convolution part

$\kappa_j$8

is then computed by exploiting the fact that the wavelet subbands of a convolution operator are circulant; each subband of $\kappa_j$9 is a rectangular circulant matrix. Consequently, $$K\psi_{jk}=\kappa_j v_{jk}, \qquad Ku_{jk}=\kappa_j \psi_{jk},$$0 can be encoded using only $$K\psi_{jk}=\kappa_j v_{jk}, \qquad Ku_{jk}=\kappa_j \psi_{jk},$$1 coefficients and computed in

$$K\psi_{jk}=\kappa_j v_{jk}, \qquad Ku_{jk}=\kappa_j \psi_{jk},$$2

operations. One then truncates $$K\psi_{jk}=\kappa_j v_{jk}, \qquad Ku_{jk}=\kappa_j \psi_{jk},$$3 to $$K\psi_{jk}=\kappa_j v_{jk}, \qquad Ku_{jk}=\kappa_j \psi_{jk},$$4 so that

$$K\psi_{jk}=\kappa_j v_{jk}, \qquad Ku_{jk}=\kappa_j \psi_{jk},$$5

The multiplier part

$$K\psi_{jk}=\kappa_j v_{jk}, \qquad Ku_{jk}=\kappa_j \psi_{jk},$$6

is sparse because

$$K\psi_{jk}=\kappa_j v_{jk}, \qquad Ku_{jk}=\kappa_j \psi_{jk},$$7

vanishes unless $$K\psi_{jk}=\kappa_j v_{jk}, \qquad Ku_{jk}=\kappa_j \psi_{jk},$$8 and $$K\psi_{jk}=\kappa_j v_{jk}, \qquad Ku_{jk}=\kappa_j \psi_{jk},$$9 overlap. If

$$\langle u_{jk},v_{j'k'}\rangle=\delta_{j-j'}\delta_{k-k'}.$$0

then

$$\langle u_{jk},v_{j'k'}\rangle=\delta_{j-j'}\delta_{k-k'}.$$1

The matrix $$\langle u_{jk},v_{j'k'}\rangle=\delta_{j-j'}\delta_{k-k'}.$$2 is computed by a two-stage sparse cascade algorithm with cost

$$\langle u_{jk},v_{j'k'}\rangle=\delta_{j-j'}\delta_{k-k'}.$$3

Finally, one forms

$$\langle u_{jk},v_{j'k'}\rangle=\delta_{j-j'}\delta_{k-k'}.$$4

The error control is spectral-norm explicit: $$\langle u_{jk},v_{j'k'}\rangle=\delta_{j-j'}\delta_{k-k'}.$$5 Because $$\langle u_{jk},v_{j'k'}\rangle=\delta_{j-j'}\delta_{k-k'}.$$6 is orthogonal and $$\langle u_{jk},v_{j'k'}\rangle=\delta_{j-j'}\delta_{k-k'}.$$7 is diagonal,

$$\langle u_{jk},v_{j'k'}\rangle=\delta_{j-j'}\delta_{k-k'}.$$8

so the choice $$\langle u_{jk},v_{j'k'}\rangle=\delta_{j-j'}\delta_{k-k'}.$$9 gives

$Y_\epsilon(dt)=(Kf)(t)\,dt+\epsilon W(dt)$0

Under the regularity assumptions, Algorithm 1 computes $Y_\epsilon(dt)=(Kf)(t)\,dt+\epsilon W(dt)$1 in

$Y_\epsilon(dt)=(Kf)(t)\,dt+\epsilon W(dt)$2

operations, and the number of nonzero coefficients in $Y_\epsilon(dt)=(Kf)(t)\,dt+\epsilon W(dt)$3 is

$Y_\epsilon(dt)=(Kf)(t)\,dt+\epsilon W(dt)$4

In many smooth applications $Y_\epsilon(dt)=(Kf)(t)\,dt+\epsilon W(dt)$5, so $Y_\epsilon(dt)=(Kf)(t)\,dt+\epsilon W(dt)$6, which is the quasi-linear regime emphasized by the authors (Escande et al., 2020).

3. Operator-adapted decompositions in inverse problems, denoising, eigenanalysis, and basis design

In inverse problems, the WVD provides an operator-adapted coordinate model in which estimation becomes levelwise penalized regression. The estimator at level $Y_\epsilon(dt)=(Kf)(t)\,dt+\epsilon W(dt)$7 is

$Y_\epsilon(dt)=(Kf)(t)\,dt+\epsilon W(dt)$8

with scale-dependent penalty calibrated by the ill-posedness of $Y_\epsilon(dt)=(Kf)(t)\,dt+\epsilon W(dt)$9. The main theorem partitions the parameter domain into dense, sparse, and critical regimes and states exact minimax rates over Besov sequence balls. In this setting, operator adaptation enters through the vaguelettes and pseudo-singular values, while adaptation to unknown regularity enters through levelwise complexity penalization (Johnstone et al., 2013).

For denoising with prior information on $$y_{jk}=\theta_{jk}+\epsilon_j z_{jk}, \qquad \epsilon_j=\epsilon \kappa_j^{-1},$$0 rather than on $$y_{jk}=\theta_{jk}+\epsilon_j z_{jk}, \qquad \epsilon_j=\epsilon \kappa_j^{-1},$$1, gamblet coefficients replace ordinary wavelet coefficients. The partial reconstruction at level $$y_{jk}=\theta_{jk}+\epsilon_j z_{jk}, \qquad \epsilon_j=\epsilon \kappa_j^{-1},$$2 is

$$y_{jk}=\theta_{jk}+\epsilon_j z_{jk}, \qquad \epsilon_j=\epsilon \kappa_j^{-1},$$3

and the main estimator is

$$y_{jk}=\theta_{jk}+\epsilon_j z_{jk}, \qquad \epsilon_j=\epsilon \kappa_j^{-1},$$4

where $$y_{jk}=\theta_{jk}+\epsilon_j z_{jk}, \qquad \epsilon_j=\epsilon \kappa_j^{-1},$$5 minimizes

$$y_{jk}=\theta_{jk}+\epsilon_j z_{jk}, \qquad \epsilon_j=\epsilon \kappa_j^{-1},$$6

The paper also gives hard and soft thresholding in gamblet coordinates. The principal guarantee is near minimax optimality, up to a multiplicative constant, over

$$y_{jk}=\theta_{jk}+\epsilon_j z_{jk}, \qquad \epsilon_j=\epsilon \kappa_j^{-1},$$7

in the energy norm induced by the operator (Yoo et al., 2018).

The same gamblet decomposition supports eigenpair computation. The multiresolution splitting

$$y_{jk}=\theta_{jk}+\epsilon_j z_{jk}, \qquad \epsilon_j=\epsilon \kappa_j^{-1},$$8

is energy-orthogonal, with level matrices

$$y_{jk}=\theta_{jk}+\epsilon_j z_{jk}, \qquad \epsilon_j=\epsilon \kappa_j^{-1},$$9

The algorithm first solves the eigenproblem on a coarse gamblet space and then performs hierarchical subspace correction. At each level $K$00, the correction solve is

$K$01

followed by a Rayleigh-Ritz step in

$K$02

The method is designed for bijective positive symmetric operators and is shown to have near-linear complexity when $K$03 is a local operator mapping $K$04 to $K$05 (Xie et al., 2018).

For direct sparse regularization, the DFD yields the explicit thresholding formula

$K$06

This is a non-iterative operator-adapted frame thresholding method. If the frame is a basis, then synthesis and analysis $K$07-regularization as well as DFD thresholding are equivalent; in the redundant case, those three approaches are pairwise different (Frikel et al., 2019).

At the level of basis construction, operator-like wavelets are generated from a Fourier multiplier operator by first defining generalized B-splines and then setting

$K$08

In the stochastic model

$K$09

the coefficient identity

$K$10

shows why sparsity of the innovation can be inherited by the wavelet coefficients (Khalidov et al., 2012).

4. Hierarchical finite-element realizations on unstructured meshes

A recent finite-element line of work constructs operator-adapted wavelet decompositions on unstructured triangular meshes with Whitney edge elements for multiscale electromagnetic problems. The continuous equation is

$K$11

and in the principal experiments the governing PDE is the vector Helmholtz equation

$K$12

with boundary conditions

$K$13

Starting from nested FEM spaces $K$14, the operator-adapted hierarchy imposes

$K$15

so that

$K$16

The transformed stiffness matrix becomes block diagonal,

$K$17

which decouples the scales. The hierarchy is built from finest to coarser levels using sparse refinement matrices $K$18, null-space matrices $K$19, and recursively defined operator-adapted transfers $K$20. The coarse and detail right-hand sides are

$K$21

and the independent systems are

$K$22

The full solution is reconstructed as

$K$23

The principal numerical claim is that the overall computational complexity is nearly linear, with sparse matrix-vector multiplications as the dominant operations (Şık et al., 23 Jul 2025).

A related polygonal-element formulation begins from a finest triangular mesh and generates coarser meshes by agglomerating triangles into convex polygons. The operator-adapted hierarchy is again

$K$24

with operator orthogonality

$K$25

On a convex polygon $K$26, generalized barycentric coordinates $K$27 are used to define Whitney one-form edge bases

$K$28

The operator-agnostic refinement matrices $K$29 and null-space matrices $K$30 are assembled first, after which the operator-adapted matrices $K$31 and compressed operators $K$32 are defined recursively. As in the triangular-mesh formulation, the detail levels can be added without recomputing the coarser levels (Şık et al., 17 Dec 2025).

5. Computational properties and demonstrated applications

The product-convolution-to-wavelet algorithm was validated on two space-varying image blurs on images up to $K$33: a vertically varying isotropic Gaussian blur and a rotated/skewed Gaussian blur varying with location. The product-convolution ranks were chosen as $K$34 and $K$35, and the decomposition precision was $K$36. For $K$37 million pixels, the reported speed-up over naive decomposition using product-convolution was $K$38 for blur 1 and $K$39 for blur 2. In wavelet-$K$40 deblurring on $K$41 images, wavelet-domain FISTA using the sparse wavelet operator was much faster than spatial-domain or product-convolution FISTA. For blur 1, FISTA-W gave $K$42 speed-up over exact FISTA, preconditioned FISTA-WP gave $K$43, and GPU $K$44; for blur 2, the corresponding figures were $K$45, $K$46, and $K$47. The loss in image quality was less than $K$48 dB relative to the exact model in that setup (Escande et al., 2020).

In gamblet denoising, the central guarantee is not an empirical benchmark but a statistical one: $K$49 The paper also proves a high-probability energy-norm bound stating that, if $K$50, then with probability at least $K$51,

$K$52

The transform itself is computable to precision $K$53 in

$K$54

complexity, with an incomplete Cholesky version stated as

$K$55

This places gamblet thresholding among near-linear operator-adapted decomposition methods (Yoo et al., 2018).

For low eigenpairs of rough operators, gamblet multilevel correction was tested on the SPE10 elliptic benchmark, random checkerboard coefficients, and Anderson localization or disordered Schrödinger-type problems. The reported observations were that gamblet-based multilevel correction converges much faster than geometric multigrid-based correction, that the first 12 eigenvalues improve dramatically level by level, and that in the SPE10 case the gamblet multilevel correction method showed about a ten-fold speedup in online CPU time versus ARPACK in the reported implementation (Xie et al., 2018).

In multiscale electromagnetics on unstructured triangular meshes, the L-shaped and U-shaped waveguide problems were used as validation. For the L-shaped waveguide, the relative $K$56 error between the finest-level FEM solution and the hierarchical operator-adapted solution with coarse plus three detail levels was about

$K$57

and the relative $K$58 error versus numerical mode matching was

$K$59

For the U-shaped waveguide, the corresponding errors were approximately

$K$60

For finest-level DoFs ranging from $K$61 to $K$62, the reported cost was about $K$63 per iteration in a two-level scenario, about $K$64 excluding precomputation in a six-level scenario, and about $K$65 including precomputation (Şık et al., 23 Jul 2025).

In the polygonal-element formulation, the wedge singularity problem and a microporous Si slab were the main examples. The reported fitted slopes from the wedge example were approximately $K$66 for the coarsest-level calculation, $K$67 for the first detail, $K$68 for the second detail, $K$69 for the third detail, and $K$70 for the fourth detail; including precomputation, the global fitted slope was about $K$71. With four detail levels, the relative $K$72-error in the near-tip region dropped below $K$73 relative to finest-level FEM, while the error versus the analytic series solution was about $K$74 or less in the same region. The paper also reports peak memory savings up to about $K$75 (Şık et al., 17 Dec 2025).

6. Adjacent methods, common misconceptions, and structural limitations

A recurrent misconception is that any use of wavelets inside an operator-learning or PDE pipeline constitutes an operator-adapted wavelet decomposition. The Wavelet Neural Operator is explicitly described as not a classical operator-adapted wavelet decomposition algorithm in the numerical-analysis sense. It uses a fixed multilevel discrete wavelet transform, with the operator-specific part learned through trainable kernels acting on wavelet coefficients. The basis is fixed in advance, with wavelet families such as $K$76 and $K$77, and there is no algorithm that derives wavelets from the PDE operator. What is adapted is the learned mapping in wavelet space, not the decomposition itself (Tripura et al., 2022).

A different nearby construction is the monogenic synchrosqueezed wavelet transform for AM-FM image decomposition. It is adapted to the monogenic/Riesz operator and to a structured class of oscillatory image models, with local wave-vector estimators

$K$78

and a reassigned representation

$K$79

The principal theorem states that significant wavelet coefficients occur only near the scale matching one mode and that each separated mode can be reconstructed by integrating the synchrosqueezed representation over a small neighborhood of its estimated local wave vector. This is an operator-adapted decomposition in the sense of analysis geometry and reassignment, but not in the sense of building wavelets from an arbitrary elliptic, convolution, or inverse operator (Clausel et al., 2012).

Across the numerical-analysis literature, the main limitations are structural. The product-convolution route assumes that the operator either is already in product-convolution form or can be approximated well by one; if the SVIR is not low-rank or local translation invariance is weak, $K$80 may need to be large. The guarantees in that setting concern decomposition of $K$81, not the original operator $K$82, unless $K$83. The WVD approach applies only to operators admitting a wavelet-vaguelette decomposition, with scale-only pseudo-singular values and frame bounds. Gamblet constructions assume a symmetric positive bijection and rely on hierarchical measurements, operator locality, and exponentially localized bases for near-linear computation. DFD thresholding requires a diagonal frame decomposition and yields a direct thresholding method, but in the redundant case it is distinct from both analysis and synthesis $K$84-regularization. The finite-element electromagnetic formulations depend on preserving sparsity under recursive operator-adapted transformations and presently focus on 2D mesh hierarchies with Whitney edge elements or convex polygonal elements (Escande et al., 2020, Johnstone et al., 2013, Yoo et al., 2018, Frikel et al., 2019, Şık et al., 23 Jul 2025, Şık et al., 17 Dec 2025).

Taken together, these works show that “operator-adapted wavelet decomposition algorithm” does not denote a single canonical procedure. It denotes a family of constructions in which the multiscale representation is tied to an operator through vaguelettes, energy-orthogonal wavelets, diagonal frame relations, Fourier-multiplier design, product-convolution structure, or operator-orthogonal finite-element detail spaces. This suggests that the most precise use of the term depends on which object is being adapted: the coefficient-extraction map, the basis functions, the transformed stiffness matrix, or the representation itself.