Undecimated Discrete Curvelet Coefficients

• The paper introduces an FFT-based UDCT that enforces sparsity through an ℓ1 penalty, capturing fracture-aligned anisotropic features.
• UDCT is a multiscale, orientation-sensitive transform that omits spatial subsampling to preserve shift-invariance and localize curvilinear structures.
• Numerical benchmarks demonstrate a 37% reduction in RMSE, validating the effectiveness of UDCT regularization over traditional isotropic priors.

Undecimated discrete curvelet coefficients (UDCT) are a multiscale, orientation-sensitive representation of fields over discretized domains, constructed via an FFT-based transform that preserves shift-invariance by omitting spatial subsampling. In the context of fractured media, as introduced in Segura’s convex SPDE inversion framework, UDCT analysis and synthesis operators play a central role in regularizing reconstructions to capture anisotropic, fracture-aligned structures that elude isotropic priors (Segura, 24 Jan 2026).

1. Mathematical Definition of the Undecimated Discrete Curvelet Transform

On each piecewise-planar fracture support, the underlying scalar field is embedded on an $n \times n$ uniform grid ($M = n^2$). The UDCT analysis operator is

$$C: \mathbb{R}^M \longrightarrow \mathbb{C}^N, \quad z \mapsto d = C z,$$

where typically $N \approx 8M$, reflecting the transform's redundancy. The transform is constructed with reference to radial scale windows $\{ W_j(r) \}_{j=0}^{J}$ and angular windows $\{ V_{j,\ell}(\theta) \}_{\ell=0}^{L_j-1}$, with $j$ denoting the scale and $\ell$ the orientation. The generating curvelet in continuous frequency is

$$\varphi_{j,\ell,m}(u,v) = \frac{1}{(2\pi)^2} \int_{\mathbb{R}^2} W_j(\|\omega\|)\, V_{j,\ell}(\arg \omega)\, e^{i\,\omega \cdot [(u,v)-m]} \, d\omega,$$

where $(u,v)$ enumerates the spatial grid and $m$ the offset. In the undecimated version, there is no subsampling after frequency masking; all locations are processed. For admissibility, each $L_j$ is a multiple of 3, and windows are smoothly overlapped and summed, following the constructions of Candès et al. (2006).

2. Enforcement of Curvelet-Sparsity via $\ell_1$ Penalty

To encourage representation of fracture-aligned ridges and edges, the framework introduces a Laplace prior on the curvelet coefficients: $$\lambda\,\|C\,z\|_1 = \lambda\sum_{j=0}^J \sum_{\ell=0}^{L_j-1} \sum_{m=0}^{n^2-1} |d_{j,\ell}[m]|$$ with $$d_{j,\ell}[m] = \langle z, \varphi_{j,\ell,m} \rangle$$. The parabolic scaling and directional selectivity of $\varphi_{j,\ell,m}$ ensure that the $\ell_1$ penalty drives most coefficients toward zero except at the multiscale ridges and edges, enforcing sparsity preferentially along dominant fracture directions.

3. Algorithmic Computation of Forward and Inverse UDCT

Implementation comprises:

1. Padding: Fields $z$ are padded to $\tilde n \times \tilde n$ ($\approx$15% expansion per dimension) to suppress FFT wrap-around.
2. Forward UDCT:
   • Compute 2D FFT of padded $z$.
   • For each $(j, \ell)$, mask spectrum: $$\hat d_{j,\ell}(\omega) = W_j(\|\omega\|)\, V_{j,\ell}(\arg \omega)\, \hat z(\omega)$$.
   • Inverse FFT to obtain $d_{j,\ell}(u,v)$.
   • Crop to $n \times n$ if needed; vectorize and stack.
3. Inverse UDCT:
   • Pad coefficient slices and FFT.
   • Multiply by conjugated window masks, sum over $(j, \ell)$, inverse FFT, crop.

Memory use is dominated by window storage and the coefficient array ($N \approx 8 M$). The transform is $O(M \log M)$ per plane due to FFT efficiency and sharing of buffers.

4. Integration with Convex MAP Reconstruction and ADMM Splitting

The total per-plane optimization problem is

$$\hat z = \arg\min_{z} \frac12\|W^{1/2}(H z-y)\|_2^2 + \frac12 z^\top Q z + \lambda \|C z\|_1,$$

where $H$ is the interpolation matrix for observed points, $W$ encodes noise weights, and $Q$ is the Matérn-type GMRF precision. Introducing $d = C z$ and recasting as constrained minimization, the problem is solved by ADMM:

• $z$-update: Solve $$(H^\top W H + Q + \rho I)\,z^{k+1} = H^\top W y + \rho C^*(d^k-u^k)$$ via sparse Cholesky or CG.
• $d$-update: Apply soft-thresholding shrinkage to $C z^{k+1} + u^k$ with $\tau = \lambda/\rho$.
• Dual update: $u^{k+1} = u^k + (C z^{k+1} - d^{k+1})$.

ADMM convergence is guaranteed under convexity, with termination upon reduction of primal and dual residuals below set tolerances. Parseval-tightness ($C^*C \approx I$) ensures efficient block decoupling.

5. Soft-Thresholding (Shrinkage) Operator for Complex Coefficients

For coefficient $w$ and threshold $\tau$, the operator is

$$\mathrm{shrink}(w, \tau) = \max\left(0, 1 - \frac{\tau}{|w|}\right) w$$

yielding zero when $|w| \le \tau$ and $w - \tau w/|w|$ otherwise. For real $w$, this becomes $\mathrm{sign}(w)\max(|w|-\tau,0)$. $\tau$ is tied to the curvelet sparsity parameter and the ADMM penalty, with its selection guided by discrepancy principles or cross-validation.

6. Numerical Benchmarks Demonstrating Effect of UDCT Regularization

Section 9 presents synthetic “ridge + blob” experiments isolating the curvelet penalty, with summary metrics:

Method	RMSE (hold-out)	$R^2$
GMRF/SPDE (no curvelets)	$1.83 \pm 0.33$	$-2.97 \pm 4.28$
+ Curvelets via ADMM	$1.15 \pm 0.30$	$-0.12 \pm 0.35$

Thus, inclusion of $\ell_1$ UDCT penalty yields a $37\%$ decrease in RMSE and pushes $R^2$ near zero, demonstrating directional sparsity value. L-curve analysis substantiates classical residual–sparsity trade-offs and justifies $\lambda$ selection strategies. Robustness is demonstrated across grid resolutions $(n = 64, 128, 256)$, and additional studies detail sensitivity to smoothness and sparsity weights.

7. Context, Implications, and Significance

In the referenced framework, the UDCT block introduces anisotropic, multiscale sparse regularization that is efficiently split from conventional quadratic blocks via ADMM. The transform’s undecimated structure ensures shift-invariance and effective localization of ridge/edge features typical of fractured media. Empirical studies establish the superiority of curvelet sparsity over isotropic priors alone, particularly in reconstructing curvilinear anomalies. A plausible implication is that similar undecimated, directional transforms may serve as essential regularization tools in other inverse problems where anisotropic structures dominate (Segura, 24 Jan 2026).