Generator-Guided Inverse Sampling for L\'evy-Driven Generative Models

# Generator-Guided Inverse Sampling for Lévy-Driven Generative Models
Source: [https://arxiv.org/html/2608.10384](https://arxiv.org/html/2608.10384)
Tianfu Qi,*Graduate Student Member, IEEE*, Jun Wang,*Senior Member, IEEE*, Jun Zhang,*Fellow, IEEE*Tianfu Qi, Jun Wang are with the National Key Laboratory of Wireless Communications, University of Electronic Science and Technology of China, Chengdu 611731, China \(e\-mail: 202311220634@std\.uestc\.edu\.cn; junwang@uestc\.edu\.cn\)\.J\. Zhang is with the Department of Electronic and Computer Engineering, Hong Kong University of Science and Technology, Hong Kong, China \(e\-mail: eejzhang@ust\.hk\)\.

###### Abstract

This paper studies inverse sampling for Lévy\-driven generative models from the perspective of Markov generators\. Unlike conventional diffusion models, Lévy\-driven dynamics involve infinite jump activities, which makes their reverse process nonlocal and difficult to characterize using score information alone\. We address this challenge by analyzing the forward and reversed generators\. It is derived that the reversed jump component generally becomes a state\-dependent Markov jump process governed by a nonlocal density ratio\. This observation motivates a structured reverse sampler that decomposes the dynamics into diffusion, small jump, and large jump components\. Based on this characterization, we develop a computationally tractable sampler for a class of isotropic linear Lévy SDEs with symmetricα\\alpha\-stable jump components\. For the jump component, the neural network is used only to amortize the rate of large jump activities, while jump amplitudes are generated from analytically derived conditional distributions, which improves interpretability and controllability\. Efficient implementation techniques are further introduced under this setting to avoid expensive high\-dimensional integration and sampling\. The sampler is further adapted to approximate observation\-guided sampling and applied to OFDM\-SISO channel estimation under mixed Gaussian and impulsive noise\. Simulations show robust estimation performance with a favorable tradeoff between complexity and performance\.

## IIntroduction

Diffusion models based on score matching have emerged as a powerful class of generative models\[[21](https://arxiv.org/html/2608.10384#bib.bib4),[9](https://arxiv.org/html/2608.10384#bib.bib5),[22](https://arxiv.org/html/2608.10384#bib.bib6),[24](https://arxiv.org/html/2608.10384#bib.bib7)\]\. They have also attracted increasing attention in signal processing applications, including inverse problem solving\[[13](https://arxiv.org/html/2608.10384#bib.bib8),[12](https://arxiv.org/html/2608.10384#bib.bib9),[5](https://arxiv.org/html/2608.10384#bib.bib10),[20](https://arxiv.org/html/2608.10384#bib.bib11),[30](https://arxiv.org/html/2608.10384#bib.bib12),[10](https://arxiv.org/html/2608.10384#bib.bib13)\], image restoration\[[26](https://arxiv.org/html/2608.10384#bib.bib14),[7](https://arxiv.org/html/2608.10384#bib.bib15),[3](https://arxiv.org/html/2608.10384#bib.bib16),[16](https://arxiv.org/html/2608.10384#bib.bib17)\], and wireless channel estimation\[[29](https://arxiv.org/html/2608.10384#bib.bib18),[2](https://arxiv.org/html/2608.10384#bib.bib19),[14](https://arxiv.org/html/2608.10384#bib.bib20)\]\. By gradually perturbing data with Gaussian noise and learning the reverse dynamics through score matching, these models provide a stochastic framework for sampling from complex data distributions\.

Despite their success, conventional diffusion models are mainly built on the Wiener process\. Their reverse process is characterized by local score information and infinitesimal Gaussian perturbations\[[21](https://arxiv.org/html/2608.10384#bib.bib4),[22](https://arxiv.org/html/2608.10384#bib.bib6)\]\. As a result, long\-range probability transport is usually achieved through a sequence of small reverse steps\. This local update mechanism may affect sampling efficiency and modeling flexibility\[[11](https://arxiv.org/html/2608.10384#bib.bib21)\], especially when the target distribution has heavy tails, impulsive components, or highly multimodal structures\[[27](https://arxiv.org/html/2608.10384#bib.bib22),[19](https://arxiv.org/html/2608.10384#bib.bib23)\]\.

Lévy\-driven generative models provide a natural extension by incorporating jump components into the forward dynamics\[[1](https://arxiv.org/html/2608.10384#bib.bib1)\]\. The resulting nonlocal transitions are well suited for modeling heavy\-tailed perturbations\. They can also provide more flexible distribution transport than standard diffusion dynamics\[[27](https://arxiv.org/html/2608.10384#bib.bib22)\]\. However, these advantages bring new challenges for inverse sampling\. Unlike the case without jumps, the time reversal of a Lévy\-driven process cannot be fully characterized by local score information alone\[[27](https://arxiv.org/html/2608.10384#bib.bib22),[19](https://arxiv.org/html/2608.10384#bib.bib23)\]\. Its reversed jump component depends on a nonlocal density ratio\. It generally becomes a Markov jump process with state\-dependent kernels, which makes the inverse sampler difficult to design\.

It has been shown in\[[27](https://arxiv.org/html/2608.10384#bib.bib22)\]that the denoising score matching \(DSM\) framework can still be applied to inverse sampling for Lévy processes\. Under this formulation, the neural network is required to learn denoising targets associated with impulsive perturbations, such as those induced byα\\alpha\-stable distributions\. However, different from white Gaussian noise \(WGN\), impulse noise \(IN\) contains multiple outliers with large amplitudes\. Compared with WGN, which has a relatively stable envelope over a long time interval, IN is much more impulsive and chaotic\. Therefore, it is difficult to learn the statistical characteristics of impulse noise at different moments\. In this case, full neural parameterization makes the sampler less interpretable and controllable\.

Motivated by these observations, we aim to establish a generator\-guided inverse sampling approach for Lévy\-driven generative models\. The generator provides a fundamental characterization of the infinitesimal statistics of Markov processes and naturally reveals the nonlocal structure of the reversed jump dynamics\. Based on this structural characterization, we design a practical sampler under isotropic linearα\\alpha\-stable assumptions, so that the main components of the reverse dynamics can be implemented in a statistically interpretable and computationally tractable manner\. The main contributions of this paper are summarized as follows\.

- •First, we provide a generator\-based characterization of the time\-reversed dynamics of Lévy\-driven Markov processes\.The analysis shows that the reversed jump component is generally a state\-dependent Markov jump process whose kernel involves a nonlocal density ratio\. This result clarifies why reverse sampling for Lévy\-driven processes cannot, in general, be reduced to score\-based local diffusion updates\.
- •Second, motivated by this characterization, we develop a practical inverse sampler for isotropic linear Lévy SDEs with symmetricα\\alpha\-stable jump components\.The reverse dynamics are decomposed into diffusion, small jump, and large jump parts\. The small jump contribution is approximated through a Gaussian surrogate, while large jumps are handled explicitly through a rate\-and\-amplitude decomposition\. The neural network is used to amortize the large jump rate, whereas jump amplitudes are sampled from analytically derived conditional distributions\.
- •Third, we adapt the proposed sampler to approximate observation\-guided sampling under observation models and demonstrate its use in OFDM\-SISO channel estimation with mixed Gaussian and impulsive noise\.In particular, the prior\-trained jump rate network is retained as a proposal mechanism and the observation likelihood is incorporated through reweighted distributions\. The resulting algorithm provides a robust estimator under heavy\-tailed perturbations and offers an interpretable alternative to fully neural reverse\-transition parameterizations\.

The remainder of this paper is organized as follows\. In Section[II](https://arxiv.org/html/2608.10384#S2), we describe the Lévy process considered throughout this paper and present the corresponding assumptions\. Statistical analyses from the generator perspective are given in Section[III](https://arxiv.org/html/2608.10384#S3)\. Based on the above theoretical results, an efficient inverse sampling algorithm is designed in Section[IV](https://arxiv.org/html/2608.10384#S4)\. Section[V](https://arxiv.org/html/2608.10384#S5)provides simulation results for channel estimation under mixed channel noise to validate the proposed inverse sampling framework\. Section[VI](https://arxiv.org/html/2608.10384#S6)concludes the paper\.

*Notations:*Vectors are denoted by bold letters and are assumed to be column vectors by default\. Superscripts and subscripts are also used to denote vectors and vector sequences\. For example, given anmm\-dimensional vector𝒙\\bm\{x\}, we have𝒙=\[x1,x2,⋯,xm\]⊤\\bm\{x\}=\[x\_\{1\},x\_\{2\},\\cdots,x\_\{m\}\]^\{\\top\}and𝒙1p=\[𝒙1,⋯,𝒙p\]⊤\\bm\{x\}\_\{1\}^\{p\}=\[\\bm\{x\}\_\{1\},\\cdots,\\bm\{x\}\_\{p\}\]^\{\\top\}\. Uppercase letters denote random variables \(RVs\), and lowercase letters denote their realizations, e\.g\.,XXandxx\. The operatorsΔ\\Delta,∇\\nabla, and∇⋅\\nabla\\cdotdenote the Laplace operator, gradient operator, and divergence operator, respectively\. The notationdiag​\(a1,a2,⋯,aL\)\\text\{diag\}\(a\_\{1\},a\_\{2\},\\cdots,a\_\{L\}\)represents a diagonal matrix with diagonal elementsa1,a2,⋯,aLa\_\{1\},a\_\{2\},\\cdots,a\_\{L\}\.⟨𝒙,𝒚⟩\\langle\\bm\{x\},\\bm\{y\}\\rangledenotes the inner product of vectors𝒙\\bm\{x\}and𝒚\\bm\{y\}\.1𝒜​\(x\)1\_\{\\mathcal\{A\}\}\(x\)denotes the indicator function, which equals 1 ifx∈𝒜x\\in\\mathcal\{A\}\. The real number domain is denoted byℝ\\mathbb\{R\}\. The notationX∼⋅X\\sim\\cdotindicates thatXXfollows a certain distribution\.

## IIPreliminaries

In this paper, we consider a Lévy process with drift and diffusion components, which can be generally written as

d​𝑿t=𝒃​\(𝑿t,t\)​d​t\+ΦG​\(t\)​d​𝑾t\\displaystyle d\{\\bm\{X\}\_\{t\}\}=\\bm\{b\}\(\\bm\{X\}\_\{t\},t\)dt\+\\Phi\_\{G\}\(t\)d\{\\bm\{W\}\_\{t\}\}\+ΦS​\(t\)​d​𝑳t,\\displaystyle\+\\Phi\_\{S\}\(t\)d\\bm\{L\}\_\{t\},t:0→T,𝑿t∈ℝD\\displaystyle t:0\\rightarrow T,\\bm\{X\}\_\{t\}\\in\\mathbb\{R\}^\{D\}\(1\)whereDDdenotes the dimension of the SDE\. The term𝒃​\(𝑿t,t\)\\bm\{b\}\(\\bm\{X\}\_\{t\},t\)is the drift coefficient of the dynamics, which depends on both the time indexttand the system state𝑿t\\bm\{X\}\_\{t\}\. The processes𝑾t\\bm\{W\}\_\{t\}and𝑳t\\bm\{L\}\_\{t\}denote the Wiener process and theα\\alpha\-stable process, respectively\. The matricesΦG​\(t\)\\Phi\_\{G\}\(t\)andΦS​\(t\)\\Phi\_\{S\}\(t\)denote the scaling matrices for𝑾t\\bm\{W\}\_\{t\}and𝑳t\\bm\{L\}\_\{t\}, respectively\. For isotropic cases,ΦG​\(t\)\\Phi\_\{G\}\(t\)andΦS​\(t\)\\Phi\_\{S\}\(t\)reduce to scaled identity matrices\. In the sequel, we focus on the SDE under the following assumptions\.

###### Assumption 1

The Wiener process𝐖t\\bm\{W\}\_\{t\}and the Lévy process𝐋t\\bm\{L\}\_\{t\}are mutually independent\.

###### Assumption 2

𝑳t\\bm\{L\}\_\{t\}is aDD\-dimensional symmetric isotropicα\\alpha\-stable Lévy process with stability indexα∈\(0,2\)\\alpha\\in\(0,2\)\.

Our ultimate goal is to investigate the characteristics of the reverse process corresponding to the SDE in \([II](https://arxiv.org/html/2608.10384#S2.Ex1)\) and design an efficient reverse sampler\. Unfortunately, unlike a diffusion process driven only by Brownian motion, \([II](https://arxiv.org/html/2608.10384#S2.Ex1)\) is a Lévy process with infinite jump activities\. Its reverse process cannot be described using only local information\. In other words, if we setΦS​\(t\)≡𝟎\\Phi\_\{S\}\(t\)\\equiv\\bm\{0\}, the backward process of \([II](https://arxiv.org/html/2608.10384#S2.Ex1)\) can be determined by the score information, i\.e\.,∇log⁡p​\(𝒙t\)\\nabla\\log p\(\\bm\{x\}\_\{t\}\)\. This information is local since it only requires the gradient of the density\. This property can be obtained by analyzing the forward and backward Kolmogorov equations, which are special cases of the Kramers\-Moyal \(KM\) expansion\. The core idea of the KM expansion is to express the time evolution of the probability density as an infinite series of state\-space changes\.

For diffusion processes driven by Brownian motion, this series can be truncated after the first two terms, which leads to the forward Kolmogorov equation\. However, ifΦS​\(t\)≠𝟎\\Phi\_\{S\}\(t\)\\neq\\bm\{0\}, large jumps may appear in \([II](https://arxiv.org/html/2608.10384#S2.Ex1)\) due tod​𝑳td\\bm\{L\}\_\{t\}\. As a result,‖𝑿t\+Δ​t−𝑿t‖\\\|\\bm\{X\}\_\{t\+\\Delta t\}\-\\bm\{X\}\_\{t\}\\\|can be quite large\. Moreover, the KM expansion involves thejj\-th moment ofXtX\_\{t\}, wherej∈ℤ\+​⋃\{0\}j\\in\\mathbb\{Z\}^\{\+\}\\bigcup\\\{0\\\}\. For a Lévy process driven by anα\\alpha\-stable process, finitepp\-th moments withp≥αp\\geq\\alphado not exist\. Therefore, the reverse process becomes much more complicated\. In the following, we use the generator to analyze the statistical properties of \([II](https://arxiv.org/html/2608.10384#S2.Ex1)\)\.

Given the forward dynamics in \([II](https://arxiv.org/html/2608.10384#S2.Ex1)\), its generator is given by\[[1](https://arxiv.org/html/2608.10384#bib.bib1)\]

\(ℒF​f\)​\(𝒙\)=\\displaystyle\\left\(\{\\mathcal\{L\}\_\{F\}f\}\\right\)\(\\bm\{x\}\)=𝒃​\(𝒙,t\)⋅∇f​\(𝒙\)\+12​Tr​\(Σ​\(t\)​∇2f​\(𝒙\)\)\\displaystyle\\bm\{b\}\(\\bm\{x\},t\)\\cdot\\nabla f\(\\bm\{x\}\)\+\\frac\{1\}\{2\}\{\\rm\{Tr\}\}\\left\(\{\\Sigma\(t\)\\nabla^\{2\}f\(\\bm\{x\}\)\}\\right\)\+\(ℒJ,F​f\)​\(𝒙\)\\displaystyle\+\(\\mathcal\{L\}\_\{J,F\}f\)\(\\bm\{x\}\)\(2\)whereΣ​\(t\)=ΦG​\(t\)​ΦG​\(t\)⊤\\Sigma\(t\)=\\Phi\_\{G\}\(t\)\\Phi\_\{G\}\(t\)^\{\\top\}and\(ℒJ,F​f\)​\(𝒙\)\(\\mathcal\{L\}\_\{J,F\}f\)\(\\bm\{x\}\)denotes the forward generator of the processΦS​\(t\)​d​𝑳t\\Phi\_\{S\}\(t\)d\\bm\{L\}\_\{t\}\. The backward generator of the reversed version of \([II](https://arxiv.org/html/2608.10384#S2.Ex1)\) will be characterized in the next section\. Before ending this section, we introduce the following assumptions, which ensure the validity of the subsequent analysis\.

###### Assumption 3

The forward SDE admits a non\-explosive Markov solution with transition densityp​\(𝐱t\)p\(\\bm\{x\}\_\{t\}\)with respect to the Lebesgue measure\.

###### Assumption 4

The densityp​\(𝐱t\)p\(\\bm\{x\}\_\{t\}\)is strictly positive and sufficiently smooth so that∇log⁡p​\(𝐱t\)\\nabla\\log p\(\\bm\{x\}\_\{t\}\)is well defined\.

## IIIGenerator\-Based Reverse Dynamics

In this section, we first derive the backward generator corresponding to the forward process in \([II](https://arxiv.org/html/2608.10384#S2.Ex1)\)\. Then, we theoretically analyze how to decompose the whole inverse sampling process and generate each component\.

### III\-ABackward generator

Throughout this paper, the term “backward generator” refers to the infinitesimal generator of the time\-reversed process\. It should not be confused with the backward Kolmogorov operator associated with the forward process\. For notational clarity, we writep​\(𝒙t\)p\(\\bm\{x\}\_\{t\}\)for the marginal density of𝑿t\\bm\{X\}\_\{t\}\.

Based on \([II](https://arxiv.org/html/2608.10384#S2.Ex1)\) and its forward generator \([II](https://arxiv.org/html/2608.10384#S2.Ex2)\), the backward generator is given in the following proposition\.

###### Proposition 1\(Backward generator\)

Consider the forward SDE in \([II](https://arxiv.org/html/2608.10384#S2.Ex1)\) with Lévy measureν\\nu\. Letνt=\(ΦS​\(t\)\)\#​ν\\nu\_\{t\}=\(\\Phi\_\{S\}\(t\)\)\_\{\\\#\}\\nube the push\-forward Lévy measure induced by the jump scaling matrixΦS​\(t\)\\Phi\_\{S\}\(t\)\. Suppose that the forward process admits a strictly positive marginal densityp​\(𝐱t\)p\(\\bm\{x\}\_\{t\}\), and thatp​\(𝐱t\)p\(\\bm\{x\}\_\{t\}\)and𝐛​\(𝐱,t\)\\bm\{b\}\(\\bm\{x\},t\)are sufficiently smooth and integrable so that the following integrations by parts and principal\-value integrals are well defined\. If the scaled Lévy measureνt\\nu\_\{t\}is symmetric, then for any test functionf∈Cc2​\(ℝD\)f\\in C\_\{c\}^\{2\}\(\\mathbb\{R\}^\{D\}\), the generator of the reversed process can be written as

\(LB​f\)​\(𝒙t\)=\\displaystyle\(L\_\{B\}f\)\(\\bm\{x\}\_\{t\}\)=\[−𝒃​\(𝒙t,t\)\+Σ​\(t\)​∇log⁡p​\(𝒙t\)\]⋅∇f​\(𝒙t\)\\displaystyle\\left\[\-\\bm\{b\}\(\\bm\{x\}\_\{t\},t\)\+\\Sigma\(t\)\\nabla\\log p\(\\bm\{x\}\_\{t\}\)\\right\]\\cdot\\nabla f\(\\bm\{x\}\_\{t\}\)\+12​Tr⁡\(Σ​\(t\)​∇2f​\(𝒙t\)\)\+\(LJ,B​f\)​\(𝒙t\)\\displaystyle\+\\frac\{1\}\{2\}\\operatorname\{Tr\}\\left\(\\Sigma\(t\)\\nabla^\{2\}f\(\\bm\{x\}\_\{t\}\)\\right\)\+\(L\_\{J,B\}f\)\(\\bm\{x\}\_\{t\}\)\(3\)\(LJ,B​f\)​\(𝒙t\)=\\displaystyle\(L\_\{J,B\}f\)\(\\bm\{x\}\_\{t\}\)=p\.v\.​∫ℝD\[f​\(𝒙t\+𝒗\)−f​\(𝒙t\)\]\\displaystyle\\operatorname\{p\.v\.\}\\int\_\{\\mathbb\{R\}^\{D\}\}\\left\[f\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\-f\(\\bm\{x\}\_\{t\}\)\\right\]×p​\(𝒙t\+𝒗\)p​\(𝒙t\)​νt​\(d​𝒗\)\.\\displaystyle\\times\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}\\nu\_\{t\}\(d\\bm\{v\}\)\.\(4\)wherep\.v\.\\operatorname\{p\.v\.\}denotes the Cauchy principal\-value interpretation\.

###### Proof:

The proof is relegated to Appendix[A](https://arxiv.org/html/2608.10384#A1)\. ∎

Remark 1:The jump integral in Proposition[1](https://arxiv.org/html/2608.10384#Thmproposition1)is understood as

p\.v\.​∫ℝD\[f​\(𝒙t\+𝒗\)−f​\(𝒙t\)\]​p​\(𝒙t\+𝒗\)p​\(𝒙t\)​νt​\(d​𝒗\)\\displaystyle\\operatorname\{p\.v\.\}\\int\_\{\\mathbb\{R\}^\{D\}\}\\left\[f\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\-f\(\\bm\{x\}\_\{t\}\)\\right\]\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}\\nu\_\{t\}\(d\\bm\{v\}\)=\\displaystyle=limϵ→0∫‖𝒗‖\>ϵ\[f​\(𝒙t\+𝒗\)−f​\(𝒙t\)\]​p​\(𝒙t\+𝒗\)p​\(𝒙t\)​νt​\(d​𝒗\)\.\\displaystyle\\lim\_\{\\epsilon\\rightarrow 0\}\\int\_\{\\\|\\bm\{v\}\\\|\>\\epsilon\}\\left\[f\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\-f\(\\bm\{x\}\_\{t\}\)\\right\]\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}\\nu\_\{t\}\(d\\bm\{v\}\)\.Forα<1\\alpha<1, the integral is locally absolutely integrable under standard smoothness assumptions\. Forα≥1\\alpha\\geq 1, the uncompensated integral is generally not absolutely integrable near the origin and should be interpreted in the Cauchy principal\-value sense\. To see this, for smoothffandp​\(𝒙t\)p\(\\bm\{x\}\_\{t\}\), we havef​\(𝒙t\+𝒗\)−f​\(𝒙t\)=∇f​\(𝒙t\)⊤​𝒗\+O​\(‖𝒗‖2\)f\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\-f\(\\bm\{x\}\_\{t\}\)=\\nabla f\(\\bm\{x\}\_\{t\}\)^\{\\top\}\\bm\{v\}\+O\(\\\|\\bm\{v\}\\\|^\{2\}\)andp​\(𝒙t\+𝒗\)p​\(𝒙t\)=1\+𝒗⊤​∇log⁡p​\(𝒙t\)\+O​\(‖𝒗‖2\)\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}=1\+\\bm\{v\}^\{\\top\}\\nabla\\log p\(\\bm\{x\}\_\{t\}\)\+O\(\\\|\\bm\{v\}\\\|^\{2\}\)\. The leading first\-order term is odd invvand is canceled by the symmetry ofνt\\nu\_\{t\}in the principal\-value sense\. The remaining second\-order terms are locally integrable forα<2\\alpha<2\. This justifies the principal\-value form used in Proposition[1](https://arxiv.org/html/2608.10384#Thmproposition1)\.

Remark 2:We specify the explicit expression of the scaled Lévy measureνt​\(d​𝒗\)\\nu\_\{t\}\(d\\bm\{v\}\)to facilitate the following analysis\. For theDD\-dimensional case, the standard Lévy measure is defined asν​\(d​𝒗\)=CD,α​‖𝒗‖−D−α\\nu\(d\\bm\{v\}\)=C\_\{D,\\alpha\}\\\|\\bm\{v\}\\\|^\{\-D\-\\alpha\}\[[1](https://arxiv.org/html/2608.10384#bib.bib1)\], whereCD,α=α​2α−1​Γ​\(D\+α2\)πD2​Γ​\(1−α2\)C\_\{D,\\alpha\}=\\frac\{\\alpha 2^\{\\alpha\-1\}\\Gamma\(\\frac\{D\+\\alpha\}\{2\}\)\}\{\\pi^\{\\frac\{D\}\{2\}\}\\Gamma\(1\-\\frac\{\\alpha\}\{2\}\)\}\. Combined with the effect ofΦS​\(t\)=σS​\(t\)​𝑰D\\Phi\_\{S\}\(t\)=\\sigma\_\{S\}\(t\)\\bm\{I\}\_\{D\}, we haveνt​\(d​𝒗\)=\|σS​\(t\)\|α​ν​\(d​𝒗\)\\nu\_\{t\}\(d\\bm\{v\}\)=\|\\sigma\_\{S\}\(t\)\|^\{\\alpha\}\\nu\(d\\bm\{v\}\)\. In the following, we denoteCD,α,σS≜\|σS​\(t\)\|α​CD,αC\_\{D,\\alpha,\\sigma\_\{S\}\}\\triangleq\|\\sigma\_\{S\}\(t\)\|^\{\\alpha\}C\_\{D,\\alpha\}to avoid redundancy\.

From Proposition[1](https://arxiv.org/html/2608.10384#Thmproposition1), we can see that the backward kernel is related to the density ratiop​\(𝒙t\+𝒗\)p​\(𝒙t\)\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}, which contains global information\. Therefore, the jump part only belongs to a Markov\-type jump process rather than a Lévy process\. This process is much more general and is not as mathematically tractable as the forward process\. Consequently, the total reversed SDE corresponding to \([II](https://arxiv.org/html/2608.10384#S2.Ex1)\) cannot be expressed as the superposition of a modified drift term, a diffusion term, and a stable process\. This makes it challenging to design an efficient sampler directly, as in the diffusion case\.

Besides guiding the design of the reverse sampler, the backward generator is also useful for understanding discrete diffusion processes\. For example, in a standard DDPM driven by a Gaussian transition kernel,p​\(𝒙t−Δ​t\|𝒙t\)p\(\\bm\{x\}\_\{t\-\\Delta t\}\|\\bm\{x\}\_\{t\}\)can be modeled by a Gaussian distribution for the following reason\. First, with sufficient diffusion,𝒙T\\bm\{x\}\_\{T\}can be approximated by a Gaussian distribution\. Next, based on the backward Kolmogorov equation, the drift of the reversed SDE is−\(𝒃​\(𝒙t,t\)−g​\(t\)2​∇log⁡p​\(𝒙t,t\)\)​d​t\-\(\\bm\{b\}\(\\bm\{x\}\_\{t\},t\)\-g\(t\)^\{2\}\\nabla\\log p\(\\bm\{x\}\_\{t\},t\)\)dt, whereg​\(t\)g\(t\)is the coefficient of the Wiener process in the forward SDE\. Att=Tt=T, we have∇log⁡p​\(𝒙T,T\)∝−𝒙T\\nabla\\log p\(\\bm\{x\}\_\{T\},T\)\\propto\-\\bm\{x\}\_\{T\}\. Then, we discretize the continuous time interval with step sizeΔ​t\\Delta t\. In this case,𝒙T−Δ​t=𝒙T−\(𝒃​\(𝒙T,T\)\+𝒙T\)​Δ​t\+g​\(t\)​Δ​t​𝒘\\bm\{x\}\_\{T\-\\Delta t\}=\\bm\{x\}\_\{T\}\-\(\\bm\{b\}\(\\bm\{x\}\_\{T\},T\)\+\\bm\{x\}\_\{T\}\)\\Delta t\+g\(t\)\\sqrt\{\\Delta t\}\\bm\{w\}, where𝒘∼𝒩​\(𝟎,𝑰D\)\\bm\{w\}\\sim\\mathcal\{N\}\(\\bm\{0\},\\bm\{I\}\_\{D\}\)and𝒃​\(𝒙t,t\)\\bm\{b\}\(\\bm\{x\}\_\{t\},t\)is usually set as an affine function of𝒙t\\bm\{x\}\_\{t\}\. This recursive formula is equivalent to the summation of Gaussian\-distributed RVs\. Therefore,p​\(𝒙t−Δ​t\|𝒙t\)p\(\\bm\{x\}\_\{t\-\\Delta t\}\|\\bm\{x\}\_\{t\}\)can be treated as a Gaussian distribution, and the objective loss function can be significantly simplified by regressing the expectation of𝒙t,t:T→0\\bm\{x\}\_\{t\},t:T\\rightarrow 0\.

Unfortunately, according to Proposition[1](https://arxiv.org/html/2608.10384#Thmproposition1), this simplification does not hold for the process in \([II](https://arxiv.org/html/2608.10384#S2.Ex1)\)\. In the sequel, we will design a simple loss function that is more suitable for neural network learning\.

### III\-BReverse sampling decomposition

From the definition of the generator, given the jump kernelK​\(d​𝒗\)K\(d\\bm\{v\}\), the term∫ℝD\(f​\(𝒙\+𝒗\)−f​\(𝒙\)\)​K​\(d​𝒗\)\\int\_\{\\mathbb\{R\}^\{D\}\}\(f\(\\bm\{x\}\+\\bm\{v\}\)\-f\(\\bm\{x\}\)\)K\(d\\bm\{v\}\)indicates that the jump amplitude follows the distributionK​\(d​𝒗\)/∫ℝDK​\(d​𝒗\)K\(d\\bm\{v\}\)/\\int\_\{\\mathbb\{R\}^\{D\}\}K\(d\\bm\{v\}\)\. However, the number of small jumps with‖𝒗‖≤ϵ\\\|\\bm\{v\}\\\|\\leq\\epsilonis usually infinite due to the singularity of the Lévy measure at the origin\. Therefore, we handle small and large jumps separately\. Specifically, we decompose the original integral by truncation with thresholdϵ\\epsilon:

∫ℝD\(f​\(𝒙\+𝒗\)−f​\(𝒙\)\)​K​\(d​𝒗\)\\displaystyle\\int\_\{\\mathbb\{R\}^\{D\}\}\(f\(\\bm\{x\}\+\\bm\{v\}\)\-f\(\\bm\{x\}\)\)K\(d\\bm\{v\}\)=\\displaystyle=∫‖𝒗‖≤ϵ\(f​\(𝒙\+𝒗\)−f​\(𝒙\)\)​K​\(d​𝒗\)\\displaystyle\\int\_\{\\\|\\bm\{v\}\\\|\\leq\\epsilon\}\(f\(\\bm\{x\}\+\\bm\{v\}\)\-f\(\\bm\{x\}\)\)K\(d\\bm\{v\}\)\+∫‖𝒗‖\>ϵ\(f​\(𝒙\+𝒗\)−f​\(𝒙\)\)​K​\(d​𝒗\)\\displaystyle\+\\int\_\{\\\|\\bm\{v\}\\\|\>\\epsilon\}\(f\(\\bm\{x\}\+\\bm\{v\}\)\-f\(\\bm\{x\}\)\)K\(d\\bm\{v\}\)\(5\)
In other words, jumps with amplitude larger thanϵ\\epsilonare treated as large jumps\. For large jumps, the occurrence frequency is finite\. The number of large jumps follows a Poisson distribution with rateλ≜∫‖𝒗‖\>ϵK​\(d​𝒗\)\\lambda\\triangleq\\int\_\{\\\|\\bm\{v\}\\\|\>\\epsilon\}K\(d\\bm\{v\}\)\. For a sufficiently small time intervalΔ​t\\Delta t, the occurrence of a large jump can also be approximated by a Bernoulli distribution with probabilityλ​Δ​t\\lambda\\Delta t\. If a large jump occurs, we need to obtain one jump sample from the normalized distributionK​\(d​𝒗\)λ\\frac\{K\(d\\bm\{v\}\)\}\{\\lambda\}\. Otherwise, no large jump occurs within the intervalΔ​t\\Delta t, and the reversed process reduces to the standard diffusion process\.

For the case with‖𝒗‖≤ϵ\\\|\\bm\{v\}\\\|\\leq\\epsilon, the Lévy measure is not integrable, and there are infinitely many small jumps\. Unlike large jump activities, these small jumps cannot be enumerated explicitly\. The distribution of small jumps is also complicated, which does not follow a Gaussian or stable distribution due to amplitude truncation\. A simple method is to directly ignore small jump activities when the thresholdϵ\\epsilonis small enough\. However, this may introduce considerable approximation error, especially when the small jumps have a nonzero mean\. Meanwhile, this effect accumulates whenΔ​t\\Delta tis small\.

Here, we adopt a more reasonable approach\. After truncation, the small jumps have a finite second moment\. Thus, we replace the small jump generator by a Gaussian generator matching the first two local moments\. This is a weak approximation whose generator error is controlled by the third local moment of the truncated Lévy measure\. Namely, the small jump part within\[t,t\+Δ​t\]\[t,t\+\\Delta t\], denoted by𝑿t,Δ​t≤ϵ\\bm\{X\}\_\{t,\\Delta t\}^\{\\leq\\epsilon\}, can be written as

𝑿t,Δ​t≤ϵ≈𝒃t,Δ​t≤ϵ​Δ​t\+Δ​t​\(Σt,Δ​t≤ϵ\)12​𝑾\\bm\{X\}\_\{t,\\Delta t\}^\{\\leq\\epsilon\}\\approx\\bm\{b\}\_\{t,\\Delta t\}^\{\\leq\\epsilon\}\\Delta t\+\\sqrt\{\\Delta t\}\(\\Sigma\_\{t,\\Delta t\}^\{\\leq\\epsilon\}\)^\{\\frac\{1\}\{2\}\}\\bm\{W\}\(6\)where𝑾\\bm\{W\}is the standard Gaussian RV and

𝒃t,Δ​t≤ϵ=\\displaystyle\\bm\{b\}\_\{t,\\Delta t\}^\{\\leq\\epsilon\}=p\.v\.​∫‖𝒗‖≤ϵ𝒗​K​\(d​𝒗\)\\displaystyle\\operatorname\{p\.v\.\}\\int\_\{\\\|\\bm\{v\}\\\|\\leq\\epsilon\}\\bm\{v\}K\(d\\bm\{v\}\)\(7\)
Substituting the state\-dependent backward jump kernel into the small jump drift, we have

𝒃t,Δ​t≤ϵ\\displaystyle\\bm\{b\}\_\{t,\\Delta t\}^\{\\leq\\epsilon\}=p\.v\.​∫‖𝒗‖≤ϵ𝒗​p​\(𝒙t\+𝒗\)p​\(𝒙t\)​νt​\(d​𝒗\)\\displaystyle=\\operatorname\{p\.v\.\}\\int\_\{\\\|\\bm\{v\}\\\|\\leq\\epsilon\}\\bm\{v\}\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}\\nu\_\{t\}\(d\\bm\{v\}\)=∫‖𝒗‖≤ϵ𝒗​\(p​\(𝒙t\+𝒗\)p​\(𝒙t\)−1\)​νt​\(d​𝒗\)\\displaystyle=\\int\_\{\\\|\\bm\{v\}\\\|\\leq\\epsilon\}\\bm\{v\}\\left\(\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}\-1\\right\)\\nu\_\{t\}\(d\\bm\{v\}\)=∫‖𝒗‖≤ϵ𝒗​\(𝒗⊤​∇log⁡p​\(𝒙t\)\+Rp​\(𝒙t,𝒗,t\)\)​νt​\(d​𝒗\)\\displaystyle=\\int\_\{\\\|\\bm\{v\}\\\|\\leq\\epsilon\}\\bm\{v\}\\left\(\\bm\{v\}^\{\\top\}\\nabla\\log p\(\\bm\{x\}\_\{t\}\)\+R\_\{p\}\(\\bm\{x\}\_\{t\},\\bm\{v\},t\)\\right\)\\nu\_\{t\}\(d\\bm\{v\}\)=Aνt​∇log⁡p​\(𝒙t\)\+Rb,ϵ​\(𝒙t,t\),\\displaystyle=A\_\{\\nu\_\{t\}\}\\nabla\\log p\(\\bm\{x\}\_\{t\}\)\+R\_\{b,\\epsilon\}\(\\bm\{x\}\_\{t\},t\),\(8\)where the truncated second\-order moment matrix is defined asAνt≜∫‖𝒗‖≤ϵ𝒗​𝒗⊤​νt​\(d​𝒗\)A\_\{\\nu\_\{t\}\}\\triangleq\\int\_\{\\\|\\bm\{v\}\\\|\\leq\\epsilon\}\\bm\{v\}\\bm\{v\}^\{\\top\}\\nu\_\{t\}\(d\\bm\{v\}\)and the Taylor remainder isRb,ϵ​\(𝒙t,t\)≜∫‖𝒗‖≤ϵ𝒗​Rp​\(𝒙t,𝒗,t\)​νt​\(d​𝒗\)R\_\{b,\\epsilon\}\(\\bm\{x\}\_\{t\},t\)\\triangleq\\int\_\{\\\|\\bm\{v\}\\\|\\leq\\epsilon\}\\bm\{v\}R\_\{p\}\(\\bm\{x\}\_\{t\},\\bm\{v\},t\)\\nu\_\{t\}\(d\\bm\{v\}\)\. Here, the second equality follows from the symmetry ofνt\\nu\_\{t\}, which gives

p\.v\.​∫‖𝒗‖≤ϵ𝒗​νt​\(d​𝒗\)=𝟎\.\\displaystyle\\operatorname\{p\.v\.\}\\int\_\{\\\|\\bm\{v\}\\\|\\leq\\epsilon\}\\bm\{v\}\\nu\_\{t\}\(d\\bm\{v\}\)=\\bm\{0\}\.\(9\)
Ifp​\(𝒙t\)p\(\\bm\{x\}\_\{t\}\)is locally smooth and strictly positive, then

p​\(𝒙t\+𝒗\)p​\(𝒙t\)=1\+𝒗⊤\\displaystyle\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}=1\+\\bm\{v\}^\{\\top\}∇log⁡p​\(𝒙t\)\+Rp​\(𝒙t,𝒗,t\),\\displaystyle\\nabla\\log p\(\\bm\{x\}\_\{t\}\)\+R\_\{p\}\(\\bm\{x\}\_\{t\},\\bm\{v\},t\),\|Rp​\(𝒙t,𝒗,t\)\|≤C​‖𝒗‖2\.\\displaystyle\|R\_\{p\}\(\\bm\{x\}\_\{t\},\\bm\{v\},t\)\|\\leq C\\\|\\bm\{v\}\\\|^\{2\}\.\(10\)
For the isotropicα\\alpha\-stable Lévy measureνt​\(d​𝒗\)=C​‖𝒗‖−D−α​d​𝒗\\nu\_\{t\}\(d\\bm\{v\}\)=C\\\|\\bm\{v\}\\\|^\{\-D\-\\alpha\}d\\bm\{v\}, the drift remainder satisfies

‖Rb,ϵ​\(𝒙t,t\)‖\\displaystyle\\\|R\_\{b,\\epsilon\}\(\\bm\{x\}\_\{t\},t\)\\\|≤C​∫‖𝒗‖≤ϵ‖𝒗‖3​νt​\(d​𝒗\)\\displaystyle\\leq C\\int\_\{\\\|\\bm\{v\}\\\|\\leq\\epsilon\}\\\|\\bm\{v\}\\\|^\{3\}\\nu\_\{t\}\(d\\bm\{v\}\)=C​ϵ3−α\.\\displaystyle=C\\epsilon^\{3\-\\alpha\}\.\(11\)
Therefore, by replacing the unknown score∇log⁡p​\(𝒙t\)\\nabla\\log p\(\\bm\{x\}\_\{t\}\)with the score networksθ1​\(𝒙t,t\)s^\{\\theta\_\{1\}\}\(\\bm\{x\}\_\{t\},t\), the practical approximation becomes

𝒃t,Δ​t≤ϵ≈Aνt​sθ1​\(𝒙t,t\)\.\\displaystyle\\bm\{b\}\_\{t,\\Delta t\}^\{\\leq\\epsilon\}\\approx A\_\{\\nu\_\{t\}\}s^\{\\theta\_\{1\}\}\(\\bm\{x\}\_\{t\},t\)\.\(12\)
Similarly, the covariance matrix of the small jump component is given by

𝚺t,Δ​t≤ϵ\\displaystyle\\bm\{\\Sigma\}\_\{t,\\Delta t\}^\{\\leq\\epsilon\}=∫‖𝒗‖≤ϵ𝒗​𝒗⊤​\(1\+𝒗⊤​∇log⁡p​\(𝒙t\)\+Rp​\(𝒙t,𝒗,t\)\)​νt​\(d​𝒗\)\\displaystyle=\\int\_\{\\\|\\bm\{v\}\\\|\\leq\\epsilon\}\\bm\{v\}\\bm\{v\}^\{\\top\}\\left\(1\+\\bm\{v\}^\{\\top\}\\nabla\\log p\(\\bm\{x\}\_\{t\}\)\+R\_\{p\}\(\\bm\{x\}\_\{t\},\\bm\{v\},t\)\\right\)\\nu\_\{t\}\(d\\bm\{v\}\)=Aνt\+RΣ,ϵ​\(𝒙t,t\)\.\\displaystyle=A\_\{\\nu\_\{t\}\}\+R\_\{\\Sigma,\\epsilon\}\(\\bm\{x\}\_\{t\},t\)\.\(13\)where the remaining Taylor term satisfies

‖RΣ,ϵ​\(𝒙t,t\)‖\\displaystyle\\\|R\_\{\\Sigma,\\epsilon\}\(\\bm\{x\}\_\{t\},t\)\\\|≤C​∫‖𝒗‖≤ϵ‖𝒗‖4​νt​\(d​𝒗\)=C​ϵ4−α\.\\displaystyle\\leq C\\int\_\{\\\|\\bm\{v\}\\\|\\leq\\epsilon\}\\\|\\bm\{v\}\\\|^\{4\}\\nu\_\{t\}\(d\\bm\{v\}\)=C\\epsilon^\{4\-\\alpha\}\.\(14\)
Thus, the covariance matrix can be approximated as

𝚺t,Δ​t≤ϵ≈Aνt\.\\displaystyle\\bm\{\\Sigma\}\_\{t,\\Delta t\}^\{\\leq\\epsilon\}\\approx A\_\{\\nu\_\{t\}\}\.\(15\)
Both \([III\-B](https://arxiv.org/html/2608.10384#S3.Ex9)\) and \([15](https://arxiv.org/html/2608.10384#S3.E15)\) admit closed\-form expressions based on the polar coordinate transform\. For instance,

Aνt=\\displaystyle A\_\{\\nu\_\{t\}\}=∫‖𝒗‖≤ϵ𝒗​𝒗⊤​νt​\(d​𝒗\)\\displaystyle\\int\_\{\\\|\\bm\{v\}\\\|\\leq\\epsilon\}\\bm\{v\}\\bm\{v\}^\{\\top\}\\nu\_\{t\}\(d\\bm\{v\}\)=\\displaystyle=\|𝕊D−1\|​CD,α,σS​𝑰DD​∫0ϵr1−α​𝑑r\\displaystyle\\frac\{\|\\mathbb\{S\}\_\{D\-1\}\|C\_\{D,\\alpha,\\sigma\_\{S\}\}\\bm\{I\}\_\{D\}\}\{D\}\\int\_\{0\}^\{\\epsilon\}r^\{1\-\\alpha\}dr=\\displaystyle=2​πD2​CD,α,σS​ϵ2−αD​Γ​\(D2\)​\(2−α\)​𝑰D,\\displaystyle\\frac\{2\\pi^\{\\frac\{D\}\{2\}\}C\_\{D,\\alpha,\\sigma\_\{S\}\}\\epsilon^\{2\-\\alpha\}\}\{D\\Gamma\(\\frac\{D\}\{2\}\)\(2\-\\alpha\)\}\\bm\{I\}\_\{D\},\(16\)where𝑰D\\bm\{I\}\_\{D\}denotes theDD\-dimentional identity matrix and𝕊D−1\\mathbb\{S\}\_\{D\-1\}represents theD−1D\-1\-dimensional sphere\. Note that this is a generator\-level weak approximation rather than an exact Gaussian representation of the small jump sum\. Denote the effect of large jumps on𝑿t\\bm\{X\}\_\{t\}by𝑿t,Δ​t\>ϵ\\bm\{X\}\_\{t,\\Delta t\}^\{\>\\epsilon\}\. Then, the total jump influence during the reverse sampling process is𝑿t,Δ​t≤ϵ\+𝑿t,Δ​t\>ϵ\\bm\{X\}\_\{t,\\Delta t\}^\{\\leq\\epsilon\}\+\\bm\{X\}\_\{t,\\Delta t\}^\{\>\\epsilon\}\. Finally, by combining the externally applied drift term and the reversed diffusion part, the recursive relation between𝑿t−Δ​t\\bm\{X\}\_\{t\-\\Delta t\}and𝑿t\\bm\{X\}\_\{t\}can be constructed\.

## IVInverse Sampling

According to the above analysis, the main challenge of inverse sampling comes from large jump sampling\. Therefore, this section mainly focuses on how to generate large jump samples based on the current system state and time instant\. First, we explain how to obtain large jump samples\. We also show that, from a theoretical perspective, a neural network is not necessary for reverse sampling\. These observations help determine the learning target and design the corresponding network structure\. Then, two techniques are introduced to effectively reduce the complexity\. Finally, since additional conditions are often required to control the generation process in practical applications, we also discuss an approximate observation\-guided version for inverse problems\.

Before proceeding, we first describe a special family of functions that is important for efficient learning\.

In many scenarios, we may want to train a network to regress a targetf​\(𝒙\)f\(\\bm\{x\}\)that is not directly accessible\. However, we can efficiently evaluate only its conditional versionf​\(𝒙\|𝒙0\)f\(\\bm\{x\}\|\\bm\{x\}\_\{0\}\)\. Thus, we aim to establish the condition under which optimizing the neural network based onf​\(𝒙\|𝒙0\)f\(\\bm\{x\}\|\\bm\{x\}\_\{0\}\)is equivalent to optimizing it based onf​\(𝒙\)f\(\\bm\{x\}\)\. This is the core of the following theorem\.

###### Theorem 1\(Marginal trainable functions\)

Define the accurate target and conditional target asl​\(𝐱\):ℝD→ℝD~l\(\\bm\{x\}\):\\mathbb\{R\}^\{D\}\\rightarrow\\mathbb\{R\}^\{\\tilde\{D\}\}andl​\(𝐱\|𝐱0\):ℝD→ℝD~l\(\\bm\{x\}\|\\bm\{x\}\_\{0\}\):\\mathbb\{R\}^\{D\}\\rightarrow\\mathbb\{R\}^\{\\tilde\{D\}\}, respectively\. Denote the neural network byfθ​\(𝐱\)f^\{\\theta\}\(\\bm\{x\}\), which is used to regressl​\(𝐱\)l\(\\bm\{x\}\)\. Then, the optimization problemarg⁡minθ⁡𝔼𝐱∼p​\(𝐱\)​\[‖fθ​\(𝐱\)−l​\(𝐱\)‖2\]\\arg\\min\_\{\\theta\}\\mathbb\{E\}\_\{\\bm\{x\}\\sim p\(\\bm\{x\}\)\}\[\\\|f^\{\\theta\}\(\\bm\{x\}\)\-l\(\\bm\{x\}\)\\\|^\{2\}\]is equivalent toargminθ𝔼𝐱,𝐱0∼p​\(𝐱,𝐱0\)\[∥fθ\(𝐱\)−l\(𝐱\|𝐱0\)∥2\]\\arg\\min\_\{\\theta\}\\mathbb\{E\}\_\{\\bm\{x\},\\bm\{x\}\_\{0\}\\sim p\(\\bm\{x\},\\bm\{x\}\_\{0\}\)\}\[\\\|f^\{\\theta\}\(\\bm\{x\}\)\-l\(\\bm\{x\}\|\\bm\{x\}\_\{0\}\)\\\|^\{2\}\]if

l​\(𝒙\)=∫l​\(𝒙\|𝒙0\)​p​\(𝒙0\|𝒙\)​𝑑𝒙0l\(\\bm\{x\}\)=\\int l\(\\bm\{x\}\|\\bm\{x\}\_\{0\}\)p\(\\bm\{x\}\_\{0\}\|\\bm\{x\}\)d\\bm\{x\}\_\{0\}\(17\)and the functionsl​\(𝐱\)l\(\\bm\{x\}\)satisfying \([17](https://arxiv.org/html/2608.10384#S4.E17)\) are referred to as marginally trainable functions\.

###### Proof:

The original optimization problem can be reformulated as follows:

𝔼𝒙∼p​\(𝒙\)​\[‖fθ​\(𝒙\)−l​\(𝒙\)‖2\]\\displaystyle\\mathbb\{E\}\_\{\\bm\{x\}\\sim p\(\\bm\{x\}\)\}\[\\\|f^\{\\theta\}\(\\bm\{x\}\)\-l\(\\bm\{x\}\)\\\|^\{2\}\]=\\displaystyle=𝔼𝒙∼p​\(𝒙\)​\[‖fθ​\(𝒙\)‖2\]−2​𝔼𝒙∼p​\(𝒙\)​\[fθ​\(𝒙\)​l​\(𝒙\)\]\\displaystyle\\mathbb\{E\}\_\{\\bm\{x\}\\sim p\(\\bm\{x\}\)\}\[\\\|f^\{\\theta\}\(\\bm\{x\}\)\\\|^\{2\}\]\-2\\mathbb\{E\}\_\{\\bm\{x\}\\sim p\(\\bm\{x\}\)\}\[f^\{\\theta\}\(\\bm\{x\}\)l\(\\bm\{x\}\)\]\+𝔼𝒙∼p​\(𝒙\)​\[‖l​\(𝒙\)‖2\]\\displaystyle\+\\mathbb\{E\}\_\{\\bm\{x\}\\sim p\(\\bm\{x\}\)\}\[\\\|l\(\\bm\{x\}\)\\\|^\{2\}\]\(18\)
The third term in \([IV](https://arxiv.org/html/2608.10384#S4.Ex18)\) is independent ofθ\\thetaand can be ignored\. For the second term,

𝔼𝒙∼p​\(𝒙\)​\[fθ​\(𝒙\)​l​\(𝒙\)\]=\\displaystyle\\mathbb\{E\}\_\{\\bm\{x\}\\sim p\(\\bm\{x\}\)\}\[f^\{\\theta\}\(\\bm\{x\}\)l\(\\bm\{x\}\)\]=∫fθ​\(𝒙\)​∫l​\(𝒙\|𝒙0\)​p​\(𝒙0\|𝒙\)​𝑑𝒙0​p​\(𝒙\)​𝑑𝒙\\displaystyle\\int f^\{\\theta\}\(\\bm\{x\}\)\\int l\(\\bm\{x\}\|\\bm\{x\}\_\{0\}\)p\(\\bm\{x\}\_\{0\}\|\\bm\{x\}\)d\\bm\{x\}\_\{0\}p\(\\bm\{x\}\)d\\bm\{x\}=\\displaystyle=∫∫fθ​\(𝒙\)​l​\(𝒙\|𝒙0\)​p​\(𝒙0\|𝒙\)​p​\(𝒙\)​𝑑𝒙0​𝑑𝒙\\displaystyle\\int\\int f^\{\\theta\}\(\\bm\{x\}\)l\(\\bm\{x\}\|\\bm\{x\}\_\{0\}\)p\(\\bm\{x\}\_\{0\}\|\\bm\{x\}\)p\(\\bm\{x\}\)d\\bm\{x\}\_\{0\}d\\bm\{x\}=\\displaystyle=𝔼𝒙,𝒙0∼p​\(𝒙,𝒙0\)​\[fθ​\(𝒙\)​l​\(𝒙\|𝒙0\)\]\\displaystyle\\mathbb\{E\}\_\{\\bm\{x\},\\bm\{x\}\_\{0\}\\sim p\(\\bm\{x\},\\bm\{x\}\_\{0\}\)\}\[f^\{\\theta\}\(\\bm\{x\}\)l\(\\bm\{x\}\|\\bm\{x\}\_\{0\}\)\]\(19\)
Plugging \([IV](https://arxiv.org/html/2608.10384#S4.Ex20)\) into \([IV](https://arxiv.org/html/2608.10384#S4.Ex18)\) completes the proof\. ∎

Theorem[1](https://arxiv.org/html/2608.10384#Thmtheorem1)provides sufficient conditions under which a wide range of functions can be used for network training\. In fact, Theorem[1](https://arxiv.org/html/2608.10384#Thmtheorem1)is general since it includes conventional training frameworks such as score matching and flow matching\. For example, the target of score matching is the score function∇log⁡p​\(𝒙t\)\\nabla\\log p\(\\bm\{x\}\_\{t\}\), which is not accessible in practice\. Therefore, the conditional score∇log⁡p​\(𝒙t\|𝒙0\)\\nabla\\log p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)is usually used, together with the simple underlying relationship between𝒙t\\bm\{x\}\_\{t\}and𝒙0\\bm\{x\}\_\{0\}\. It can be verified that

∇log⁡p​\(𝒙t\)=\\displaystyle\\nabla\\log p\(\\bm\{x\}\_\{t\}\)=1p​\(𝒙t\)​∇​∫p​\(𝒙t\|𝒙0\)​p​\(𝒙0\)​𝑑𝒙0\\displaystyle\\frac\{1\}\{p\(\\bm\{x\}\_\{t\}\)\}\\nabla\\int p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)p\(\\bm\{x\}\_\{0\}\)d\\bm\{x\}\_\{0\}=\\displaystyle=∫∇p​\(𝒙t\|𝒙0\)p​\(𝒙t\)​p​\(𝒙t\|𝒙0\)​p​\(𝒙0\)​p​\(𝒙t\|𝒙0\)​𝑑𝒙0\\displaystyle\\int\\frac\{\\nabla p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\}\{p\(\\bm\{x\}\_\{t\}\)p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\}p\(\\bm\{x\}\_\{0\}\)p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)d\\bm\{x\}\_\{0\}=\\displaystyle=∫∇p​\(𝒙t\|𝒙0\)p​\(𝒙t\|𝒙0\)​p​\(𝒙0\|𝒙t\)​𝑑𝒙0\\displaystyle\\int\\frac\{\\nabla p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\}\{p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\}p\(\\bm\{x\}\_\{0\}\|\\bm\{x\}\_\{t\}\)d\\bm\{x\}\_\{0\}=\\displaystyle=∫∇log⁡p​\(𝒙t\|𝒙0\)​p​\(𝒙0\|𝒙t\)​𝑑𝒙0\\displaystyle\\int\\nabla\\log p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)p\(\\bm\{x\}\_\{0\}\|\\bm\{x\}\_\{t\}\)d\\bm\{x\}\_\{0\}\(20\)
Hence,∇log⁡p​\(𝒙t\)\\nabla\\log p\(\\bm\{x\}\_\{t\}\)belongs to the class of marginally trainable functions\. Similar derivations apply to the training of the conditional vector field in flow matching\[[15](https://arxiv.org/html/2608.10384#bib.bib25)\]\.

### IV\-ASampling of large jumps

Suppose that the dataset𝒙0,j,j=1,⋯,N\\bm\{x\}\_\{0,j\},j=1,\\cdots,Nis available\. According to the forward SDE, the expressions ofp​\(𝒙t\|𝒙0\)p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)andp​\(𝒙t\+𝒗\|𝒙0\)p​\(𝒙t\|𝒙0\)\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{x\}\_\{0\}\)\}\{p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\}can usually be efficiently obtained\. Consequently,

p​\(𝒙t\+𝒗\)p​\(𝒙t\)=\\displaystyle\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}=∫p​\(𝒙t\+𝒗\|𝒙0\)p​\(𝒙t\|𝒙0\)​p​\(𝒙0\|𝒙t\)​𝑑𝒙0\\displaystyle\\int\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{x\}\_\{0\}\)\}\{p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\}p\(\\bm\{x\}\_\{0\}\|\\bm\{x\}\_\{t\}\)d\\bm\{x\}\_\{0\}=\\displaystyle=∫p​\(𝒙t\+𝒗\|𝒙0\)p​\(𝒙t\|𝒙0\)​p​\(𝒙0\)​p​\(𝒙t\|𝒙0\)p​\(𝒙t\)​𝑑𝒙0\\displaystyle\\int\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{x\}\_\{0\}\)\}\{p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\}\\frac\{p\(\\bm\{x\}\_\{0\}\)p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}d\\bm\{x\}\_\{0\}=\\displaystyle=∫p​\(𝒙t\+𝒗\|𝒙0\)p​\(𝒙t\|𝒙0\)​w​\(𝒙t,𝒙0\)​p​\(𝒙0\)​𝑑𝒙0\\displaystyle\\int\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{x\}\_\{0\}\)\}\{p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\}w\(\\bm\{x\}\_\{t\},\\bm\{x\}\_\{0\}\)p\(\\bm\{x\}\_\{0\}\)d\\bm\{x\}\_\{0\}≈\\displaystyle\\approx1N​∑j=1Np​\(𝒙t\+𝒗\|𝒙0,j\)p​\(𝒙t\|𝒙0,j\)​w​\(𝒙t,𝒙0,j\)\\displaystyle\\frac\{1\}\{N\}\\sum\_\{j=1\}^\{N\}\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{x\}\_\{0,j\}\)\}\{p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0,j\}\)\}w\(\\bm\{x\}\_\{t\},\\bm\{x\}\_\{0,j\}\)\(21\)where we define

w​\(𝒙t,𝒙0\)≜\\displaystyle w\(\\bm\{x\}\_\{t\},\\bm\{x\}\_\{0\}\)\\triangleqp​\(𝒙t\|𝒙0\)p​\(𝒙t\)=p​\(𝒙t\|𝒙0\)∫p​\(𝒙t\|𝒙0\)​p​\(𝒙0\)​𝑑𝒙0\\displaystyle\\frac\{p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}=\\frac\{p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\}\{\\int p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)p\(\\bm\{x\}\_\{0\}\)d\\bm\{x\}\_\{0\}\}≈\\displaystyle\\approxp​\(𝒙t\|𝒙0\)1N​∑j=1Np​\(𝒙t\|𝒙0,j\)\\displaystyle\\frac\{p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\}\{\\frac\{1\}\{N\}\\sum\_\{j=1\}^\{N\}p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0,j\}\)\}\(22\)
As for the overall jump rate, we have

λ​\(𝒙t\)≜\\displaystyle\\lambda\(\\bm\{x\}\_\{t\}\)\\triangleq∫‖𝒗‖\>ϵp​\(𝒙t\+𝒗\)p​\(𝒙t\)​νt​\(d​𝒗\)\\displaystyle\\int\_\{\\\|\\bm\{v\}\\\|\>\\epsilon\}\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}\\nu\_\{t\}\(d\\bm\{v\}\)=\\displaystyle=1∑j=1Np​\(𝒙t\|𝒙0,j\)​∑j=1N∫‖𝒗‖\>ϵp​\(𝒙t\+𝒗\|𝒙0,j\)​νt​\(d​𝒗\)\\displaystyle\\frac\{1\}\{\\sum\_\{j=1\}^\{N\}p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0,j\}\)\}\\sum\_\{j=1\}^\{N\}\\int\_\{\\\|\\bm\{v\}\\\|\>\\epsilon\}p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{x\}\_\{0,j\}\)\\nu\_\{t\}\(d\\bm\{v\}\)\(23\)
Finally, the unnormalized jump distribution is written as follows,

q​\(𝒙t,𝒗\)≜\\displaystyle q\(\\bm\{x\}\_\{t\},\\bm\{v\}\)\\triangleqp​\(𝒙t\+𝒗\)p​\(𝒙t\)​νt​\(d​𝒗\)\\displaystyle\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}\\nu\_\{t\}\(d\\bm\{v\}\)≈\\displaystyle\\approx∑j=1Np​\(𝒙t\+𝒗\|𝒙0,j\)∑j=1Np​\(𝒙t\|𝒙0,j\)​νt​\(d​𝒗\),‖𝒗‖\>ϵ\\displaystyle\\frac\{\\sum\_\{j=1\}^\{N\}p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{x\}\_\{0,j\}\)\}\{\\sum\_\{j=1\}^\{N\}p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0,j\}\)\}\\nu\_\{t\}\(d\\bm\{v\}\),\\\|\\bm\{v\}\\\|\>\\epsilon\(24\)
The standard score function can also be directly approximated from the dataset based on \([IV](https://arxiv.org/html/2608.10384#S4.Ex22)\) as follows:

∇log⁡p​\(𝒙t\)=\\displaystyle\\nabla\\log p\(\\bm\{x\}\_\{t\}\)=∫∇log⁡p​\(𝒙t\|𝒙0\)​p​\(𝒙0\|𝒙t\)​𝑑𝒙0\\displaystyle\\int\\nabla\\log p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)p\(\\bm\{x\}\_\{0\}\|\\bm\{x\}\_\{t\}\)d\\bm\{x\}\_\{0\}≈\\displaystyle\\approx1N​∑j=1Nw​\(𝒙t,𝒙0,j\)​∇p​\(𝒙t\|𝒙0,j\)p​\(𝒙t\|𝒙0,j\)\\displaystyle\\frac\{1\}\{N\}\\sum\_\{j=1\}^\{N\}\\frac\{w\(\\bm\{x\}\_\{t\},\\bm\{x\}\_\{0,j\}\)\\nabla p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0,j\}\)\}\{p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0,j\}\)\}\(25\)
With the above procedure, reverse sampling can be performed using the distributionp​\(𝒙t\|𝒙0\)p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\. However, this introduces an unacceptable computational burden during inference\. For example,w​\(𝒙t,𝒙0\)w\(\\bm\{x\}\_\{t\},\\bm\{x\}\_\{0\}\)andλ​\(𝒙t\)\\lambda\(\\bm\{x\}\_\{t\}\)need to be calculated at every iteration\. This remains time\-consuming, even though the explicit expression ofp​\(𝒙t\|𝒙0\)p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)and the integral in \([IV\-A](https://arxiv.org/html/2608.10384#S4.Ex29)\) are available\.

### IV\-BNetwork amortization

To make large jump sampling more efficient, we use a neural network to amortize part of the computational burden and accelerate inference\.

Based on the previous analysis, there are several possible learnable targets, including the density ratiop​\(𝒙t\+𝒗\)p​\(𝒙t\)\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}, the prior distributionp​\(𝒙t\)p\(\\bm\{x\}\_\{t\}\), and the logarithmic density ratiolog⁡p​\(𝒙t\+𝒗\)p​\(𝒙t\)\\log\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}\. However, these choices lead to numerical challenges during both learning and sampling\. For example, consider learningp​\(𝒙t\+𝒗\)p​\(𝒙t\)\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}\. Its advantage is that it is marginally trainable, so the learning target can be easily computed from its conditional version\. However, large jump sampling requires numerical integration with respect tod​𝒗d\\bm\{v\}to obtain the overall jump rate\. This is infeasible for high\-dimensional data, since the integrand value must be obtained from the network output\. Even though analytical expressions ofp​\(𝒙t\|𝒙0,j\)p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0,j\}\)and∫‖𝒗‖\>ϵp​\(𝒙t\+𝒗\|𝒙0,j\)​νt​\(d​𝒗\)\\int\_\{\\\|\\bm\{v\}\\\|\>\\epsilon\}p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{x\}\_\{0,j\}\)\\nu\_\{t\}\(d\\bm\{v\}\)can be derived,2​N2Nadditions are still required at each iteration\. In addition, numerical stability is not guaranteed whenp​\(𝒙t\)p\(\\bm\{x\}\_\{t\}\)is very small butp​\(𝒙t\+𝒗\)p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)is relatively large\. In this case,p​\(𝒙t\+𝒗\)p​\(𝒙t\)\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}behaves like an outlier with a large amplitude\.

Although this issue can be alleviated by applying the logarithm,log⁡p​\(𝒙t\+𝒗\)p​\(𝒙t\)\\log\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}does not satisfy the conditions in Theorem[1](https://arxiv.org/html/2608.10384#Thmtheorem1)\. Moreover, the integration difficulty still remains\. As forp​\(𝒙t\)p\(\\bm\{x\}\_\{t\}\), it requires the network to learn the absolute values of the global probability distribution, which is difficult, especially whenDDis large\. Note that the method in\[[27](https://arxiv.org/html/2608.10384#bib.bib22)\]is equivalent to directly learning reverse sampling, and its large jump part is generated by the target distribution\. This approach places a heavier burden on the network, making it more like a black box and difficult to control\.

Therefore, we choose to learn the large jump rate, i\.e\.,λ​\(𝒙t\)\\lambda\(\\bm\{x\}\_\{t\}\)\. This task is similar to parameter estimation for a specialized distribution\. From this perspective, it is more learnable thanp​\(𝒙t\)p\(\\bm\{x\}\_\{t\}\)andp​\(𝒙t\+𝒗\)p​\(𝒙t\)\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}\. Moreover, althoughλ​\(𝒙t\)\\lambda\(\\bm\{x\}\_\{t\}\)is the integral of the density ratio with respect to the Lévy measure, it still belongs to the class of marginally trainable functions\. Therefore, it can be efficiently trained usingλ​\(𝒙t\|𝒙0\)≜∫‖𝒗‖\>ϵp​\(𝒙t\+𝒗\|𝒙0\)p​\(𝒙t\|𝒙0\)​νt​\(d​𝒗\)\\lambda\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\\triangleq\\int\_\{\\\|\\bm\{v\}\\\|\>\\epsilon\}\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{x\}\_\{0\}\)\}\{p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\}\\nu\_\{t\}\(d\\bm\{v\}\)\.

Based on the above discussion, the next proposition provides the loss function for the jump part\. For the diffusion part, the standard score matching framework can be directly used\.

###### Proposition 2\(Loss function for jump activities\)

The loss function of the neural network for learning the overall jump rate is given as follows:

ℒ​\(θ\)=\\displaystyle\\mathcal\{L\}\(\\theta\)=𝔼t∼𝒰​\[0,1\],𝒙t∼p​\(𝒙t\|𝒙0\),𝒙0∼p​\(𝒙0\)\[∥gθ2\(𝒙t,t\)\\displaystyle\\mathbb\{E\}\_\{t\\sim\\mathcal\{U\}\[0,1\],\\bm\{x\}\_\{t\}\\sim p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\),\\bm\{x\}\_\{0\}\\sim p\(\\bm\{x\}\_\{0\}\)\}\\big\[\\\|g^\{\\theta\_\{2\}\}\(\\bm\{x\}\_\{t\},t\)−λ\(𝒙t\|𝒙0\)∥2\]\\displaystyle\-\\lambda\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\\\|^\{2\}\\big\]\(26\)

### IV\-CEfficient calculation of large jump rate

The remaining two main challenges in bridging the gap between theoretical analysis and practical applications are how to calculate the target value in \([2](https://arxiv.org/html/2608.10384#S4.Ex32)\) and how to rapidly sample a large jump from the distributionq​\(𝒙t,𝒗\)q\(\\bm\{x\}\_\{t\},\\bm\{v\}\)\. Here, we first focus on the first problem\.

So far, we have not imposed specific assumptions on the coefficients in \([II](https://arxiv.org/html/2608.10384#S2.Ex1)\), except that the drift coefficient depends on both the system state and the time index, while the scaling matrices of the Wiener process and the stable process depend only on the time index\. For the design of a specific reverse sampling algorithm, its expressions need to be further specified\.The principle is to make the learning process as simple as possible without sacrificing much generality\.Similar to conventional diffusion models and SDEs driven only by the Wiener process\[[21](https://arxiv.org/html/2608.10384#bib.bib4),[22](https://arxiv.org/html/2608.10384#bib.bib6)\], we make the following assumptions\.

###### Assumption 5

𝒃​\(𝑿t,t\)\\bm\{b\}\(\\bm\{X\}\_\{t\},t\)is an affine function of𝐗t\\bm\{X\}\_\{t\}, i\.e\.,𝐛​\(𝐗t,t\)=R​𝐗t\+𝐬\\bm\{b\}\(\\bm\{X\}\_\{t\},t\)=R\\bm\{X\}\_\{t\}\+\\bm\{s\}, which is similar to the Ornstein\-Uhlenbeck \(OU\) process\. Since𝐛​\(𝐗t,t\)\\bm\{b\}\(\\bm\{X\}\_\{t\},t\)is a time\-homogeneous drift, we rewrite it as𝐛​\(𝐗t\)\\bm\{b\}\(\\bm\{X\}\_\{t\}\)in the sequel\.

###### Assumption 6

ForΦG​\(t\)\\Phi\_\{G\}\(t\)andΦS​\(t\)\\Phi\_\{S\}\(t\), we restrict them to diagonal matrices, i\.e\.,ΦG​\(t\)=diag​\(σG,1​\(t\),⋯,σG,D​\(t\)\)\\Phi\_\{G\}\(t\)=\\text\{diag\}\(\\sigma\_\{G,1\}\(t\),\\cdots,\\sigma\_\{G,D\}\(t\)\)andΦS​\(t\)=diag​\(σS,1​\(t\),⋯,σS,D​\(t\)\)\\Phi\_\{S\}\(t\)=\\text\{diag\}\(\\sigma\_\{S,1\}\(t\),\\cdots,\\sigma\_\{S,D\}\(t\)\)\.

In the following, we rewrite𝒃​\(𝑿t,t\)\\bm\{b\}\(\\bm\{X\}\_\{t\},t\)as𝒃​\(𝑿t\)\\bm\{b\}\(\\bm\{X\}\_\{t\}\)because it is directly related only to𝑿t\\bm\{X\}\_\{t\}\. Under these configurations, the distribution of𝑿t\\bm\{X\}\_\{t\}generated from a given𝑿0\\bm\{X\}\_\{0\}can be explicitly determined\.

###### Proposition 3\(Distribution ofXt\\bm\{X\}\_\{t\}\)

Consider the forward SDE in \([II](https://arxiv.org/html/2608.10384#S2.Ex1)\)\. Let𝐛​\(𝐗t\)\\bm\{b\}\(\\bm\{X\}\_\{t\}\),ΦG​\(t\)\\Phi\_\{G\}\(t\), andΦS​\(t\)\\Phi\_\{S\}\(t\)satisfy Assumptions[5](https://arxiv.org/html/2608.10384#Thmassumption5)and[6](https://arxiv.org/html/2608.10384#Thmassumption6)\. Then,𝐗t=𝛍​\(t,𝐱0\)\+𝐆​\(t\)\+𝐒​\(t\)\\bm\{X\}\_\{t\}=\\bm\{\\mu\}\(t,\\bm\{x\}\_\{0\}\)\+\\bm\{G\}\(t\)\+\\bm\{S\}\(t\), where𝐆​\(t\)∼𝒩​\(𝟎,ΣG\)\\bm\{G\}\(t\)\\sim\\mathcal\{N\}\(\\bm\{0\},\\Sigma\_\{G\}\)and𝐒​\(t\)\\bm\{S\}\(t\)follows anα\\alpha\-stable distribution with characteristic functionϕS​\(t\)\\phi\_\{S\}\(t\)\. Moreover,

𝝁​\(t,𝒙0\)=exp⁡\(R​t\)​𝒙0\+exp⁡\(R​t\)​∫0texp⁡\(−R​l\)​𝒔​𝑑l\\bm\{\\mu\}\(t,\\bm\{x\}\_\{0\}\)=\\exp\(Rt\)\\bm\{x\}\_\{0\}\+\\exp\(Rt\)\\int\_\{0\}^\{t\}\\exp\(\-Rl\)\\bm\{s\}dl\(27\)ΣG​\(t\)=∫0texp⁡\(\(t−l\)​R\)​ΦG​\(l\)​ΦG​\(l\)⊤​exp⁡\(\(t−l\)​R⊤\)​𝑑l\\Sigma\_\{G\}\(t\)=\\int\_\{0\}^\{t\}\\exp\(\(t\-l\)R\)\\Phi\_\{G\}\(l\)\\Phi\_\{G\}\(l\)^\{\\top\}\\exp\(\(t\-l\)R^\{\\top\}\)dl\(28\)ϕS​\(𝒖,t\)=exp⁡\(−∫0t‖ΦS​\(l\)⊤​exp⁡\(\(t−l\)​R⊤\)​𝒖‖2α​𝑑l\)\\phi\_\{S\}\(\\bm\{u\},t\)=\\exp\\bigg\(\-\\int\_\{0\}^\{t\}\\big\\\|\\Phi\_\{S\}\(l\)^\{\\top\}\\exp\(\(t\-l\)R^\{\\top\}\)\\bm\{u\}\\big\\\|\_\{2\}^\{\\alpha\}dl\\bigg\)\(29\)

###### Proof:

The proof is relegated to Appendix[B](https://arxiv.org/html/2608.10384#A2)\. ∎

###### Corollary 1\(Distribution ofXt\\bm\{X\}\_\{t\}for special cases\)

Assume thatRR,ΦG​\(t\)\\Phi\_\{G\}\(t\), andΦS​\(t\)\\Phi\_\{S\}\(t\)are all scaled identity matrices, i\.e\.,R=R0​𝐈DR=R\_\{0\}\\bm\{I\}\_\{D\},ΦG​\(t\)=σG​\(t\)​𝐈D\\Phi\_\{G\}\(t\)=\\sigma\_\{G\}\(t\)\\bm\{I\}\_\{D\}, andΦS​\(t\)=σS​\(t\)​𝐈D\\Phi\_\{S\}\(t\)=\\sigma\_\{S\}\(t\)\\bm\{I\}\_\{D\}\. The other settings are the same as those in Proposition[3](https://arxiv.org/html/2608.10384#Thmproposition3)\. Then,ΣG​\(t\)=γG​\(t\)2​𝐈D\\Sigma\_\{G\}\(t\)=\\gamma\_\{G\}\(t\)^\{2\}\\bm\{I\}\_\{D\}, whereγG​\(t\)=∫0texp⁡\(2​R0​\(t−l\)\)​σG​\(l\)2​𝑑l\\gamma\_\{G\}\(t\)=\\sqrt\{\\int\_\{0\}^\{t\}\\exp\(2R\_\{0\}\(t\-l\)\)\\sigma\_\{G\}\(l\)^\{2\}dl\}\. Moreover,ϕS​\(𝐮,t\)=exp⁡\(−γα​\(t\)α​‖𝐮‖2α\)\\phi\_\{S\}\(\\bm\{u\},t\)=\\exp\(\-\\gamma\_\{\\alpha\}\(t\)^\{\\alpha\}\\\|\\bm\{u\}\\\|\_\{2\}^\{\\alpha\}\), whereγα​\(t\)=\(∫0t\|σS​\(l\)\|α​exp⁡\(α​R0​\(t−l\)\)​𝑑l\)1α\\gamma\_\{\\alpha\}\(t\)=\(\\int\_\{0\}^\{t\}\|\\sigma\_\{S\}\(l\)\|^\{\\alpha\}\\exp\(\\alpha R\_\{0\}\(t\-l\)\)dl\)^\{\\frac\{1\}\{\\alpha\}\}\.

Remark 3:From Proposition[3](https://arxiv.org/html/2608.10384#Thmproposition3), it can be observed that the final distribution of𝑿T\\bm\{X\}\_\{T\}is still related to𝑿0\\bm\{X\}\_\{0\}\. In practice, we can chooseTTandR<0R<0such thatexp⁡\(R​T\)​𝑿0≈0\\exp\(RT\)\\bm\{X\}\_\{0\}\\approx 0\. In this case, the original distribution of the inverse sampling can be set to the mixed distribution composed of Gaussian distribution andα\\alpha\-stable distribution\.

Proposition[3](https://arxiv.org/html/2608.10384#Thmproposition3)and Corollary[1](https://arxiv.org/html/2608.10384#Thmcorollary1)imply that ifRR,ΦG​\(t\)\\Phi\_\{G\}\(t\), andΦS​\(t\)\\Phi\_\{S\}\(t\)are diagonal matrices with different diagonal elements, the isotropic property of the distribution of𝑿t\\bm\{X\}\_\{t\}will be destroyed\. This makes the following derivations much more cumbersome\. In the remainder of this paper, we use the assumptions in Corollary[1](https://arxiv.org/html/2608.10384#Thmcorollary1)to facilitate the design of the reverse sampler\.

Recall that the first challenge comes from the integral operation\. During training, performing high\-dimensional numerical integration for every data sample is time\-consuming\. Our goal is to derive an analytical expression for∫‖𝒗‖\>ϵp​\(𝒙t\+𝒗\|𝒙0\)​νt​\(d​𝒗\)\\int\_\{\\\|\\bm\{v\}\\\|\>\\epsilon\}p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{x\}\_\{0\}\)\\nu\_\{t\}\(d\\bm\{v\}\)\. According to Proposition[3](https://arxiv.org/html/2608.10384#Thmproposition3),p​\(𝒙t\|𝒙0\)p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)is the distribution of an RV composed of a bias term, a Gaussian RV, and a stable RV, which are mutually independent\. Due to the general lack of an analytical PDF for theα\\alpha\-stable distribution, it is challenging to obtain an explicit expression forp​\(𝒙t\|𝒙0\)p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\. Therefore, PDF approximation is needed to facilitate the following analysis\.

LetXB∼𝒩​\(0,2​γg2\)X\_\{B\}\\sim\\mathcal\{N\}\(0,2\\gamma\_\{g\}^\{2\}\)andXS∼𝒮​\(α,0,γs,0\)X\_\{S\}\\sim\\mathcal\{S\}\(\\alpha,0,\\gamma\_\{s\},0\)be mutually independent\. Denote the PDFs ofXBX\_\{B\}andXSX\_\{S\}byfXB​\(x\)f\_\{X\_\{B\}\}\(x\)andfXS​\(x\)f\_\{X\_\{S\}\}\(x\), respectively\. Then, the exact PDF ofXM≜XB\+XSX\_\{M\}\\triangleq X\_\{B\}\+X\_\{S\}equals the convolution offXB​\(x\)f\_\{X\_\{B\}\}\(x\)andfXS​\(x\)f\_\{X\_\{S\}\}\(x\)\. Unfortunately, there is no general explicit expression forfXM​\(x\)f\_\{X\_\{M\}\}\(x\)\. In\[[23](https://arxiv.org/html/2608.10384#bib.bib26)\], an approximate PDF is proposed to describe the statistical model ofXMX\_\{M\}\. It can be regarded as a weighted sum of a Gaussian kernel and a tail component\. Specifically,

fXM​\(x\)≈f~XM​\(x\)=1k1​\[c1​g0​exp⁡\(−x24​γg2\)\+α​γs​Cαc2\+\|x\|α\+1\],f\_\{X\_\{M\}\}\(x\)\\approx\\tilde\{f\}\_\{X\_\{M\}\}\(x\)=\\frac\{1\}\{k\_\{1\}\}\\bigg\[c\_\{1\}g\_\{0\}\\exp\\bigg\(\-\\frac\{x^\{2\}\}\{4\\gamma\_\{g\}^\{2\}\}\\bigg\)\+\\frac\{\\alpha\\gamma\_\{s\}C\_\{\\alpha\}\}\{c\_\{2\}\+\|x\|^\{\\alpha\+1\}\}\\bigg\],\(30\)wherek1k\_\{1\}is the normalization factor\. The parameterc2≥0c\_\{2\}\\geq 0is used to avoid singularity at the origin and control the shape of the tail component\. The parametersγg\\gamma\_\{g\}andγs\\gamma\_\{s\}represent the scaling parameters of the Gaussian kernel and the tail component, respectively\. The quantitiesg0g\_\{0\}andCαC\_\{\\alpha\}are expressed by the above parameters in closed form\. It has been shown in\[[23](https://arxiv.org/html/2608.10384#bib.bib26)\]thatf~XM​\(𝒙\)\\tilde\{f\}\_\{X\_\{M\}\}\(\\bm\{x\}\)can accurately modelXMX\_\{M\}under various scenarios\.

The tail component in \([30](https://arxiv.org/html/2608.10384#S4.E30)\) is derived from the asymptotic behavior of theα\\alpha\-stable distribution\[[23](https://arxiv.org/html/2608.10384#bib.bib26),[18](https://arxiv.org/html/2608.10384#bib.bib3)\]\. Following a similar procedure, \([30](https://arxiv.org/html/2608.10384#S4.E30)\) can also be extended to high\-dimensional cases\[[17](https://arxiv.org/html/2608.10384#bib.bib27)\]\. For example, let𝑿B∼𝒩​\(𝟎,σ2​𝑰\)\\bm\{X\}\_\{B\}\\sim\\mathcal\{N\}\(\\bm\{0\},\\sigma^\{2\}\\bm\{I\}\),𝑿S=A​𝑿G\\bm\{X\}\_\{S\}=\\sqrt\{A\}\\bm\{X\}\_\{G\}, and𝑿M=𝑿B\+𝑿S\\bm\{X\}\_\{M\}=\\bm\{X\}\_\{B\}\+\\bm\{X\}\_\{S\}\. Here, we consider the sub\-Gaussian representation of the multivariate isotropic stable distribution, whereAAfollows a right\-skewedα\\alpha\-stable distribution\[[18](https://arxiv.org/html/2608.10384#bib.bib3)\]\. Note that the characteristic function of𝑿M\\bm\{X\}\_\{M\}isϕ𝑿M​\(𝒍\)=exp⁡\(−γg2​‖𝒍‖2−γsα​‖𝒍‖α\)\\phi\_\{\\bm\{X\}\_\{M\}\}\(\\bm\{l\}\)=\\exp\(\-\\gamma\_\{g\}^\{2\}\\\|\\bm\{l\}\\\|^\{2\}\-\\gamma\_\{s\}^\{\\alpha\}\\\|\\bm\{l\}\\\|^\{\\alpha\}\)\. Then, the PDF of𝑿M\\bm\{X\}\_\{M\}can be similarly approximated as\[[17](https://arxiv.org/html/2608.10384#bib.bib27)\]

f~𝑿M​\(𝒙\)=\\displaystyle\\tilde\{f\}\_\{\\bm\{X\}\_\{M\}\}\(\\bm\{x\}\)=1k2\[c1g0exp\(−‖𝒙‖24​γg2\)\\displaystyle\\frac\{1\}\{k\_\{2\}\}\\bigg\[c\_\{1\}g\_\{0\}\\exp\\bigg\(\-\\frac\{\\\|\\bm\{x\}\\\|^\{2\}\}\{4\\gamma\_\{g\}^\{2\}\}\\bigg\)\+α​2α−1​Γ​\(α\+D2\)​γsαπD/2​Γ​\(1−α2\)\(c2\+∥𝒙∥2\)−α\+D2\]\\displaystyle\+\\frac\{\\alpha 2^\{\\alpha\-1\}\\Gamma\(\\frac\{\\alpha\+D\}\{2\}\)\\gamma\_\{s\}^\{\\alpha\}\}\{\\pi^\{D/2\}\\Gamma\(1\-\\frac\{\\alpha\}\{2\}\)\}\(c\_\{2\}\+\\\|\\bm\{x\}\\\|^\{2\}\)^\{\-\\frac\{\\alpha\+D\}\{2\}\}\\bigg\]\(31\)wherek2k\_\{2\}denotes the normalization factor\.

With \([IV\-C](https://arxiv.org/html/2608.10384#S4.Ex33)\), the target in \([2](https://arxiv.org/html/2608.10384#S4.Ex32)\) can be further derived\. Sincep​\(𝒙t\|𝒙0,j\)p\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0,j\}\)is independent ofd​𝒗d\\bm\{v\}, we focus only on∫‖𝒗‖\>ϵp​\(𝒙t\+𝒗\|𝒙0\)​νt​\(d​𝒗\)\\int\_\{\\\|\\bm\{v\}\\\|\>\\epsilon\}p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{x\}\_\{0\}\)\\nu\_\{t\}\(d\\bm\{v\}\)in the sequel\. For notational simplicity, we derive explicit expressions for the following integrals:

P1​\(ϵ\)≜∫‖𝒗‖\>ϵexp⁡\(−‖𝒗\+𝝁‖2\)​νt​\(d​𝒗\)P\_\{1\}\(\\epsilon\)\\triangleq\\int\_\{\\\|\\bm\{v\}\\\|\>\\epsilon\}\\exp\(\-\\\|\\bm\{v\}\+\\bm\{\\mu\}\\\|^\{2\}\)\\nu\_\{t\}\(d\\bm\{v\}\)\(32\)P2​\(ϵ\)≜∫‖𝒗‖\>ϵ\(c2\+‖𝒗\+𝝁‖2\)−α\+D2​νt​\(d​𝒗\)P\_\{2\}\(\\epsilon\)\\triangleq\\int\_\{\\\|\\bm\{v\}\\\|\>\\epsilon\}\(c\_\{2\}\+\\\|\\bm\{v\}\+\\bm\{\\mu\}\\\|^\{2\}\)^\{\-\\frac\{\\alpha\+D\}\{2\}\}\\nu\_\{t\}\(d\\bm\{v\}\)\(33\)
Note that𝝁\\bm\{\\mu\}is introduced only to characterize the mismatch betweenp​\(𝒙t\+𝒗\|𝒙0\)p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{x\}\_\{0\}\)and the Lévy measure\. Moreover, the scaling parameterγg\\gamma\_\{g\}is omitted in \([32](https://arxiv.org/html/2608.10384#S4.E32)\) for brevity\. It does not denote the mean of𝑿t\\bm\{X\}\_\{t\}\. Direct numerical integration is infeasible in high\-dimensional settings\. To address this issue, we transformP1​\(ϵ\)P\_\{1\}\(\\epsilon\)andP2​\(ϵ\)P\_\{2\}\(\\epsilon\)into polar coordinates as follows:

P1​\(ϵ\)=\\displaystyle P\_\{1\}\(\\epsilon\)=2​πD−12​CD,α,σSΓ​\(D−12\)∫ϵ\+∞∫−11exp\(−\(r2\+∥𝝁∥2\\displaystyle\\frac\{2\\pi^\{\\frac\{D\-1\}\{2\}\}C\_\{D,\\alpha,\\sigma\_\{S\}\}\}\{\\Gamma\(\\frac\{D\-1\}\{2\}\)\}\\int\_\{\\epsilon\}^\{\+\\infty\}\\int\_\{\-1\}^\{1\}\\exp\(\-\(r^\{2\}\+\\\|\\bm\{\\mu\}\\\|^\{2\}\+2rz∥𝝁∥\)\)r−α−1\(1−z2\)D−32dzdr\\displaystyle\+2rz\\\|\\bm\{\\mu\}\\\|\)\)r^\{\-\\alpha\-1\}\(1\-z^\{2\}\)^\{\\frac\{D\-3\}\{2\}\}dzdr\(34\)P2​\(ϵ\)=\\displaystyle P\_\{2\}\(\\epsilon\)=2​πD−12​CD,α,σSΓ​\(D−12\)∫ϵ\+∞∫−11\(c2\+r2\+∥𝝁∥2\\displaystyle\\frac\{2\\pi^\{\\frac\{D\-1\}\{2\}\}C\_\{D,\\alpha,\\sigma\_\{S\}\}\}\{\\Gamma\(\\frac\{D\-1\}\{2\}\)\}\\int\_\{\\epsilon\}^\{\+\\infty\}\\int\_\{\-1\}^\{1\}\(c\_\{2\}\+r^\{2\}\+\\\|\\bm\{\\mu\}\\\|^\{2\}\+2rz∥𝝁∥\)−α\+D2r−α−1\(1−z2\)D−32dzdr\\displaystyle\+2rz\\\|\\bm\{\\mu\}\\\|\)^\{\-\\frac\{\\alpha\+D\}\{2\}\}r^\{\-\\alpha\-1\}\(1\-z^\{2\}\)^\{\\frac\{D\-3\}\{2\}\}dzdr\(35\)which does not admit a closed\-form expression in general\. However, the original integral has been significantly simplified into a two\-dimensional integral\. In fact, the inner integral with respect tozzin \([32](https://arxiv.org/html/2608.10384#S4.E32)\) and \([33](https://arxiv.org/html/2608.10384#S4.E33)\) can be further simplified\. First, we need the following two integral identities\[[31](https://arxiv.org/html/2608.10384#bib.bib2)\]:

∫−11exp⁡\(a​x\)\\displaystyle\\int\_\{\-1\}^\{1\}\\exp\(ax\)\(1−x2\)p​d​x\\displaystyle\(1\-x^\{2\}\)^\{p\}dx=\\displaystyle=π​Γ​\(p\+1\)​\(2​a−1\)p\+12​Ip\+12​\(a\),\\displaystyle\\sqrt\{\\pi\}\\Gamma\(p\+1\)\(2a^\{\-1\}\)^\{\\frac\{p\+1\}\{2\}\}I\_\{\\frac\{p\+1\}\{2\}\}\(a\),\(36\)∫01xa−1\\displaystyle\\int\_\{0\}^\{1\}x^\{a\-1\}\(1−x\)b−a−1​\(1−z​x\)−c​d​x\\displaystyle\(1\-x\)^\{b\-a\-1\}\(1\-zx\)^\{\-c\}dx=\\displaystyle=B​\(a,b−a\)​F12​\(c,a;b;z\),b\>a\>0\\displaystyle B\(a,b\-a\)\\prescript\{\}\{2\}\{F\}\_\{1\}\(c,a;b;z\),b\>a\>0\(37\)whereIa​\(x\)I\_\{a\}\(x\)denotes the modified Bessel function of the first kind andB​\(⋅\)B\(\\cdot\)denotes the Beta function\. ForP1​\(ϵ\)P\_\{1\}\(\\epsilon\), applying \([IV\-C](https://arxiv.org/html/2608.10384#S4.Ex36)\) yields

P1​\(ϵ\)=\\displaystyle P\_\{1\}\(\\epsilon\)=2​πD2​CD,α,σS​‖𝝁‖2−D2​exp⁡\(−‖𝝁‖2\)\\displaystyle 2\\pi^\{\\frac\{D\}\{2\}\}C\_\{D,\\alpha,\\sigma\_\{S\}\}\\\|\\bm\{\\mu\}\\\|^\{\\frac\{2\-D\}\{2\}\}\\exp\(\-\\\|\\bm\{\\mu\}\\\|^\{2\}\)×∫ϵ\+∞exp\(−r2\)r−α−D2ID−22\(2r∥𝝁∥\)dr\\displaystyle\\times\\int\_\{\\epsilon\}^\{\+\\infty\}\\exp\(\-r^\{2\}\)r^\{\-\\alpha\-\\frac\{D\}\{2\}\}I\_\{\\frac\{D\-2\}\{2\}\}\(2r\\\|\\bm\{\\mu\}\\\|\)dr\(38\)
P2​\(ϵ\)=\\displaystyle P\_\{2\}\(\\epsilon\)=CD,α,σS​πD−12​2D−1​B​\(D−12,D−12\)Γ​\(D−12\)​∫ϵ\+∞r−α−1​F12​\(α\+D2,D−12;D−1;−4​r​‖𝝁‖c2\+r2\+‖𝝁‖2−2​r​‖𝝁‖\)\(c2\+r2\+‖𝝁‖2−2​r​‖𝝁‖\)α\+D2​𝑑r\\displaystyle\\frac\{C\_\{D,\\alpha,\\sigma\_\{S\}\}\\pi^\{\\frac\{D\-1\}\{2\}\}2^\{D\-1\}B\(\\frac\{D\-1\}\{2\},\\frac\{D\-1\}\{2\}\)\}\{\\Gamma\(\\frac\{D\-1\}\{2\}\)\}\\int\_\{\\epsilon\}^\{\+\\infty\}\\frac\{r^\{\-\\alpha\-1\}\\prescript\{\}\{2\}\{F\}\_\{1\}\(\\frac\{\\alpha\+D\}\{2\},\\frac\{D\-1\}\{2\};D\-1;\-\\frac\{4r\\\|\\bm\{\\mu\}\\\|\}\{c\_\{2\}\+r^\{2\}\+\\\|\\bm\{\\mu\}\\\|^\{2\}\-2r\\\|\\bm\{\\mu\}\\\|\}\)\}\{\(c\_\{2\}\+r^\{2\}\+\\\|\\bm\{\\mu\}\\\|^\{2\}\-2r\\\|\\bm\{\\mu\}\\\|\)^\{\\frac\{\\alpha\+D\}\{2\}\}\}dr\(39\)
ForP2​\(ϵ\)P\_\{2\}\(\\epsilon\), we use the variable transformationx=2​t−1x=2t\-1\. After some manipulations, \([39](https://arxiv.org/html/2608.10384#S4.E39)\) can be obtained\. In this way, calculating the total jump rateλ​\(𝒙t\|𝒙0,j\)\\lambda\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0,j\}\)reduces to evaluating a one\-dimensional integral, which significantly decreases the computational complexity compared toDD\-dimensional integrals\.

### IV\-DHigh\-dimensional sampling

As for the second problem, numerical sampling algorithms such as Hamiltonian Monte Carlo \(HMC\) can be applied, as the PDF and its gradient information are available\. However, reverse sampling of the SDE is still slow due to the cold\-start problem\. Specifically, for each reverse sampling iteration with sufficiently smallΔ​t\\Delta t, only one sample is needed for a large jump if it occurs\. However, the target sampling distribution varies withtt\. In this case, the HMC sampler should be “retrained” repeatedly to reach the steady state of the dynamic system, which is quite inefficient\.

Then, we design the fast sampling method for the large jump\. Recall that the normalized target distribution is

q~​\(𝒙t,𝒗\)≜\\displaystyle\\tilde\{q\}\(\\bm\{x\}\_\{t\},\\bm\{v\}\)\\triangleqp​\(𝒙t\+𝒗\)​νt​\(d​𝒗\)Qt,ϵ\\displaystyle\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\\nu\_\{t\}\(d\\bm\{v\}\)\}\{Q\_\{t,\\epsilon\}\}≈\\displaystyle\\approx∑j=1Np​\(𝒙t\+𝒗\|𝒙0,j\)​νt​\(d​𝒗\)N​Qt,ϵ\\displaystyle\\frac\{\\sum\_\{j=1\}^\{N\}p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{x\}\_\{0,j\}\)\\nu\_\{t\}\(d\\bm\{v\}\)\}\{NQ\_\{t,\\epsilon\}\}\(40\)where we defineQt,ϵ=∫‖𝒗‖\>ϵp​\(𝒙t\+𝒗\)​νt​\(d​𝒗\)Q\_\{t,\\epsilon\}=\\int\_\{\\\|\\bm\{v\}\\\|\>\\epsilon\}p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\\nu\_\{t\}\(d\\bm\{v\}\)\. Note thatp​\(𝒙t\)p\(\\bm\{x\}\_\{t\}\)is omitted because it is not related to𝒗\\bm\{v\}, and our objective is to generate a sample from the distribution with respect to𝒗\\bm\{v\}\.

Different from the one\-dimensional case, we sample fromq~​\(𝒙t,𝒗\)\\tilde\{q\}\(\\bm\{x\}\_\{t\},\\bm\{v\}\)under polar coordinates\. First, we rewriteq~​\(𝒙t,𝒗\)\\tilde\{q\}\(\\bm\{x\}\_\{t\},\\bm\{v\}\)asq~​\(𝒙t,r,z,θ\)\\tilde\{q\}\(\\bm\{x\}\_\{t\},r,z,\\theta\), wherer≜‖𝒗‖r\\triangleq\\\|\\bm\{v\}\\\|,z≜⟨𝝁,𝒗⟩‖𝝁‖​‖𝒗‖z\\triangleq\\frac\{\\langle\\bm\{\\mu\},\\bm\{v\}\\rangle\}\{\\\|\\bm\{\\mu\}\\\|\\\|\\bm\{v\}\\\|\}, andθ∈ℛD−2\\theta\\in\\mathcal\{R\}^\{D\-2\}represents the remaining angular variables\. Then,q~​\(𝒙t,𝒗\)\\tilde\{q\}\(\\bm\{x\}\_\{t\},\\bm\{v\}\)can be rewritten as in \([IV\-D](https://arxiv.org/html/2608.10384#S4.Ex42)\), where𝝁\\bm\{\\mu\}is defined as before and incorporates the information of𝒙t\\bm\{x\}\_\{t\}\. We then need to samplerr,zz, andθ\\thetafrom the joint distributionq~​\(𝒙t,r,z,θ\)\\tilde\{q\}\(\\bm\{x\}\_\{t\},r,z,\\theta\)\. Givenrrandzz,θ\\thetacan be generated uniformly on the sphere orthogonal to the direction represented byzz\. Therefore,θ\\thetadoes not explicitly appear in the PDF\. Based on this observation, we omitθ\\thetainq~​\(𝒙t,r,z,θ\)\\tilde\{q\}\(\\bm\{x\}\_\{t\},r,z,\\theta\)for simplicity\.

Note that the derivations of the marginal distributions are the same as those ofP1​\(ϵ\)P\_\{1\}\(\\epsilon\)andP2​\(ϵ\)P\_\{2\}\(\\epsilon\)\. Hence, the previous lookup tables can be directly used for fast sampling based on inverse CDF methods\. In other words, we can first samplerrfrom the marginal distributionq~​\(𝒙t,r\)\\tilde\{q\}\(\\bm\{x\}\_\{t\},r\)and then samplezzfromq~​\(z\|𝒙t,r\)\\tilde\{q\}\(z\|\\bm\{x\}\_\{t\},r\)\. To distinguish the distribution under polar coordinates fromq~​\(𝒙t,𝒗\)\\tilde\{q\}\(\\bm\{x\}\_\{t\},\\bm\{v\}\), we denote the CDFs ofq~​\(𝒙t,r\)\\tilde\{q\}\(\\bm\{x\}\_\{t\},r\)andq~​\(z\|𝒙t,r\)\\tilde\{q\}\(z\|\\bm\{x\}\_\{t\},r\)byFq~,\(r,t\)​\(⋅\)F\_\{\\tilde\{q\},\(r,t\)\}\(\\cdot\)andFq~,z\|\(r,t\)​\(⋅\)F\_\{\\tilde\{q\},z\|\(r,t\)\}\(\\cdot\), respectively, in the sequel\.

However, even ifP1​\(x\)P\_\{1\}\(x\)andP2​\(x\)P\_\{2\}\(x\)withx∈\[ϵ,M\]x\\in\[\\epsilon,M\]and a large constantMMcan be precomputed and stored in the lookup table, the inference stage is still time\-consuming\. This is because, for each pair\(r,Fq~,\(r,t\)​\(⋅\)\)\(r,F\_\{\\tilde\{q\},\(r,t\)\}\(\\cdot\)\), the corresponding conditional pair\(r,Fq~\(\|\),\(r,t\)​\(⋅\)\)\(r,F\_\{\\tilde\{q\}\(\|\),\(r,t\)\}\(\\cdot\)\)needs to be obtainedNNtimes\. This complexity can be considerably reduced by another simple sampling step\. Indeed,Qt,ϵQ\_\{t,\\epsilon\}can be decomposed as

Qt,ϵ=1N​∑j=1N∫‖𝒗‖\>ϵp​\(𝒙t\+𝒗\|𝒙0,j\)​νt​\(d​𝒗\)≜1N​∑j=1NQt,ϵ,j\\displaystyle Q\_\{t,\\epsilon\}=\\frac\{1\}\{N\}\\sum\_\{j=1\}^\{N\}\\int\_\{\\\|\\bm\{v\}\\\|\>\\epsilon\}p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{x\}\_\{0,j\}\)\\nu\_\{t\}\(d\\bm\{v\}\)\\triangleq\\frac\{1\}\{N\}\\sum\_\{j=1\}^\{N\}Q\_\{t,\\epsilon,j\}\(41\)
Based on \([IV\-D](https://arxiv.org/html/2608.10384#S4.Ex39)\),q~​\(𝒙t,𝒗\)\\tilde\{q\}\(\\bm\{x\}\_\{t\},\\bm\{v\}\)can be rewritten as follows:

q~​\(𝒙t,𝒗\)=\\displaystyle\\tilde\{q\}\(\\bm\{x\}\_\{t\},\\bm\{v\}\)=∑j=1Np​\(𝒙t\+𝒗\|𝒙0,j\)​νt​\(d​𝒗\)∑j=1NQt,ϵ,j\\displaystyle\\frac\{\\sum\_\{j=1\}^\{N\}p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{x\}\_\{0,j\}\)\\nu\_\{t\}\(d\\bm\{v\}\)\}\{\\sum\_\{j=1\}^\{N\}Q\_\{t,\\epsilon,j\}\}=\\displaystyle=∑j=1NQt,ϵ,j​q~​\(𝒙t,𝒗\|𝒙0,j\)∑j=1NQt,ϵ,j\\displaystyle\\frac\{\\sum\_\{j=1\}^\{N\}Q\_\{t,\\epsilon,j\}\\tilde\{q\}\(\\bm\{x\}\_\{t\},\\bm\{v\}\|\\bm\{x\}\_\{0,j\}\)\}\{\\sum\_\{j=1\}^\{N\}Q\_\{t,\\epsilon,j\}\}≜\\displaystyle\\triangleq∑j=1NgQ​\(𝒙0,j\|𝒙t,t\)​q~​\(𝒙t,𝒗\|𝒙0,j\)\\displaystyle\\sum\_\{j=1\}^\{N\}g\_\{Q\}\(\\bm\{x\}\_\{0,j\}\|\\bm\{x\}\_\{t\},t\)\\tilde\{q\}\(\\bm\{x\}\_\{t\},\\bm\{v\}\|\\bm\{x\}\_\{0,j\}\)\(42\)where we define the weighted distributiongQ​\(𝒙\|𝒙t,t\)g\_\{Q\}\(\\bm\{x\}\|\\bm\{x\}\_\{t\},t\)as

gQ​\(𝒙\|𝒙t,t\)=∑j=1N1\{𝒙0,j\}​\(𝒙\)​Qt,ϵ,j∑j=1NQt,ϵ,jg\_\{Q\}\(\\bm\{x\}\|\\bm\{x\}\_\{t\},t\)=\\sum\_\{j=1\}^\{N\}\\frac\{1\_\{\\\{\\bm\{x\}\_\{0,j\}\\\}\}\(\\bm\{x\}\)Q\_\{t,\\epsilon,j\}\}\{\\sum\_\{j=1\}^\{N\}Q\_\{t,\\epsilon,j\}\}\(43\)
In this case,q~​\(𝒙t,𝒗\)\\tilde\{q\}\(\\bm\{x\}\_\{t\},\\bm\{v\}\)can be regarded as a mixture distribution of the conditional distributionsq~​\(𝒙t,𝒗\|𝒙0,j\)\\tilde\{q\}\(\\bm\{x\}\_\{t\},\\bm\{v\}\|\\bm\{x\}\_\{0,j\}\)\. Then,q~​\(𝒙t,𝒗\)\\tilde\{q\}\(\\bm\{x\}\_\{t\},\\bm\{v\}\)is equivalent to a distribution controlled by an underlying Markov process\. Based on the above observation, we first randomly choose a data sample𝒙0∗∈\{𝒙0,j,j=1,⋯,N\}\\bm\{x\}^\{\*\}\_\{0\}\\in\\\{\\bm\{x\}\_\{0,j\},j=1,\\cdots,N\\\}according to the distributiongQ​\(𝒙\|𝒙t,t\)g\_\{Q\}\(\\bm\{x\}\|\\bm\{x\}\_\{t\},t\)\. Then, we generate a sample from the distribution conditioned on𝒙0∗\\bm\{x\}^\{\*\}\_\{0\}\. In this way, we only need to focus on the conditional distributionq~​\(𝒙t\|𝒙0∗\)\\tilde\{q\}\(\\bm\{x\}\_\{t\}\|\\bm\{x\}^\{\*\}\_\{0\}\)and obtain the pair\(r,Fq~\(\|\),\(r,t\)​\(⋅\)\)\(r,F\_\{\\tilde\{q\}\(\|\),\(r,t\)\}\(\\cdot\)\)only once\. This is much faster than directly sampling fromq~​\(𝒙t,𝒗\)\\tilde\{q\}\(\\bm\{x\}\_\{t\},\\bm\{v\}\)\.

To obtain the completegQ​\(𝒙\|𝒙t,t\)g\_\{Q\}\(\\bm\{x\}\|\\bm\{x\}\_\{t\},t\), allQt,ϵ,jQ\_\{t,\\epsilon,j\}withj=1,⋯,Nj=1,\\cdots,Nshould be computed\. This is equivalent to searching the lookup tableNNtimes\. In fact, manyQt,ϵ,jQ\_\{t,\\epsilon,j\}may be rather small, and their influence can be neglected\. Therefore, instead of constructing the completegQ​\(𝒙\|𝒙t,t\)g\_\{Q\}\(\\bm\{x\}\|\\bm\{x\}\_\{t\},t\), we truncategQ​\(𝒙\|𝒙t,t\)g\_\{Q\}\(\\bm\{x\}\|\\bm\{x\}\_\{t\},t\)by choosing theKKsamples with the largest approximate mixture weights, or equivalently the largestQt,ϵ,jQ\_\{t,\\epsilon,j\}in the prior sampler\. This operation is reasonable because the contribution of𝒙0,j\\bm\{x\}\_\{0,j\}is determined by the overlap between the conditional densityp​\(𝒙t\+𝒗\|𝒙0,j\)p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{x\}\_\{0,j\}\)and the truncated Lévy measure\. In the isotropic case,νt​\(d​𝒗\)∝‖𝒗‖−D−α​d​𝒗\\nu\_\{t\}\(d\\bm\{v\}\)\\propto\\\|\\bm\{v\}\\\|^\{\-D\-\\alpha\}d\\bm\{v\}\. After truncation, it favors jumps with‖𝒗‖\\\|\\bm\{v\}\\\|close toϵ\\epsilon\. Therefore,Qt,ϵ,jQ\_\{t,\\epsilon,j\}is large only when the high\-density region ofp​\(𝒙t\+𝒗\|𝒙0,j\)p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{x\}\_\{0,j\}\)intersects the Lévy\-favored shell around𝒙t\\bm\{x\}\_\{t\}\. This motivates a top\-KKapproximation that keeps only the candidates with the largest approximateQt,ϵ,jQ\_\{t,\\epsilon,j\}and renormalizes their weights\. For convenience, we denote the index set by𝒦\\mathcal\{K\}and the truncated version ofgQ​\(𝒙\|𝒙t,t\)g\_\{Q\}\(\\bm\{x\}\|\\bm\{x\}\_\{t\},t\)byg~Q​\(𝒙\)\\tilde\{g\}\_\{Q\}\(\\bm\{x\}\)\.

Remark 4:If we choose the largestKKelements from allQt,ϵ,jQ\_\{t,\\epsilon,j\}via the trivial traversal, the complexity still remainsO​\(N\)O\(N\)\. To solve this problem, more advanced and efficient ordering algorithms can be utilized, such as approximate nearest\-neighbor \(ANN\) methods\.

q~​\(𝒙t,r,z,θ\)∝\\displaystyle\\tilde\{q\}\(\\bm\{x\}\_\{t\},r,z,\\theta\)\\propto\(c1​g0k2​exp⁡\(−14​γg2​\(r2\+‖𝝁‖2\+2​r​z​‖𝝁‖\)\)\+α​2α−1​Γ​\(α\+D2\)​γsαk2​πD/2​Γ​\(1−α2\)​\(c2\+r2\+‖𝝁‖2\+2​r​z​‖𝝁‖\)−α\+D2\)\\displaystyle\\bigg\(\\frac\{c\_\{1\}g\_\{0\}\}\{k\_\{2\}\}\\exp\\bigg\(\-\\frac\{1\}\{4\\gamma\_\{g\}^\{2\}\}\(r^\{2\}\+\\\|\\bm\{\\mu\}\\\|^\{2\}\+2rz\\\|\\bm\{\\mu\}\\\|\)\\bigg\)\+\\frac\{\\alpha 2^\{\\alpha\-1\}\\Gamma\(\\frac\{\\alpha\+D\}\{2\}\)\\gamma\_\{s\}^\{\\alpha\}\}\{k\_\{2\}\\pi^\{D/2\}\\Gamma\(1\-\\frac\{\\alpha\}\{2\}\)\}\(c\_\{2\}\+r^\{2\}\+\\\|\\bm\{\\mu\}\\\|^\{2\}\+2rz\\\|\\bm\{\\mu\}\\\|\)^\{\-\\frac\{\\alpha\+D\}\{2\}\}\\bigg\)×r−α−1\(1−z2\)D−32,r\>ϵ,z∈\[−1,1\],\\displaystyle\\times r^\{\-\\alpha\-1\}\(1\-z^\{2\}\)^\{\\frac\{D\-3\}\{2\}\},r\>\\epsilon,z\\in\[\-1,1\],\(44\)
In addition to the inverse CDF method, rejection sampling can also be efficient for sampling the radial variablerrbased on \([IV\-D](https://arxiv.org/html/2608.10384#S4.Ex42)\)\. This converts the originalDD\-dimensional sampling problem into two one\-dimensional sampling steps and one remaining direction sampling step\. Since \([IV\-D](https://arxiv.org/html/2608.10384#S4.Ex42)\) has an analytical PDF expression with explicit asymptotic behavior, the proposal distribution for samplingrrcan be chosen as a Student’s t\-distribution with2​α\+D2\\alpha\+Ddegrees of freedom\. This choice matches the tail behavior of the marginal distribution induced by \([IV\-D](https://arxiv.org/html/2608.10384#S4.Ex42)\)\. Note that marginalization with respect tozzdoes not change the tail order of the distribution in terms ofrr\. Its closed\-form expression has been analyzed in \([IV\-C](https://arxiv.org/html/2608.10384#S4.Ex34)\)∼\\sim\([39](https://arxiv.org/html/2608.10384#S4.E39)\)\. Forzz, rejection sampling is not necessary, since the proposal distribution is not straightforward to design due to the bounded domain, and the corresponding lookup table is quite small\.

Finally, the overall training procedure, the large jump sampling method, and the complete reverse sampling algorithm are summarized in Algorithm[1](https://arxiv.org/html/2608.10384#alg1), Algorithm[2](https://arxiv.org/html/2608.10384#alg2), and Algorithm[3](https://arxiv.org/html/2608.10384#alg3), respectively\. Note that complete reverse sampling starts from samples generated from the terminal prior distributionpt​p​\(𝒙\)p\_\{tp\}\(\\bm\{x\}\)\. In our setting,pt​p​\(𝒙\)p\_\{tp\}\(\\bm\{x\}\)is the convolution of a multivariate Gaussian distribution and a stable distribution\. The corresponding samples can be generated separately from these two distributions\. Here, if we set𝒔≠𝟎\\bm\{s\}\\neq\\bm\{0\}, the mean ofpt​p​\(𝒙\)p\_\{tp\}\(\\bm\{x\}\)is−𝒔R0\-\\frac\{\\bm\{s\}\}\{R\_\{0\}\}\. Without loss of generality, we can also set𝒔=𝟎\\bm\{s\}=\\bm\{0\}so thatpt​p​\(𝒙\)p\_\{tp\}\(\\bm\{x\}\)is symmetric about the origin\.

Remark 5:In Algorithm[2](https://arxiv.org/html/2608.10384#alg2), we only consider the case where𝝁t≠𝟎\\bm\{\\mu\}\_\{t\}\\neq\\bm\{0\}andD≥3D\\geq 3\. If𝝁t=𝟎\\bm\{\\mu\}\_\{t\}=\\bm\{0\}, the conditional density is radially symmetric\. Therefore, we only need to samplerr\. WhenD=1,2D=1,2, the sampling procedure can be significantly simplified\. For example, whenD=1D=1, the orthogonal direction does not exist, and onlyrrneeds to be sampled\.

Remark 6:We emphasize that the generator\-based reverse characterization is stated at the level of Lévy\-driven Markov processes, whereas the computational sampler developed in this paper is derived under additional structural assumptions\. In particular, the practical algorithm focuses on isotropic linear Lévy SDEs with symmetricα\\alpha\-stable jump measures, for which the conditional transition density, jump rate approximation, and polar coordinate sampling procedure can be implemented efficiently\. This distinction allows us to retain the general insight provided by the reversed generator while avoiding an overstatement of the scope of the proposed numerical sampler\.

Algorithm 1Training algorithm1:Dataset

𝒙0,j,j=1,⋯,N\\bm\{x\}\_\{0,j\},j=1,\\cdots,Nand the training epoch

NepochN\_\{\\text\{epoch\}\}\.

2:Well\-trained networks

sθ1​\(⋅\)s^\{\\theta\_\{1\}\}\(\\cdot\)and

gθ2​\(⋅\)g^\{\\theta\_\{2\}\}\(\\cdot\)\.

3:\# Training ofsθ1​\(⋅\)s^\{\\theta\_\{1\}\}\(\\cdot\):

4:Use the dataset to optimize

θ1\\theta\_\{1\}via standard denoising score matching\.

5:\# Training ofgθ2​\(⋅\)g^\{\\theta\_\{2\}\}\(\\cdot\):

6:

k←1k\\leftarrow 1\.

7:while

k≤Nepochk\\leq N\_\{\\text\{epoch\}\}do

8:Sample

t∼𝒰​\(0,1\)t\\sim\\mathcal\{U\}\(0,1\)\.

9:Calculate the bias and scaling parameters of the Gaussian and stable RVs used to describe

𝒙t\\bm\{x\}\_\{t\}based on \([27](https://arxiv.org/html/2608.10384#S4.E27)\)

∼\\sim\([29](https://arxiv.org/html/2608.10384#S4.E29)\)\.

10:Calculate the target based on \([2](https://arxiv.org/html/2608.10384#S4.Ex32)\), \([IV\-C](https://arxiv.org/html/2608.10384#S4.Ex38)\), and \([39](https://arxiv.org/html/2608.10384#S4.E39)\)\.

11:Optimize

θ2\\theta\_\{2\}by stochastic gradient descent based on \([2](https://arxiv.org/html/2608.10384#S4.Ex32)\)\.

12:

k←k\+1k\\leftarrow k\+1\.

13:end while

Algorithm 2Sampling algorithm for large jumps1:

𝒙t\\bm\{x\}\_\{t\},

𝒙0∗\\bm\{x\}^\{\*\}\_\{0\}, lookup tables containing

\(r,Fq~,\(r,t\)​\(⋅\)\)\(r,F\_\{\\tilde\{q\},\(r,t\)\}\(\\cdot\)\)and

\(z,Fq~,z\|\(r,t\)​\(⋅\)\)\(z,F\_\{\\tilde\{q\},z\|\(r,t\)\}\(\\cdot\)\)\.

2:Large jump sample

𝝎\\bm\{\\omega\}\.

3:

𝝁t←𝒙t−𝝁​\(t,𝒙0∗\)\\bm\{\\mu\}\_\{t\}\\leftarrow\\bm\{x\}\_\{t\}\-\\bm\{\\mu\}\(t,\\bm\{x\}\_\{0\}^\{\*\}\), where

𝝁​\(t,𝒙0∗\)\\bm\{\\mu\}\(t,\\bm\{x\}\_\{0\}^\{\*\}\)is given in \([27](https://arxiv.org/html/2608.10384#S4.E27)\)\.

4:

𝝁˘t←𝝁t‖𝝁t‖\\breve\{\\bm\{\\mu\}\}\_\{t\}\\leftarrow\\frac\{\\bm\{\\mu\}\_\{t\}\}\{\\\|\\bm\{\\mu\}\_\{t\}\\\|\}\.

5:Sample

u1,u2∼𝒰​\(0,1\)u\_\{1\},u\_\{2\}\\sim\\mathcal\{U\}\(0,1\)and

𝒖3∼𝒩​\(𝟎,𝑰D\)\\bm\{u\}\_\{3\}\\sim\\mathcal\{N\}\(\\bm\{0\},\\bm\{I\}\_\{D\}\)\.

6:

r←Fq~,\(r,t\)−1​\(u1\)r\\leftarrow F^\{\-1\}\_\{\\tilde\{q\},\(r,t\)\}\(u\_\{1\}\), or equivalently, apply rejection sampling based on

q~​\(𝒙t,r\)\\tilde\{q\}\(\\bm\{x\}\_\{t\},r\)with Student’s t\-distribution as the proposal distribution\.

7:

z←Fq~,z\|\(r,t\)−1​\(u2\)z\\leftarrow F^\{\-1\}\_\{\\tilde\{q\},z\|\(r,t\)\}\(u\_\{2\}\)\.

8:

𝒖3←𝒖3−⟨𝝁˘t,𝒖3⟩​𝝁˘t\\bm\{u\}\_\{3\}\\leftarrow\\bm\{u\}\_\{3\}\-\\langle\\breve\{\\bm\{\\mu\}\}\_\{t\},\\bm\{u\}\_\{3\}\\rangle\\breve\{\\bm\{\\mu\}\}\_\{t\}\.

9:

𝒖˘3←𝒖3‖𝒖3‖\\breve\{\\bm\{u\}\}\_\{3\}\\leftarrow\\frac\{\\bm\{u\}\_\{3\}\}\{\\\|\\bm\{u\}\_\{3\}\\\|\}\.

10:

𝝎←r​z​𝝁˘t\+r​1−z2​𝒖˘3\\bm\{\\omega\}\\leftarrow rz\\breve\{\\bm\{\\mu\}\}\_\{t\}\+r\\sqrt\{1\-z^\{2\}\}\\breve\{\\bm\{u\}\}\_\{3\}\.

Algorithm 3Total reverse sampling algorithm1:

𝒙T,j∼pt​p​\(𝒙\),j=1,⋯,N\\bm\{x\}\_\{T,j\}\\sim p\_\{tp\}\(\\bm\{x\}\),j=1,\\cdots,N, well\-trained networks

sθ1​\(⋅\)s^\{\\theta\_\{1\}\}\(\\cdot\)and

gθ2​\(⋅\)g^\{\\theta\_\{2\}\}\(\\cdot\), weighted distribution

gQ​\(𝒙\|𝒙t,t\)g\_\{Q\}\(\\bm\{x\}\|\\bm\{x\}\_\{t\},t\), lookup tables containing

\(r,Fq~,\(r,t\)​\(⋅\)\)\(r,F\_\{\\tilde\{q\},\(r,t\)\}\(\\cdot\)\)and

\(z,Fq~,z\|\(r,t\)​\(⋅\)\)\(z,F\_\{\\tilde\{q\},z\|\(r,t\)\}\(\\cdot\)\), and time step size

Δ​t\\Delta t\.

2:Reversed samples

𝒙0,j,j=1,⋯,N\\bm\{x\}\_\{0,j\},j=1,\\cdots,N\.

3:

j←1j\\leftarrow 1\.

4:while

j≤Nj\\leq Ndo

5:

t←Tt\\leftarrow T\.

6:while

t\>0t\>0do

7:\# Score\-based update:

8:Sample

𝝃1,𝝃2∼𝒩​\(𝟎,𝑰D\)\\bm\{\\xi\}\_\{1\},\\bm\{\\xi\}\_\{2\}\\sim\\mathcal\{N\}\(\\bm\{0\},\\bm\{I\}\_\{D\}\)\.

9:

𝒙~t,j←𝒙t,j\+\(−𝒃​\(𝒙t,j\)\+Σ​\(t\)​sθ1​\(𝒙t,j,t\)\)​Δ​t\+ΦG​\(t\)​Δ​t​𝝃1\\tilde\{\\bm\{x\}\}\_\{t,j\}\\leftarrow\\bm\{x\}\_\{t,j\}\+\(\-\\bm\{b\}\(\\bm\{x\}\_\{t,j\}\)\+\\Sigma\(t\)s^\{\\theta\_\{1\}\}\(\\bm\{x\}\_\{t,j\},t\)\)\\Delta t\+\\Phi\_\{G\}\(t\)\\sqrt\{\\Delta t\}\\bm\{\\xi\}\_\{1\}\.

10:\# Large jump generation:

11:

λ​\(𝒙t,j\)←gθ2​\(𝒙t,j,t\)\\lambda\(\\bm\{x\}\_\{t,j\}\)\\leftarrow g^\{\\theta\_\{2\}\}\(\\bm\{x\}\_\{t,j\},t\)\.

12:Sample

β∼𝒰​\(0,1\)\\beta\\sim\\mathcal\{U\}\(0,1\)\.

13:if

β<1−exp⁡\(−λ​\(𝒙t,j\)​Δ​t\)\\beta<1\-\\exp\(\-\\lambda\(\\bm\{x\}\_\{t,j\}\)\\Delta t\)then

14:Randomly choose a sample

𝒙0∗\\bm\{x\}^\{\*\}\_\{0\}from the dataset according to the truncated weighted distribution

g~Q​\(𝒙\)\\tilde\{g\}\_\{Q\}\(\\bm\{x\}\)\.

15:Obtain a jump sample

𝝎\\bm\{\\omega\}from Algorithm[2](https://arxiv.org/html/2608.10384#alg2)using

𝒙t,j\\bm\{x\}\_\{t,j\},

𝒙0∗\\bm\{x\}^\{\*\}\_\{0\}, and the lookup tables\.

16:else

17:

𝝎←𝟎\\bm\{\\omega\}\\leftarrow\\bm\{0\}\.

18:end if

19:

𝒙~t,j←𝒙~t,j\+𝝎\\tilde\{\\bm\{x\}\}\_\{t,j\}\\leftarrow\\tilde\{\\bm\{x\}\}\_\{t,j\}\+\\bm\{\\omega\}\.

20:\# Small jump generation:

21:Calculate the drift

𝒃t,Δ​t≤ϵ\\bm\{b\}\_\{t,\\Delta t\}^\{\\leq\\epsilon\}in \([III\-B](https://arxiv.org/html/2608.10384#S3.Ex9)\) and the scaling parameter

Σt,Δ​t≤ϵ\\Sigma\_\{t,\\Delta t\}^\{\\leq\\epsilon\}in \([15](https://arxiv.org/html/2608.10384#S3.E15)\)\.

22:

𝒙t−Δ​t,j←𝒙~t,j\+𝒃t,Δ​t≤ϵ​Δ​t\+Δ​t​\(Σt,Δ​t≤ϵ\)12​𝝃2\\bm\{x\}\_\{t\-\\Delta t,j\}\\leftarrow\\tilde\{\\bm\{x\}\}\_\{t,j\}\+\\bm\{b\}\_\{t,\\Delta t\}^\{\\leq\\epsilon\}\\Delta t\+\\sqrt\{\\Delta t\}\(\\Sigma\_\{t,\\Delta t\}^\{\\leq\\epsilon\}\)^\{\\frac\{1\}\{2\}\}\\bm\{\\xi\}\_\{2\}\.

23:

t←t−Δ​tt\\leftarrow t\-\\Delta t\.

24:end while

25:

j←j\+1j\\leftarrow j\+1\.

26:end while

It should be remarked that for the large jump component, we should theoretically sample the number of jumps fromNt\>ϵ∼Poisson⁡\(gθ2​\(𝒙t,t\)​Δ​t\)N\_\{t\}^\{\>\\epsilon\}\\sim\\operatorname\{Poisson\}\\left\(g^\{\\theta\_\{2\}\}\(\\bm\{x\}\_\{t\},t\)\\Delta t\\right\)\. IfNt\>ϵ\>0N\_\{t\}^\{\>\\epsilon\}\>0, independently generate jump amplitudes𝑽t,ℓ\>ϵ,ℓ=1,…,Nt\>ϵ\\bm\{V\}\_\{t,\\ell\}^\{\>\\epsilon\},\\ell=1,\\ldots,N\_\{t\}^\{\>\\epsilon\}\. Then the large jump contribution is𝑿t,Δ​t\>ϵ=∑ℓ=1Nt\>ϵ𝑽t,ℓ\>ϵ\\bm\{X\}\_\{t,\\Delta t\}^\{\>\\epsilon\}=\\sum\_\{\\ell=1\}^\{N\_\{t\}^\{\>\\epsilon\}\}\\bm\{V\}\_\{t,\\ell\}^\{\>\\epsilon\}\. In practice, one may further approximate the compound Poisson update by allowing at most one large jump within each time step\. LetΛ​\(𝒙t,t\)=λ​\(𝒙t,t\)​Δ​t\\Lambda\(\\bm\{x\}\_\{t\},t\)=\\lambda\(\\bm\{x\}\_\{t\},t\)\\Delta t\. The probability of two or more large jumps in the same interval is

Pr⁡\(Nt\>ϵ≥2\)=\\displaystyle\\Pr\\left\(N\_\{t\}^\{\>\\epsilon\}\\geq 2\\right\)=1−e−Λ​\(𝒙t,t\)​\(1\+Λ​\(𝒙t,t\)\)\\displaystyle 1\-e^\{\-\\Lambda\(\\bm\{x\}\_\{t\},t\)\}\\left\(1\+\\Lambda\(\\bm\{x\}\_\{t\},t\)\\right\)≤\\displaystyle\\leqΛ​\(𝒙t,t\)22\\displaystyle\\frac\{\\Lambda\(\\bm\{x\}\_\{t\},t\)^\{2\}\}\{2\}\(45\)
Therefore, the single jump approximation introduces an additional error of order

Emulti=O​\(supt,𝒙tλ​\(𝒙t,t\)2​Δ​t2\),\\displaystyle E\_\{\\mathrm\{multi\}\}=O\\left\(\\sup\_\{t,\\bm\{x\}\_\{t\}\}\\lambda\(\\bm\{x\}\_\{t\},t\)^\{2\}\\Delta t^\{2\}\\right\),\(46\)provided thatλ​\(𝒙t,t\)​Δ​t\\lambda\(\\bm\{x\}\_\{t\},t\)\\Delta tis uniformly small\.

### IV\-EApproximate observation\-guided sampler

Posterior sampling for the diffusion process can be derived straightforwardly based on Bayes’ rule, i\.e\.,∇log⁡p​\(𝒙t\|𝒚\)=∇log⁡p​\(𝒚\|𝒙t\)\+∇log⁡p​\(𝒙t\)\\nabla\\log p\(\\bm\{x\}\_\{t\}\|\\bm\{y\}\)=\\nabla\\log p\(\\bm\{y\}\|\\bm\{x\}\_\{t\}\)\+\\nabla\\log p\(\\bm\{x\}\_\{t\}\)\. The score function is approximated by the output of the neural network, and the likelihood is obtained by approximatingp​\(𝒚\|𝒙t\)p\(\\bm\{y\}\|\\bm\{x\}\_\{t\}\)withp\(𝒚\|E\[𝒙^0\|𝒙t\]\)p\(\\bm\{y\}\|E\[\\hat\{\\bm\{x\}\}\_\{0\}\|\\bm\{x\}\_\{t\}\]\)\[[4](https://arxiv.org/html/2608.10384#bib.bib28)\]\.

However, a similar strategy is not directly applicable to the jump component because the reversed jump kernel involves a nonlocal density ratio\. Exact posterior sampling for Lévy\-driven SDEs is therefore considerably more challenging\. In principle, the large jump rate should be replaced by the posterior large jump rate, which requires likelihood\-reweighted integration over the empirical mixture at each reverse step\. This computation is rather expensive, especially because the integral cannot be simplified in general\.

In this work, we therefore use the prior\-trained rate network as an amortized proposal rate for large jump events in the posterior sampler\. This approximation means that the prior rate network is not interpreted as the exact posterior rate\. Instead, it provides a computationally efficient proposal mechanism for deciding when nonlocal transitions are proposed\. The observation information is then incorporated into the large jump amplitude and direction through likelihood\-reweighted empirical mixture weights\. This design prioritizes the posterior correction of the jump destination, which has a direct impact on the nonlocal transition once a large jump is triggered, while avoiding the high cost of recomputing the exact posterior jump rate at every reverse step\.

In the experiments, we will evaluate its effect by comparing the prior proposal rate with the likelihood\-reweighted posterior rate approximation and by examining the sensitivity of the final estimation performance to the jump\-rate choice\. In the following proposition, we analyze how the observation modifies the empirical mixture weights used for large jump amplitude sampling\.

###### Proposition 4

Let𝐲\\bm\{y\}be the observations used for posterior sampling of the forward SDE \([II](https://arxiv.org/html/2608.10384#S2.Ex1)\)\. Assume that the observation depends on the clean stateX0X\_\{0\}only, i\.e\.,p𝐘\|𝐗0,𝐗t​\(𝐲\|𝐱0,𝐱t\)=p𝐘\|𝐗0​\(𝐲\|𝐱0\)p\_\{\\bm\{Y\}\|\\bm\{X\}\_\{0\},\\bm\{X\}\_\{t\}\}\(\\bm\{y\}\|\\bm\{x\}\_\{0\},\\bm\{x\}\_\{t\}\)=p\_\{\\bm\{Y\}\|\\bm\{X\}\_\{0\}\}\(\\bm\{y\}\|\\bm\{x\}\_\{0\}\)\. Then, the weighted distributiongQ​\(𝐱\|𝐱t,t\)g\_\{Q\}\(\\bm\{x\}\|\\bm\{x\}\_\{t\},t\)used for generating large jump samples should be replaced by

gposterior,Q​\(𝒙\)=∑j=1N1\{𝒙0,j\}​\(𝒙\)​p𝒀\|𝑿0​\(𝒚\|𝒙0,j\)​Qt,ϵ,j∑j=1Np𝒀\|𝑿0​\(𝒚\|𝒙0,j\)​Qt,ϵ,jg\_\{\\text\{posterior\},Q\}\(\\bm\{x\}\)=\\sum\_\{j=1\}^\{N\}\\frac\{1\_\{\\\{\\bm\{x\}\_\{0,j\}\\\}\}\(\\bm\{x\}\)p\_\{\\bm\{Y\}\|\\bm\{X\}\_\{0\}\}\(\\bm\{y\}\|\\bm\{x\}\_\{0,j\}\)Q\_\{t,\\epsilon,j\}\}\{\\sum\_\{j=1\}^\{N\}p\_\{\\bm\{Y\}\|\\bm\{X\}\_\{0\}\}\(\\bm\{y\}\|\\bm\{x\}\_\{0,j\}\)Q\_\{t,\\epsilon,j\}\}\(47\)

###### Proof:

The proof is relegated to Appendix[C](https://arxiv.org/html/2608.10384#A3)\. ∎

Note that the likelihood can usually be efficiently obtained from the transition function between𝒙0\\bm\{x\}\_\{0\}and𝒚\\bm\{y\}\. For instance, this transition function is the channel model for communication tasks\. Proposition[4](https://arxiv.org/html/2608.10384#Thmproposition4)shows that the large jump amplitude sampler can be adapted to the observation\-guided setting by modifying the empirical mixture weights\. Similarly,gposterior,Q​\(𝒙\)g\_\{\\text\{posterior\},Q\}\(\\bm\{x\}\)can also be truncated to reduce the inference latency\.

## VNumerical results

In this section, several numerical experiments are provided to validate the advantages of the proposed inverse sampler\. We consider channel estimation for OFDM\-SISO systems under mixed noise\.

![Refer to caption](https://arxiv.org/html/2608.10384v1/x1.png)\(a\)
![Refer to caption](https://arxiv.org/html/2608.10384v1/x2.png)\(b\)
![Refer to caption](https://arxiv.org/html/2608.10384v1/x3.png)\(c\)

Figure 1:Training loss, test loss, and accuracy validation of jump rate prediction using the proposed lightweight neural network under different values ofα\\alpha\. \(a\) Comparison of training and test losses; \(b\) Comparison of jump rates forα=1\.2\\alpha=1\.2; \(c\) Comparison of jump rates forα=1\.8\\alpha=1\.8\.For the channel model, we consider a frequency\-selective SISO\-OFDM channel modeled by a tapped\-delay\-line \(TDL\) response\. The path gains are modeled as independent circularly symmetric complex Gaussian random variables with prescribed average powers\. Moreover, the delay\-power profiles from 3GPP TR 38\.901 are used\. Due to space limitations, we choose the TDL\-A, TDL\-C, and TDL\-D delay\-power profiles as representative scenarios\. The mixed channel noise consists of WGN and IN described by the Sα\\alphaS distribution\. They are mutually independent and memoryless\. Let their scaling parameters beγg\\gamma\_\{g\}andγs\\gamma\_\{s\}, respectively\. Then, the generalized signal\-to\-noise ratio \(GSNR\) is defined as follows:

GSNR\(dB\)=10​log10⁡Ps2​\(γg2\+γs2\),\\text\{GSNR\(dB\)\}=10\\log\_\{10\}\\frac\{P\_\{s\}\}\{2\(\\gamma\_\{g\}^\{2\}\+\\gamma\_\{s\}^\{2\}\)\},\(48\)wherePsP\_\{s\}is the transmit signal power\.

For the baselines, we consider existing model\-based algorithms and generative methods, including:

- •Linear minimum mean square error \(LMMSE\): LMMSE is an efficient channel estimation method under WGN\. However, its performance and stability may deteriorate considerably in mixed\-noise scenarios due to impulsive outliers\.
- •Clipped LMMSE: Clipped LMMSE first suppresses large\-amplitude pilot observations through clipping and then applies the conventional LMMSE estimator\. This simple preprocessing improves the robustness of LMMSE against IN while retaining low computational complexity\.
- •Clipped orthogonal matching pursuit \(OMP\)\[[25](https://arxiv.org/html/2608.10384#bib.bib30)\]: OMP iteratively selects dominant delay atoms from a common candidate dictionary and reconstructs the frequency\-domain channel\. Clipped OMP is a clipped version of OMP for scenarios with mixed noise\.
- •Clipped sparse Bayesian learning \(SBL\): SBL estimates the delay\-domain sparse channel by learning the hyperparameters that reflect the sparsity\. Clipped SBL is a clipped version of conventional SBL\.
- •Outlier\-aware SBL\[[8](https://arxiv.org/html/2608.10384#bib.bib29),[28](https://arxiv.org/html/2608.10384#bib.bib31)\]: Different from Clipped SBL, outlier\-aware SBL explicitly models the received pilot observations as the sum of a sparse delay\-domain channel component, a sparse outlier component, and Gaussian background noise\. By jointly estimating the channel and the impulsive outliers, this method is designed for mixed Gaussian and IN environments\.
- •Diffusion method\[[29](https://arxiv.org/html/2608.10384#bib.bib18)\]: This method uses a generative model based on the diffusion forward process and posterior inverse inference\. Different from the proposed framework, it does not introduce the jump process and uses only the local score function to recover the target distribution\.
- •Lévy\-DSM\[[27](https://arxiv.org/html/2608.10384#bib.bib22)\]: This method is based on the inverse sampler designed for the Lévy process under the DSM framework\. It shows that the network can be trained to regress noise followingα\\alpha\-stable distributions to achieve inverse sampling\.
- •Benchmark: The benchmark is a genie\-aided LMMSE estimator under WGN with known true delay support\. It uses ideal prior delay information that is unavailable to practical estimators\. Therefore, it serves as an NMSE reference for the achievable estimation performance\.

Finally, the other parameter configurations are as follows: the number of subcarriers isNc=64N\_\{c\}=64; the number of OFDM symbols per frame isNt=16N\_\{t\}=16; the pilot pattern is “comb”; the pilot spacing is 4 or 8; the modulation scheme is QPSK; and the number of paths, path delays, and path powers follow TDL\-A/C/D\. The pilot spacing is set to either 4 or 8 to validate the performance under different pilot densities\.

### V\-AStructure and training of networkgθ2g^\{\\theta\_\{2\}\}

Since the large jump rate is a scalar statistical quantity, a lightweight neural network is sufficient for its amortized estimation\. Specifically, we parameterizegθ2​\(𝒙t,t\)g\_\{\\theta\_\{2\}\}\(\\bm\{x\}\_\{t\},t\)by a convolutional neural network \(CNN\) and a multilayer perceptron \(MLP\)\. Its input is the normalized state𝒙t\\bm\{x\}\_\{t\}concatenated with a time embedding oftt\. To ensure the nonnegativity of the estimated jump rate, the final output is passed through a softplus activation\. Compared with ReLU, softplus provides a smoother parameterization and avoids permanently inactive zero\-rate outputs\. In our implementation, the network is trained by minimizing the mean square error \(MSE\) between its output andλ​\(𝒙t\|𝒙0\)\\lambda\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\. Note that the original sample𝒙0\\bm\{x\}\_\{0\}is only used to construct the conditional training target\. It is not fed into the network as input because an accurate𝒙0\\bm\{x\}\_\{0\}is not available at the inference stage\.

For the OFDM channel estimation task, the complex\-valued state is first represented by its real and imaginary parts\. For example, when one frame containsNtN\_\{t\}OFDM symbols andNcN\_\{c\}subcarriers, we reshape the state into a two\-channel tensor𝐗t∈ℝ2×Nt×Nc\\mathbf\{X\}\_\{t\}\\in\\mathbb\{R\}^\{2\\times N\_\{t\}\\times N\_\{c\}\}\. The rate network adopts a small convolutional architecture to exploit the local time\-frequency structure while keeping the number of trainable parameters low\. Specifically, it consists of four3×33\\times 3convolutional layers with channel widths\{16,16,32,32\}\\\{16,16,32,32\\\}, respectively, followed by a global average pooling layer\. The pooled feature is concatenated with the time embedding and then passed to a two\-layer MLP head to produce a scalar output\. The estimated large jump rate is given by

gθ2​\(𝐱t,t\)=\\displaystyle g^\{\\theta\_\{2\}\}\(\\mathbf\{x\}\_\{t\},t\)=softplus\(MLPθ2\(\[GAP\(CNNθ2\(𝐗t\)\),Emb\(t\)\)\)\\displaystyle\\mathrm\{softplus\}\(\\mathrm\{MLP\}^\{\\theta\_\{2\}\}\(\[\\mathrm\{GAP\}\(\\mathrm\{CNN\}^\{\\theta\_\{2\}\}\(\\mathbf\{X\}\_\{t\}\)\),\\mathrm\{Emb\}\(t\)\)\)\+λϵ,\\displaystyle\+\\lambda\_\{\\epsilon\},\(49\)whereGAP​\(⋅\)\\mathrm\{GAP\}\(\\cdot\)denotes global average pooling,Emb​\(t\)\\mathrm\{Emb\}\(t\)is the time embedding, andλϵ\>0\\lambda\_\{\\epsilon\}\>0is a small constant for numerical stability\. ForNt=16N\_\{t\}=16andNc=64N\_\{c\}=64, this small CNN contains approximately2×1042\\times 10^\{4\}trainable parameters\. This is much smaller than directly parameterizing the whole reverse transition with a neural network\.

The MSE loss between the accurate rate and the network output is shown in Fig\.[1](https://arxiv.org/html/2608.10384#S5.F1)\. Fig\.[1b](https://arxiv.org/html/2608.10384#S5.F1.sf2)and Fig\.[1c](https://arxiv.org/html/2608.10384#S5.F1.sf3)present the results of 1000 experiments, where the time instant and data sample𝒙0\\bm\{x\}\_\{0\}are chosen randomly\. These 1000 tests are then sorted according to their accurate jump rates\. For validation, the reference large jump rate is computed using the empirical marginal expression in \([IV\-A](https://arxiv.org/html/2608.10384#S4.Ex29)\), rather than a single conditional targetλ​\(𝒙t\|𝒙0\)\\lambda\(\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0\}\)\. From Fig\.[1a](https://arxiv.org/html/2608.10384#S5.F1.sf1), it can be observed that the neural network converges very fast and only requires about 60 epochs with different values ofα\\alpha\. Moreover, the network can accurately estimate the jump rate under different impulsiveness levels\.

![Refer to caption](https://arxiv.org/html/2608.10384v1/x4.png)\(a\)NMSE, TDL\-A,α=1\.2\\alpha=1\.2
![Refer to caption](https://arxiv.org/html/2608.10384v1/x5.png)\(b\)BER, TDL\-A,α=1\.2\\alpha=1\.2
![Refer to caption](https://arxiv.org/html/2608.10384v1/x6.png)\(c\)NMSE, TDL\-C,α=1\.2\\alpha=1\.2
![Refer to caption](https://arxiv.org/html/2608.10384v1/x7.png)\(d\)BER, TDL\-C,α=1\.2\\alpha=1\.2
![Refer to caption](https://arxiv.org/html/2608.10384v1/x8.png)\(e\)NMSE, TDL\-D,α=1\.2\\alpha=1\.2
![Refer to caption](https://arxiv.org/html/2608.10384v1/x9.png)\(f\)BER, TDL\-D,α=1\.2\\alpha=1\.2

Figure 2:NMSE and BER comparisons of channel estimation under different scenarios \(TDL\-A/C/D\) withα=1\.2\\alpha=1\.2\.
### V\-BVarious channel configurations

Here, we investigate the channel estimation performance under variousα\\alpha\. First, the channel estimation accuracy and the corresponding detection performance are compared\. We use normalized mean square error \(NMSE\) and bit error rate \(BER\) as the criteria\. The experimental results are presented in Fig\.[2](https://arxiv.org/html/2608.10384#S5.F2)and Fig\.[3](https://arxiv.org/html/2608.10384#S5.F3)\.

In general, the proposed approach achieves the best performance among the baselines in terms of both NMSE and BER\. For “Clipped LMMSE”, “Clipped SBL”, and “Clipped OMP”, the clipping operation suppresses most IN samples and thus improves robustness\. However, clipping also distorts the transmitted signals, especially in the high\-GSNR regime\. Therefore, these three methods show a lower NMSE decay rate as the GSNR continues to increase\. “Outlier\-aware SBL” has better robustness under mixed channel noise because it explicitly considers the influence of IN samples\. However, it can only achieve an accurate estimation at the pilot locations\. For learning\-based methods, both “Score\-based” and “Proposed” learn the statistical structure of the channel model and thus achieve better estimation accuracy in various scenarios\. Compared with the “Score\-based” algorithm, the proposed method is more likely to reach a better performance because it allows large jumps to enable nonlocal transitions between distant high\-probability regions\. In contrast, “Score\-based” only uses gradient information and small step sizes for optimization\. As a result, it may get trapped in local optima, especially when the channel model is more complicated\. Finally, the proposed method is also more robust and stable compared to the “Lévy\-DSM”, which is consistent with the above analyses\.

Fig\.[3](https://arxiv.org/html/2608.10384#S5.F3)presents the estimation results under different impulsiveness levels\. Whenα=1\.2\\alpha=1\.2, there are multiple outliers with large amplitudes, and the channel becomes much more complex\. In this case, the proposed approach is not only robust but also achieves a larger performance gain compared with the cases with largerα\\alpha\. However, whenα\\alphais close to 2 and the impulsiveness is very weak, the existing baselines also show comparable performance\. This is consistent with the above analysis\.

![Refer to caption](https://arxiv.org/html/2608.10384v1/x10.png)\(a\)NMSE, TDL\-C,α=1\.8\\alpha=1\.8
![Refer to caption](https://arxiv.org/html/2608.10384v1/x11.png)\(b\)BER, TDL\-C,α=1\.8\\alpha=1\.8

Figure 3:NMSE and BER comparisons of channel estimation with TDL\-C as the delay\-power profile andα=1\.8\\alpha=1\.8\.
### V\-CAblations

Ablation studies are conducted to explore the influence from the approximation of the forward noise distribution, the top\-KKtruncation and the approximated posterior jump rate\.

During the design of the inverse sampler, the distribution of the forward mixed noise at arbitrary time instants is required\. However, there is no exact closed\-form PDF expression for the impulsive component or the mixed noise\. To avoid numerical integration, we use the approximated PDF of the mixed noise given in \([30](https://arxiv.org/html/2608.10384#S4.E30)\) and \([IV\-C](https://arxiv.org/html/2608.10384#S4.Ex33)\)\. The following experiments demonstrate the effectiveness of this approximation\. We still use channel estimation as the application\.

In Fig\.[4](https://arxiv.org/html/2608.10384#S5.F4), Fig\.[5](https://arxiv.org/html/2608.10384#S5.F5)and Fig\.[6](https://arxiv.org/html/2608.10384#S5.F6), we separately focus on the effects of noise distribution approximation, top\-KKtruncation and large jump rate approximation\. The time consumption is shown in Table[I](https://arxiv.org/html/2608.10384#S5.T1), using a CPU configuration of ‘12th Gen Intel\(R\) Core\(TM\) i9\-12900H’\. It can be observed that using the accurate PDF does not bring a considerable performance gain, while the complexity of “Accurate PDF” is unacceptable\. This is because the distribution needs to be calculated by numerical integration at every step\. More importantly, even though the high\-dimensional PDF can be transformed into polar coordinates, the integral is still challenging to further simplify\. This remains time\-consuming when a high resolution is required to ensure accuracy\.

Then, the approximation error introduced by the top\-KKtruncation is examined\. For the benchmark, we setKKequal to the total number of samples used for training\. Here, for the truncation scheme, we setKKas10%10\\%of the total number of training samples\. Similar to the noise model approximation, the truncation error is also not significant\. This is because, for a given time instant and sample𝒙t\\bm\{x\}\_\{t\}, the probability that most original data samples generate𝒙t\\bm\{x\}\_\{t\}is quite small\. Therefore, the influence of these samples on jump generation is negligible\. As a result,KKcan be much smaller than the total number of samples, which reduces redundancy and improves efficiency\.

Finally, Fig\.[6](https://arxiv.org/html/2608.10384#S5.F6)implies that the posterior inverse sampling based on the exact large jump rate achieves a better NMSE\. However, its computational cost significantly increases since the jump rate calculation needs to rely on the numerical integral\. Moreover, in terms of detection error probability, the performance gain is less considerable than that with respect to NMSE\. This result supports that using prior large jump rate is an efficient scheme\.

In summary, the approximation operations in the inverse sampling algorithm achieve a tradeoff between performance and complexity\.

![Refer to caption](https://arxiv.org/html/2608.10384v1/x12.png)\(a\)NMSE
![Refer to caption](https://arxiv.org/html/2608.10384v1/x13.png)\(b\)BER

Figure 4:Ablation results compared with the accurate mixed noise PDF in terms of NMSE and BER\.![Refer to caption](https://arxiv.org/html/2608.10384v1/x14.png)\(a\)NMSE
![Refer to caption](https://arxiv.org/html/2608.10384v1/x15.png)\(b\)BER

Figure 5:Ablation results compared with the no\-truncation case in terms of NMSE and BER\.![Refer to caption](https://arxiv.org/html/2608.10384v1/x16.png)\(a\)NMSE
![Refer to caption](https://arxiv.org/html/2608.10384v1/x17.png)\(b\)BER

Figure 6:Ablation results compared with the accurate posterior large jump rate in terms of NMSE and BER\.TABLE I:Time consumption comparison of simulations in Fig\.[4](https://arxiv.org/html/2608.10384#S5.F4)∼\\sim[6](https://arxiv.org/html/2608.10384#S5.F6)under different values ofα\\alpha\.Algorithms/Scenariosα=1\.2\\alpha=1\.2α=1\.8\\alpha=1\.8Proposed5\.6s5\.5sAccurate PDF1717\.6s1698\.4sNo truncation52\.7s54\.1sAccurate rate6125\.3s6379\.6s

## VIConclusions

In this paper, we proposed a generator\-guided inverse sampling algorithm for a class of isotropic linear Lévy\-driven SDEs with symmetricα\\alpha\-stable jump components\. We first derived the expression of the backward generator and used it to decompose the full reverse sampling process into the diffusion part and the jump components\. Then, by theoretically analyzing the inverse sampling process from a statistical perspective, we determined the objective that can be learned by a lightweight network\. We also designed an efficient reverse sampling method by introducing several efficient approximation operations\. As an application, the proposed inverse sampling algorithm outperformed the considered baselines in channel estimation under mixed channel noise across different scenarios\. Moreover, numerical results showed that the proposed approach achieves a good balance between performance and complexity\.

## Appendix AProof of Proposition[1](https://arxiv.org/html/2608.10384#Thmproposition1)

For a Lévy process with infinite activity, the generator associated with the test functionffis given by

\(ℒJ,F​f\)​\(𝒙t\)=𝒃J,Fϵ​\(𝒙t,t\)⋅∇f​\(𝒙t\)\\displaystyle\(\\mathcal\{L\}\_\{J,F\}f\)\(\\bm\{x\}\_\{t\}\)=\\bm\{b\}\_\{J,F\}^\{\\epsilon\}\(\\bm\{x\}\_\{t\},t\)\\cdot\\nabla f\(\\bm\{x\}\_\{t\}\)\+∫RD\(f​\(𝒙t\+𝒗\)−f​\(𝒙t\)−1\{𝒗:‖𝒗‖≤ϵ\}​𝒗⋅∇f​\(𝒙t\)\)​νt​\(d​𝒗\)\\displaystyle\+\\int\_\{R^\{D\}\}\{\\left\(\{f\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\-f\(\\bm\{x\}\_\{t\}\)\-1\_\{\\\{\\bm\{v\}:\\\|\\bm\{v\}\\\|\\leq\\epsilon\\\}\}\\bm\{v\}\\cdot\\nabla f\(\\bm\{x\}\_\{t\}\)\}\\right\)\\nu\_\{t\}\(d\\bm\{v\}\)\}\(A\.1\)where𝒃J,Fϵ​\(𝒙t,t\)\\bm\{b\}\_\{J,F\}^\{\\epsilon\}\(\\bm\{x\}\_\{t\},t\)is the canonical drift induced byd​𝑳td\\bm\{L\}\_\{t\}\. For the isotropic case, we have𝒃J,Fϵ​\(𝒙t,t\)=0\\bm\{b\}\_\{J,F\}^\{\\epsilon\}\(\\bm\{x\}\_\{t\},t\)=0\. By substituting \([A](https://arxiv.org/html/2608.10384#A1.Ex45)\) into \([II](https://arxiv.org/html/2608.10384#S2.Ex2)\), the complete forward generator is obtained as

\(ℒF​f\)​\(𝒙t\)=𝒃Fϵ​\(𝒙t,t\)⋅∇f​\(𝒙t\)\+12​Tr​\(Σ​\(t\)​∇2f​\(𝒙t\)\)\\displaystyle\\left\(\{\\mathcal\{L\}\_\{F\}f\}\\right\)\(\\bm\{x\}\_\{t\}\)=\\bm\{b\}\_\{F\}^\{\\epsilon\}\(\\bm\{x\}\_\{t\},t\)\\cdot\\nabla f\(\\bm\{x\}\_\{t\}\)\+\\frac\{1\}\{2\}\{\\rm\{Tr\}\}\\left\(\{\\Sigma\(t\)\\nabla^\{2\}f\(\\bm\{x\}\_\{t\}\)\}\\right\)\+∫RD\(f​\(𝒙t\+𝒗\)−f​\(𝒙t\)−1\{𝒗:‖𝒗‖≤ϵ\}​𝒗⋅∇f​\(𝒙t\)\)​νt​\(d​𝒗\)\\displaystyle\+\\int\_\{R^\{D\}\}\{\\left\(\{f\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\-f\(\\bm\{x\}\_\{t\}\)\-1\_\{\\\{\\bm\{v\}:\\\|\\bm\{v\}\\\|\\leq\\epsilon\\\}\}\\bm\{v\}\\cdot\\nabla f\(\\bm\{x\}\_\{t\}\)\}\\right\)\\nu\_\{t\}\(d\\bm\{v\}\)\}\(A\.2\)with𝒃Fϵ​\(𝒙t,t\)=𝒃​\(𝒙t,t\)\+𝒃J,Fϵ​\(𝒙t,t\)\\bm\{b\}\_\{F\}^\{\\epsilon\}\(\\bm\{x\}\_\{t\},t\)=\\bm\{b\}\(\\bm\{x\}\_\{t\},t\)\+\\bm\{b\}\_\{J,F\}^\{\\epsilon\}\(\\bm\{x\}\_\{t\},t\)\. Next, time reversal is applied to derive the backward generator\. Based on\[[6](https://arxiv.org/html/2608.10384#bib.bib24)\], the forward and backward jump kernels satisfy the flux equation

p​\(𝒙t\)​𝐾←t,𝒙t​\(d​𝒗\)=p​\(𝒙t\+𝒗\)​𝐾→t,𝒙t\+𝒗​\(d​\(−𝒗\)\)\\displaystyle p\(\\bm\{x\}\_\{t\}\)\{\{\\mathord\{\\mathrel\{\\mathop\{\\kern 0\.0ptK\}\\limits^\{\{\\lower 3\.0pt\\hbox\{$\\scriptscriptstyle\\leftarrow$\}\}\}\}\}\}\_\{t,\\bm\{x\}\_\{t\}\}\}\(d\\bm\{v\}\)=p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\{\{\\mathord\{\\mathrel\{\\mathop\{\\kern 0\.0ptK\}\\limits^\{\{\\lower 3\.0pt\\hbox\{$\\scriptscriptstyle\\rightarrow$\}\}\}\}\}\}\_\{t,\\bm\{x\}\_\{t\}\+\\bm\{v\}\}\}\(d\(\-\\bm\{v\}\)\)\(A\.3\)which means that the jump density from𝒙t\\bm\{x\}\_\{t\}to𝒙t\+𝒗\\bm\{x\}\_\{t\}\+\\bm\{v\}should be equal to that from𝒙t\+𝒗\\bm\{x\}\_\{t\}\+\\bm\{v\}back to𝒙t\\bm\{x\}\_\{t\}\. Moreover, the backward canonical drift satisfies\[[6](https://arxiv.org/html/2608.10384#bib.bib24)\]

𝒃J,Fϵ​\(𝒙t,t\)\\displaystyle\\bm\{b\}\_\{J,F\}^\{\\epsilon\}\(\\bm\{x\}\_\{t\},t\)\+𝒃J,Bϵ​\(𝒙t,t\)\\displaystyle\+\\bm\{b\}\_\{J,B\}^\{\\epsilon\}\(\\bm\{x\}\_\{t\},t\)=\\displaystyle=∫RD1\{𝒗:‖𝒗‖≤ϵ\}​𝒗​\(𝐾→t,𝒙t​\(d​𝒗\)\+𝐾←t,𝒙t​\(d​𝒗\)\)\\displaystyle\\int\_\{R^\{D\}\}\{1\_\{\\\{\\bm\{v\}:\\\|\\bm\{v\}\\\|\\leq\\epsilon\\\}\}\\bm\{v\}\\left\(\{\{\\mathord\{\\mathrel\{\\mathop\{\\kern 0\.0ptK\}\\limits^\{\{\\lower 3\.0pt\\hbox\{$\\scriptscriptstyle\\rightarrow$\}\}\}\}\}\}\_\{t,\\bm\{x\}\_\{t\}\}\(d\\bm\{v\}\)\+\{\\mathord\{\\mathrel\{\\mathop\{\\kern 0\.0ptK\}\\limits^\{\{\\lower 3\.0pt\\hbox\{$\\scriptscriptstyle\\leftarrow$\}\}\}\}\}\}\_\{t,\\bm\{x\}\_\{t\}\}\(d\\bm\{v\}\)\}\\right\)\}\(A\.4\)
Therefore, we have

𝐾←t,𝒙t​\(d​𝒗\)=\\displaystyle\{\{\\mathord\{\\mathrel\{\\mathop\{\\kern 0\.0ptK\}\\limits^\{\{\\lower 3\.0pt\\hbox\{$\\scriptscriptstyle\\leftarrow$\}\}\}\}\}\}\_\{t,\\bm\{x\}\_\{t\}\}\}\(d\\bm\{v\}\)=p​\(𝒙t\+𝒗\)p​\(𝒙t\)​K→t,𝒙t\+𝒗​\(d​\(−𝒗\)\)\\displaystyle\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}\{\{\\vec\{K\}\}\_\{t,\\bm\{x\}\_\{t\}\+\\bm\{v\}\}\}\\left\(\{d\\left\(\{\-\\bm\{v\}\}\\right\)\}\\right\)=\\displaystyle=p​\(𝒙t\+𝒗\)p​\(𝒙t\)​νt​\(d​\(−𝒗\)\)\\displaystyle\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\{p\(\\bm\{x\}\_\{t\}\)\}\\nu\_\{t\}\(d\(\-\\bm\{v\}\)\)\(A\.5\)
In this case, the backward generator of the Lévy process is given by

\(ℒJ,B​f\)​\(𝒙t\)\\displaystyle\(\\mathcal\{L\}\_\{J,B\}f\)\(\\bm\{x\}\_\{t\}\)=\\displaystyle=𝒃J,Bϵ\(𝒙t,t\)⋅∇f\(𝒙t\)\+∫RD\(f\(𝒙t\+𝒗\)−f\(𝒙t\)\\displaystyle\\bm\{b\}\_\{J,B\}^\{\\epsilon\}\(\\bm\{x\}\_\{t\},t\)\\cdot\\nabla f\(\\bm\{x\}\_\{t\}\)\+\\int\_\{R^\{D\}\}\\big\(f\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\-f\(\\bm\{x\}\_\{t\}\)−1\{𝒗:‖𝒗‖≤ϵ\}𝒗⋅∇f\(𝒙t\)\)𝐾←t,𝒙t\(d𝒗\)\\displaystyle\-1\_\{\\\{\\bm\{v\}:\\\|\\bm\{v\}\\\|\\leq\\epsilon\\\}\}\\bm\{v\}\\cdot\\nabla f\(\\bm\{x\}\_\{t\}\)\\big\)\{\{\\mathord\{\\mathrel\{\\mathop\{\\kern 0\.0ptK\}\\limits^\{\{\\lower 3\.0pt\\hbox\{$\\scriptscriptstyle\\leftarrow$\}\}\}\}\}\}\_\{t,\\bm\{x\}\_\{t\}\}\}\(d\\bm\{v\}\)=\\displaystyle=\(−𝒃J,Fϵ\(𝒙t,t\)\+∫RD1\{𝒗:‖𝒗‖≤ϵ\}𝒗⋅\(νt\(d𝒗\)\\displaystyle\\bigg\(\-\\bm\{b\}\_\{J,F\}^\{\\epsilon\}\(\\bm\{x\}\_\{t\},t\)\+\\int\_\{R^\{D\}\}1\_\{\\\{\\bm\{v\}:\\\|\\bm\{v\}\\\|\\leq\\epsilon\\\}\}\\bm\{v\}\\cdot\\Big\(\\nu\_\{t\}\(d\\bm\{v\}\)\+𝐾←t,𝒙t\(d𝒗\)\)\)⋅∇f\(𝒙t\)\+∫RD\(f\(𝒙t\+𝒗\)−f\(𝒙t\)\\displaystyle\+\{\{\\mathord\{\\mathrel\{\\mathop\{\\kern 0\.0ptK\}\\limits^\{\{\\lower 3\.0pt\\hbox\{$\\scriptscriptstyle\\leftarrow$\}\}\}\}\}\}\_\{t,\\bm\{x\}\_\{t\}\}\}\(d\\bm\{v\}\)\\Big\)\\bigg\)\\cdot\\nabla f\(\\bm\{x\}\_\{t\}\)\+\\int\_\{R^\{D\}\}\\big\(f\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\-f\(\\bm\{x\}\_\{t\}\)−1\{𝒗:‖𝒗‖≤ϵ\}𝒗⋅∇f\(𝒙t\)\)𝐾←t,𝒙t\(d𝒗\)\\displaystyle\-1\_\{\\\{\\bm\{v\}:\\\|\\bm\{v\}\\\|\\leq\\epsilon\\\}\}\\bm\{v\}\\cdot\\nabla f\(\\bm\{x\}\_\{t\}\)\\big\)\\mathord\{\\mathrel\{\\mathop\{\\kern 0\.0ptK\}\\limits^\{\{\\lower 3\.0pt\\hbox\{$\\scriptscriptstyle\\leftarrow$\}\}\}\}\}\_\{t,\\bm\{x\}\_\{t\}\}\(d\\bm\{v\}\)=\\displaystyle=\(−𝒃J,Fϵ​\(𝒙t,t\)\+p\.v\.​∫RD1\{𝒗:‖𝒗‖≤ϵ\}​𝒗​νt​\(d​𝒗\)\)⋅∇f​\(𝒙t\)\\displaystyle\\bigg\(\-\\bm\{b\}\_\{J,F\}^\{\\epsilon\}\(\\bm\{x\}\_\{t\},t\)\+\\operatorname\{p\.v\.\}\\int\_\{R^\{D\}\}1\_\{\\\{\\bm\{v\}:\\\|\\bm\{v\}\\\|\\leq\\epsilon\\\}\}\\bm\{v\}\\nu\_\{t\}\(d\\bm\{v\}\)\\bigg\)\\cdot\\nabla f\(\\bm\{x\}\_\{t\}\)\+p\.v\.​∫RD\(f​\(𝒙t\+𝒗\)−f​\(𝒙t\)\)​p​\(𝒙t\+𝒗\)p​\(𝒙t\)​νt​\(d​𝒗\)\\displaystyle\+\\operatorname\{p\.v\.\}\\int\_\{R^\{D\}\}\\left\(f\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\-f\(\\bm\{x\}\_\{t\}\)\\right\)\\frac\{\{\{p\}\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\}\}\{\{\{p\}\(\\bm\{x\}\_\{t\}\)\}\}\\nu\_\{t\}\(d\\bm\{v\}\)\(A\.6\)
Combining \([A](https://arxiv.org/html/2608.10384#A1.Ex46)\) with the time reversal of the diffusion part, the complete backward generator is given by

\(ℒB​f\)​\(𝒙t\)=\\displaystyle\\left\(\{\\mathcal\{L\}\_\{B\}f\}\\right\)\(\\bm\{x\}\_\{t\}\)=\(−𝒃​\(𝒙t,t\)\+Σ​\(t\)​∇log⁡p​\(𝒙t\)\)⋅∇f​\(𝒙t\)\\displaystyle\\big\(\-\\bm\{b\}\(\\bm\{x\}\_\{t\},t\)\+\\Sigma\(t\)\\nabla\\log p\(\\bm\{x\}\_\{t\}\)\\big\)\\cdot\\nabla f\(\\bm\{x\}\_\{t\}\)\+12​Tr​\(Σ​\(t\)​∇2f​\(𝒙t\)\)\+\(ℒJ,B​f\)​\(𝒙t\)\\displaystyle\+\\frac\{1\}\{2\}\{\\rm\{Tr\}\}\\left\(\{\\Sigma\(t\)\\nabla^\{2\}f\(\\bm\{x\}\_\{t\}\)\}\\right\)\+\(\\mathcal\{L\}\_\{J,B\}f\)\(\\bm\{x\}\_\{t\}\)\(A\.7\)
Finally, if the scaled Lévy measureνt\\nu\_\{t\}is symmetric, we haveνt​\(d​\(−𝒗\)\)=νt​\(d​𝒗\)\\nu\_\{t\}\(d\(\-\\bm\{v\}\)\)=\\nu\_\{t\}\(d\\bm\{v\}\)\. For the symmetricα\\alpha\-stable case, the canonical jump drift vanishes\. This further simplifies the backward generator and completes the proof\.

## Appendix BProof of Proposition[3](https://arxiv.org/html/2608.10384#Thmproposition3)

We first use a variable transformation to absorb the drift term\. Let𝒀t=exp⁡\(−R​t\)​𝑿t\\bm\{Y\}\_\{t\}=\\exp\(\-Rt\)\\bm\{X\}\_\{t\}\. Based on Ito’s formula,

d​𝒀t=\\displaystyle d\\bm\{Y\}\_\{t\}=exp⁡\(−R​t\)​d​𝑿t\+d​\(exp⁡\(−R​t\)\)​𝑿t\\displaystyle\\exp\(\-Rt\)d\\bm\{X\}\_\{t\}\+d\(\\exp\(\-Rt\)\)\\bm\{X\}\_\{t\}=\\displaystyle=exp⁡\(−R​t\)​d​𝑿t−R​exp⁡\(−R​t\)​𝑿t\\displaystyle\\exp\(\-Rt\)d\\bm\{X\}\_\{t\}\-R\\exp\(\-Rt\)\\bm\{X\}\_\{t\}=\\displaystyle=exp⁡\(−R​t\)​\(𝒔​d​t\+ΦG​\(t\)​d​𝑾t\+ΦS​\(t\)​d​𝑳t\)\\displaystyle\\exp\(\-Rt\)\(\\bm\{s\}dt\+\\Phi\_\{G\}\(t\)d\\bm\{W\}\_\{t\}\+\\Phi\_\{S\}\(t\)d\\bm\{L\}\_\{t\}\)\(B\.1\)
According to𝒀t=exp⁡\(−R​t\)​𝑿t\\bm\{Y\}\_\{t\}=\\exp\(\-Rt\)\\bm\{X\}\_\{t\},

𝑿t=\\displaystyle\\bm\{X\}\_\{t\}=exp⁡\(R​t\)​𝑿0\+exp⁡\(R​t\)​∫0texp⁡\(−R​l\)​𝒔​𝑑l⏟𝝁​\(t,𝒙0\)\\displaystyle\\underbrace\{\\exp\(Rt\)\\bm\{X\}\_\{0\}\+\\exp\(Rt\)\\int\_\{0\}^\{t\}\\exp\(\-Rl\)\\bm\{s\}dl\}\_\{\\bm\{\\mu\}\(t,\\bm\{x\}\_\{0\}\)\}\+∫0texp⁡\(\(t−l\)​R\)​ΦG​\(l\)​𝑑𝑾l⏟𝑮​\(t\)\\displaystyle\+\\underbrace\{\\int\_\{0\}^\{t\}\\exp\(\(t\-l\)R\)\\Phi\_\{G\}\(l\)d\\bm\{W\}\_\{l\}\}\_\{\\bm\{G\}\(t\)\}\+∫0texp⁡\(\(t−l\)​R\)​ΦS​\(l\)​𝑑𝑳l⏟𝑺​\(t\)\\displaystyle\+\\underbrace\{\\int\_\{0\}^\{t\}\\exp\(\(t\-l\)R\)\\Phi\_\{S\}\(l\)d\\bm\{L\}\_\{l\}\}\_\{\\bm\{S\}\(t\)\}\(B\.2\)
Then, based on the linearity of Gaussian RV andα\\alpha\-stable RV, the proof is completed\.

## Appendix CProof of Proposition[4](https://arxiv.org/html/2608.10384#Thmproposition4)

First, we havep​\(𝒙t\+𝒗\|𝒚\)p​\(𝒙t\|𝒚\)​νt​\(d​𝒗\)=p​\(𝒙t\+𝒗,𝒚\)p​\(𝒙t,𝒚\)​νt​\(d​𝒗\)\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\|\\bm\{y\}\)\}\{p\(\\bm\{x\}\_\{t\}\|\\bm\{y\}\)\}\\nu\_\{t\}\(d\\bm\{v\}\)=\\frac\{p\(\\bm\{x\}\_\{t\}\+\\bm\{v\},\\bm\{y\}\)\}\{p\(\\bm\{x\}\_\{t\},\\bm\{y\}\)\}\\nu\_\{t\}\(d\\bm\{v\}\)\. Since the denominator is independent of𝒗\\bm\{v\}, it can be omitted during sampling\. For clarity, we add subscripts to indicate the arguments of the distributions\. The target sampling distribution is given by

∫‖𝒗‖\>ϵp𝑿t,𝒀​\(𝒙t\+𝒗,𝒚\)​νt​\(d​𝒗\)\\displaystyle\\int\_\{\\\|\\bm\{v\}\\\|\>\\epsilon\}p\_\{\\bm\{X\}\_\{t\},\\bm\{Y\}\}\(\\bm\{x\}\_\{t\}\+\\bm\{v\},\\bm\{y\}\)\\nu\_\{t\}\(d\\bm\{v\}\)∝\\displaystyle\\propto∫‖𝒗‖\>ϵp𝒀\|𝑿t​\(𝒚\|𝒙t\+𝒗\)​p𝑿t​\(𝒙t\+𝒗\)​‖𝒗‖−α−D​𝑑𝒗\\displaystyle\\int\_\{\\\|\\bm\{v\}\\\|\>\\epsilon\}p\_\{\\bm\{Y\}\|\\bm\{X\}\_\{t\}\}\(\\bm\{y\}\|\\bm\{x\}\_\{t\}\+\\bm\{v\}\)p\_\{\\bm\{X\}\_\{t\}\}\(\\bm\{x\}\_\{t\}\+\\bm\{v\}\)\\\|\\bm\{v\}\\\|^\{\-\\alpha\-D\}d\\bm\{v\}=\\displaystyle=∫‖𝒖−𝒙t‖\>ϵp𝒀\|𝑿t​\(𝒚\|𝒖\)​p𝑿t​\(𝒖\)​‖𝒖−𝒙t‖−α−D​𝑑𝒖\\displaystyle\\int\_\{\\\|\\bm\{u\}\-\\bm\{x\}\_\{t\}\\\|\>\\epsilon\}p\_\{\\bm\{Y\}\|\\bm\{X\}\_\{t\}\}\(\\bm\{y\}\|\\bm\{u\}\)p\_\{\\bm\{X\}\_\{t\}\}\(\\bm\{u\}\)\\\|\\bm\{u\}\-\\bm\{x\}\_\{t\}\\\|^\{\-\\alpha\-D\}d\\bm\{u\}=\\displaystyle=∫‖𝒖−𝒙t‖\>ϵ∫ℛDp𝒀\|𝑿0​\(𝒚\|𝒙0\)​p𝑿0\|𝑿t​\(𝒙0\|𝒖\)​𝑑𝒙0​p𝑿t​\(𝒖\)\\displaystyle\\int\_\{\\\|\\bm\{u\}\-\\bm\{x\}\_\{t\}\\\|\>\\epsilon\}\\int\_\{\\mathcal\{R\}^\{D\}\}p\_\{\\bm\{Y\}\|\\bm\{X\}\_\{0\}\}\(\\bm\{y\}\|\\bm\{x\}\_\{0\}\)p\_\{\\bm\{X\}\_\{0\}\|\\bm\{X\}\_\{t\}\}\(\\bm\{x\}\_\{0\}\|\\bm\{u\}\)d\\bm\{x\}\_\{0\}p\_\{\\bm\{X\}\_\{t\}\}\(\\bm\{u\}\)×‖𝒖−𝒙t‖−α−D​d​𝒖\\displaystyle\\times\\\|\\bm\{u\}\-\\bm\{x\}\_\{t\}\\\|^\{\-\\alpha\-D\}d\\bm\{u\}=\\displaystyle=∫‖𝒖−𝒙t‖\>ϵ∫ℛDp𝒀\|𝑿0​\(𝒚\|𝒙0\)​p𝑿t\|𝑿0​\(𝒖\|𝒙0\)​p𝑿0​\(𝒙0\)​𝑑𝒙0\\displaystyle\\int\_\{\\\|\\bm\{u\}\-\\bm\{x\}\_\{t\}\\\|\>\\epsilon\}\\int\_\{\\mathcal\{R\}^\{D\}\}p\_\{\\bm\{Y\}\|\\bm\{X\}\_\{0\}\}\(\\bm\{y\}\|\\bm\{x\}\_\{0\}\)p\_\{\\bm\{X\}\_\{t\}\|\\bm\{X\}\_\{0\}\}\(\\bm\{u\}\|\\bm\{x\}\_\{0\}\)p\_\{\\bm\{X\}\_\{0\}\}\(\\bm\{x\}\_\{0\}\)d\\bm\{x\}\_\{0\}×‖𝒖−𝒙t‖−α−D​d​𝒖\\displaystyle\\times\\\|\\bm\{u\}\-\\bm\{x\}\_\{t\}\\\|^\{\-\\alpha\-D\}d\\bm\{u\}≈\\displaystyle\\approx∫‖𝒖−𝒙t‖\>ϵ1N​∑j=1Np𝒀\|𝑿0​\(𝒚\|𝒙0,j\)​p𝑿t\|𝑿0​\(𝒖\|𝒙0,j\)\\displaystyle\\int\_\{\\\|\\bm\{u\}\-\\bm\{x\}\_\{t\}\\\|\>\\epsilon\}\\frac\{1\}\{N\}\\sum\_\{j=1\}^\{N\}p\_\{\\bm\{Y\}\|\\bm\{X\}\_\{0\}\}\(\\bm\{y\}\|\\bm\{x\}\_\{0,j\}\)p\_\{\\bm\{X\}\_\{t\}\|\\bm\{X\}\_\{0\}\}\(\\bm\{u\}\|\\bm\{x\}\_\{0,j\}\)×‖𝒖−𝒙t‖−α−D​d​𝒖\\displaystyle\\times\\\|\\bm\{u\}\-\\bm\{x\}\_\{t\}\\\|^\{\-\\alpha\-D\}d\\bm\{u\}∝\\displaystyle\\propto∫‖𝒗‖\>ϵ1N​∑j=1Np𝒀\|𝑿0​\(𝒚\|𝒙0,j\)​p𝑿t\|𝑿0​\(𝒗\+𝒙t\|𝒙0,j\)​νt​\(d​𝒗\)\\displaystyle\\int\_\{\\\|\\bm\{v\}\\\|\>\\epsilon\}\\frac\{1\}\{N\}\\sum\_\{j=1\}^\{N\}p\_\{\\bm\{Y\}\|\\bm\{X\}\_\{0\}\}\(\\bm\{y\}\|\\bm\{x\}\_\{0,j\}\)p\_\{\\bm\{X\}\_\{t\}\|\\bm\{X\}\_\{0\}\}\(\\bm\{v\}\+\\bm\{x\}\_\{t\}\|\\bm\{x\}\_\{0,j\}\)\\nu\_\{t\}\(d\\bm\{v\}\)\(C\.1\)
By introducing the normalization factors and following procedures similar to those in \([41](https://arxiv.org/html/2608.10384#S4.E41)\) and \([IV\-D](https://arxiv.org/html/2608.10384#S4.Ex40)\), we can deduce that the target sampling distribution can also be decomposed into the conditional distribution and the weighted distributiongposterior,Q​\(𝒙\)g\_\{\\text\{posterior\},Q\}\(\\bm\{x\}\), where

gposterior,Q​\(𝒙\)=∑j=1N1𝒙0,j​\(𝒙\)​p𝒀\|𝑿0​\(𝒚\|𝒙0,j\)​Qt,ϵ,j∑j=1Np𝒀\|𝑿0​\(𝒚\|𝒙0,j\)​Qt,ϵ,jg\_\{\\text\{posterior\},Q\}\(\\bm\{x\}\)=\\sum\_\{j=1\}^\{N\}\\frac\{1\_\{\\bm\{x\}\_\{0,j\}\}\(\\bm\{x\}\)p\_\{\\bm\{Y\}\|\\bm\{X\}\_\{0\}\}\(\\bm\{y\}\|\\bm\{x\}\_\{0,j\}\)Q\_\{t,\\epsilon,j\}\}\{\\sum\_\{j=1\}^\{N\}p\_\{\\bm\{Y\}\|\\bm\{X\}\_\{0\}\}\(\\bm\{y\}\|\\bm\{x\}\_\{0,j\}\)Q\_\{t,\\epsilon,j\}\}\(C\.2\)

## References

- \[1\]\(2004\)Lévy processes and stochastic calculus\.Cambridge Studies in Advanced Mathematics,Cambridge University Press,Cambridge\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p3.1),[§II](https://arxiv.org/html/2608.10384#S2.p4.4),[§III\-A](https://arxiv.org/html/2608.10384#S3.SS1.p4.7)\.
- \[2\]M\. Arvinte and J\. I\. Tamir\(2023\-Jun\.\)MIMO channel estimation using score\-based generative models\.IEEE Transactions on Wireless Communications22\(6\),pp\. 3698–3713\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p1.1)\.
- \[3\]Y\. Chen and T\. Pock\(2017\-Jun\.\)Trainable nonlinear reaction diffusion: a flexible framework for fast and effective image restoration\.IEEE Transactions on Pattern Analysis and Machine Intelligence39\(6\),pp\. 1256–1272\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p1.1)\.
- \[4\]H\. Chung, J\. Kim, M\. T\. McCann, M\. L\. Klasky, and J\. C\. Ye\(2023\)Diffusion posterior sampling for general noisy inverse problems\.InThe Eleventh International Conference on Learning Representations \(ICLR 2023\),Kigali, Rwanda\.Cited by:[§IV\-E](https://arxiv.org/html/2608.10384#S4.SS5.p1.3)\.
- \[5\]H\. Chung, B\. Sim, D\. Ryu, and J\. C\. Ye\(2022\)Improving diffusion models for inverse problems using manifold constraints\.Advances in Neural Information Processing Systems35,pp\. 25683–25696\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p1.1)\.
- \[6\]G\. Conforti and C\. Léonard\(2022\-Feb\.\)Time reversal of markov processes with jumps under a finite entropy condition\.Stochastic Processes and their Applications144,pp\. 85–124\.Cited by:[Appendix A](https://arxiv.org/html/2608.10384#A1.p1.5),[Appendix A](https://arxiv.org/html/2608.10384#A1.p1.9)\.
- \[7\]B\. Fei, Z\. Lyu, L\. Pan, J\. Zhang, W\. Yang, T\. Luo, B\. Zhang, and B\. Dai\(2023\)Generative diffusion prior for unified image restoration and enhancement\.In2023 IEEE/CVF Conference on Computer Vision and Pattern Recognition \(CVPR\),Vol\.,pp\. 9935–9946\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p1.1)\.
- \[8\]X\. Feng, J\. Wang, X\. Kuai, M\. Zhou, H\. Sun, and J\. Li\(2022\-Jan\.\)Message passing\-based impulsive noise mitigation and channel estimation for underwater acoustic ofdm communications\.IEEE Transactions on Vehicular Technology71\(1\),pp\. 611–625\.Cited by:[5th item](https://arxiv.org/html/2608.10384#S5.I1.i5.p1.1.1)\.
- \[9\]J\. Ho, A\. Jain, and P\. Abbeel\(2020\)Denoising diffusion probabilistic models\.Advances in neural information processing systems33,pp\. 6840–6851\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p1.1)\.
- \[10\]Y\. Hu, E\. Bell, G\. Wang, and Y\. Sun\(2026\)PRISM: probabilistic and robust inverse solver with measurement\-conditioned diffusion prior for blind inverse problems\.InICASSP 2026 \- 2026 IEEE International Conference on Acoustics, Speech and Signal Processing \(ICASSP\),Vol\.,pp\. 11432–11436\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p1.1)\.
- \[11\]T\. Karras, M\. Aittala, T\. Aila, and S\. Laine\(2022\)Elucidating the design space of diffusion\-based generative models\.Advances in neural information processing systems35,pp\. 26565–26577\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p2.1)\.
- \[12\]B\. Kawar, M\. Elad, S\. Ermon, and J\. Song\(2022\)Denoising diffusion restoration models\.Advances in neural information processing systems35,pp\. 23593–23606\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p1.1)\.
- \[13\]J\. Li and C\. Wang\(2025\)Efficient diffusion posterior sampling for noisy inverse problems\.SIAM Journal on Imaging Sciences18\(2\),pp\. 1468–1492\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p1.1)\.
- \[14\]R\. Li, J\. Sun, and J\. Xue\(2026\-Dec\.\)Generative diffusion\-based bayesian modeling for universal channel estimation\.IEEE Journal on Selected Areas in Communications44\(\),pp\. 3104–3119\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p1.1)\.
- \[15\]Y\. Lipman, R\. T\. Q\. Chen, H\. Ben\-Hamu, M\. Nickel, and M\. Le\(2023\)Flow matching for generative modeling\.InThe Eleventh International Conference on Learning Representations,Kigali, Rwanda\.Cited by:[§IV](https://arxiv.org/html/2608.10384#S4.p5.1)\.
- \[16\]O\. Özdenizci and R\. Legenstein\(2023\-Aug\.\)Restoring vision in adverse weather conditions with patch\-based denoising diffusion models\.IEEE Transactions on Pattern Analysis and Machine Intelligence45\(8\),pp\. 10346–10357\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p1.1)\.
- \[17\]T\. Qi, J\. Zhang, J\. Wang, and Y\. Zhu\(2025\-Sept\.\)Bursty mixed gaussian\-impulsive noise model and parameter estimation\.IEEE Trans\. Commun\.73\(9\),pp\. 8274–8288\.Cited by:[§IV\-C](https://arxiv.org/html/2608.10384#S4.SS3.p8.9)\.
- \[18\]G\. Samorodnitsky, M\. S\. T\. Chapman, and Hall\(1994\)Stable non\-gaussian random processes: stochastic models with infinite variance\.Chapman & Hall,New York\.Cited by:[§IV\-C](https://arxiv.org/html/2608.10384#S4.SS3.p8.9)\.
- \[19\]D\. Shariatian, U\. Simsekli, and A\. Oliviero Durmus\(2025\)HEAVY\-tailed diffusion with denoising lévy probabilistic models\.InInternational Conference on Learning Representations,Vol\.2025,pp\. 96991–97024\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p2.1),[§I](https://arxiv.org/html/2608.10384#S1.p3.1)\.
- \[20\]B\. Song, S\. M\. Kwon, Z\. Zhang, X\. Hu, Q\. Qu, and L\. Shen\(2024\)Solving inverse problems with latent diffusion models via hard data consistency\.InInternational Conference on Learning Representations,Vol\.2024,pp\. 7624–7654\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p1.1)\.
- \[21\]Y\. Song and S\. Ermon\(2019\)Generative modeling by estimating gradients of the data distribution\.Advances in neural information processing systems32\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p1.1),[§I](https://arxiv.org/html/2608.10384#S1.p2.1),[§IV\-C](https://arxiv.org/html/2608.10384#S4.SS3.p2.1)\.
- \[22\]Y\. Song, J\. Sohl\-Dickstein, D\. P\. Kingma, A\. Kumar, S\. Ermon, and B\. Poole\(2021\)Score\-based generative modeling through stochastic differential equations\.International Conference on Learning Representations9\.External Links:[Link](https://mlanthology.org/iclr/2021/song2021iclr-scorebased/)Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p1.1),[§I](https://arxiv.org/html/2608.10384#S1.p2.1),[§IV\-C](https://arxiv.org/html/2608.10384#S4.SS3.p2.1)\.
- \[23\]G\. Sureka and K\. Kiasaleh\(2013\-May\.\)Sub\-optimum receiver architecture for awgn channel with symmetric alpha\-stable interference\.IEEE Trans\. Commun\.61\(5\),pp\. 1926–1935\.Cited by:[§IV\-C](https://arxiv.org/html/2608.10384#S4.SS3.p7.11),[§IV\-C](https://arxiv.org/html/2608.10384#S4.SS3.p7.19),[§IV\-C](https://arxiv.org/html/2608.10384#S4.SS3.p8.9)\.
- \[24\]P\. Vincent\(2011\-Jul\.\)A connection between score matching and denoising autoencoders\.Neural Computation23\(7\),pp\. 1661–1674\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p1.1)\.
- \[25\]Z\. Wang, Y\. Li, C\. Wang, D\. Ouyang, and Y\. Huang\(2021\-Aug\.\)A\-omp: an adaptive omp algorithm for underwater acoustic ofdm channel estimation\.IEEE Wireless Communications Letters10\(8\),pp\. 1761–1765\.Cited by:[3rd item](https://arxiv.org/html/2608.10384#S5.I1.i3.p1.1.1)\.
- \[26\]B\. Xia, Y\. Zhang, S\. Wang, Y\. Wang, X\. Wu, Y\. Tian, W\. Yang, and L\. Van Gool\(2023\)Diffir: efficient diffusion model for image restoration\.InProceedings of the IEEE/CVF international conference on computer vision,pp\. 13095–13105\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p1.1)\.
- \[27\]E\. B\. Yoon, K\. Park, S\. Kim, and S\. Lim\(2023\)Score\-based generative models with lévy processes\.Advances in neural information processing systems36,pp\. 40694–40707\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p2.1),[§I](https://arxiv.org/html/2608.10384#S1.p3.1),[§I](https://arxiv.org/html/2608.10384#S1.p4.1),[§IV\-B](https://arxiv.org/html/2608.10384#S4.SS2.p3.3),[7th item](https://arxiv.org/html/2608.10384#S5.I1.i7.p1.1.1)\.
- \[28\]Z\. Zhang, X\. Han, W\. Li, L\. Wei, and J\. Yin\(2026\-May\.\)Robust time\-correlated sparse bayesian channel estimation for short\-block underwater acoustic communications under impulsive noise\.IEEE Wireless Communications Letters15\(\),pp\. 3194–3198\.Cited by:[5th item](https://arxiv.org/html/2608.10384#S5.I1.i5.p1.1.1)\.
- \[29\]X\. Zhou, L\. Liang, J\. Zhang, P\. Jiang, Y\. Li, and S\. Jin\(2025\-Jul\.\)Generative diffusion models for high dimensional channel estimation\.IEEE Transactions on Wireless Communications24\(7\),pp\. 5840–5854\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p1.1),[6th item](https://arxiv.org/html/2608.10384#S5.I1.i6.p1.1.1)\.
- \[30\]N\. Zilberstein, A\. Sabharwal, and S\. Segarra\(2024\-Jan\.\)Solving linear inverse problems using higher\-order annealed langevin diffusion\.IEEE Transactions on Signal Processing72\(\),pp\. 492–505\.Cited by:[§I](https://arxiv.org/html/2608.10384#S1.p1.1)\.
- \[31\]D\. Zwillinger, V\. Moll, I\.S\. Gradshteyn, and I\.M\. Ryzhik \(Eds\.\)\(2014\)Table of integrals, series, and products \(eighth edition\)\.Academic Press,Boston\.Cited by:[§IV\-C](https://arxiv.org/html/2608.10384#S4.SS3.p10.7)\.