Learning Normalized Energy Models for Linear Inverse Problems

# Learning Normalized Energy Models for Linear Inverse Problems
Source: [https://arxiv.org/html/2605.15487](https://arxiv.org/html/2605.15487)
###### Abstract

Generative diffusion models can provide powerful prior probability models for inverse problems in imaging, but existing implementations suffer from two key limitations:\(i\)\(i\)the prior density is represented implicitly, and\(i​i\)\(ii\)they rely on likelihood approximations that introduce sampling biases\. We address these challenges by introducing a new energy\-based model trained for denoising with a covariance\-based regularization term that enforces consistency across different measurement conditions\. The trained model can compute normalized posterior densities for diverse linear inverse problems, without additional retraining or fine tuning\. In addition to preserving the sampling capabilities of diffusion models, this enables previously unavailable capabilities: energy\-guided adaptive sampling that adjusts schedules on\-the\-fly, unbiased Metropolis\-Hastings correction steps, and blind estimation of the degradation operator via Bayes rule\. We validate the method on multiple datasets \(ImageNet, CelebA, AFHQ\) and tasks \(inpainting, deblurring\), demonstrating competitive or superior performance to established baselines\. Code is available at[https://github\.com/nzilberstein/Anisotropic\-energy\-Model](https://github.com/nzilberstein/Anisotropic-energy-Model)\.

Machine Learning, ICML

## 1Introduction

Generative diffusion models\(Sohl\-Dicksteinet al\.,[2015](https://arxiv.org/html/2605.15487#bib.bib36); Hoet al\.,[2020](https://arxiv.org/html/2605.15487#bib.bib23); Song and Ermon,[2019](https://arxiv.org/html/2605.15487#bib.bib151); Songet al\.,[2021b](https://arxiv.org/html/2605.15487#bib.bib37)\)have achieved remarkable success in image generative modeling\. Beyond sampling, the primary use of such models is to provide prior probabilities for inverse problems\(Daraset al\.,[2024](https://arxiv.org/html/2605.15487#bib.bib92)\), such as deblurring or inpainting\. For this purpose, current solutions fall into two categories:*Bayesian*methods\(Kadkhodaie and Simoncelli,[2021](https://arxiv.org/html/2605.15487#bib.bib46); Kawaret al\.,[2022](https://arxiv.org/html/2605.15487#bib.bib49); Chunget al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib54)\)combine the prior embedded in a pre\-trained unconditional diffusion model, with measurement likelihoods via Bayes rule\. While this separation offers flexibility, the likelihood term is intractable, and thus computing solutions relies on approximations that can lead to biased sampling\(Chunget al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib54); Bruna and Han,[2024](https://arxiv.org/html/2605.15487#bib.bib129)\)\.*Regression*methods\(Sahariaet al\.,[2022](https://arxiv.org/html/2605.15487#bib.bib79); Liuet al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib122); Delbracio and Milanfar,[2023](https://arxiv.org/html/2605.15487#bib.bib124); Negrelet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib104); Elataet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib120)\)learn conditional diffusion models that bridge measurements and clean signals directly\. Such methods avoid likelihood approximations, but this comes at the cost of having a unique prior that is independent of the considered degradation\. In addition, both approaches are score\-based, rendering posterior density evaluation computationally expensive\.

We address these challenges by developing an energy\-based model \(EBM\) that learns explicit and normalized prior and posterior densities while providing efficient access to both posterior means and posterior samples\. A single trained model can be used to recover images corrupted by a diverse range of degradations in a consistent manner\. This unified formulation, which leverages the connection between linear inverse problems and anisotropic denoising, also unlocks novel capabilities that we detail below\. A number of previous publications have combined unnormalized EBMs with unconditional diffusion models\(Duet al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib128); Thorntonet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib125)\), leveraging explicit energy values for compositional generation\. While recent proposals use regularization terms that enable proper normalization\(Guthet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib114); Yuet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib115); Plaineret al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib126)\), these methods have so far been limited to isotropic noise and cannot properly handle the correlated noise that commonly arises when solving linear inverse problems\.

Here, we introduce Anisotropic Covariance Score Matching \(A\-CSM\), an extension of the “dual score\-matching” framework of\(Guthet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib114)\)that can handle anisotropic \(colored\) noise\. We develop a novel regularizer that enforces consistency by constraining the gradient of the energy with respect to the noise covariance\. Our primary contributions are:

- •a novel training objective derived from the Fokker\-Planck equation that enforces consistency of energies across noise covariance matrices;
- •a complementary architecture that supports anisotropic noise scheduling via a novel noise embedding;
- •experiments demonstrating that the resulting learned normalized densities enable new capabilities: energy\-guided adaptive sampling, unbiased MCMC correction steps, and blind inverse problem solutions;
- •experiments validating the model across multiple datasets \(ImageNet, CelebA, AFHQ\) and tasks \(inpainting, deblurring\), demonstrating that our single unified model can achieve performance competitive with or superior to established diffusion\-based solvers\.

## 2Background and related works

### 2\.1Diffusion models

Diffusion models\(Sohl\-Dicksteinet al\.,[2015](https://arxiv.org/html/2605.15487#bib.bib36); Hoet al\.,[2020](https://arxiv.org/html/2605.15487#bib.bib23); Songet al\.,[2021b](https://arxiv.org/html/2605.15487#bib.bib37)\)are derived from two continuous\-time Markov processes:\(i\)\(i\)a forward process that gradually corrupts clean data with noise, and\(i​i\)\(ii\)a reverse process that generates samples by iterative denoising\. Throughout this work, we assume a variance\-exploding \(VE\) formulation of the forward process\(Songet al\.,[2021b](https://arxiv.org/html/2605.15487#bib.bib37)\):

𝐱t=𝐱\+σt​𝐯,𝐯∼𝒩​\(0,𝐈\)\.\{\\mathbf\{x\}\}\_\{t\}=\{\\mathbf\{x\}\}\+\\sigma\_\{t\}\{\\mathbf\{v\}\},\\quad\{\\mathbf\{v\}\}\\sim\{\\mathcal\{N\}\}\(0,\{\\mathbf\{I\}\}\)\.\(1\)The corresponding discrete reverse process is

𝐱t−1=𝐱t\+δ​σt2​∇log⁡p​\(𝐱t\)\+δ​σt​𝐯t,\{\\mathbf\{x\}\}\_\{t\-1\}=\{\\mathbf\{x\}\}\_\{t\}\+\\delta\\sigma\_\{t\}^\{2\}\\nabla\\log p\(\{\\mathbf\{x\}\}\_\{t\}\)\+\\delta\\sigma\_\{t\}\{\\mathbf\{v\}\}\_\{t\},\(2\)whereδ​σt2=σt2−σt−12\\delta\\sigma\_\{t\}^\{2\}=\\sigma\_\{t\}^\{2\}\-\\sigma\_\{t\-1\}^\{2\}and the unknown score function∇𝐱tlog⁡p​\(𝐱t\)\\nabla\_\{\{\\mathbf\{x\}\}\_\{t\}\}\\log p\(\{\\mathbf\{x\}\}\_\{t\}\)is approximated by a neural networks𝜽​\(𝐱t,t\)s\_\{\\boldsymbol\{\\theta\}\}\(\{\\mathbf\{x\}\}\_\{t\},t\)trained with denoising score matching \(DSM\)\(Vincent,[2011](https://arxiv.org/html/2605.15487#bib.bib21); Raphan and Simoncelli,[2011](https://arxiv.org/html/2605.15487#bib.bib141)\)\. The DSM loss provides an upper\-bound onKL​\(p​\(𝐱\)∥p𝜽​\(𝐱\)\)\\mathrm\{KL\}\(p\(\{\\mathbf\{x\}\}\)\\,\\\|\\,p\_\{\\boldsymbol\{\\theta\}\}\(\{\\mathbf\{x\}\}\)\), and is therefore equivalent to maximum\-likelihood training\(Songet al\.,[2021a](https://arxiv.org/html/2605.15487#bib.bib68)\)for a particular weighting of the loss\. Notably, the model provides an*implicit*representation ofp​\(𝐱\)p\(\{\\mathbf\{x\}\}\), represented using a family of score functions learned for all noise levels\(Kadkhodaie and Simoncelli,[2021](https://arxiv.org/html/2605.15487#bib.bib46)\)\.

##### Beyond isotropic noise\.

While the majority of diffusion models assume an isotropic forward process, several recent studies have explored more general processes, including anisotropic processes with vector\-valued parameterizations\(Daraset al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib105); Hoogeboom and Salimans,[2023](https://arxiv.org/html/2605.15487#bib.bib106); Bansalet al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib107); Gerdeset al\.,[2024](https://arxiv.org/html/2605.15487#bib.bib109); Chenet al\.,[2024](https://arxiv.org/html/2605.15487#bib.bib108); Asthanaet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib143)\)and learnable\(Bartoshet al\.,[2024](https://arxiv.org/html/2605.15487#bib.bib110); Sahooet al\.,[2024](https://arxiv.org/html/2605.15487#bib.bib111)\)and higher\-order processes\(Dockhornet al\.,[2021](https://arxiv.org/html/2605.15487#bib.bib25); Singhalet al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib113)\)\. Here, we make use of anisotropic forward diffusion process with covariance schedule\(𝚺t\)t=0T\(\\boldsymbol\{\\Sigma\}\_\{t\}\)\_\{t=0\}^\{T\}:

𝐱t=𝐱\+𝚺t12​𝐯,𝐯∼𝒩​\(0,𝐈\),\{\\mathbf\{x\}\}\_\{t\}=\{\\mathbf\{x\}\}\+\\boldsymbol\{\\Sigma\}\_\{t\}^\{\\frac\{1\}\{2\}\}\{\\mathbf\{v\}\},\\quad\{\\mathbf\{v\}\}\\sim\{\\mathcal\{N\}\}\(0,\{\\mathbf\{I\}\}\),\(3\)with the corresponding reverse process

𝐱t−1=𝐱t\+𝜹​𝚺t​∇log⁡p​\(𝐱t\)\+𝜹​𝚺t12​𝐯t\.\{\\mathbf\{x\}\}\_\{t\-1\}=\{\\mathbf\{x\}\}\_\{t\}\+\\boldsymbol\{\\delta\}\\boldsymbol\{\\Sigma\}\_\{t\}\\nabla\\log p\(\{\\mathbf\{x\}\}\_\{t\}\)\+\\boldsymbol\{\\delta\}\\boldsymbol\{\\Sigma\}\_\{t\}^\{\\frac\{1\}\{2\}\}\{\\mathbf\{v\}\}\_\{t\}\.\(4\)

### 2\.2Linear inverse problems

Linear inverse problems consist of estimating an unknown signal𝐱\{\\mathbf\{x\}\}from a noisy measurement𝐲\{\\mathbf\{y\}\}, obtained via a degradation operator𝐇\{\\mathbf\{H\}\}\(assumed linear\):

𝐲=𝐇𝐱\+σ​𝐯,𝐯∼𝒩​\(0,𝐈\)\.\{\\mathbf\{y\}\}=\{\\mathbf\{H\}\}\{\\mathbf\{x\}\}\+\\sigma\{\\mathbf\{v\}\},\\quad\{\\mathbf\{v\}\}\\sim\{\\mathcal\{N\}\}\(0,\{\\mathbf\{I\}\}\)\.\(5\)While most solutions aim for a point estimate \(typically the posterior mean\), one can also consider sampling from the posteriorp​\(𝐱\|𝐲\)∝p​\(𝐲\|𝐱\)​p​\(𝐱\)p\(\{\\mathbf\{x\}\}\|\{\\mathbf\{y\}\}\)\\propto p\(\{\\mathbf\{y\}\}\|\{\\mathbf\{x\}\}\)p\(\{\\mathbf\{x\}\}\), wherep​\(𝐲\|𝐱\)p\(\{\\mathbf\{y\}\}\|\{\\mathbf\{x\}\}\)is the measurement model \([5](https://arxiv.org/html/2605.15487#S2.E5)\) andp​\(𝐱\)p\(\{\\mathbf\{x\}\}\)is the prior implicitly embedded in a diffusion model\. Two classes of solutions exist in the literature:*Bayesian approaches*, which approximate the priorp​\(𝐱\)p\(\{\\mathbf\{x\}\}\)via a diffusion model and combine it with the measurement model at inference time; and*regression approaches*, which approximate the posteriorp​\(𝐱\|𝐲\)p\(\{\\mathbf\{x\}\}\|\{\\mathbf\{y\}\}\)directly through supervised regression on paired data\(𝐱,𝐲\)\(\{\\mathbf\{x\}\},\{\\mathbf\{y\}\}\)\.

##### Bayesian approaches\.

These methods generate a sample from the posterior by conditioning the reverse process \([2](https://arxiv.org/html/2605.15487#S2.E2)\) on the observation𝐲\{\\mathbf\{y\}\}, through use of Bayes’ rule:

∇𝐱tlog⁡p​\(𝐱t\|𝐲\)=∇𝐱tlog⁡p​\(𝐱t\)\+∇𝐱tlog⁡p​\(𝐲\|𝐱t\)\.\\nabla\_\{\{\\mathbf\{x\}\}\_\{t\}\}\\log p\(\{\\mathbf\{x\}\}\_\{t\}\|\{\\mathbf\{y\}\}\)=\\nabla\_\{\{\\mathbf\{x\}\}\_\{t\}\}\\log p\(\{\\mathbf\{x\}\}\_\{t\}\)\+\\nabla\_\{\{\\mathbf\{x\}\}\_\{t\}\}\\log p\(\{\\mathbf\{y\}\}\|\{\\mathbf\{x\}\}\_\{t\}\)\.\(6\)While the first term is typically computed using a pre\-trainedisotropicdiffusion model, the second term \(the so\-called “guidance”\) is intractable, due to the need for high\-dimensional integration:p​\(𝐲\|𝐱t\)=∫p​\(𝐲\|𝐱\)​p​\(𝐱\|𝐱t\)​d𝐱p\(\{\\mathbf\{y\}\}\|\{\\mathbf\{x\}\}\_\{t\}\)=\\int p\(\{\\mathbf\{y\}\}\|\{\\mathbf\{x\}\}\)p\(\{\\mathbf\{x\}\}\|\{\\mathbf\{x\}\}\_\{t\}\)\\mathrm\{d\}\{\\mathbf\{x\}\}\. Previous methods\(Kadkhodaie and Simoncelli,[2021](https://arxiv.org/html/2605.15487#bib.bib46); Kawaret al\.,[2022](https://arxiv.org/html/2605.15487#bib.bib49); Chunget al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib54)\)have resorted to approximations\. But computational costs are high even for simple Gaussian approximations\(Chunget al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib54); Songet al\.,[2022](https://arxiv.org/html/2605.15487#bib.bib55)\), as the guidance term involves the*Jacobian*of the learned score\. Beyond these guidance approaches, more recent studies have explored optimization\-based approaches\(Mardaniet al\.,[2024](https://arxiv.org/html/2605.15487#bib.bib6); Zilbersteinet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib135); Zhuet al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib41); Fenget al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib95)\), which frame posterior sampling as stochastic optimization; and sequential Monte Carlo\-based\(Wuet al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib136)\)approaches, which leverage particle methods\.

##### Regression\-based approaches\.

A second category focuses on directly learning the conditional score \([6](https://arxiv.org/html/2605.15487#S2.E6)\)\. Such works learn a diffusion model between the noisy measurement𝐲\{\\mathbf\{y\}\}and the clean signal𝐱\{\\mathbf\{x\}\}, effectively treating the inverse problem as conditional generation by absorbing the degradation into the training data\(Liuet al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib122); Delbracio and Milanfar,[2023](https://arxiv.org/html/2605.15487#bib.bib124); Negrelet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib104); Huet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib119)\)\. While these methods do not explicitly input the degradation operator𝐇\{\\mathbf\{H\}\}into the score network, a few recent papers\(Elataet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib120); Terriset al\.,[2026](https://arxiv.org/html/2605.15487#bib.bib121)\)have proposed degradation\-aware parameterizations of the conditional denoiser𝔼​\[𝐱\|𝐲,𝐇\]\\mathbb\{E\}\[\{\\mathbf\{x\}\}\\,\|\\,\{\\mathbf\{y\}\},\{\\mathbf\{H\}\}\]\. However, they rely on a conditioning mechanism that requires repeated mappings between the ambient and feature spaces at each level of the architecture, increasing the computational complexity\.

## 3Learning EBMs through anisotropic denoising

We aim to learn a*single*modelpθ​\(𝐱\|𝐲\)p\_\{\\theta\}\(\{\\mathbf\{x\}\}\|\{\\mathbf\{y\}\}\)that explicitly approximates the posterior*density*for measurements𝐲\{\\mathbf\{y\}\}arising from a variety of different degradation operators, in contrast to previous score\-based approaches\. Learning the posterior density gives access not only to the conditional score \([6](https://arxiv.org/html/2605.15487#S2.E6)\) through differentiation, but also unlocks several new applications that we describe in Section[4](https://arxiv.org/html/2605.15487#S4)\. We first describe in Section[3\.1](https://arxiv.org/html/2605.15487#S3.SS1)the connection between linear inverse problems and anisotropic denoising, and how this can be leveraged to learn an energy model\. We introduce our generalized dual score matching loss in Section[3\.2](https://arxiv.org/html/2605.15487#S3.SS2)and our anisotropic energy architecture in Section[3\.3](https://arxiv.org/html/2605.15487#S3.SS3)\. Finally, we show in Section[3\.4](https://arxiv.org/html/2605.15487#S3.SS4)how our energy model can be used both as a denoiser and as a posterior sampler for linear inverse problems\.

### 3\.1Energy\-based linear inverse problem solvers

We leverage the equivalence between linear degradations and colored noise to show how an energy\-based model conditioned on different noise covariances enables access to normalized posterior densities, posterior means, and posterior samples\.

![Refer to caption](https://arxiv.org/html/2605.15487v1/x1.png)Figure 1:Illustration of possible paths in image space \(left\) and covariance space \(right\)\. Isotropic models are limited to covariance schedules that are reparameterizations of the diagonal \(orange\), while anisotropic models can explore paths where different signal components are noise\-corrupted at different rates \(red, black, and purple\)\.##### From linear inverse problems to anisotropic denoising\.

Consider a linear inverse problem of the form𝐲=𝐇𝐱\+σ​𝐯\{\\mathbf\{y\}\}=\{\\mathbf\{H\}\}\{\\mathbf\{x\}\}\+\\sigma\{\\mathbf\{v\}\}\. The posterior is unchanged by redefining𝐲\{\\mathbf\{y\}\}as𝐇−1​𝐲\{\\mathbf\{H\}\}^\{\-1\}\{\\mathbf\{y\}\}, which corresponds to observing

𝐲=𝐱\+𝚺12​𝐯′,𝐯′∼𝒩​\(0,𝐈\),\{\\mathbf\{y\}\}=\{\\mathbf\{x\}\}\+\\boldsymbol\{\\Sigma\}^\{\\frac\{1\}\{2\}\}\{\\mathbf\{v\}\}^\{\\prime\},\\quad\{\\mathbf\{v\}\}^\{\\prime\}\\sim\\mathcal\{N\}\(0,\{\\mathbf\{I\}\}\),\(7\)where𝚺=σ2​𝐇−1​\(𝐇−1\)⊤\\boldsymbol\{\\Sigma\}=\\sigma^\{2\}\{\\mathbf\{H\}\}^\{\-1\}\(\{\\mathbf\{H\}\}^\{\-1\}\)^\{\\top\}\. Linear inverse problems are thus equivalent to denoising problems for images contaminated with correlated additive noise, with the noise covariances𝚺\\boldsymbol\{\\Sigma\}depending on the linear measurement𝐇\{\\mathbf\{H\}\}\. If𝐇\{\\mathbf\{H\}\}is not invertible, it can be stabilized by addition of a small multiple of the identity, which corresponds to adding large noise in the nullspace of𝐇\{\\mathbf\{H\}\}\. In the following,𝐲\{\\mathbf\{y\}\}refers to the observation model \([7](https://arxiv.org/html/2605.15487#S3.E7)\)\.

##### Anisotropic EBMs\.

Using Bayes rule, the posterior distribution can be expressed asp​\(𝐱\|𝐲\)=p​\(𝐱\)​p​\(𝐲\|𝐱\)/p​\(𝐲\)p\(\{\\mathbf\{x\}\}\|\{\\mathbf\{y\}\}\)=p\(\{\\mathbf\{x\}\}\)p\(\{\\mathbf\{y\}\}\|\{\\mathbf\{x\}\}\)/p\(\{\\mathbf\{y\}\}\)\. Its normalization requires having a model ofp​\(𝐲\)p\(\{\\mathbf\{y\}\}\), which depends on𝚺\\boldsymbol\{\\Sigma\}\. To this end, we propose to learn an energy modelUθ​\(𝐲,𝚺\)U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)that approximatesU​\(𝐲,𝚺\)=−log⁡p​\(𝐲\|𝚺\)U\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)=\-\\log p\(\{\\mathbf\{y\}\}\|\\boldsymbol\{\\Sigma\}\)\. Notice that we can recoverp​\(𝐱\)p\(\{\\mathbf\{x\}\}\)by setting𝚺=0\\boldsymbol\{\\Sigma\}=0\.

The posterior mean can be computed using an anisotropic generalization of Tweedie’s identity:\(Robbins,[1956](https://arxiv.org/html/2605.15487#bib.bib139); Miyasawa and others,[1961](https://arxiv.org/html/2605.15487#bib.bib140); Raphan and Simoncelli,[2011](https://arxiv.org/html/2605.15487#bib.bib141)\):

𝔼​\[𝐱\|𝐲,𝚺\]=𝐲−𝚺​∇𝐲U​\(𝐲,𝚺\),\\mathbb\{E\}\[\{\\mathbf\{x\}\}\\,\|\\,\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\]=\{\\mathbf\{y\}\}\-\\boldsymbol\{\\Sigma\}\\nabla\_\{\{\\mathbf\{y\}\}\}U\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\),\(8\)where the score can be recovered from the energy model as∇𝐲U​\(𝐲,𝚺\)≈∇𝐲Uθ​\(𝐲,𝚺\)\\nabla\_\{\{\\mathbf\{y\}\}\}U\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\\approx\\nabla\_\{\{\\mathbf\{y\}\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\. Therefore, by*learning an anisotropic denoiser*, we gain access to the posterior mean for different linear inverse problems\.

Learning a score model over a family of covariances𝚺\\boldsymbol\{\\Sigma\}enables posterior sampling\. This is achieved by defining a noise schedule that bridges between measurements and the corresponding clean data\. In particular, we need to specify a sequence of covariance matrices\(𝚺t\)t=0T\(\\boldsymbol\{\\Sigma\}\_\{t\}\)\_\{t=0\}^\{T\}that interpolate between the measurement covariance𝚺\\boldsymbol\{\\Sigma\}att=Tt=Tand the zero covariance att=0t=0, as in \([4](https://arxiv.org/html/2605.15487#S2.E4)\)\. Several potential paths are illustrated in Figure[1](https://arxiv.org/html/2605.15487#S3.F1)\. Note that unlike the isotropic setting, these schedules do not simply correspond to different discretizations of a correponsding continuous\-time process\.

### 3\.2Learning with anisotropic dual score matching

We learn a model of the energyU​\(𝐲,𝚺\)=−log⁡p​\(𝐲\|𝚺\)U\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)=\-\\log p\(\{\\mathbf\{y\}\}\|\\boldsymbol\{\\Sigma\}\)by developing an anisotropic generalization of the \(isotropic\) dual score matching approach ofGuthet al\.\([2025](https://arxiv.org/html/2605.15487#bib.bib114)\)\. In a nutshell, dual score matching consists of matching*both*derivatives of an energy modelUθ​\(𝐲,𝚺\)U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)to the data, i\.e\., the gradients of the energy with respect to*\(i\)*the noisy input𝐲\{\\mathbf\{y\}\}\(the “data score” used in diffusion models\) and*\(ii\)*the covariance𝚺\\boldsymbol\{\\Sigma\}\(which we dub the “covariance score”\)\. Each of these derivatives can be trained by minimizing a corresponding loss function\. As we explain below, the dual score matching formulation enables learning of*normalized*densities\.

##### Data score\.

The data score can be learned using the anisotropic generalization of Tweedie’s formula \([8](https://arxiv.org/html/2605.15487#S3.E8)\), which leads to the anisotropic denoising score matching objective:

ℓA−DSM\(θ\)=𝔼\[∥𝚺12∇𝐲Uθ\(𝐲,𝚺\)−𝚺−12\(𝐲−𝐱\)∥2\]\.\\ell\_\{\\mathrm\{A\-DSM\}\}\(\\theta\)=\\mathbb\{E\}\\mathopen\{\}\\mathclose\{\{\\left\[\\lVert\\boldsymbol\{\\Sigma\}^\{\\frac\{1\}\{2\}\}\\nabla\_\{\{\\mathbf\{y\}\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\-\\boldsymbol\{\\Sigma\}^\{\-\\frac\{1\}\{2\}\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)\\rVert^\{2\}\}\}\\right\]\.\(9\)The reweighting by𝚺12\\boldsymbol\{\\Sigma\}^\{\\frac\{1\}\{2\}\}leads to a scale\-invariant loss, a generalization of the maximum\-likelihood weighting used in\(Songet al\.,[2021a](https://arxiv.org/html/2605.15487#bib.bib68); Guthet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib114)\)\.

##### Covariance score\.

The covariance score can be shown to satisfy an analogue of the Tweedie formula \(see Appendix[A](https://arxiv.org/html/2605.15487#A1)\):

∇𝚺U\(𝐲,𝚺\)=𝔼\[12𝚺−1−12𝚺−1\(𝐲−𝐱\)\(𝐲−𝐱\)⊤𝚺−1\]\.\\nabla\_\{\\boldsymbol\{\\Sigma\}\}U\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)=\\mathbb\{E\}\\mathopen\{\}\\mathclose\{\{\\left\[\\frac\{1\}\{2\}\\boldsymbol\{\\Sigma\}^\{\-1\}\-\\frac\{1\}\{2\}\\boldsymbol\{\\Sigma\}^\{\-1\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)^\{\\top\}\\boldsymbol\{\\Sigma\}^\{\-1\}\}\}\\right\]\.This expression generalizes the time score identity ofGuthet al\.\([2025](https://arxiv.org/html/2605.15487#bib.bib114)\); Yuet al\.\([2025](https://arxiv.org/html/2605.15487#bib.bib115)\); Plaineret al\.\([2025](https://arxiv.org/html/2605.15487#bib.bib126)\), and leads to a covariance score matching objective:

ℓA−CSM\(θ\)=𝔼𝐱,𝐲,𝚺\[∥𝚺12∇𝚺Uθ\(𝐲,𝚺\)𝚺12−12𝐈\+12𝚺−12\(𝐲−𝐱\)\(𝐲−𝐱\)⊤𝚺−12∥22\],\\ell\_\{\\mathrm\{A\-CSM\}\}\(\\theta\)=\\mathbb\{E\}\_\{\{\\mathbf\{x\}\},\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\}\\Big\[\\big\\lVert\\boldsymbol\{\\Sigma\}^\{\\frac\{1\}\{2\}\}\\nabla\_\{\\boldsymbol\{\\Sigma\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\\boldsymbol\{\\Sigma\}^\{\\frac\{1\}\{2\}\}\-\\\\ \\frac\{1\}\{2\}\{\\mathbf\{I\}\}\+\\frac\{1\}\{2\}\\boldsymbol\{\\Sigma\}^\{\-\\frac\{1\}\{2\}\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)^\{\\top\}\\boldsymbol\{\\Sigma\}^\{\-\\frac\{1\}\{2\}\}\\big\\rVert\_\{2\}^\{2\}\\Big\],\(10\)that has been weighted by𝚺12\\boldsymbol\{\\Sigma\}^\{\\frac\{1\}\{2\}\}on both sides to be scale\-invariant, and where∥⋅∥2\\\|\\cdot\\\|\_\{2\}is the Frobenius norm\. The A\-CSM objective \([10](https://arxiv.org/html/2605.15487#S3.E10)\) acts as a regularizer that ties all the different marginalsp​\(𝐲\|𝚺\)p\(\{\\mathbf\{y\}\}\|\\boldsymbol\{\\Sigma\}\)together in a meaningful way\. It can be thought of as a proxy to enforce the continuity equation given by the Fokker\-Planck equation\(Pavliotis,[2014](https://arxiv.org/html/2605.15487#bib.bib127)\)\. This allows learning of a normalized density \(as we show below\), but also improves the denoising performance of the model, as illustrated in Appendix[D\.2](https://arxiv.org/html/2605.15487#A4.SS2)\.

##### Normalization of the learned density model\.

We define an overall objective as a weighted sum of the two objectives defined above:1d​ℓA−DSM\+1d2​ℓA−CSM\\frac\{1\}\{d\}\\ell\_\{\\mathrm\{A\-DSM\}\}\+\\frac\{1\}\{d^\{2\}\}\\ell\_\{\\mathrm\{A\-CSM\}\}\(see Appendix[A](https://arxiv.org/html/2605.15487#A1)\)\. After training, the learned energy modelUθ​\(𝐲,𝚺\)U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)approximates the true energy up to a constant:Uθ​\(𝐲,𝚺\)≈U​\(𝐲,𝚺\)\+cstU\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\\approx U\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\+\\mathrm\{cst\}\. Importantly, the A\-CSM objective ensures mass conservation through the space of covariances, and thus this constant is independent of𝚺\\boldsymbol\{\\Sigma\}\. For large𝚺\\boldsymbol\{\\Sigma\},𝐲\{\\mathbf\{y\}\}is approximately distributed as𝒩​\(0,𝚺\)\\mathcal\{N\}\(0,\\boldsymbol\{\\Sigma\}\), allowing to normalize the trained energy model as follows:

Uθ​\(𝐲,𝚺\)→Uθ​\(𝐲,𝚺\)−𝔼𝐲​\[Uθ​\(𝐲,𝚺\)∣𝚺\]\+12​log​det\(2​π​e​𝚺\)\.U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\\rightarrow U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\\\\ \-\\mathbb\{E\}\_\{\{\\mathbf\{y\}\}\}\[U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\\mid\\boldsymbol\{\\Sigma\}\]\+\\frac\{1\}\{2\}\\log\\det\(2\\pi\\mathrm\{e\}\\boldsymbol\{\\Sigma\}\)\.\(11\)

### 3\.3Architecture and covariance conditioning

We computeUθ​\(𝐲,𝚺\)U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)with a neural network, with a novel architecture that is suitable for the problem\. Since our proposed energy model is trained as an anisotropic denoiser, the first requirement is that the score∇𝐲Uθ​\(𝐲,𝚺\)\\nabla\_\{\{\\mathbf\{y\}\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)derived from the energy preserves the inductive biases of score architectures used in previous literature to achieve high\-quality denoising \(and thus, score estimation\)\. FollowingRomanoet al\.\([2017](https://arxiv.org/html/2605.15487#bib.bib45)\); Guthet al\.\([2025](https://arxiv.org/html/2605.15487#bib.bib114)\); Thorntonet al\.\([2025](https://arxiv.org/html/2605.15487#bib.bib125)\), we define

Uθ​\(𝐲,𝚺\)=12​⟨𝐲,𝐬𝜽​\(𝐲,𝚺\)⟩,U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)=\\frac\{1\}\{2\}\\langle\{\\mathbf\{y\}\},\{\\mathbf\{s\}\}\_\{\\boldsymbol\{\\theta\}\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\\rangle,\(12\)where𝐬θ\{\\mathbf\{s\}\}\_\{\\theta\}is an existing score network\. Specifically, we implement𝐬θ\{\\mathbf\{s\}\}\_\{\\theta\}using the UNet architecture ofSonget al\.\([2021b](https://arxiv.org/html/2605.15487#bib.bib37)\); Karraset al\.\([2022](https://arxiv.org/html/2605.15487#bib.bib118)\)\.

The second requirement is that the architecture can operate under conditioning with a wide range of covariances𝚺\\boldsymbol\{\\Sigma\}arising from different linear inverse problems\. However, an arbitrary covariance hasd​\(d−1\)/2d\(d\-1\)/2degrees of freedom, which is problematic in terms of both memory and computational cost\. To alleviate this, we limit ourselves to covariances that are diagonal in either the spatial or spatial frequency domain\. This reduces the parameterization toddcoefficients, while still covering many standard inverse problems: inpainting \(block\-diagonal covariances in pixel space\), deblurring \(diagonal in frequency domain\), and super\-resolution \(approximately diagonal in frequency domain\)\.

Finally, we must define a mechanism by which the conditioning covariance matrix𝚺\\boldsymbol\{\\Sigma\}is incorporated into the score network\. Our design builds on the multiplicative conditioning mechanism used in isotropic score architectures via gain control\(Karraset al\.,[2024](https://arxiv.org/html/2605.15487#bib.bib154)\)\. At each layerℓ\\ell, consisting ofcℓc\_\{\\ell\}feature channels, an embedding vector𝐞ℓ∈ℝcℓ\{\\mathbf\{e\}\}\_\{\\ell\}\\in\{\\mathbb\{R\}\}^\{c\_\{\\ell\}\}of the input noise varianceσ2\\sigma^\{2\}is computed and used to modulate \(multiply\) the channels\. We represent spatial covariance matrices as spatially varying noise maps inℝd\{\\mathbb\{R\}\}^\{d\}, yielding*spatially varying embeddings*𝐞ℓ∈ℝcℓ×dℓ\{\\mathbf\{e\}\}\_\{\\ell\}\\in\{\\mathbb\{R\}\}^\{c\_\{\\ell\}\\times d\_\{\\ell\}\}, wheredℓd\_\{\\ell\}denotes the spatial resolution at layerℓ\\ell\. Spectral covariance matrices are incorporated through an analogous mechanism in the frequency domain, producing embeddings that modulate the corresponding feature channels \(but not the spatial dimensions\)\. At each layer, conditioning is applied through a gain modulation of the form𝐱ℓ←SiLU\(𝐱ℓ⊙\(1\+𝐞ℓ\)\)\{\\mathbf\{x\}\}\_\{\\ell\}\\leftarrow\\mathrm\{SiLU\}\\mathopen\{\}\\mathclose\{\{\\left\(\{\\mathbf\{x\}\}\_\{\\ell\}\\odot\\mathopen\{\}\\mathclose\{\{\\left\(1\+\{\\mathbf\{e\}\}\_\{\\ell\}\}\}\\right\)\}\}\\right\), whereSiLU\(\.\)\\mathrm\{SiLU\}\(\.\)is the swish function\(Elfwinget al\.,[2018](https://arxiv.org/html/2605.15487#bib.bib148)\)\. The resulting architecture is schematically illustrated in Fig\.[2](https://arxiv.org/html/2605.15487#S3.F2)\. Additional implementation details are provided in Appendix[C\.1](https://arxiv.org/html/2605.15487#A3.SS1)\. We emphasize that this new embedding module does not introduce significant additional computational overhead relative to those used for isotropic noiseKarraset al\.\([2022](https://arxiv.org/html/2605.15487#bib.bib118)\)\.

![Refer to caption](https://arxiv.org/html/2605.15487v1/x2.png)Figure 2:Proposed architecture based on UNet\. We incorporate the covariance information through the embedding network, with two dedicated branches for the two covariance domains \(spatial and spectral\)\.##### Training and sampling\.

We now describe our implementation and training pipeline\. We assume a covariance distributionp​\(𝚺\)p\(\\boldsymbol\{\\Sigma\}\)for training in which the two classes \(spatial and frequency covariances\) are selected with probability0\.50\.5\. In the spatial domain, we consider center\-box and horizontal\-box masks with sizes from 1 tod2=642d^\{2\}=64^\{2\}\. In the spectral domain, we consider Gaussian deblurring and×4\\times 4super\-resolution\. Additional details regarding hyperparameters and examples of the different covariances are provided in Appendix[C](https://arxiv.org/html/2605.15487#A3), and the full training and sampling algorithms are provided in Appendix[A\.1](https://arxiv.org/html/2605.15487#A1.SS1)\.

### 3\.4Validation: comparison to baselines

##### Experimental setting\.

We compare our energy\-based approach with several well\-established methods\(Daraset al\.,[2024](https://arxiv.org/html/2605.15487#bib.bib92)\)\. Among Bayesian baselines, we include the guidance methods DPS\(Chunget al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib54)\)and DAPS\(Zhanget al\.,[2025a](https://arxiv.org/html/2605.15487#bib.bib130)\)and the variational inference\-based RED\-Diff\(Mardaniet al\.,[2024](https://arxiv.org/html/2605.15487#bib.bib6)\)\. We also include Palette\(Sahariaet al\.,[2022](https://arxiv.org/html/2605.15487#bib.bib79)\)as a conditional model trained to solve a particular inverse problem via regression\. Details of each method are provided Appendix[C\.4](https://arxiv.org/html/2605.15487#A3.SS4)\. We implement all methods, including our own, using the base architecture illustrated in Fig\.[2](https://arxiv.org/html/2605.15487#S3.F2)\. The differences lie inp​\(𝚺\)p\(\\boldsymbol\{\\Sigma\}\)—Bayesian models have𝚺=σ2​𝐈\\boldsymbol\{\\Sigma\}=\\sigma^\{2\}\{\\mathbf\{I\}\}—and in the input provided to the model—conditional models stack the measurements with the input noisy image\. Models were trained on CelebA \(64×6464\\times 64\)\(Liuet al\.,[2015](https://arxiv.org/html/2605.15487#bib.bib144)\), ImageNet64 \(64×6464\\times 64\)\(Russakovskyet al\.,[2015](https://arxiv.org/html/2605.15487#bib.bib81)\), AFHQ\-Cat \(192×192192\\times 192\)\(Choiet al\.,[2020](https://arxiv.org/html/2605.15487#bib.bib146)\), and MNIST \(28×2828\\times 28\)\(LeCunet al\.,[2010](https://arxiv.org/html/2605.15487#bib.bib145)\)\.

We assess the reconstruction quality of the different samplers with mean squared error \(in terms ofPSNR=−10​log10⁡MSE\\mathrm\{PSNR\}=\-10\\log\_\{10\}\\mathrm\{MSE\}\) and two perceptual metrics, LPIPS\(Zhanget al\.,[2018](https://arxiv.org/html/2605.15487#bib.bib131)\)and DISTS\(Dinget al\.,[2020](https://arxiv.org/html/2605.15487#bib.bib147)\)\. Each metric is averaged over 400 examples with ground truth images from the test set\. In addition, we compute the FID\(Heuselet al\.,[2017](https://arxiv.org/html/2605.15487#bib.bib132)\)between ensembles of test set images and posterior samples\. All methods use up to1,0001,000neural function evaluations \(1,2001,200for inpainting on CelebA\)\. For each method, the number of steps is tuned via ablation and then fixed for all experiments; no ground\-truth information is used at test time\. This is especially relevant for RED\-Diff, which performs better with fewer noise levels, consistent withMardaniet al\.\([2024](https://arxiv.org/html/2605.15487#bib.bib6)\)\.

Table 1:Quantitative results for inpainting and Gaussian deblurring across the CelebA and ImageNet64 datasets\.
##### Experimental results\.

We consider the two linear inverse problems that were used to jointly train the model\. To test the spatial covariance configuration, we perform inpainting on a square of size45×4545\\times 45in the center of the image, withσ=10−4\\sigma=10^\{\-4\}\. For the spectral configuration, we perform deblurring with a Gaussian kernel of size8×88\\times 8and a standard deviation of0\.80\.8, withσ=10−2\\sigma=10^\{\-2\}\. Results are shown in Table[1](https://arxiv.org/html/2605.15487#S3.T1)for both CelebA and ImageNet, and in Fig\.[3](https://arxiv.org/html/2605.15487#S3.F3)we illustrate a few inpainting examples on CelebA\. We observe that our proposed sampler is competitive with, and in most cases outperforms, the baseline samplers\. We include additional comparisons on AFHQ\-Cat \(192×192192\\times 192\)\(Choiet al\.,[2020](https://arxiv.org/html/2605.15487#bib.bib146)\)in Appendix[D\.1\.3](https://arxiv.org/html/2605.15487#A4.SS1.SSS3)\. Overall, these experiments demonstrate that energy models are fully competitive with or even improves on score models for solving linear inverse problems, and that a single trained model can be used to tackle a variety of inverse problems\.

![Refer to caption](https://arxiv.org/html/2605.15487v1/x3.png)Figure 3:Inpainting results on CelebA\. DPS and our method generate sharper images while RED\-Diff’s solutions are smoother, which is consistent with our analysis of their posterior density values in Section[4\.1](https://arxiv.org/html/2605.15487#S4.SS1)\.

## 4Applications of anisotropic energy models

Our energy model offers benefits beyond its use in solving linear inverse problems\. We demonstrate in Section[4\.1](https://arxiv.org/html/2605.15487#S4.SS1)that having access to the learned energy enables more refined comparisons between posterior sampling algorithms\. In Sections[4\.2\.1](https://arxiv.org/html/2605.15487#S4.SS2.SSS1)and[4\.2\.2](https://arxiv.org/html/2605.15487#S4.SS2.SSS2), we develop samplers that leverage access to the energy to automatically define adaptive anisotropic schedules and use*unbiased*corrector steps via Metropolis\-Hastings proposals\. Finally, we show in Section[4\.3](https://arxiv.org/html/2605.15487#S4.SS3)that our energy model can estimate unknown degradation operators and noise levels, enabling blind inverse problem solvers\.

### 4\.1The probabilities of posterior samples

In this section, we compare the prior and posterior probabilities of samples generated by three solvers, which include a Bayesian model that relies on likelihood approximation \(DPS\), a variational inference procedure that approximates maximum\-a\-posterior optimization \(RED\-Diff\), and a learned posterior sampler \(ours\)\. We fix one particular observation𝐲\{\\mathbf\{y\}\}and generate several potential solutions\{𝐱^i\}i=1N\\\{\{\\hat\{\\mathbf\{x\}\}\}\_\{i\}\\\}\_\{i=1\}^\{N\}with each sampler\. We then compute the probabilities of these samples under the priorp​\(𝐱\)p\(\{\\mathbf\{x\}\}\)and the posteriorp​\(𝐱\|𝐲\)p\(\{\\mathbf\{x\}\}\|\{\\mathbf\{y\}\}\)\(via Bayes rule\) using our learned energy model\. The resulting distributions of log probabilities are shown in Fig\.[4](https://arxiv.org/html/2605.15487#S4.F4)\. Examples of solutions at different prior probabilities are shown in the right panel: generally, smoother images tend to have higher probability, while more detailed, textured images have lower probability\. In particular, we see that DPS generates samples that are lower\-probability under the prior due to their likelihood approximations, whereas RED\-Diff generates higher\-probability samples due to its a maximum\-a\-posterior behavior\. Both methods produce samples of low posterior probability\. Our energy\-based sampler achieves a better balance: notably, the prior probabilities of solutions lie near the ground truth, and their posterior probability is only slightly lower than that of the ground truth\. These results demonstrate that our anisotropic energy\-based approach leads to accurate posterior samples that are consistent with both the prior and the measurements\. A similar analysis across multiple observations𝐲\{\\mathbf\{y\}\}pairs is provided in Appendix[D\.3](https://arxiv.org/html/2605.15487#A4.SS3)\.

![Refer to caption](https://arxiv.org/html/2605.15487v1/x4.png)Figure 4:Left and middle:Histograms oflog⁡pθ​\(𝐱^\)\\log p\_\{\\theta\}\(\{\\hat\{\\mathbf\{x\}\}\}\)andlog⁡pθ​\(𝐱^\|𝐲\)\\log p\_\{\\theta\}\(\{\\hat\{\\mathbf\{x\}\}\}\|\{\\mathbf\{y\}\}\)for inpainting solutions𝐱^\{\\hat\{\\mathbf\{x\}\}\}generated by DPS, RED\-Diff, and our energy model from a given measurement𝐲\{\\mathbf\{y\}\}, along with the ground truth𝐱\{\\mathbf\{x\}\}\. Our energy model is well\-calibrated with respect to both prior and posterior probabilities\.Right:Examples of generated images𝐱^\{\\hat\{\\mathbf\{x\}\}\}sorted from lowest to highest prior probability \(in reading order\)\. The colored border indicates the sampler used to produce the image \(green for DPS, red for RED\-Diff, blue for ours\) and its prior probability \(darker shades for lower probability images\), matching the arrows in the left panel\. Note that higher\-probability images are smoother and less detailed\. The ground truth and measurements are shown on the left with blue/black borders, respectively\.
### 4\.2Energy\-based schedules and corrector steps

#### 4\.2\.1Energy\-guided generation

An important feature of anisotropic diffusion is the freedom it allows in designing sampling paths\(Negrelet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib104); Gerdeset al\.,[2024](https://arxiv.org/html/2605.15487#bib.bib109)\), as illustrated in Fig\.[1](https://arxiv.org/html/2605.15487#S3.F1)\. For any sequence of covariance steps𝜹​𝚺t=𝚺t−𝚺t−1\\boldsymbol\{\\delta\}\\boldsymbol\{\\Sigma\}\_\{t\}=\\boldsymbol\{\\Sigma\}\_\{t\}\-\\boldsymbol\{\\Sigma\}\_\{t\-1\}, a sample can be generated with the iterations

\{𝐱t−1=𝐱t−𝜹​𝚺t​∇𝐲Uθ​\(𝐱t,𝚺t\)\+\(𝜹​𝚺t\)12​𝐯t𝚺t−1=𝚺t−𝜹​𝚺t\\begin\{cases\}\{\\mathbf\{x\}\}\_\{t\-1\}=\{\\mathbf\{x\}\}\_\{t\}\-\\boldsymbol\{\\delta\}\\boldsymbol\{\\Sigma\}\_\{t\}\\nabla\_\{\{\\mathbf\{y\}\}\}U\_\{\\theta\}\(\{\\mathbf\{x\}\}\_\{t\},\\boldsymbol\{\\Sigma\}\_\{t\}\)\+\(\\boldsymbol\{\\delta\}\\boldsymbol\{\\Sigma\}\_\{t\}\)^\{\\frac\{1\}\{2\}\}\{\\mathbf\{v\}\}\_\{t\}\\\\ \\boldsymbol\{\\Sigma\}\_\{t\-1\}=\\boldsymbol\{\\Sigma\}\_\{t\}\-\\boldsymbol\{\\delta\}\\boldsymbol\{\\Sigma\}\_\{t\}\\end\{cases\}\(13\)Designing samplers therefore reduces to selecting the*covariance step*𝜹​𝚺t\\boldsymbol\{\\delta\}\\boldsymbol\{\\Sigma\}\_\{t\}at each iteration\. After demonstrating this flexibility in the setting of any\-order generation, we introduce an energy\-guided sampler that selectsδ​𝚺t\\delta\\boldsymbol\{\\Sigma\}\_\{t\}automatically and adaptively by exploiting the variations of the energy with the noise covariance𝚺\\boldsymbol\{\\Sigma\}, leading to superior sampling quality\.

##### Any\-order generation\.

We consider a family of block\-diagonal covariances in pixel space\. Each covariance from the family has a fixed variance on eachb×bb\\times bimage patch, wherebbcan take any value in\{1,2,4,7,14,28\}\\\{1,2,4,7,14,28\\\}\(see Fig\.[9](https://arxiv.org/html/2605.15487#A3.F9)in Appendix[C\.2](https://arxiv.org/html/2605.15487#A3.SS2)for an illustration\)\. We trained a model on MNIST with this set of covariances, where the noise variance in each patch was independent and log uniformly distributed in\[10−9,103\]\[10^\{\-9\},10^\{3\}\]\. This allows generating image patches in any order in an autoregressive manner\. We illustrate three different orders with different values ofbbin Fig\.[5](https://arxiv.org/html/2605.15487#S4.F5)\. Each path generates a different solution even with identical random seeds\.

![Refer to caption](https://arxiv.org/html/2605.15487v1/x5.png)
![Refer to caption](https://arxiv.org/html/2605.15487v1/x6.png)
![Refer to caption](https://arxiv.org/html/2605.15487v1/x7.png)

Figure 5:Three autoregressive generations for a model trained on MNIST, with identical initial and injected noise\. The first two rows generate the four quadrants of the image respectively in reading and reversed reading order\. The third row generates the image in 16 patches, ordered from the center to the outer border\.
##### Energy\-guided adaptive sampler\.

In theory, with a perfect score and an infinite number of sampling steps, all covariance schedules are equivalent\. In practice, different schedules yield different approximations to the score and discretization errors\. Which covariance paths lead to higher sampling quality, and how can we define them? We introduce a heuristically motivated*energy\-guided*schedule that sets the covariance step𝜹​𝚺t\\boldsymbol\{\\delta\}\\boldsymbol\{\\Sigma\}\_\{t\}adaptively by doing a descent step on the energy with respect to𝚺t\\boldsymbol\{\\Sigma\}\_\{t\}, using the covariance score∇𝚺Uθ​\(𝐱t,𝚺t\)\\nabla\_\{\\boldsymbol\{\\Sigma\}\}U\_\{\\theta\}\(\{\\mathbf\{x\}\}\_\{t\},\\boldsymbol\{\\Sigma\}\_\{t\}\)\. Specifically, we set

𝜹​𝚺t∝𝚺t​∇𝚺Uθ​\(𝐱t,𝚺t\)​𝚺t\.\\boldsymbol\{\\delta\}\\boldsymbol\{\\Sigma\}\_\{t\}\\propto\\boldsymbol\{\\Sigma\}\_\{t\}\\nabla\_\{\\boldsymbol\{\\Sigma\}\}U\_\{\\theta\}\(\{\\mathbf\{x\}\}\_\{t\},\\boldsymbol\{\\Sigma\}\_\{t\}\)\\boldsymbol\{\\Sigma\}\_\{t\}\.\(14\)The multiplication by𝚺t\\boldsymbol\{\\Sigma\}\_\{t\}on both sides can be thought of as a preconditioning of the gradient \(notice that𝜹​𝚺t\\boldsymbol\{\\delta\}\\boldsymbol\{\\Sigma\}\_\{t\}has the same units as𝚺t\\boldsymbol\{\\Sigma\}\_\{t\}\) and can be motivated as a steepest descent in the Bregman geometry generated by−log​det𝚺t\-\\log\\det\\boldsymbol\{\\Sigma\}\_\{t\}\(see Appendix[B](https://arxiv.org/html/2605.15487#A2)\)\.[Equation14](https://arxiv.org/html/2605.15487#S4.E14)removes the need for choosing a discretization schedule during sampling and instead defines one automatically and adaptively to the current iterate𝐱t\{\\mathbf\{x\}\}\_\{t\}, a property shared by the algorithm ofKadkhodaie and Simoncelli \([2021](https://arxiv.org/html/2605.15487#bib.bib46)\)but in a different setting\.

To illustrate the behavior of this sampler, we consider the task of reconstructing MNIST digits from a subset ofkkrandomly selected pixels \(unobserved pixels are replaced by high\-variance noise\)\. Fig\.[6](https://arxiv.org/html/2605.15487#S4.F6)shows the classification error over1,0001,000samples as a function of the number of measurementskk\(experimental details provided in Appendix[C\.1](https://arxiv.org/html/2605.15487#A3.SS1)\)\. Both samplers converge to clean\-image baseline performance for largekk, but the energy\-guided sampler \([14](https://arxiv.org/html/2605.15487#S4.E14)\) achieves a better classification error than a geometric schedule𝜹​𝚺t∝𝚺t\\boldsymbol\{\\delta\}\\boldsymbol\{\\Sigma\}\_\{t\}\\propto\\boldsymbol\{\\Sigma\}\_\{t\}for smallerkk\. This illustrates the advantage of allowing the model to unmask different pixels at different rates in challenging inverse problems\.

Preliminary experiments on a CelebA inpainting task indicate that the energy\-guided sampler underperforms predefined schedules in this setting\. We hypothesize that on more challenging tasks with more spatial variability, diagonal covariances in the spatial domain are not flexible enough to define meaningful paths\. We believe that more general covariance spaces that are not restricted to a fixed basis should lead to gains similar to those we observe on MNIST, and leave exploration of those questions to future research\.

![Refer to caption](https://arxiv.org/html/2605.15487v1/x8.png)Figure 6:Classification error as a function of the number of measurements\. The energy\-guided scheduling \([14](https://arxiv.org/html/2605.15487#S4.E14), orange\) consistently achieves a lower error rate than the fixed geometric scheduling𝜹​𝚺∝𝚺\\boldsymbol\{\\delta\}\\boldsymbol\{\\Sigma\}\\propto\\boldsymbol\{\\Sigma\}\(blue\) with few measurements, with both methods reaching the baseline classifier error at300300measurements\.

#### 4\.2\.2Unbiased corrector steps with MALA

Typically, diffusion\-based posterior samplers iterate between prediction and correction steps\(Songet al\.,[2021b](https://arxiv.org/html/2605.15487#bib.bib37); Bradley and Nakkiran,[2025](https://arxiv.org/html/2605.15487#bib.bib137)\), where the correction is given by a version of the*Unadjusted Langevin Algorithm*\(ULA\):

𝐱t′=𝐱t−η2\[∇𝐲Uθ\(𝐱t,𝚺t\)−∇𝐱tlogp\(𝐱T\|𝐱t\)\]\+η𝐯′t\.\{\\mathbf\{x\}\}\_\{t\}^\{\\prime\}=\{\\mathbf\{x\}\}\_\{t\}\-\\frac\{\\eta\}\{2\}\\mathopen\{\}\\mathclose\{\{\\left\[\\nabla\_\{\{\\mathbf\{y\}\}\}U\_\{\\theta\}\(\{\\mathbf\{x\}\}\_\{t\},\\boldsymbol\{\\Sigma\}\_\{t\}\)\-\\nabla\_\{\{\\mathbf\{x\}\}\_\{t\}\}\\log p\(\{\\mathbf\{x\}\}\_\{T\}\|\{\\mathbf\{x\}\}\_\{t\}\)\}\}\\right\]\+\\sqrt\{\\eta\}\{\\mathbf\{v\}\}^\{\\prime\}\_\{t\}\.\(15\)where𝐱T=𝐲\{\\mathbf\{x\}\}\_\{T\}=\{\\mathbf\{y\}\}andp​\(𝐱T\|𝐱t\)∼𝒩​\(𝐱T;𝐱t,ΣT−Σt\)p\(\{\\mathbf\{x\}\}\_\{T\}\|\{\\mathbf\{x\}\}\_\{t\}\)\\sim\{\\mathcal\{N\}\}\(\{\\mathbf\{x\}\}\_\{T\};\{\\mathbf\{x\}\}\_\{t\},\\Sigma\_\{T\}\-\\Sigma\_\{t\}\)\. Notice that we need an additional term in the correction step as we aim to sample fromp​\(𝐱t\|𝚺t,𝐲\)p\(\{\\mathbf\{x\}\}\_\{t\}\|\\boldsymbol\{\\Sigma\}\_\{t\},\{\\mathbf\{y\}\}\)rather than justp​\(𝐱t\|𝚺t\)p\(\{\\mathbf\{x\}\}\_\{t\}\|\\boldsymbol\{\\Sigma\}\_\{t\}\)\. Since the energy model providesUθ​\(𝐱t,𝚺t\)≈−log⁡p​\(𝐱t\|𝚺t\)U\_\{\\theta\}\(\{\\mathbf\{x\}\}\_\{t\},\\boldsymbol\{\\Sigma\}\_\{t\}\)\\approx\-\\log p\(\{\\mathbf\{x\}\}\_\{t\}\|\\boldsymbol\{\\Sigma\}\_\{t\}\), the likelihood term∇𝐱tlog⁡p​\(𝐱T=𝐲\|𝐱t\)\\nabla\_\{\{\\mathbf\{x\}\}\_\{t\}\}\\log p\(\{\\mathbf\{x\}\}\_\{T\}=\{\\mathbf\{y\}\}\|\{\\mathbf\{x\}\}\_\{t\}\)acts as guidance to keep the iterates consistent with the observation𝐲\{\\mathbf\{y\}\}at covariance𝚺T\\boldsymbol\{\\Sigma\}\_\{T\}\. It has been shown that using better MCMC samplers than ULA can boost the performance of diffusion\-based samplers\(Duet al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib128)\)\. In this context, an advantage of our energy formulation is that it provides access to the densityp​\(𝐱t\|𝚺t\)p\(\{\\mathbf\{x\}\}\_\{t\}\|\\boldsymbol\{\\Sigma\}\_\{t\}\), in addition to the score, at every step\. This enables the use of*unbiased*correctors such as the*Metropolis\-Adjusted Langevin Algorithm*\(MALA\)\(Roberts and Tweedie,[1996](https://arxiv.org/html/2605.15487#bib.bib138)\), which treats \([15](https://arxiv.org/html/2605.15487#S4.E15)\) as a proposalq​\(𝐱t′\|𝐱\)q\(\{\\mathbf\{x\}\}\_\{t\}^\{\\prime\}\|\{\\mathbf\{x\}\}\)that is accepted with probabilitymin\(1,p​\(𝐱t′\|𝚺T,𝐲\)​q​\(𝐱t\|𝐱t′\)p​\(𝐱t\|𝚺T,𝐲\)​q​\(𝐱t′\|𝐱t\)\)\\min\\mathopen\{\}\\mathclose\{\{\\left\(1,\\frac\{p\(\{\\mathbf\{x\}\}^\{\\prime\}\_\{t\}\|\\boldsymbol\{\\Sigma\}\_\{T\},\{\\mathbf\{y\}\}\)\\,q\(\{\\mathbf\{x\}\}\_\{t\}\|\{\\mathbf\{x\}\}^\{\\prime\}\_\{t\}\)\}\{p\(\{\\mathbf\{x\}\}\_\{t\}\|\\boldsymbol\{\\Sigma\}\_\{T\},\{\\mathbf\{y\}\}\)\\,q\(\{\\mathbf\{x\}\}^\{\\prime\}\_\{t\}\|\{\\mathbf\{x\}\}\_\{t\}\)\}\}\}\\right\)\. We emphasize that the ability to compute this acceptance probability during*posterior*sampling is a distinctive feature of our anisotropic energy model, which explicitly provides access to the densityp​\(𝐲\|𝚺\)p\(\{\\mathbf\{y\}\}\|\\boldsymbol\{\\Sigma\}\), setting it apart from prior isotropic energy\-based diffusion approaches\(Duet al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib128); Thorntonet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib125)\)\.

We report in Table[2](https://arxiv.org/html/2605.15487#S4.T2)a comparison between ULA and MALA as a function of the number of corrector steps\. We observe that increasing the number of MALA steps consistently improves reconstruction quality as measured by LPIPS with the ground truth, whereas ULA does not yield further gains\. A visual comparison is provided in Fig\.[29](https://arxiv.org/html/2605.15487#A4.F29)\.

Table 2:LPIPS distance between ground truth𝐱\{\\mathbf\{x\}\}and reconstruction𝐱^\{\\hat\{\\mathbf\{x\}\}\}for different correction schemes on a CelebA inpainting task\.

### 4\.3Blind inverse problems

We show that our normalized energy model can be used to solve*blind inverse problems*, where the noise covariance𝚺\\boldsymbol\{\\Sigma\}\(or equivalently, the measurement operator𝐇\{\\mathbf\{H\}\}\) is unknown\. Our proposed strategy is straightforward: we estimate𝚺\\boldsymbol\{\\Sigma\}by maximizing the posterior probabilitylog⁡p​\(𝚺\|𝐲\)=log⁡p​\(𝐲\|𝚺\)\+log⁡p​\(𝚺\)\+cst\\log p\(\\boldsymbol\{\\Sigma\}\|\{\\mathbf\{y\}\}\)=\\log p\(\{\\mathbf\{y\}\}\|\\boldsymbol\{\\Sigma\}\)\+\\log p\(\\boldsymbol\{\\Sigma\}\)\+\\mathrm\{cst\}, given a prior on covariance matricesp​\(𝚺\)p\(\\boldsymbol\{\\Sigma\}\)\. Here, we choose a uniform prior over a set𝒮\\mathcal\{S\}of covariances, leading to

𝚺^=argmax𝚺∈𝒮log⁡pθ​\(𝐲\|𝚺\)\.\\hat\{\\boldsymbol\{\\Sigma\}\}=\\operatornamewithlimits\{argmax\}\_\{\\boldsymbol\{\\Sigma\}\\in\\mathcal\{S\}\}\\,\\log p\_\{\\theta\}\(\{\\mathbf\{y\}\}\|\\boldsymbol\{\\Sigma\}\)\.\(16\)We illustrate this procedure on an inpainting task where the box size is unknown and is estimated through \([16](https://arxiv.org/html/2605.15487#S4.E16)\)\. More precisely, we observe𝐲=𝐱\+𝚺s1/2​𝐯\{\\mathbf\{y\}\}=\{\\mathbf\{x\}\}\+\\boldsymbol\{\\Sigma\}\_\{s\}^\{1/2\}\{\\mathbf\{v\}\}, where𝚺s\\boldsymbol\{\\Sigma\}\_\{s\}is a block\-diagonal covariance with noise standard deviationσ1\\sigma\_\{1\}inside a centrals×ss\\times sbox andσ2=10−4\\sigma\_\{2\}=10^\{\-4\}outside of it\. Fig\.[7](https://arxiv.org/html/2605.15487#S4.F7)shows that our energy model estimates accurately both the box sizessand the noise levelσ1\\sigma\_\{1\}\. This success critically relies on the normalization of the energy values across covariances𝚺\\boldsymbol\{\\Sigma\}: an energy model trained without the covariance score matching objective fails at this task\.

![Refer to caption](https://arxiv.org/html/2605.15487v1/x9.png)
![Refer to caption](https://arxiv.org/html/2605.15487v1/x10.png)
![Refer to caption](https://arxiv.org/html/2605.15487v1/x11.png)

Figure 7:Blind reconstruction experiment\.Top:two observations corrupted by noise of standard deviationσ1\\sigma\_\{1\}inside a centrals×ss\\times ssquare, with\(σ1,s\)∈\{\(0\.1,14\),\(2,20\)\}\(\\sigma\_\{1\},s\)\\in\\\{\(0\.1,14\),\(2,20\)\\\}\)\.Middle and bottom:Log probabilitylog⁡pθ​\(𝐲\|𝚺s\)\\log p\_\{\\theta\}\(\{\\mathbf\{y\}\}\|\\boldsymbol\{\\Sigma\}\_\{s\}\)as a function of box sizessandσ1\\sigma\_\{1\}, respectively\. The noise covariance parameters estimated by the energy model trained with dual score matching \(red dashed vertical lines\) accurately track their ground truth values \(solid green vertical lines\)\.

## 5Discussion

We have introduced a framework for learning normalized conditional densities that can be used to solve linear inverse problems\. At the core of our method is Anisotropic Covariance Score Matching \(A\-CSM\), a regularization that enforces consistency across covariance matrices via the Fokker\-Planck equation, allowing accurate normalization of the learned density\. We validate our approach on multiple datasets and inverse problems, demonstrating that covariance\-conditioned energy\-based models share the strong performance of current diffusion solvers while enabling new capabilities that score\-based models do not have\. In particular, explicit access to normalized densities enables blind estimation, energy\-guided adaptive sampling, unbiased MCMC correction steps, and principled comparisons between posterior sampling algorithms\.

It is useful to position our method within the literature on linear inverse problems, building on the distinction between Bayesian and regression\-based methods introduced in Section[2\.2](https://arxiv.org/html/2605.15487#S2.SS2)\. Bayes’ rule,p​\(𝐱\|𝐲\)∝p​\(𝐱\)​p​\(𝐲\|𝐱\)p\(\{\\mathbf\{x\}\}\|\{\\mathbf\{y\}\}\)\\propto p\(\{\\mathbf\{x\}\}\)\\,p\(\{\\mathbf\{y\}\}\|\{\\mathbf\{x\}\}\), implies an appealing separation between a prior modelp​\(𝐱\)p\(\{\\mathbf\{x\}\}\)and a likelihood termp​\(𝐲\|𝐱\)p\(\{\\mathbf\{y\}\}\|\{\\mathbf\{x\}\}\): once a prior is learned, it can in principle be applied to any inverse problem with an explicit forward model without additional training\. This approach however faces a fundamental challenge: computing posterior means𝔼​\[𝐱\|𝐲\]\\mathbb\{E\}\[\{\\mathbf\{x\}\}\\,\|\\,\{\\mathbf\{y\}\}\]or drawing posterior samples is intractable in high dimensions\. For instance, for diffusion priors, posterior sampling requires the likelihoodp​\(𝐲\|𝐱t\)p\(\{\\mathbf\{y\}\}\|\{\\mathbf\{x\}\}\_\{t\}\)at every noise level, which involves an intractable marginalization over the clean signal\. The opposite approach—directly learning𝔼​\[𝐱\|𝐲\]\\mathbb\{E\}\[\{\\mathbf\{x\}\}\\,\|\\,\{\\mathbf\{y\}\}\], or more generally the conditional score∇log⁡p​\(𝐱t\|𝐲\)\\nabla\\log p\(\{\\mathbf\{x\}\}\_\{t\}\|\{\\mathbf\{y\}\}\), via regression—can be both efficient and accurate, but it breaks the separation between an explicit prior shared across degradations and a problem\-specific likelihood\. As a result, modifying the likelihood term requires in principle additional training, and the relationships between the learned conditional expectation operators across different degradations are no longer explicit\.

Within this context, our anisotropic energy model can be viewed as an intermediate approach that aims to reconcile these two viewpoints and combine their respective advantages of flexibility and tractability\. Because it yields a learned normalized posterior density—and therefore the corresponding score function—it can be regarded as a generalization of the standard score\-based methods discussed in Section[2\.2](https://arxiv.org/html/2605.15487#S2.SS2)\. It also amounts to a change of perspective: rather than forming the posteriorp​\(𝐱\|𝐲\)p\(\{\\mathbf\{x\}\}\|\{\\mathbf\{y\}\}\)from a priorp​\(𝐱\)p\(\{\\mathbf\{x\}\}\)and a likelihoodp​\(𝐲\|𝐱\)p\(\{\\mathbf\{y\}\}\|\{\\mathbf\{x\}\}\), we rely on the anisotropic diffusion \(Fokker–Planck\) equation∇𝚺p​\(𝐲\|𝚺\)=12​∇𝐲2p​\(𝐲\|𝚺\)\\nabla\_\{\\boldsymbol\{\\Sigma\}\}\\,p\(\{\\mathbf\{y\}\}\|\\boldsymbol\{\\Sigma\}\)=\\tfrac\{1\}\{2\}\\nabla^\{2\}\_\{\{\\mathbf\{y\}\}\}\\,p\(\{\\mathbf\{y\}\}\|\\boldsymbol\{\\Sigma\}\)to characterize the measurement densityp​\(𝐲\|𝚺\)p\(\{\\mathbf\{y\}\}\|\\boldsymbol\{\\Sigma\}\), whose score gives efficient access to the posterior mean\. The diffusion equation is enforced implicitly by the A\-CSM objective, which guarantees consistency between a common prior and the posterior means associated with each𝚺\\boldsymbol\{\\Sigma\}\. In this formulation, adding a novel degradation \(i\.e\., covariance matrix𝚺\\boldsymbol\{\\Sigma\}\) not seen during training requires additional learning to extend the range oflog⁡p​\(𝐲\|𝚺\)\\log p\(\{\\mathbf\{y\}\}\|\\boldsymbol\{\\Sigma\}\), but regularized by the A\-CSM objective, leveraging what has already been learned\.

Our method has several limitations\. First, training energy\-based models via \(dual\) score matching is more computationally intensive than training score\-based models\. It could be accelerated using sliced score matching\(Songet al\.,[2020](https://arxiv.org/html/2605.15487#bib.bib152)\), which enables replacing the extra backpropagation step with more memory\-efficient \(forward\-mode\) Jacobian vector products\. Second, the spatial and spectral diagonal covariance parameterization covers a variety of standard inverse problems\. Extending it to more general covariances would broaden the applications of our model and unlock the full potential of energy\-guided sampling \(at the cost of additional memory and computation\)\.

We envision several applications of energy\-based models that should be explored\. First, they can provide estimates of posterior entropy and mutual information between measurements and ground truth signals, which could be used to optimize*measurement design*, offering a more direct method thanZhanget al\.\([2025b](https://arxiv.org/html/2605.15487#bib.bib142)\)\. Second, the dependence of the energyUθ​\(𝐲,𝚺\)U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)on the noise covariance𝚺\\boldsymbol\{\\Sigma\}provides valuable information about the local geometry of the priorp​\(𝐱\)p\(\{\\mathbf\{x\}\}\), allowing for example estimation of local curvature and “tangent subspaces” in the vicinity of individual images\.

## Acknowledgments

The majority of this work was conducted during Nicolas Zilberstein’s summer internship at Center for Computational Neuroscience, Flatiron Institute, a division of the Simons Foundation\. The authors thank Jona Bruna, Florentin Coeurdoux, Pierre\-Etienne Fiquet, Zahra Kadkhodaie, and Guy Ohayon for useful discussions\. The authors thank the Scientific Computing Core, Flatiron Institute, for computing facilities and support\. This research was partially sponsored by the Army Research Office under Grant Number W911NF\-17\-S\-0002 and by the National Science Foundation under awards CCF\-2340481 and EF\-2126387\. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the Army Research Office, the U\.S\. Army, or the U\.S\. Government\. The U\.S\. Government is authorized to reproduce and distribute reprints for Government purposes, notwithstanding any copyright notation herein\.

## References

- X\. Ai, Y\. He, A\. Gu, R\. Salakhutdinov, J\. Z\. Kolter, N\. M\. Boffi, and M\. Simchowitz \(2026\)Joint distillation for fast likelihood evaluation and sampling in flow\-based models\.Intl\. Conf\. Learn\. Repr\. \(ICLR\)\.Cited by:[3rd item](https://arxiv.org/html/2605.15487#A3.I1.i2.I1.i3.p1.1),[§C\.3](https://arxiv.org/html/2605.15487#A3.SS3.p1.1)\.
- T\. Asthana, Y\. Bao, A\. Awad, and H\. Krim \(2025\)Accelerated image\-aware diffusion modeling\.InEuropean Signal Process\. Conf\. \(EUSIPCO\),pp\. 1862–1866\.Cited by:[§2\.1](https://arxiv.org/html/2605.15487#S2.SS1.SSS0.Px1.p1.1)\.
- M\. Balcerak, T\. Amiranashvili, A\. Terpin, S\. Shit, L\. Bogensperger, S\. Kaltenbach, P\. Koumoutsakos, and B\. Menze \(2025\)Energy matching: unifying flow matching and energy\-based models for generative modeling\.InAdvances in Neural Inf\. Process\. Syst\. \(NIPS\),Cited by:[2nd item](https://arxiv.org/html/2605.15487#A3.I1.i2.I1.i2.p1.1),[§C\.3](https://arxiv.org/html/2605.15487#A3.SS3.p1.1)\.
- A\. Bansal, E\. Borgnia, H\. Chu, J\. Li, H\. Kazemi, F\. Huang, M\. Goldblum, J\. Geiping, and T\. Goldstein \(2023\)Cold diffusion: inverting arbitrary image transforms without noise\.Advances in Neural Inf\. Process\. Syst\. \(NIPS\)36,pp\. 41259–41282\.Cited by:[§2\.1](https://arxiv.org/html/2605.15487#S2.SS1.SSS0.Px1.p1.1)\.
- G\. Bartosh, D\. Vetrov, and C\. A\. Naesseth \(2024\)Neural flow diffusion models: learnable forward process for improved diffusion modelling\.Advances in Neural Inf\. Process\. Syst\. \(NIPS\)37,pp\. 73952–73985\.Cited by:[§2\.1](https://arxiv.org/html/2605.15487#S2.SS1.SSS0.Px1.p1.1)\.
- A\. Bradley and P\. Nakkiran \(2025\)Classifier\-free guidance is a predictor\-corrector\.Trans\. Mach\. Learn\. Res\.\.Cited by:[§4\.2\.2](https://arxiv.org/html/2605.15487#S4.SS2.SSS2.p1.13)\.
- J\. Bruna and J\. Han \(2024\)Provable posterior sampling with denoising oracles via tilted transport\.Advances in Neural Inf\. Process\. Syst\. \(NIPS\)37,pp\. 82863–82894\.Cited by:[§1](https://arxiv.org/html/2605.15487#S1.p1.1)\.
- B\. Chen, D\. Martí Monsó, Y\. Du, M\. Simchowitz, R\. Tedrake, and V\. Sitzmann \(2024\)Diffusion forcing: next\-token prediction meets full\-sequence diffusion\.Advances in Neural Inf\. Process\. Syst\. \(NIPS\)37,pp\. 24081–24125\.Cited by:[§2\.1](https://arxiv.org/html/2605.15487#S2.SS1.SSS0.Px1.p1.1)\.
- Y\. Choi, Y\. Uh, J\. Yoo, and J\. Ha \(2020\)StarGAN v2: diverse image synthesis for multiple domains\.InProceedings of the IEEE/CVF Int\. Conf\. Comput\. Vis\. Pattern Recogn\. \(CVPR\),Cited by:[§3\.4](https://arxiv.org/html/2605.15487#S3.SS4.SSS0.Px1.p1.6),[§3\.4](https://arxiv.org/html/2605.15487#S3.SS4.SSS0.Px2.p1.6)\.
- H\. Chung, J\. Kim, M\. T\. Mccann, M\. L\. Klasky, and J\. C\. Ye \(2023\)Diffusion posterior sampling for general noisy inverse problems\.InIntl\. Conf\. Learn\. Repr\. \(ICLR\),Cited by:[§C\.4](https://arxiv.org/html/2605.15487#A3.SS4.SSS0.Px2.p1.1),[Table 7](https://arxiv.org/html/2605.15487#A4.T7.11.8.1.1),[Table 8](https://arxiv.org/html/2605.15487#A4.T8.8.8.1.1),[§1](https://arxiv.org/html/2605.15487#S1.p1.1),[§2\.2](https://arxiv.org/html/2605.15487#S2.SS2.SSS0.Px1.p1.2),[§3\.4](https://arxiv.org/html/2605.15487#S3.SS4.SSS0.Px1.p1.6),[Table 1](https://arxiv.org/html/2605.15487#S3.T1.18.18.20.1.1)\.
- G\. Daras, H\. Chung, C\. Lai, Y\. Mitsufuji, J\. C\. Ye, P\. Milanfar, A\. G\. Dimakis, and M\. Delbracio \(2024\)A survey on diffusion models for inverse problems\.arXiv preprint arXiv:2410\.00083\.Cited by:[§1](https://arxiv.org/html/2605.15487#S1.p1.1),[§3\.4](https://arxiv.org/html/2605.15487#S3.SS4.SSS0.Px1.p1.6)\.
- G\. Daras, M\. Delbracio, H\. Talebi, A\. Dimakis, and P\. Milanfar \(2023\)Soft diffusion: score matching with general corruptions\.Trans\. Mach\. Learn\. Res\.\.Cited by:[§2\.1](https://arxiv.org/html/2605.15487#S2.SS1.SSS0.Px1.p1.1)\.
- M\. Delbracio and P\. Milanfar \(2023\)Inversion by direct iteration: an alternative to denoising diffusion for image restoration\.Trans\. Mach\. Learn\. Res\.\.Cited by:[§1](https://arxiv.org/html/2605.15487#S1.p1.1),[§2\.2](https://arxiv.org/html/2605.15487#S2.SS2.SSS0.Px2.p1.4)\.
- K\. Ding, K\. Ma, S\. Wang, and E\. P\. Simoncelli \(2020\)Image quality assessment: unifying structure and texture similarity\.IEEE Trans\. Pattern Anal\. Mach\. Intell\.44\(5\),pp\. 2567–2581\.Cited by:[§3\.4](https://arxiv.org/html/2605.15487#S3.SS4.SSS0.Px1.p2.3)\.
- T\. Dockhorn, A\. Vahdat, and K\. Kreis \(2021\)Score\-based generative modeling with critically\-damped Langevin diffusion\.InIntl\. Conf\. Learn\. Repr\. \(ICLR\),Cited by:[§2\.1](https://arxiv.org/html/2605.15487#S2.SS1.SSS0.Px1.p1.1)\.
- Y\. Du, C\. Durkan, R\. Strudel, J\. B\. Tenenbaum, S\. Dieleman, R\. Fergus, J\. Sohl\-Dickstein, A\. Doucet, and W\. S\. Grathwohl \(2023\)Reduce, reuse, recycle: compositional generation with energy\-based diffusion models and MCMC\.InIntl\. Conf\. on Mach\. Learn\. \(ICML\),pp\. 8489–8510\.Cited by:[4th item](https://arxiv.org/html/2605.15487#A3.I1.i2.I1.i4.p1.1),[§C\.3](https://arxiv.org/html/2605.15487#A3.SS3.p1.1),[§1](https://arxiv.org/html/2605.15487#S1.p2.1),[§4\.2\.2](https://arxiv.org/html/2605.15487#S4.SS2.SSS2.p1.12)\.
- N\. Elata, H\. Chung, J\. C\. Ye, T\. Michaeli, and M\. Elad \(2025\)InvFussion: bridging supervised and zero\-shot diffusion for inverse problems\.Advances in Neural Inf\. Process\. Syst\. \(NIPS\)\.Cited by:[§1](https://arxiv.org/html/2605.15487#S1.p1.1),[§2\.2](https://arxiv.org/html/2605.15487#S2.SS2.SSS0.Px2.p1.4)\.
- S\. Elfwing, E\. Uchibe, and K\. Doya \(2018\)Sigmoid\-weighted linear units for neural network function approximation in reinforcement learning\.Neural Netw\.107,pp\. 3–11\.Cited by:[§3\.3](https://arxiv.org/html/2605.15487#S3.SS3.p3.11)\.
- B\. T\. Feng, J\. Smith, M\. Rubinstein, H\. Chang, K\. L\. Bouman, and W\. T\. Freeman \(2023\)Score\-based diffusion models as principled priors for inverse imaging\.InProceedings of the IEEE/CVF International Conference on Computer Vision,pp\. 10520–10531\.Cited by:[§2\.2](https://arxiv.org/html/2605.15487#S2.SS2.SSS0.Px1.p1.2)\.
- M\. Gerdes, M\. Welling, and M\. C\. Cheng \(2024\)GUD: generation with unified diffusion\.arXiv preprint arXiv:2410\.02667\.Cited by:[§2\.1](https://arxiv.org/html/2605.15487#S2.SS1.SSS0.Px1.p1.1),[§4\.2\.1](https://arxiv.org/html/2605.15487#S4.SS2.SSS1.p1.1)\.
- F\. Guth, Z\. Kadkhodaie, and E\. P\. Simoncelli \(2025\)Learning normalized image densities via dual score matching\.Advances in Neural Inf\. Process\. Syst\. \(NIPS\)\.Cited by:[item 1](https://arxiv.org/html/2605.15487#A3.I1.i1.p1.6),[§D\.3](https://arxiv.org/html/2605.15487#A4.SS3.SSS0.Px2.p1.1),[§D\.4](https://arxiv.org/html/2605.15487#A4.SS4.p1.1),[§D\.5](https://arxiv.org/html/2605.15487#A4.SS5.p1.14),[§1](https://arxiv.org/html/2605.15487#S1.p2.1),[§1](https://arxiv.org/html/2605.15487#S1.p3.1),[§3\.2](https://arxiv.org/html/2605.15487#S3.SS2.SSS0.Px1.p1.1),[§3\.2](https://arxiv.org/html/2605.15487#S3.SS2.SSS0.Px2.p1.5),[§3\.2](https://arxiv.org/html/2605.15487#S3.SS2.p1.4),[§3\.3](https://arxiv.org/html/2605.15487#S3.SS3.p1.2)\.
- K\. He, X\. Zhang, S\. Ren, and J\. Sun \(2015\)Deep residual learning for image recognition\. 2016 ieee conf\.InProceedings of the IEEE/CVF Int\. Conf\. Comput\. Vis\. Pattern Recogn\. \(CVPR\),pp\. 770–778\.Cited by:[§C\.1](https://arxiv.org/html/2605.15487#A3.SS1.SSS0.Px1.p1.6)\.
- M\. Heusel, H\. Ramsauer, T\. Unterthiner, B\. Nessler, and S\. Hochreiter \(2017\)GANs trained by a two time\-scale update rule converge to a local Nash equilibrium\.Advances in Neural Inf\. Process\. Syst\. \(NIPS\)30\.Cited by:[§3\.4](https://arxiv.org/html/2605.15487#S3.SS4.SSS0.Px1.p2.3)\.
- J\. Ho, A\. Jain, and P\. Abbeel \(2020\)Denoising diffusion probabilistic models\.Advances in Neural Inf\. Process\. Syst\. \(NIPS\)33,pp\. 6840–6851\.Cited by:[§1](https://arxiv.org/html/2605.15487#S1.p1.1),[§2\.1](https://arxiv.org/html/2605.15487#S2.SS1.p1.2)\.
- E\. Hoogeboom and T\. Salimans \(2023\)Blurring diffusion models\.InIntl\. Conf\. Learn\. Repr\. \(ICLR\),Cited by:[§2\.1](https://arxiv.org/html/2605.15487#S2.SS1.SSS0.Px1.p1.1)\.
- Y\. Hu, A\. Peng, W\. Gan, P\. Milanfar, M\. Delbracio, and U\. S\. Kamilov \(2025\)Stochastic deep restoration priors for imaging inverse problems\.InIntl\. Conf\. on Mach\. Learn\. \(ICML\),Cited by:[§2\.2](https://arxiv.org/html/2605.15487#S2.SS2.SSS0.Px2.p1.4)\.
- Z\. Kadkhodaie and E\. Simoncelli \(2021\)Stochastic solutions for linear inverse problems using the prior implicit in a denoiser\.Advances in Neural Inf\. Process\. Syst\. \(NIPS\)34,pp\. 13242–13254\.Cited by:[§1](https://arxiv.org/html/2605.15487#S1.p1.1),[§2\.1](https://arxiv.org/html/2605.15487#S2.SS1.p1.7),[§2\.2](https://arxiv.org/html/2605.15487#S2.SS2.SSS0.Px1.p1.2),[§4\.2\.1](https://arxiv.org/html/2605.15487#S4.SS2.SSS1.Px2.p1.8)\.
- T\. Karras, M\. Aittala, T\. Aila, and S\. Laine \(2022\)Elucidating the design space of diffusion\-based generative models\.Advances in Neural Inf\. Process\. Syst\. \(NIPS\)35,pp\. 26565–26577\.Cited by:[§C\.1](https://arxiv.org/html/2605.15487#A3.SS1.p1.1),[§3\.3](https://arxiv.org/html/2605.15487#S3.SS3.p1.4),[§3\.3](https://arxiv.org/html/2605.15487#S3.SS3.p3.11)\.
- T\. Karras, M\. Aittala, J\. Lehtinen, J\. Hellsten, T\. Aila, and S\. Laine \(2024\)Analyzing and improving the training dynamics of diffusion models\.InProceedings of the IEEE/CVF Int\. Conf\. Comput\. Vis\. Pattern Recogn\. \(CVPR\),pp\. 24174–24184\.Cited by:[§3\.3](https://arxiv.org/html/2605.15487#S3.SS3.p3.11)\.
- B\. Kawar, M\. Elad, S\. Ermon, and J\. Song \(2022\)Denoising diffusion restoration models\.Advances in Neural Inf\. Process\. Syst\. \(NIPS\)35,pp\. 23593–23606\.Cited by:[§1](https://arxiv.org/html/2605.15487#S1.p1.1),[§2\.2](https://arxiv.org/html/2605.15487#S2.SS2.SSS0.Px1.p1.2)\.
- Y\. LeCun, C\. Cortes, and C\. Burges \(2010\)MNIST handwritten digit database\.ATT Labs \[Online\]\. Available: http://yann\.lecun\.com/exdb/mnist2\.Cited by:[§3\.4](https://arxiv.org/html/2605.15487#S3.SS4.SSS0.Px1.p1.6)\.
- G\. Liu, A\. Vahdat, D\. Huang, E\. Theodorou, W\. Nie, and A\. Anandkumar \(2023\)I2I^\{2\}SB: image\-to\-image Schrödinger bridge\.InIntl\. Conf\. on Mach\. Learn\. \(ICML\),pp\. 22042–22062\.Cited by:[§1](https://arxiv.org/html/2605.15487#S1.p1.1),[§2\.2](https://arxiv.org/html/2605.15487#S2.SS2.SSS0.Px2.p1.4)\.
- Z\. Liu, P\. Luo, X\. Wang, and X\. Tang \(2015\)Deep learning face attributes in the wild\.InProceedings of the IEEE/CVF Int\. Conf\. Comput\. Vis\. \(ICCV\),Cited by:[§3\.4](https://arxiv.org/html/2605.15487#S3.SS4.SSS0.Px1.p1.6)\.
- M\. Mardani, J\. Song, J\. Kautz, and A\. Vahdat \(2024\)A variational perspective on solving inverse problems with diffusion models\.Intl\. Conf\. Learn\. Repr\. \(ICLR\)\.Cited by:[§C\.4](https://arxiv.org/html/2605.15487#A3.SS4.SSS0.Px1.p1.7),[Table 7](https://arxiv.org/html/2605.15487#A4.T7.11.9.2.1),[Table 8](https://arxiv.org/html/2605.15487#A4.T8.8.9.2.1),[§2\.2](https://arxiv.org/html/2605.15487#S2.SS2.SSS0.Px1.p1.2),[§3\.4](https://arxiv.org/html/2605.15487#S3.SS4.SSS0.Px1.p1.6),[§3\.4](https://arxiv.org/html/2605.15487#S3.SS4.SSS0.Px1.p2.3),[Table 1](https://arxiv.org/html/2605.15487#S3.T1.18.18.21.2.1)\.
- K\. Miyasawaet al\.\(1961\)An empirical Bayes estimator of the mean of a normal population\.Bull\. Inst\. Internat\. Statist38\(181\-188\),pp\. 1–2\.Cited by:[Appendix A](https://arxiv.org/html/2605.15487#A1.SS0.SSS0.Px1.p1.2),[§3\.1](https://arxiv.org/html/2605.15487#S3.SS1.SSS0.Px2.p2.2)\.
- H\. Negrel, F\. Coeurdoux, M\. S\. Albergo, and E\. Vanden\-Eijnden \(2025\)Multitask learning with stochastic interpolants\.InAdvances in Neural Inf\. Process\. Syst\. \(NIPS\),Cited by:[§1](https://arxiv.org/html/2605.15487#S1.p1.1),[§2\.2](https://arxiv.org/html/2605.15487#S2.SS2.SSS0.Px2.p1.4),[§4\.2\.1](https://arxiv.org/html/2605.15487#S4.SS2.SSS1.p1.1)\.
- G\. A\. Pavliotis \(2014\)Stochastic processes and applications: diffusion processes, the fokker\-planck and langevin equations\.Springer\.Cited by:[§3\.2](https://arxiv.org/html/2605.15487#S3.SS2.SSS0.Px2.p1.3)\.
- M\. Plainer, H\. Wu, L\. Klein, S\. Günnemann, and F\. Noé \(2025\)Consistent sampling and simulation: molecular dynamics with energy\-based diffusion models\.Advances in Neural Inf\. Process\. Syst\. \(NIPS\)\.Cited by:[§1](https://arxiv.org/html/2605.15487#S1.p2.1),[§3\.2](https://arxiv.org/html/2605.15487#S3.SS2.SSS0.Px2.p1.5)\.
- M\. Raphan and E\. P\. Simoncelli \(2011\)Least squares estimation without priors or supervision\.Neural Comp\.23\(2\),pp\. 374–420\.Cited by:[Appendix A](https://arxiv.org/html/2605.15487#A1.SS0.SSS0.Px1.p1.2),[§2\.1](https://arxiv.org/html/2605.15487#S2.SS1.p1.7),[§3\.1](https://arxiv.org/html/2605.15487#S3.SS1.SSS0.Px2.p2.2)\.
- H\. E\. Robbins \(1956\)An empirical Bayes approach to statistics\.InProc\. Third Berkeley Symp\. Math\. Statist\. Prob\.,pp\. 388–394\.Cited by:[§3\.1](https://arxiv.org/html/2605.15487#S3.SS1.SSS0.Px2.p2.2)\.
- G\. O\. Roberts and R\. L\. Tweedie \(1996\)Exponential convergence of Langevin distributions and their discrete approximations\.Bernoulli2,pp\. 341–363\.Cited by:[§4\.2\.2](https://arxiv.org/html/2605.15487#S4.SS2.SSS2.p1.12)\.
- Y\. Romano, M\. Elad, and P\. Milanfar \(2017\)The little engine that could: regularization by denoising \(RED\)\.SIAM J\. Imag\. Sciences10\(4\),pp\. 1804–1844\.Cited by:[§3\.3](https://arxiv.org/html/2605.15487#S3.SS3.p1.2)\.
- O\. Russakovsky, J\. Deng, H\. Su, J\. Krause, S\. Satheesh, S\. Ma, Z\. Huang, A\. Karpathy, A\. Khosla, M\. Bernstein,et al\.\(2015\)Imagenet large scale visual recognition challenge\.Int\. J\. Comp\. Vision115,pp\. 211–252\.Cited by:[§3\.4](https://arxiv.org/html/2605.15487#S3.SS4.SSS0.Px1.p1.6)\.
- C\. Saharia, W\. Chan, H\. Chang, C\. Lee, J\. Ho, T\. Salimans, D\. Fleet, and M\. Norouzi \(2022\)Palette: image\-to\-image diffusion models\.InACM SIGGRAPH,pp\. 1–10\.Cited by:[§C\.4](https://arxiv.org/html/2605.15487#A3.SS4.SSS0.Px4.p1.2),[§1](https://arxiv.org/html/2605.15487#S1.p1.1),[§3\.4](https://arxiv.org/html/2605.15487#S3.SS4.SSS0.Px1.p1.6),[Table 1](https://arxiv.org/html/2605.15487#S3.T1.18.18.23.4.1)\.
- S\. Sahoo, A\. Gokaslan, C\. M\. De Sa, and V\. Kuleshov \(2024\)Diffusion models with learned adaptive noise\.Advances in Neural Inf\. Process\. Syst\. \(NIPS\)37,pp\. 105730–105779\.Cited by:[§2\.1](https://arxiv.org/html/2605.15487#S2.SS1.SSS0.Px1.p1.1)\.
- R\. Singhal, M\. Goldstein, and R\. Ranganath \(2023\)Where to diffuse, how to diffuse, and how to get back: automated learning for multivariate diffusions\.InIntl\. Conf\. Learn\. Repr\. \(ICLR\),Cited by:[§2\.1](https://arxiv.org/html/2605.15487#S2.SS1.SSS0.Px1.p1.1)\.
- J\. Sohl\-Dickstein, E\. Weiss, N\. Maheswaranathan, and S\. Ganguli \(2015\)Deep unsupervised learning using nonequilibrium thermodynamics\.InIntl\. Conf\. on Mach\. Learn\. \(ICML\),Cited by:[§1](https://arxiv.org/html/2605.15487#S1.p1.1),[§2\.1](https://arxiv.org/html/2605.15487#S2.SS1.p1.2)\.
- J\. Song, A\. Vahdat, M\. Mardani, and J\. Kautz \(2022\)Pseudoinverse\-guided diffusion models for inverse problems\.InIntl\. Conf\. Learn\. Repr\. \(ICLR\),Cited by:[§2\.2](https://arxiv.org/html/2605.15487#S2.SS2.SSS0.Px1.p1.2)\.
- Y\. Song, C\. Durkan, I\. Murray, and S\. Ermon \(2021a\)Maximum likelihood training of score\-based diffusion models\.Advances in Neural Inf\. Process\. Syst\. \(NIPS\)34,pp\. 1415–1428\.Cited by:[§2\.1](https://arxiv.org/html/2605.15487#S2.SS1.p1.7),[§3\.2](https://arxiv.org/html/2605.15487#S3.SS2.SSS0.Px1.p1.1)\.
- Y\. Song and S\. Ermon \(2019\)Generative modeling by estimating gradients of the data distribution\.Advances in neural information processing systems32\.Cited by:[§1](https://arxiv.org/html/2605.15487#S1.p1.1)\.
- Y\. Song, S\. Garg, J\. Shi, and S\. Ermon \(2020\)Sliced score matching: a scalable approach to density and score estimation\.InUncertainty in artificial intelligence,pp\. 574–584\.Cited by:[§5](https://arxiv.org/html/2605.15487#S5.p4.1)\.
- Y\. Song, J\. Sohl\-Dickstein, D\. P\. Kingma, A\. Kumar, S\. Ermon, and B\. Poole \(2021b\)Score\-based generative modeling through stochastic differential equations\.InIntl\. Conf\. Learn\. Repr\. \(ICLR\),Cited by:[§C\.1](https://arxiv.org/html/2605.15487#A3.SS1.p1.1),[§1](https://arxiv.org/html/2605.15487#S1.p1.1),[§2\.1](https://arxiv.org/html/2605.15487#S2.SS1.p1.2),[§3\.3](https://arxiv.org/html/2605.15487#S3.SS3.p1.4),[§4\.2\.2](https://arxiv.org/html/2605.15487#S4.SS2.SSS2.p1.13)\.
- M\. Terris, S\. Hurault, M\. Song, and J\. Tachella \(2026\)Reconstruct anything model: a lightweight foundation model for computational imaging\.Intl\. Conf\. Learn\. Repr\. \(ICLR\)\.Cited by:[§2\.2](https://arxiv.org/html/2605.15487#S2.SS2.SSS0.Px2.p1.4)\.
- J\. Thornton, L\. Béthune, R\. Zhang, A\. Bradley, P\. Nakkiran, and S\. Zhai \(2025\)Composition and control with distilled energy diffusion models and sequential monte carlo\.InInt\. Conf\. on Artif\. Intell\. and Stat\.,Cited by:[1st item](https://arxiv.org/html/2605.15487#A3.I1.i2.I1.i1.p1.1),[2nd item](https://arxiv.org/html/2605.15487#A3.I1.i2.I1.i2.p1.1),[§C\.3](https://arxiv.org/html/2605.15487#A3.SS3.p1.1),[§1](https://arxiv.org/html/2605.15487#S1.p2.1),[§3\.3](https://arxiv.org/html/2605.15487#S3.SS3.p1.2),[§4\.2\.2](https://arxiv.org/html/2605.15487#S4.SS2.SSS2.p1.12)\.
- P\. Vincent \(2011\)A connection between score matching and denoising autoencoders\.Neural computation23\(7\),pp\. 1661–1674\.Cited by:[§2\.1](https://arxiv.org/html/2605.15487#S2.SS1.p1.7)\.
- L\. Wu, B\. Trippe, C\. Naesseth, D\. Blei, and J\. P\. Cunningham \(2023\)Practical and asymptotically exact conditional sampling in diffusion models\.Advances in Neural Inf\. Process\. Syst\. \(NIPS\)\.Cited by:[§2\.2](https://arxiv.org/html/2605.15487#S2.SS2.SSS0.Px1.p1.2)\.
- H\. Yu, A\. Klami, A\. Hyvarinen, A\. Korba, and O\. Chehab \(2025\)Density ratio estimation with conditional probability paths\.InIntl\. Conf\. on Mach\. Learn\. \(ICML\),Cited by:[§1](https://arxiv.org/html/2605.15487#S1.p2.1),[§3\.2](https://arxiv.org/html/2605.15487#S3.SS2.SSS0.Px2.p1.5)\.
- B\. Zhang, W\. Chu, J\. Berner, C\. Meng, A\. Anandkumar, and Y\. Song \(2025a\)Improving diffusion inverse problem solving with decoupled noise annealing\.InProceedings of the IEEE/CVF Int\. Conf\. Comput\. Vis\. Pattern Recogn\. \(CVPR\),pp\. 20895–20905\.Cited by:[§C\.4](https://arxiv.org/html/2605.15487#A3.SS4.SSS0.Px3.p1.5),[§3\.4](https://arxiv.org/html/2605.15487#S3.SS4.SSS0.Px1.p1.6),[Table 1](https://arxiv.org/html/2605.15487#S3.T1.18.18.22.3.1)\.
- L\. Zhang, Z\. Kadkhodaie, E\. P\. Simoncelli, and D\. H\. Brainard \(2025b\)Generalized compressed sensing for image reconstruction with diffusion probabilistic models\.Trans\. Mach\. Learn\. Res\.\.Cited by:[§5](https://arxiv.org/html/2605.15487#S5.p5.3)\.
- R\. Zhang, P\. Isola, A\. A\. Efros, E\. Shechtman, and O\. Wang \(2018\)The unreasonable effectiveness of deep features as a perceptual metric\.InProceedings of the IEEE/CVF Int\. Conf\. Comput\. Vis\. \(ICCV\),pp\. 586–595\.Cited by:[§3\.4](https://arxiv.org/html/2605.15487#S3.SS4.SSS0.Px1.p2.3)\.
- Y\. Zhu, K\. Zhang, J\. Liang, J\. Cao, B\. Wen, R\. Timofte, and L\. Van Gool \(2023\)Denoising diffusion models for plug\-and\-play image restoration\.InProceedings of the IEEE/CVF Int\. Conf\. Comput\. Vis\. Pattern Recogn\. \(CVPR\),pp\. 1219–1229\.Cited by:[§2\.2](https://arxiv.org/html/2605.15487#S2.SS2.SSS0.Px1.p1.2)\.
- N\. Zilberstein, M\. Mardani, and S\. Segarra \(2025\)Repulsive latent score distillation for solving inverse problems\.InIntl\. Conf\. Learn\. Repr\. \(ICLR\),Cited by:[§2\.2](https://arxiv.org/html/2605.15487#S2.SS2.SSS0.Px1.p1.2)\.

## Appendices

## Appendix ADual\-score matching for anisotropic noise

We expand on the dual\-score matching formulation for anisotropic noise introduced in Section[3\.2](https://arxiv.org/html/2605.15487#S3.SS2)\. Recall that we want to learn a covariance\-dependent energy function,Uθ​\(𝐲,𝚺\)=−log⁡p​\(𝐲\|𝚺\)U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)=\-\\log p\(\{\\mathbf\{y\}\}\|\\boldsymbol\{\\Sigma\}\), to approximate the negative log\-likelihood \(NLL\) of the noisy image distribution, for a family of covariance matrices𝚺\\boldsymbol\{\\Sigma\}\. We model this energy as follows

Uθ\(𝐲,𝚺\)=−logp\(𝐲\|𝚺\)=−log\(∫p\(𝐱\)p\(𝐲\|𝐱,𝚺\)d𝐱\)U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)=\-\\log p\(\{\\mathbf\{y\}\}\|\\boldsymbol\{\\Sigma\}\)=\-\\log\\mathopen\{\}\\mathclose\{\{\\left\(\\int p\(\{\\mathbf\{x\}\}\)p\(\{\\mathbf\{y\}\}\|\{\\mathbf\{x\}\},\\boldsymbol\{\\Sigma\}\)\\mathrm\{d\}\{\\mathbf\{x\}\}\}\}\\right\)\(17\)To learn this model, we use dual score matching, which requires computing the gradients w\.r\.t\.𝐲\{\\mathbf\{y\}\}and𝚺\\boldsymbol\{\\Sigma\}\.

##### Data score\.

We can use the Miyasawa\-Tweedie identity\(Miyasawa and others,[1961](https://arxiv.org/html/2605.15487#bib.bib140); Raphan and Simoncelli,[2011](https://arxiv.org/html/2605.15487#bib.bib141)\)given by

∇𝐲Uθ\(𝐲,𝚺\)=𝔼𝐱\[𝚺−1\(𝐲−𝐱\)\|𝐲\]\.\\nabla\_\{\{\\mathbf\{y\}\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)=\\mathbb\{E\}\_\{\{\\mathbf\{x\}\}\}\\mathopen\{\}\\mathclose\{\{\\left\[\\boldsymbol\{\\Sigma\}^\{\-1\}\\mathopen\{\}\\mathclose\{\{\\left\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\}\}\\right\)\|\{\\mathbf\{y\}\}\}\}\\right\]\.\(18\)We define the anisotropic denoising score matching loss as

ℓ~A−DSM=𝔼𝐱,𝐲\[∥∇𝐲Uθ\(𝐲,𝚺\)−𝚺−1\(𝐲−𝐱\)∥22\]\\tilde\{\\ell\}\_\{\\mathrm\{A\-DSM\}\}=\\mathbb\{E\}\_\{\{\\mathbf\{x\}\},\{\\mathbf\{y\}\}\}\\mathopen\{\}\\mathclose\{\{\\left\[\\lVert\\nabla\_\{\{\\mathbf\{y\}\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\-\\boldsymbol\{\\Sigma\}^\{\-1\}\\mathopen\{\}\\mathclose\{\{\\left\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\}\}\\right\)\\rVert\_\{2\}^\{2\}\}\}\\right\]\(19\)To have a scale\-invariant loss, we reweight by𝚺12\\boldsymbol\{\\Sigma\}^\{\\frac\{1\}\{2\}\}:

ℓA−DSM=𝔼\[∥𝚺12∇𝐲Uθ\(𝐲,𝚺\)−𝚺−12\(𝐲−𝐱\)∥2\]\.\\ell\_\{\\mathrm\{A\-DSM\}\}=\\mathbb\{E\}\\mathopen\{\}\\mathclose\{\{\\left\[\\lVert\\boldsymbol\{\\Sigma\}^\{\\frac\{1\}\{2\}\}\\nabla\_\{\{\\mathbf\{y\}\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\-\\boldsymbol\{\\Sigma\}^\{\-\\frac\{1\}\{2\}\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)\\rVert^\{2\}\}\}\\right\]\.\(20\)

##### Covariance score\.

To obtain this second score, we differentiate the energy in \([17](https://arxiv.org/html/2605.15487#A1.E17)\) w\.r\.t\.𝚺\\boldsymbol\{\\Sigma\}\. Given thatp​\(𝐲\|𝐱,𝚺\)=e−12​\(𝐲−𝐱\)⊤​𝚺−1​\(𝐲−𝐱\)−d2​log⁡\(2​π\)−12​log​det\(𝚺\)p\(\{\\mathbf\{y\}\}\|\{\\mathbf\{x\}\},\\boldsymbol\{\\Sigma\}\)=e^\{\-\\frac\{1\}\{2\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)^\{\\top\}\\boldsymbol\{\\Sigma\}^\{\-1\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)\-\\frac\{d\}\{2\}\\log\(2\\pi\)\-\\frac\{1\}\{2\}\\log\\det\(\\boldsymbol\{\\Sigma\}\)\}, we obtain

∇𝚺Uθ​\(𝐲,𝚺\)\\displaystyle\\nabla\_\{\\boldsymbol\{\\Sigma\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)=𝔼𝐱\[−∇𝚺logp\(𝐲\|𝐱,𝚺\)\|𝐲,𝚺\]\\displaystyle=\\mathbb\{E\}\_\{\{\\mathbf\{x\}\}\}\\mathopen\{\}\\mathclose\{\{\\left\[\-\\nabla\_\{\\boldsymbol\{\\Sigma\}\}\\log p\(\{\\mathbf\{y\}\}\|\{\\mathbf\{x\}\},\\boldsymbol\{\\Sigma\}\)\|\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\}\}\\right\]\(21\)=𝔼𝐱\[−∇𝚺\(−12\(𝐲−𝐱\)⊤𝚺−1\(𝐲−𝐱\)−d2log\(2π\)−12logdet\(𝚺\)\)\|𝐲,𝚺\]\.\\displaystyle=\\mathbb\{E\}\_\{\{\\mathbf\{x\}\}\}\\mathopen\{\}\\mathclose\{\{\\left\[\-\\nabla\_\{\\boldsymbol\{\\Sigma\}\}\\mathopen\{\}\\mathclose\{\{\\left\(\-\\frac\{1\}\{2\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)^\{\\top\}\\boldsymbol\{\\Sigma\}^\{\-1\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)\-\\frac\{d\}\{2\}\\log\(2\\pi\)\-\\frac\{1\}\{2\}\\log\\det\(\\boldsymbol\{\\Sigma\}\)\}\}\\right\)\\bigg\|\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\}\}\\right\]\.The two gradient terms can be simplified as follows:

∇𝚺\[−12\(𝐲−𝐱\)⊤𝚺−1\(𝐲−𝐱\)\]\\displaystyle\\nabla\_\{\\boldsymbol\{\\Sigma\}\}\\mathopen\{\}\\mathclose\{\{\\left\[\-\\frac\{1\}\{2\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)^\{\\top\}\\boldsymbol\{\\Sigma\}^\{\-1\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)\}\}\\right\]=12​𝚺−1​\(𝐲−𝐱\)​\(𝐲−𝐱\)⊤​𝚺−1\\displaystyle=\\frac\{1\}\{2\}\\boldsymbol\{\\Sigma\}^\{\-1\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)^\{\\top\}\\boldsymbol\{\\Sigma\}^\{\-1\}\(22\)∇𝚺\[−12logdet𝚺\]\\displaystyle\\nabla\_\{\\boldsymbol\{\\Sigma\}\}\\mathopen\{\}\\mathclose\{\{\\left\[\-\\frac\{1\}\{2\}\\log\\det\\boldsymbol\{\\Sigma\}\}\}\\right\]=−12​𝚺−1\.\\displaystyle=\-\\frac\{1\}\{2\}\\boldsymbol\{\\Sigma\}^\{\-1\}\.\(23\)Substituting \([22](https://arxiv.org/html/2605.15487#A1.E22)\) and \([23](https://arxiv.org/html/2605.15487#A1.E23)\) into \([21](https://arxiv.org/html/2605.15487#A1.E21)\) gives

∇𝚺Uθ\(𝐲,𝚺\)=𝔼𝐱\[12𝚺−1−12𝚺−1\(𝐲−𝐱\)\(𝐲−𝐱\)⊤𝚺−1\|𝐲,𝚺\]\.\\nabla\_\{\\boldsymbol\{\\Sigma\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)=\\mathbb\{E\}\_\{\{\\mathbf\{x\}\}\}\\mathopen\{\}\\mathclose\{\{\\left\[\\frac\{1\}\{2\}\\boldsymbol\{\\Sigma\}^\{\-1\}\-\\frac\{1\}\{2\}\\boldsymbol\{\\Sigma\}^\{\-1\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)^\{\\top\}\\boldsymbol\{\\Sigma\}^\{\-1\}\\bigg\|\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\}\}\\right\]\.\(24\)Finally, we define the anisotropic covariance score matching loss as

ℓ~A−CSM=𝔼𝐱,𝐲\[∥∇𝚺Uθ\(𝐲,𝚺\)−\(12𝚺−1−12𝚺−1\(𝐲−𝐱\)\(𝐲−𝐱\)⊤𝚺−1\)∥F2\],\\tilde\{\\ell\}\_\{\\mathrm\{A\-CSM\}\}=\\mathbb\{E\}\_\{\{\\mathbf\{x\}\},\{\\mathbf\{y\}\}\}\\mathopen\{\}\\mathclose\{\{\\left\[\\mathopen\{\}\\mathclose\{\{\\left\\lVert\\nabla\_\{\\boldsymbol\{\\Sigma\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\-\\mathopen\{\}\\mathclose\{\{\\left\(\\frac\{1\}\{2\}\\boldsymbol\{\\Sigma\}^\{\-1\}\-\\frac\{1\}\{2\}\\boldsymbol\{\\Sigma\}^\{\-1\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)^\{\\top\}\\boldsymbol\{\\Sigma\}^\{\-1\}\}\}\\right\)\}\}\\right\\rVert\_\{F\}^\{2\}\}\}\\right\],\(25\)Similarly to the data score, we weight the loss by𝚺12\\boldsymbol\{\\Sigma\}^\{\\frac\{1\}\{2\}\}on both sides to be scale\-invariant:

ℓA−CSM=𝔼𝐱,𝐲,𝚺​\[∥𝚺12​∇𝚺Uθ​\(𝐲,𝚺\)​𝚺12−12​𝐈\+12​𝚺−12​\(𝐲−𝐱\)​\(𝐲−𝐱\)⊤​𝚺−12∥22\]\.\\ell\_\{\\mathrm\{A\-CSM\}\}=\\mathbb\{E\}\_\{\{\\mathbf\{x\}\},\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\}\\Big\[\\big\\lVert\\boldsymbol\{\\Sigma\}^\{\\frac\{1\}\{2\}\}\\nabla\_\{\\boldsymbol\{\\Sigma\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\\boldsymbol\{\\Sigma\}^\{\\frac\{1\}\{2\}\}\-\\frac\{1\}\{2\}\{\\mathbf\{I\}\}\+\\frac\{1\}\{2\}\\boldsymbol\{\\Sigma\}^\{\-\\frac\{1\}\{2\}\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)^\{\\top\}\\boldsymbol\{\\Sigma\}^\{\-\\frac\{1\}\{2\}\}\\big\\rVert\_\{2\}^\{2\}\\Big\]\.\(26\)

##### Final loss\.

Now that we have the dual\-score matching loss, we define an overall objective by integrating over all possible covariance matrices𝚺\\boldsymbol\{\\Sigma\}\. Directly extending the combined loss from the paper yields

ℓ\(θ,𝚺\)=𝔼𝚺\[k1ℓA−DSM\(θ,𝚺\)\+k2ℓA−CSM\(θ,𝚺\)\]\\ell\(\\theta,\\boldsymbol\{\\Sigma\}\)=\\mathbb\{E\}\_\{\\boldsymbol\{\\Sigma\}\}\\mathopen\{\}\\mathclose\{\{\\left\[k\_\{1\}\\ell\_\{\\mathrm\{A\-DSM\}\}\(\\theta,\\boldsymbol\{\\Sigma\}\)\+k\_\{2\}\\ell\_\{\\mathrm\{A\-CSM\}\}\(\\theta,\\boldsymbol\{\\Sigma\}\)\}\}\\right\]\(27\)wherek1,k2k\_\{1\},k\_\{2\}are constants that normalize each term with respect to the ambient dimensiondd\. In particular, we considerk1=1/dk\_\{1\}=1/dandk2=1/d2k\_\{2\}=1/d^\{2\}:

ℓ\(θ,𝚺\)=𝔼𝚺\[1dℓA−DSM\(θ,𝚺\)\+1d2ℓA−CSM\(θ,𝚺\)\]\\ell\(\\theta,\\boldsymbol\{\\Sigma\}\)=\\mathbb\{E\}\_\{\\boldsymbol\{\\Sigma\}\}\\mathopen\{\}\\mathclose\{\{\\left\[\\frac\{1\}\{d\}\\ell\_\{\\mathrm\{A\-DSM\}\}\(\\theta,\\boldsymbol\{\\Sigma\}\)\+\\frac\{1\}\{d^\{2\}\}\\ell\_\{\\mathrm\{A\-CSM\}\}\(\\theta,\\boldsymbol\{\\Sigma\}\)\}\}\\right\]\(28\)While this loss is appealing, it is impractical due to the size of the covariance\. As mentioned in Section[3\.2](https://arxiv.org/html/2605.15487#S3.SS2), we consider a diagonal matrix instead\.

##### Diagonal𝚺\\boldsymbol\{\\Sigma\}with independent parameters\.

We consider diagonal matrices where𝚺=diag⁡\(𝚽\)\\boldsymbol\{\\Sigma\}=\\operatorname\{diag\}\(\\boldsymbol\{\\Phi\}\)with𝚽=\[Φ1,…,Φd\]⊤\\boldsymbol\{\\Phi\}=\[\\Phi\_\{1\},\\dots,\\Phi\_\{d\}\]^\{\\top\}\. Each parameter is independent with priorp​\(𝚽\)=∏i=1dp​\(Φi\)p\(\\boldsymbol\{\\Phi\}\)=\\prod\_\{i=1\}^\{d\}p\(\\Phi\_\{i\}\), wherep​\(Φi\)∝Φi−1p\(\\Phi\_\{i\}\)\\propto\\Phi\_\{i\}^\{\-1\}forΦi∈\[Φmin,Φmax\]\\Phi\_\{i\}\\in\[\\Phi\_\{\\min\},\\Phi\_\{\\max\}\]\.

To compute the gradient with respect to𝚽\\boldsymbol\{\\Phi\}, we utilize the chain rule\. Note that∂vec⁡\(𝚺\)∂Φi\\frac\{\\partial\\operatorname\{vec\}\(\\boldsymbol\{\\Sigma\}\)\}\{\\partial\\Phi\_\{i\}\}is a vector of lengthd2d^\{2\}with a11at the position corresponding to theii\-th diagonal elementΣi​i\\Sigma\_\{ii\}, and0elsewhere\. For a generic functionUθ​\(𝐲,𝚺\)U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\), the gradient componentiiis:

\[∇𝚽Uθ\(𝐲,𝚺\(𝚽\)\)\]i\\displaystyle\\mathopen\{\}\\mathclose\{\{\\left\[\\nabla\_\{\\boldsymbol\{\\Phi\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\(\\boldsymbol\{\\Phi\}\)\)\}\}\\right\]\_\{i\}=vec\(∇𝚺Uθ\(𝐲,𝚺\)\)⊤⋅∂vec⁡\(𝚺\)∂Φi\\displaystyle=\\operatorname\{vec\}\\mathopen\{\}\\mathclose\{\{\\left\(\\nabla\_\{\\boldsymbol\{\\Sigma\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\}\}\\right\)^\{\\top\}\\cdot\\frac\{\\partial\\operatorname\{vec\}\(\\boldsymbol\{\\Sigma\}\)\}\{\\partial\\Phi\_\{i\}\}\(29\)=\[∇𝚺Uθ\(𝐲,𝚺\)\]i​i\\displaystyle=\\mathopen\{\}\\mathclose\{\{\\left\[\\nabla\_\{\\boldsymbol\{\\Sigma\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\}\}\\right\]\_\{ii\}\(30\)
Specifically, for the log\-determinant term in \([21](https://arxiv.org/html/2605.15487#A1.E21)\), the gradient is:

\[∇𝚽\(−12logdet𝚺\)\]i\\displaystyle\\mathopen\{\}\\mathclose\{\{\\left\[\\nabla\_\{\\boldsymbol\{\\Phi\}\}\\mathopen\{\}\\mathclose\{\{\\left\(\-\\frac\{1\}\{2\}\\log\\det\\boldsymbol\{\\Sigma\}\}\}\\right\)\}\}\\right\]\_\{i\}=vec\(−12𝚺−1\)⊤⋅∂vec⁡\(𝚺\)∂Φi\\displaystyle=\\operatorname\{vec\}\\mathopen\{\}\\mathclose\{\{\\left\(\-\\frac\{1\}\{2\}\\boldsymbol\{\\Sigma\}^\{\-1\}\}\}\\right\)^\{\\top\}\\cdot\\frac\{\\partial\\operatorname\{vec\}\(\\boldsymbol\{\\Sigma\}\)\}\{\\partial\\Phi\_\{i\}\}\(31\)=−12​\(𝚺−1\)i​i=−12​Φi\\displaystyle=\-\\frac\{1\}\{2\}\(\\boldsymbol\{\\Sigma\}^\{\-1\}\)\_\{ii\}=\-\\frac\{1\}\{2\\Phi\_\{i\}\}\(32\)Thus, in vector form, we obtain the simplified result:

∇𝚽\[−12logdet𝚺\]=−12𝚽−1\\nabla\_\{\\boldsymbol\{\\Phi\}\}\\mathopen\{\}\\mathclose\{\{\\left\[\-\\frac\{1\}\{2\}\\log\\det\\boldsymbol\{\\Sigma\}\}\}\\right\]=\-\\frac\{1\}\{2\}\\boldsymbol\{\\Phi\}^\{\-1\}\(33\)For the quadratic term, we have the following \(assuming𝐳=\(𝐲−𝐱\)\{\\mathbf\{z\}\}=\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)\)

∇𝚽\[−12\(𝐲−𝐱\)⊤𝚺−1\(𝐲−𝐱\)\]\\displaystyle\\nabla\_\{\\boldsymbol\{\\Phi\}\}\\mathopen\{\}\\mathclose\{\{\\left\[\-\\frac\{1\}\{2\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)^\{\\top\}\\boldsymbol\{\\Sigma\}^\{\-1\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)\}\}\\right\]=vec\(12𝚺−1𝐳𝐳⊤𝚺−1\)⋅∂vec⁡\(𝚺\)∂𝚽,\\displaystyle=\\operatorname\{vec\}\\mathopen\{\}\\mathclose\{\{\\left\(\\frac\{1\}\{2\}\\boldsymbol\{\\Sigma\}^\{\-1\}\{\\mathbf\{z\}\}\{\\mathbf\{z\}\}^\{\\top\}\\boldsymbol\{\\Sigma\}^\{\-1\}\}\}\\right\)\\cdot\\frac\{\\partial\\operatorname\{vec\}\(\\boldsymbol\{\\Sigma\}\)\}\{\\partial\\boldsymbol\{\\Phi\}\},\(34\)=12vec\(\(𝚽−1∘𝐳\)\(𝚽−1∘𝐳\)⊤\)⋅∂vec⁡\(𝚺\)∂𝚽,\\displaystyle=\\frac\{1\}\{2\}\\operatorname\{vec\}\\mathopen\{\}\\mathclose\{\{\\left\(\(\\boldsymbol\{\\Phi\}^\{\-1\}\\circ\{\\mathbf\{z\}\}\)\(\\boldsymbol\{\\Phi\}^\{\-1\}\\circ\{\\mathbf\{z\}\}\)^\{\\top\}\}\}\\right\)\\cdot\\frac\{\\partial\\operatorname\{vec\}\(\\boldsymbol\{\\Sigma\}\)\}\{\\partial\\boldsymbol\{\\Phi\}\},where∘\\circis the element\-wise product\. Theii\-th element is

\[∇𝚽\[−12\(𝐲−𝐱\)⊤𝚺−1\(𝐲−𝐱\)\]\]i=\[12vec\(\(𝚽−1∘𝐳\)\(𝚽−1∘𝐳\)⊤\)\]i​i\.\\mathopen\{\}\\mathclose\{\{\\left\[\\nabla\_\{\\boldsymbol\{\\Phi\}\}\\mathopen\{\}\\mathclose\{\{\\left\[\-\\frac\{1\}\{2\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)^\{\\top\}\\boldsymbol\{\\Sigma\}^\{\-1\}\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)\}\}\\right\]\}\}\\right\]\_\{i\}=\\mathopen\{\}\\mathclose\{\{\\left\[\\frac\{1\}\{2\}\\operatorname\{vec\}\\mathopen\{\}\\mathclose\{\{\\left\(\(\\boldsymbol\{\\Phi\}^\{\-1\}\\circ\{\\mathbf\{z\}\}\)\(\\boldsymbol\{\\Phi\}^\{\-1\}\\circ\{\\mathbf\{z\}\}\)^\{\\top\}\}\}\\right\)\}\}\\right\]\_\{ii\}\.\(35\)Combining both \([35](https://arxiv.org/html/2605.15487#A1.E35)\) and \([33](https://arxiv.org/html/2605.15487#A1.E33)\), we get

∇𝚽Uθ\(𝐲,𝚺\(𝚽\)\)=𝔼𝐱\[12𝚽−1−12\(𝚽−1∘𝐳\)2\|𝐲,𝚺\(𝚽\)\]\.\\nabla\_\{\\boldsymbol\{\\Phi\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\(\\boldsymbol\{\\Phi\}\)\)=\\mathbb\{E\}\_\{\{\\mathbf\{x\}\}\}\\mathopen\{\}\\mathclose\{\{\\left\[\\frac\{1\}\{2\}\\boldsymbol\{\\Phi\}^\{\-1\}\-\\frac\{1\}\{2\}\(\\boldsymbol\{\\Phi\}^\{\-1\}\\circ\{\\mathbf\{z\}\}\)^\{2\}\\middle\|\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\(\\boldsymbol\{\\Phi\}\)\}\}\\right\]\.\(36\)The training algorithm is described in Alg\.[1](https://arxiv.org/html/2605.15487#alg1)\. Notice we are usingλ=d−12\\lambda=d^\{\-\\frac\{1\}\{2\}\}in \([28](https://arxiv.org/html/2605.15487#A1.E28)\), which corresponds to the loss in the main paper\.

### A\.1Algorithms

Algorithm 1Trainingrepeat

𝐱∼p​\(𝐱\)\{\\mathbf\{x\}\}\\sim p\(\{\\mathbf\{x\}\}\)\(This is the dataset given beforehand\)

𝚽∼p​\(𝚽\)\\boldsymbol\{\\Phi\}\\sim p\(\\boldsymbol\{\\Phi\}\)

𝐯∼𝒩​\(0,𝐈\)\{\\mathbf\{v\}\}\\sim\{\\mathcal\{N\}\}\(0,\{\\mathbf\{I\}\}\)

𝐲=𝐱\+𝚺​\(𝚽\)1/2​𝐯\{\\mathbf\{y\}\}=\{\\mathbf\{x\}\}\+\\boldsymbol\{\\Sigma\}\(\\boldsymbol\{\\Phi\}\)^\{1/2\}\{\\mathbf\{v\}\}

Compute

ℓA−DSM\(θ,𝚺\(𝚽\)\)=1N∑𝐱,𝐲∥∇𝐲Uθ\(𝐲,𝚺\(𝚽\)\)−𝚽−1∘\(𝐲−𝐱\)∥𝚽d2\\ell\_\{\\mathrm\{A\-DSM\}\}\(\\theta,\\boldsymbol\{\\Sigma\}\(\\boldsymbol\{\\Phi\}\)\)=\\frac\{1\}\{N\}\\sum\_\{\{\\mathbf\{x\}\},\{\\mathbf\{y\}\}\}\\lVert\\nabla\_\{\{\\mathbf\{y\}\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\(\\boldsymbol\{\\Phi\}\)\)\-\\boldsymbol\{\\Phi\}^\{\-1\}\\circ\\mathopen\{\}\\mathclose\{\{\\left\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\}\}\\right\)\\rVert\_\{\\frac\{\\boldsymbol\{\\Phi\}\}\{d\}\}^\{2\}
Compute

ℓA−CSM\(θ,𝚺\(𝚽\)\)=1N∑𝐱,𝐲∥∇𝚽Uθ\(𝐲,𝚺\(𝚽\)\)−\(12𝚽−1−12\(𝚽−1∘\(𝐲−𝐱\)\)2\)∥𝚽2d2\\ell\_\{\\mathrm\{A\-CSM\}\}\(\\theta,\\boldsymbol\{\\Sigma\}\(\\boldsymbol\{\\Phi\}\)\)=\\frac\{1\}\{N\}\\sum\_\{\{\\mathbf\{x\}\},\{\\mathbf\{y\}\}\}\\mathopen\{\}\\mathclose\{\{\\left\\lVert\\nabla\_\{\\boldsymbol\{\\Phi\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\(\\boldsymbol\{\\Phi\}\)\)\-\\mathopen\{\}\\mathclose\{\{\\left\(\\frac\{1\}\{2\}\\boldsymbol\{\\Phi\}^\{\-1\}\-\\frac\{1\}\{2\}\(\\boldsymbol\{\\Phi\}^\{\-1\}\\circ\(\{\\mathbf\{y\}\}\-\{\\mathbf\{x\}\}\)\)^\{2\}\}\}\\right\)\}\}\\right\\rVert\_\{\\frac\{\\boldsymbol\{\\Phi\}^\{2\}\}\{d^\{2\}\}\}
Combine

ℓ​\(θ\)=ℓA−DSM​\(θ,𝚺​\(𝚽\)\)\+ℓA−CSM​\(θ,𝚺​\(𝚽\)\)\{\\ell\}\(\\theta\)=\\ell\_\{\\mathrm\{A\-DSM\}\}\(\\theta,\\boldsymbol\{\\Sigma\}\(\\boldsymbol\{\\Phi\}\)\)\+\\ell\_\{\\mathrm\{A\-CSM\}\}\(\\theta,\\boldsymbol\{\\Sigma\}\(\\boldsymbol\{\\Phi\}\)\)
Perform one gradient descent step on

ℓ​\(θ\)\\ell\(\\theta\)w\.r\.t\.

θ\\theta
untilConvergence

Algorithm 2Predictor–Corrector Sampling0:Noise schedule

\{𝚺0,…,𝚺T\}\\\{\\boldsymbol\{\\Sigma\}\_\{0\},\\dots,\\boldsymbol\{\\Sigma\}\_\{T\}\\\}, number of corrector steps

nstepsn\_\{\\mathrm\{steps\}\}, energy model

UθU\_\{\\theta\}, measurement

𝐲\{\\mathbf\{y\}\}, binary matrix with the structure of the degradation

𝚺i​d=clamp​\(𝐇−1,0,1\)\\boldsymbol\{\\Sigma\}\_\{id\}=\\textrm\{clamp\}\(\{\\mathbf\{H\}\}^\{\-1\},0,1\),

r\\mathrm\{r\}, temperature

η\\eta
1:Initialize

𝐱T←𝐲\{\\mathbf\{x\}\}\_\{T\}\\leftarrow\{\\mathbf\{y\}\}
2:for

c=T−1c=T\-1to

0do

3:// PREDICTOR STEP\.

4:Compute score gradient

𝐠c=−∇𝐱cUθ​\(𝐱c,𝚺c\)\{\\mathbf\{g\}\}\_\{c\}=\-\\nabla\_\{\{\\mathbf\{x\}\}\_\{c\}\}U\_\{\\theta\}\(\{\\mathbf\{x\}\}\_\{c\},\\boldsymbol\{\\Sigma\}\_\{c\}\)
5:ifSchedule = Adaptivethen

6:Compute covariance gradients

𝐠Σc=−∇𝚺Uθ​\(𝐱c,0,𝚺c\)\{\\mathbf\{g\}\}\_\{\\Sigma\_\{c\}\}=\-\\nabla\_\{\\boldsymbol\{\\Sigma\}\}U\_\{\\theta\}\(\{\\mathbf\{x\}\}\_\{c,0\},\\boldsymbol\{\\Sigma\}\_\{c\}\)
7:Covariance gradient descent \+ PSD constraint:

𝚺c\+1←max⁡\(𝚺c−ηΣ​𝚺∘𝐠Σ∘𝚺,ε​𝐈\)\\boldsymbol\{\\Sigma\}\_\{c\+1\}\\leftarrow\\max\(\\boldsymbol\{\\Sigma\}\_\{c\}\-\\eta\_\{\\Sigma\}\\,\\boldsymbol\{\\Sigma\}\\circ\{\\mathbf\{g\}\}\_\{\\Sigma\}\\circ\\boldsymbol\{\\Sigma\},\\varepsilon\{\\mathbf\{I\}\}\)
8:Define

Δ​𝚺=𝚺c∘𝐠Σc∘𝚺c\\Delta\\boldsymbol\{\\Sigma\}=\\boldsymbol\{\\Sigma\}\_\{c\}\\circ\{\\mathbf\{g\}\}\_\{\\Sigma\_\{c\}\}\\circ\\boldsymbol\{\\Sigma\}\_\{c\}
9:elseifSchedule = Fixedthen

10:

Δ​𝚺=𝚺c−𝚺c\+1\\Delta\\boldsymbol\{\\Sigma\}=\\boldsymbol\{\\Sigma\}\_\{c\}\-\\boldsymbol\{\\Sigma\}\_\{c\+1\}
11:endif

12:Sample noise

𝐳∼𝒩​\(0,𝐈\)\{\\mathbf\{z\}\}\\sim\\mathcal\{N\}\(0,\{\\mathbf\{I\}\}\)𝐱c\+1,0←𝐱c\+Δ​𝚺∘𝐠\+Δ​𝚺1/2∘𝐳\{\\mathbf\{x\}\}\_\{c\+1,0\}\\leftarrow\{\\mathbf\{x\}\}\_\{c\}\+\\Delta\\boldsymbol\{\\Sigma\}\\circ\{\\mathbf\{g\}\}\+\\Delta\\boldsymbol\{\\Sigma\}^\{1/2\}\\circ\{\\mathbf\{z\}\}
13:// CORRECTOR STEP\.

14:for

s=0s=0to

nstepsn\_\{\\mathrm\{steps\}\}do

15:Compute score gradient

𝐠c\+1,s=−∇𝐱c\+1,sUθ​\(𝐱c\+1,s,𝚺c\+1\)\{\\mathbf\{g\}\}\_\{c\+1,s\}=\-\\nabla\_\{\{\\mathbf\{x\}\}\_\{c\+1,s\}\}U\_\{\\theta\}\(\{\\mathbf\{x\}\}\_\{c\+1,s\},\\boldsymbol\{\\Sigma\}\_\{c\+1\}\)
16:Sample noise

𝐳∼𝒩​\(0,𝐈\)\{\\mathbf\{z\}\}\\sim\\mathcal\{N\}\(0,\{\\mathbf\{I\}\}\)
17:ifDomain = pixelthen

18:Compute step size

ϵ←\(r⋅‖𝚺i​d∘𝐳‖‖𝚺i​d∘𝐠‖\)2⋅2\\epsilon\\leftarrow\\mathopen\{\}\\mathclose\{\{\\left\(\\frac\{\\mathrm\{r\}\\cdot\\\|\\boldsymbol\{\\Sigma\}\_\{id\}\\circ\{\\mathbf\{z\}\}\\\|\}\{\\\|\\boldsymbol\{\\Sigma\}\_\{id\}\\circ\{\\mathbf\{g\}\}\\\|\}\}\}\\right\)^\{2\}\\cdot 2
19:Construct element\-wise step size

𝚺ϵ=𝚺i​d∘ϵ\\boldsymbol\{\\Sigma\}\_\{\\epsilon\}=\\boldsymbol\{\\Sigma\}\_\{id\}\\circ\\epsilon𝐱c\+1,s\+1←𝐱c\+1,s\+𝚺ϵ∘𝐠c\+1,s\+2​T​𝚺ϵ1/2∘𝐳\{\\mathbf\{x\}\}\_\{c\+1,s\+1\}\\leftarrow\{\\mathbf\{x\}\}\_\{c\+1,s\}\+\\boldsymbol\{\\Sigma\}\_\{\\epsilon\}\\circ\{\\mathbf\{g\}\}\_\{c\+1,s\}\+\\sqrt\{2\\,T\}\\,\\boldsymbol\{\\Sigma\}\_\{\\epsilon\}^\{1/2\}\\circ\{\\mathbf\{z\}\}
20:elseifDomain = frequencythen

21:Compute step size

ϵ←\(r⋅‖𝐳‖‖𝐠‖\)2⋅2\\epsilon\\leftarrow\\mathopen\{\}\\mathclose\{\{\\left\(\\frac\{\\mathrm\{r\}\\cdot\\\|\{\\mathbf\{z\}\}\\\|\}\{\\\|\{\\mathbf\{g\}\}\\\|\}\}\}\\right\)^\{2\}\\cdot 2𝐱c\+1,s\+1←𝐱c\+1,s\+ϵ\(𝐠c\+1,s\+∇𝐱tlogp\(𝐱T=y\|𝐱c\+1,s\)\)\+2​T​ϵ𝐳\{\\mathbf\{x\}\}\_\{c\+1,s\+1\}\\leftarrow\{\\mathbf\{x\}\}\_\{c\+1,s\}\+\\epsilon\\mathopen\{\}\\mathclose\{\{\\left\(\{\\mathbf\{g\}\}\_\{c\+1,s\}\+\\nabla\_\{\{\\mathbf\{x\}\}\_\{t\}\}\\log p\(\{\\mathbf\{x\}\}\_\{T\}=y\|\{\\mathbf\{x\}\}\_\{c\+1,s\}\)\}\}\\right\)\+\\sqrt\{2\\,T\\,\\epsilon\}\{\\mathbf\{z\}\}
22:endif

23:ifMALAthen

24:Compute Metropolis–Hastings acceptance probability

α=min\(1,exp⁡\(−Uθ​\(𝐱⋆,𝚺c\+1\)\)​q​\(𝐱c\+1,s∣𝐱⋆\)exp⁡\(−Uθ​\(𝐱c\+1,s,𝚺c\+1\)\)​q​\(𝐱⋆∣𝐱c\+1,s\)\)\\alpha=\\min\\\!\\mathopen\{\}\\mathclose\{\{\\left\(1,\\frac\{\\exp\\\!\\big\(\-U\_\{\\theta\}\(\{\\mathbf\{x\}\}^\{\\star\},\\boldsymbol\{\\Sigma\}\_\{c\+1\}\)\\big\)\\,q\(\{\\mathbf\{x\}\}\_\{c\+1,s\}\\mid\{\\mathbf\{x\}\}^\{\\star\}\)\}\{\\exp\\\!\\big\(\-U\_\{\\theta\}\(\{\\mathbf\{x\}\}\_\{c\+1,s\},\\boldsymbol\{\\Sigma\}\_\{c\+1\}\)\\big\)\\,q\(\{\\mathbf\{x\}\}^\{\\star\}\\mid\{\\mathbf\{x\}\}\_\{c\+1,s\}\)\}\}\}\\right\)
25:Sample

u∼𝒰​\(0,1\)u\\sim\\mathcal\{U\}\(0,1\)and accept/reject

26:

𝐱c\+1,s\+1←\{𝐱⋆,u<α𝐱c\+1,s,otherwise\{\\mathbf\{x\}\}\_\{c\+1,s\+1\}\\leftarrow\\begin\{cases\}\{\\mathbf\{x\}\}^\{\\star\},&u<\\alpha\\\\ \{\\mathbf\{x\}\}\_\{c\+1,s\},&\\text\{otherwise\}\\end\{cases\}
27:endif

28:endfor

29:endforReturn

𝐱0,ns​t​e​p​s\{\\mathbf\{x\}\}\_\{0,n\_\{steps\}\}

## Appendix BEnergy\-guided sampling as a Bregman mirror descent

While the gradient∇𝚺Uθ​\(𝐲,𝚺\)\\nabla\_\{\\boldsymbol\{\\Sigma\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)corresponds to the steepest descent with respect to the Euclidean norm, it is more natural here to consider the steepest descent with respect to the Bregman divergence generated by the negative log determinant:

D​\(𝚺∥𝚺′\)\\displaystyle D\(\\boldsymbol\{\\Sigma\}\\,\\\|\\,\\boldsymbol\{\\Sigma\}^\{\\prime\}\)=tr​\(𝚺′⁣−1​𝚺\)−log​det\(𝚺′⁣−1​𝚺\)−d\.\\displaystyle=\\text\{tr\}\(\\boldsymbol\{\\Sigma\}^\{\\prime\-1\}\\boldsymbol\{\\Sigma\}\)\-\\log\\det\(\\boldsymbol\{\\Sigma\}^\{\\prime\-1\}\\boldsymbol\{\\Sigma\}\)\-d\.\(37\)Note that this Bregman divergence corresponds to twice the KL divergence between two zero\-mean Gaussian distributions with covariance𝚺\\boldsymbol\{\\Sigma\}and𝚺′\\boldsymbol\{\\Sigma\}^\{\\prime\}\. For a stepsize parameterγ\>0\\gamma\>0, a Bregman mirror descent step computes

argmin𝚺′Uθ​\(𝐲,𝚺′\)\+1γ​D​\(𝚺′∥𝚺\)\\displaystyle\\operatornamewithlimits\{argmin\}\_\{\\boldsymbol\{\\Sigma\}^\{\\prime\}\}\\,U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}^\{\\prime\}\)\+\\frac\{1\}\{\\gamma\}D\(\\boldsymbol\{\\Sigma\}^\{\\prime\}\\,\\\|\\,\\boldsymbol\{\\Sigma\}\)=\(𝚺−1\+γ​∇𝚺Uθ​\(𝐲,𝚺\)\)−1\\displaystyle=\(\\boldsymbol\{\\Sigma\}^\{\-1\}\+\\gamma\\nabla\_\{\\boldsymbol\{\\Sigma\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\)^\{\-1\}\(38\)=𝚺−γ​𝚺​∇𝚺Uθ​\(𝐲,𝚺\)​𝚺\+o​\(γ\),\\displaystyle=\\boldsymbol\{\\Sigma\}\-\\gamma\\boldsymbol\{\\Sigma\}\\nabla\_\{\\boldsymbol\{\\Sigma\}\}U\_\{\\theta\}\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\\boldsymbol\{\\Sigma\}\+o\(\\gamma\),\(39\)where the last step is a first\-order Taylor approximation for small stepsizesγ→0\\gamma\\to 0\. This leads to the update

δ​𝚺=𝚺​∇𝚺U​\(𝐲,𝚺\)​𝚺\.\\delta\\boldsymbol\{\\Sigma\}=\\boldsymbol\{\\Sigma\}\\nabla\_\{\\boldsymbol\{\\Sigma\}\}U\(\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\)\\boldsymbol\{\\Sigma\}\.\(40\)

## Appendix CImplementation

### C\.1Architecture and training hyperparameters

Our architecture is a modification of the SongNet NCSNPP network\(Songet al\.,[2021b](https://arxiv.org/html/2605.15487#bib.bib37)\)implemented inKarraset al\.\([2022](https://arxiv.org/html/2605.15487#bib.bib118)\)\. We use the same hyperparameters as the original implementation, which are reported in Table[3](https://arxiv.org/html/2605.15487#A3.T3); the only modification is the noise embedding, which we describe below\. We also describe in Table[5](https://arxiv.org/html/2605.15487#A3.T5)the hyperparameters of the training\. The MNIST classifier used in[Figure6](https://arxiv.org/html/2605.15487#S4.F6)is a CNN with 3 layers and ReLU activation functions\.

Table 3:SongUNet Architecture Summary##### Noise Embedding Modulation\.

Given the diagonal noise of the covariancediag​\(𝚺t\)\\text\{diag\}\(\\boldsymbol\{\\Sigma\}\_\{t\}\), we first compute embeddings:

𝐞spatial\\displaystyle\{\\mathbf\{e\}\}\_\{\\text\{spatial\}\}←fspatial​\(diag​\(𝚺t\)\)∈ℝB×cemb×H×W,\\displaystyle\\leftarrow f\_\{\\text\{spatial\}\}\(\\text\{diag\}\(\\boldsymbol\{\\Sigma\}\_\{t\}\)\)\\in\\mathbb\{R\}^\{B\\times c\_\{\\text\{emb\}\}\\times H\\times W\},\(41\)𝐞frequency\\displaystyle\{\\mathbf\{e\}\}\_\{\\text\{frequency\}\}←ffrequency​\(diag​\(𝚺t\)\)∈ℝB×cemb×1×1,\\displaystyle\\leftarrow f\_\{\\text\{frequency\}\}\(\\text\{diag\}\(\\boldsymbol\{\\Sigma\}\_\{t\}\)\)\\in\\mathbb\{R\}^\{B\\times c\_\{\\text\{emb\}\}\\times 1\\times 1\},\(42\)wherefspatial\(\.\)f\_\{\\text\{spatial\}\}\(\.\)andffrequency​\(𝚺t\)f\_\{\\text\{frequency\}\}\(\\boldsymbol\{\\Sigma\}\_\{t\}\)are two residual architectures \(ResNet\)\(Heet al\.,[2015](https://arxiv.org/html/2605.15487#bib.bib153)\), consisting of3×33\\times 3convolutions paired with Group Normalization and SiLU activations, and a residual connection\. Given that the main architecture \(a UNet\) performs down and up\-sampling operations both in the channel and in the spatial domain, we need to accommodate the dimensions of the embedding network to each layerℓ\\ell\. For the channel dimension, we use1×11\\times 1convolutional layers as follows:

𝐞ℓ,spatial\\displaystyle\{\\mathbf\{e\}\}\_\{\\ell,\\text\{spatial\}\}←Conv1×1​\(𝐞spatial\)∈ℝB×cℓ×H×W\\displaystyle\\leftarrow\\text\{Conv\}\_\{1\\times 1\}\(\{\\mathbf\{e\}\}\_\{\\text\{spatial\}\}\)\\in\\mathbb\{R\}^\{B\\times c\_\{\\ell\}\\times H\\times W\}\(43\)𝐞ℓ,frequency\\displaystyle\{\\mathbf\{e\}\}\_\{\\ell,\\text\{frequency\}\}←Conv1×1​\(𝐞frequency\)∈ℝB×cℓ×1×1\\displaystyle\\leftarrow\\text\{Conv\}\_\{1\\times 1\}\(\{\\mathbf\{e\}\}\_\{\\text\{frequency\}\}\)\\in\\mathbb\{R\}^\{B\\times c\_\{\\ell\}\\times 1\\times 1\}\(44\)For the spatial dimensions, we use the following operations

𝐞ℓ,spatial\\displaystyle\{\\mathbf\{e\}\}\_\{\\ell,\\text\{spatial\}\}←Interpolate​\(𝐞ℓ,spatial\)∈ℝB×cℓ×Hℓ×Wℓ\\displaystyle\\leftarrow\\text\{Interpolate\}\(\{\\mathbf\{e\}\}\_\{\\ell,\\text\{spatial\}\}\)\\in\\mathbb\{R\}^\{B\\times c\_\{\\ell\}\\times H\_\{\\ell\}\\times W\_\{\\ell\}\}\(45\)Given the embedding vector for each layer, we apply a multiplicative conditioning mechanism, given by

𝐱ℓ=SiLU\(GroupNorm\(𝐱ℓ\)⊙\(1\+𝐞ℓ,spatial/frequency\)\)\{\\mathbf\{x\}\}\_\{\\ell\}=\\text\{SiLU\}\\mathopen\{\}\\mathclose\{\{\\left\(\\text\{GroupNorm\}\(\{\\mathbf\{x\}\}\_\{\\ell\}\)\\odot\(1\+\{\\mathbf\{e\}\}\_\{\\ell,\\text\{spatial/frequency\}\}\)\}\}\\right\)\(46\)whereGroupNorm\(\.\)\\text\{GroupNorm\}\(\.\)is the group normalization, with a number of channels in each group ofmin⁡\(32,cℓ/4\)\\min\(32,c\_\{\\ell\}/4\)\.

Table 4:Training hyperparameters for CelebA and ImageNet64
Table 5:Training hyperparameters for AFHQ

### C\.2Distribution covariances used for training

As described in Section[3\.3](https://arxiv.org/html/2605.15487#S3.SS3), we used four types of covariances: Gaussian deblurring, shown in Fig\.[8\(a\)](https://arxiv.org/html/2605.15487#A3.F8.sf1), super\-resolution \(a combination of Gaussian kernel with a downsampling operation\), and box and half\-mask inpainting, illustrated in Fig\.[8\(b\)](https://arxiv.org/html/2605.15487#A3.F8.sf2)\. The family of autoregressive patch\-based covariances used in the MNIST experiments in Section[4\.2\.1](https://arxiv.org/html/2605.15487#S4.SS2.SSS1)are shown in Fig\.[9](https://arxiv.org/html/2605.15487#A3.F9)\.

![Refer to caption](https://arxiv.org/html/2605.15487v1/x12.png)\(a\)Gaussian Kernel
![Refer to caption](https://arxiv.org/html/2605.15487v1/x13.png)\(b\)Box and half\-mask inpainting

Figure 8:Examples of particular𝚺t\\boldsymbol\{\\Sigma\}\_\{t\}for different kernels and masks\.![Refer to caption](https://arxiv.org/html/2605.15487v1/x14.png)Figure 9:Examples of autoregressive patch\-based covariances𝚺t\\boldsymbol\{\\Sigma\}\_\{t\}used for training\. We first sample the number of boxes, and then the variance of each box\.
### C\.3Comparison with previous energy\-based models

Our method is more computationally intensive to train than standard score\-based models\. This overhead primarily comes from double backpropagation, which is inherent to energy\-based formulations and is widely acknowledged in prior work\(Thorntonet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib125); Balceraket al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib149); Aiet al\.,[2026](https://arxiv.org/html/2605.15487#bib.bib150); Duet al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib128)\)\. In this sense, the increased cost is not specific to our method, but rather a known trade\-off when learning energy models\. That said, we further contextualize below that our model has a computational cost comparable to prior work when accounting for model size and resolution\.

1. 1\.We use the Song architecture from EDM \(≈70\\approx 70M parameters\), which is substantially larger than the network used in the original DSM paper\(Guthet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib114)\)\(≈13\\approx 13M parameters\)\. The larger architecture naturally increases training time, while also yielding improved performance\. For reference, EDM training at64×6464\\times 64takes approximately 4 days on 8 V100 GPUs \(≈4\\approx 4A100s\) with a batch size of 256 \(see[https://github\.com/NVlabs/edm](https://github.com/NVlabs/edm)\)\. Scaling to192×192192\\times 192increases the pixel count by9×9\\times, which significantly increases the computational cost, either in terms of GPU count or total training time\.
2. 2\.This computational challenge is shared by energy\-based models more broadly: - •Thornton\(Thorntonet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib125)\)used 8 A100 GPUs at64×6464\\times 64resolution \(without the dual objective\)\. - •Balcerak\(Balceraket al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib149)\)trained at64×6464\\times 64on CelebA with 4 A100s and a batch size of 32; neither\(Thorntonet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib125)\)nor\(Balceraket al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib149)\)report total training time\. - •Boffi\(Aiet al\.,[2026](https://arxiv.org/html/2605.15487#bib.bib150)\)required approximately 8 days for the flow matching teacher and 6 days for the shortcut model on ImageNet64×6464\\times 64\. - •Du\(Duet al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib128)\)trained at128×128128\\times 128using TPU clusters, making direct comparison difficult\.

### C\.4Details of baselines

##### RED\-Diff\.

It is a variational framework that frames the sampling problem as stochastic optimization by minimizing the KL divergence

q\(𝐱0\|𝐲\)=argminq​\(𝐱0\|𝐲\)KL\(q\(𝐱0\|𝐲\)\|\|p\(𝐱0\|𝐲\)\)\.q\(\{\\mathbf\{x\}\}\_\{0\}\|\{\\mathbf\{y\}\}\)=\\operatornamewithlimits\{argmin\}\_\{q\(\{\\mathbf\{x\}\}\_\{0\}\|\{\\mathbf\{y\}\}\)\}\\mathrm\{KL\}\(q\(\{\\mathbf\{x\}\}\_\{0\}\|\{\\mathbf\{y\}\}\)\|\|p\(\{\\mathbf\{x\}\}\_\{0\}\|\{\\mathbf\{y\}\}\)\)\.\(47\)Whenq​\(𝐱0\|𝐲\)∼𝒩​\(𝝁,σ2​𝐈\)q\(\{\\mathbf\{x\}\}\_\{0\}\|\{\\mathbf\{y\}\}\)\\sim\{\\mathcal\{N\}\}\(\\boldsymbol\{\\mu\},\\sigma^\{2\}\{\\mathbf\{I\}\}\), the gradient update boils down to

∇𝝁=𝔼t∼𝒰​\[0,T\]\[∇𝝁∥𝐲−𝐇𝝁∥2\+𝔼ϵ∼𝒩​\(0,I\)\[λt\(ϵ𝜽\(𝐱t,σt\)−ϵ\)\]\]\\nabla\_\{\\boldsymbol\{\\mu\}\}=\\mathbb\{E\}\_\{t\\sim\\mathcal\{U\}\[0,T\]\}\\mathopen\{\}\\mathclose\{\{\\left\[\\nabla\_\{\\boldsymbol\{\\mu\}\}\\\|\{\\mathbf\{y\}\}\-\{\\mathbf\{H\}\}\\boldsymbol\{\\mu\}\\\|^\{2\}\+\\mathbb\{E\}\_\{\\epsilon\\sim\\mathcal\{N\}\(0,I\)\}\\mathopen\{\}\\mathclose\{\{\\left\[\\lambda\_\{t\}\\mathopen\{\}\\mathclose\{\{\\left\(\\boldsymbol\{\\epsilon\}\_\{\\boldsymbol\{\\theta\}\}\(\{\\mathbf\{x\}\}\_\{t\},\\sigma\_\{t\}\)\-\\boldsymbol\{\\epsilon\}\}\}\\right\)\}\}\\right\]\}\}\\right\]whereϵ𝜽​\(𝐱t,σt\)=−σt​𝐬𝜽​\(𝐱t,σt\)\\boldsymbol\{\\epsilon\}\_\{\\boldsymbol\{\\theta\}\}\(\{\\mathbf\{x\}\}\_\{t\},\\sigma\_\{t\}\)=\-\\sigma\_\{t\}\{\\mathbf\{s\}\}\_\{\\boldsymbol\{\\theta\}\}\(\{\\mathbf\{x\}\}\_\{t\},\\sigma\_\{t\}\)and𝐱t=αt​𝝁\+σt​ϵ\{\\mathbf\{x\}\}\_\{t\}=\\alpha\_\{t\}\\boldsymbol\{\\mu\}\+\\sigma\_\{t\}\\boldsymbol\{\\epsilon\}; in the signal domain, we have\(ϵ𝜽​\(𝐱t,σt\)−ϵ\)∝\(𝝁−𝔼​\[𝐱0\|𝐱t\]\)\(\\boldsymbol\{\\epsilon\}\_\{\\boldsymbol\{\\theta\}\}\(\{\\mathbf\{x\}\}\_\{t\},\\sigma\_\{t\}\)\-\\boldsymbol\{\\epsilon\}\)\\propto\(\\boldsymbol\{\\mu\}\-\\mathbb\{E\}\[\{\\mathbf\{x\}\}\_\{0\}\|\{\\mathbf\{x\}\}\_\{t\}\]\)\. In this formulation, the full trajectory acts as multi\-scale regularizer, imposing coarse\-to\-fine details\. We adapted the implementation from the original paper\(Mardaniet al\.,[2024](https://arxiv.org/html/2605.15487#bib.bib6)\)to our code\. We follow the same weighting scheme for inpainting, and we useλ=0\.25\\lambda=0\.25andlr=0\.1l\_\{r\}=0\.1, while for Gaussina deblurring we useλ=2\\lambda=2\. As mentioned in Section[3\.4](https://arxiv.org/html/2605.15487#S3.SS4), we train from scratch an isotropic version of the proposed network\.

##### DPS\.

We adapted the implementation of the original paper\(Chunget al\.,[2023](https://arxiv.org/html/2605.15487#bib.bib54)\)to our code\. We tried using the same weight for the gradient of the likelihood, but for our implementation and a box size of 45, the sampler was diverging\. Hence, we perform a grid search to get the best performance, with 0\.05 the best one for inpainting, while 0\.01 for deblurring\. For all the experiments, we match the number of steps than our energy\-guided sampler to have a fair comparison\. We use the same base network as RED\-Diff\.

##### DAPS\.

We also implemented DAPS\(Zhanget al\.,[2025a](https://arxiv.org/html/2605.15487#bib.bib130)\), which alternates between ODE\-based denoising and Langevin dynamics for measurement consistency at each noise level\. At each step, a short ODE flow \(5 steps\) estimates the clean imagex^0\\hat\{x\}\_\{0\}, followed by Langevin sampling that minimizes both the measurement residual and a prior term anchored tox^0\\hat\{x\}\_\{0\}\. We use 100 Langevin steps per noise level with step sizeη=2×10−6\\eta=2\\times 10^\{\-6\}, and250250\(for inpainting\) and200200\(for deblurring\) number of ODE runs\. We use the same base network as RED\-Diff\.

##### Palette\.

We follow the formulation from the original paper\(Sahariaet al\.,[2022](https://arxiv.org/html/2605.15487#bib.bib79)\)\. Hence, we modified the Song architecture by stacking the noisy image𝐱t\{\\mathbf\{x\}\}\_\{t\}with the measurement𝐲\{\\mathbf\{y\}\}, and using the same base architecture that we used for the guidance methods\.

##### Ours\.

We describe the setup for our method in Table[6](https://arxiv.org/html/2605.15487#A3.T6)\.

Table 6:Inpainting \- Hyperparameters for the Predictor\-Corrector Sampler with fixed corrector\.
##### Computational cost for inference\.

We expand on the running times for Gaussian deblurring across all methods, considering the running time for sampling a batch of 16 images at64×6464\\times 64:

Importantly, the inference time of our method is comparable to DPS and standard diffusion samplers, so the computational overhead is concentrated entirely in training and does not affect sampling time\. In contrast to these methods, our energy formualtion learns a full posterior density\. We believe this makes the overhead modest relative to the added capabilities our approach unlocks, such as blind estimation and normalization constant computation, which are not possible with DPS or RED\-Diff\.

For sampling, the additional∇yUθ\\nabla\_\{y\}U\_\{\\theta\}gradient is similar in cost to DPS’s backpropagation and architectural overhead is minimal\. Thus, wall\-clock time is approximately equal at the same NFEs\. For AFHQ192×192192\\times 192: DPS 55 sec, ours 74 sec\.

## Appendix DAdditional experiments

### D\.1Additional inverse problems

#### D\.1\.1In\-distribution covariances

We include the results forsuper\-resolution×4\\times 4andhorizontal half\-mask inpaintingin Table[7](https://arxiv.org/html/2605.15487#A4.T7); visual results can be found in Appendix[D\.6](https://arxiv.org/html/2605.15487#A4.SS6)\.

For super\-resolution, we do not use corrector steps; in terms of numerical results, our method outperforms Bayesian approaches by a significant margin, similarly to the deblurring case\. As illustrated in Fig\.[26](https://arxiv.org/html/2605.15487#A4.F26), the generated images present some artifacts\. We attribute this partly to the known limitations of the frequency\-domain conditioning: our frequency\-domain degradations are defined with a Fourier transform with periodic boundary conditions, which can cause border effects\. For practical use, a more careful implementation with more advanced boundary handling \(e\.g\., a suitable discrete cosine transform\) should be preferred\.

Regarding half\-mask inpainting, our method achieves better FID than baselines\. However, in contrast to box inpainting, our method have a similar LPIPS than Bayesian baselines\.

Table 7:Experiments on CelebA64×6464\\times 64: Comparison on super\-resolution×4\\times 4and half\-mask inpainting\.
#### D\.1\.2Out\-of\-distribution covariances

In this subsection, we illustrate how our method can handle out\-of\-distribution degradations, i\.e\., those inverse problems with an associated covariance𝚺\\boldsymbol\{\\Sigma\}not used for training\. In particular, we considerrandom inpainting, where we mask out a 70% of pixels,vertical half mask inpainting, with a mask of size of64×3064\\times 30, andmotion deblurring; for the latter, we considered the kernel from[https://github\.com/LeviBorodenko/motionblur](https://github.com/LeviBorodenko/motionblur)of size15×1515\\times 15and intensity 0\.3, andσv=2×10−3\\sigma\_\{v\}=2\\times 10^\{\-3\}; and for scheduling, we use a Gaussian deblurring kernel of size15×1515\\times 15and intensity0\.80\.8, which serves as an approximation\. The results are shown in Table[8](https://arxiv.org/html/2605.15487#A4.T8)and Table[9](https://arxiv.org/html/2605.15487#A4.T9)\. Remarkably, we observe that our models achieves a similar performance than Bayesian methods\.

Table 8:Experiments on CelebA64×6464\\times 64: Comparison on random inpainting \(0\.7%\) and vertical half\-mask inpainting\.Table 9:Experiments on CelebA64×6464\\times 64: Comparison on motion deblurring\.
#### D\.1\.3Large\-scale experiments

In this subsection, we include experiments on AFHQ\-Cat192×192192\\times 192, demonstrating that our proposed model can be trained on large\-scale images\. To simplify training, we consider only inpainting with a center\-box of sizess×ss\\times s, with values froms=20s=20tos=50s=50, and half\-mask with sizess×64s\\times 64, for the same values ofss\.

##### Comparison with baselines\.

We compare with DPS for a box inpainting task; the result is shown in Tables[11](https://arxiv.org/html/2605.15487#A4.T11)and[13](https://arxiv.org/html/2605.15487#A4.T13)\.

Table 10:Box30×3030\\times 30Table 11:Experiments on CelebA64×6464\\times 64: Comparison on box inpainting30×3030\\times 30\.Table 12:Box50×5050\\times 50Table 13:Experiments on CelebA64×6464\\times 64: Comparison on box inpainting50×5050\\times 50\.
##### Blind experiment\.

We also consider a blind experiment, following the setup from Section[4\.3](https://arxiv.org/html/2605.15487#S4.SS3)\. In particular, we considerσ2=10−4\\sigma\_\{2\}=10^\{\-4\}, andσ1,s=\{1,30\}\\sigma\_\{1\},s=\\\{1,30\\\}\. The result is illustrated in Fig\.[11](https://arxiv.org/html/2605.15487#A4.F11); the value of the estimatedσ1\\sigma\_\{1\}is 0\.2\.

![Refer to caption](https://arxiv.org/html/2605.15487v1/x15.png)Figure 10:Noisy observation with a box covariance of size30×3030\\times 30andσ1=1\\sigma\_\{1\}=1\.![Refer to caption](https://arxiv.org/html/2605.15487v1/x16.png)
![Refer to caption](https://arxiv.org/html/2605.15487v1/x17.png)

Figure 11:Blind reconstruction experiment\. Log probabilitylog⁡pθ​\(𝐲\|𝚺s\)\\log p\_\{\\theta\}\(\{\\mathbf\{y\}\}\|\\boldsymbol\{\\Sigma\}\_\{s\}\)as a function of box sizessandσ1\\sigma\_\{1\}, respectively\. The energy model trained with dual score matching correctly identifies the true box size \(dashed red for estimated size and solid green vertical lines\)\.

### D\.2Comparison between dual and single score matching

##### Blind experiment\.

Validating the correctness of the dual score matching is difficult since we do not have access to the the true density\. However, we can use some experiments to assess this, like the blind experiment in Section[4\.3](https://arxiv.org/html/2605.15487#S4.SS3)\. In particular, this experiment aims to find the covariance that best explains a given measurement, relying on a proper density across different covariances\. While in Section[4\.3](https://arxiv.org/html/2605.15487#S4.SS3)we ony compared the dual score matching with the ground truth, here we add the single objective to show that, in this additional case, the estimation fails; this is illustrated in Fig\.[12](https://arxiv.org/html/2605.15487#A4.F12)\.

![Refer to caption](https://arxiv.org/html/2605.15487v1/x18.png)Figure 12:Log probabilitylog⁡pθ​\(𝐲\|𝚺s\)\\log p\_\{\\theta\}\(\{\\mathbf\{y\}\}\|\\boldsymbol\{\\Sigma\}\_\{s\}\)as a function of box sizessfors=\{5,14,20\}s=\\\{5,14,20\\\}withσ1=\{2,0\.1,2\}\\sigma\_\{1\}=\\\{2,0\.1,2\\\}respectively\. The energy model trained with dual score matching correctly identifies the true box size \(dashed red for estimated size and solid green vertical lines\), which is not the case for the model trained with single score matching \(dashed blue vertical line\)\.
##### One\-shot conditional MMSE estimator\.

One advantage of our proposed method is that we can compute𝔼​\[𝐱\|𝐲,𝚺\]\\mathbb\{E\}\[\{\\mathbf\{x\}\}\|\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\]in one shot\. We illustrate this for a denoising task with box of size20×2020\\times 20\(the outside does not have noise\), where we evaluate the one\-shot denoiser for multiple noise level in the box\. The ablation is shown in Fig\.[13](https://arxiv.org/html/2605.15487#A4.F13), where we compute the MSE as a function ofσ12\\sigma\_\{1\}^\{2\}, the noise inside the box\. In this experiment, we also consider our proposed energy model trained with the single objective \(only A\-DSM\)\. From this experiment, we observe that the energy\-model trained with the dual\-objective yields a better model in terms of denoising performance, showcasing that the regularization not only allows learning a normalized density model, but also helps for improving the denoising performance\. To complement this experiment, we illustrate a few visual examples in Fig\.[14](https://arxiv.org/html/2605.15487#A4.F14)using our energy model trained with the dual objective, and for a covariance with box of size10×1010\\times 10,σ1=0\.5\\sigma\_\{1\}=0\.5\(inside the box\) andσ2=0\.1\\sigma\_\{2\}=0\.1\(outside the box\)\.

![Refer to caption](https://arxiv.org/html/2605.15487v1/x19.png)
![Refer to caption](https://arxiv.org/html/2605.15487v1/x20.png)

Figure 13:MMSE for one\-shot denoising in box inpainting \(10×1010\\times 10\) as a function of the noise in the box\.![Refer to caption](https://arxiv.org/html/2605.15487v1/x21.png)Figure 14:MMSE estimator computed by our model, which corresponds to𝐱^=𝔼​\[𝐱\|𝐲,𝚺\]\{\\hat\{\\mathbf\{x\}\}\}=\\mathbb\{E\}\[\{\\mathbf\{x\}\}\|\{\\mathbf\{y\}\},\\boldsymbol\{\\Sigma\}\]\.

### D\.3Analysis of the posterior distribution

##### CelebA\.

We include here additional results for the analysis of the posterior distribution described in Section[4\.1](https://arxiv.org/html/2605.15487#S4.SS1)\. We start showing the distribution across multiple\(𝐱,𝐲\)\(\{\\mathbf\{x\}\},\{\\mathbf\{y\}\}\)in Fig\.[15](https://arxiv.org/html/2605.15487#A4.F15), with generated images sorted from high to low probability in Fig\.[16](https://arxiv.org/html/2605.15487#A4.F16)\. Notice that the behavior is similar to the single case, where our energy model is the one closer the posterior \(left figure\)\.

![Refer to caption](https://arxiv.org/html/2605.15487v1/x22.png)Figure 15:Histograms oflog⁡pθ​\(𝐱^\)\\log p\_\{\\theta\}\(\{\\hat\{\\mathbf\{x\}\}\}\)andlog⁡pθ​\(𝐱^\|𝐲\)\\log p\_\{\\theta\}\(\{\\hat\{\\mathbf\{x\}\}\}\|\{\\mathbf\{y\}\}\)for inpainting solutions𝐱^\{\\hat\{\\mathbf\{x\}\}\}generated by DPS, RED\-Diff, and our energy model for different measurements𝐲\{\\mathbf\{y\}\}associated to different𝐱\{\\mathbf\{x\}\}, along with the ground truths𝐱\{\\mathbf\{x\}\}\. Our energy model is well\-calibrated with respect to both prior and posterior probabilities\.![Refer to caption](https://arxiv.org/html/2605.15487v1/x23.png)Figure 16:Images sorted from high to low probability \(from left\-up corner to right\-down corner\)\. Images with highest prior probability contains less details than those associated with higer probability\.
##### ImageNet\.

We evaluate our model on ImageNet64, first validating that the learned prior aligns with the complexity\-probability relationship described inGuthet al\.\([2025](https://arxiv.org/html/2605.15487#bib.bib114)\)\. As illustrated in Fig\.[17](https://arxiv.org/html/2605.15487#A4.F17), images sorted by decreasing probability transition from sparse, low\-detail structures to high\-complexity textures; the highest\-probability samples typically feature single objects on plain backgrounds, while the lowest\-probability samples are densely textured\.

We extend this analysis to samples sorted by the posterior distributionlog⁡p​\(𝐱\|𝐲\)\\log p\(\\mathbf\{x\}\|\\mathbf\{y\}\)in Fig\.[18](https://arxiv.org/html/2605.15487#A4.F18)\. We observe that the highest posterior probability does not necessarily correlate with the best reconstruction or visual quality\. Furthermore, the complexity\-probability trend observed in the prior persists in the posterior case\.

![Refer to caption](https://arxiv.org/html/2605.15487v1/x24.png)Figure 17:ImageNet posterior samples for a21×2121\\times 21center inpainting task, sorted by decreasing prior probability\. Consistent with previous observations, samples with the highest probability feature uniform backgrounds and minimal detail, whereas lower\-probability samples exhibit significantly higher structural complexity and dense detail\.![Refer to caption](https://arxiv.org/html/2605.15487v1/x25.png)Figure 18:ImageNet posterior samples \(red border\) generated by sampling from measurements𝐲\{\\mathbf\{y\}\}generated using a box inpainting mask with a box of size21×2121\\times 21, and sorted from high to low posterior probability, next to their corresponding ground truth \(blue border\)\.

### D\.4Recovering the Gumbel distribution for the isotropic case

Our model generalizes the energy model of\(Guthet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib114)\), which is recovered as a special case by setting𝚺t=σt2​𝐈\\boldsymbol\{\\Sigma\}\_\{t\}=\\sigma\_\{t\}^\{2\}\\mathbf\{I\}\(isotropic covariance\)\. Under this reduction, the anisotropic score matching objective reduces to the standard dual score matching of \[2\], and the Gumbel distribution result therefore carries over directly\. We include a curve computed over 50k images alongside with some of the images, which can be found in Figs\.[19](https://arxiv.org/html/2605.15487#A4.F19)and[20](https://arxiv.org/html/2605.15487#A4.F20)

![Refer to caption](https://arxiv.org/html/2605.15487v1/x26.png)Figure 19:Histogram of log probabilitieslog⁡pθ​\(𝐱\)​\[dB/dim\]\\log p\_\{\\theta\}\(\{\\mathbf\{x\}\}\)\[\\text\{dB\}/\\text\{dim\}\]for 50k in the ImageNet test set\. The distribution is well\-fit by a Gumbel distribution \(red line\), but there is a second mode \(centered near 38 dB/dim\) in contrast to the model in \[1\]\. This additional mode is caused by grayscale images within ImageNet; these samples possess significantly higher probability under the modelUθU\_\{\\theta\}compared to standard RGB samples, creating an additional mode\.![Refer to caption](https://arxiv.org/html/2605.15487v1/x27.png)Figure 20:Images from the test set of ImageNet64×6464\\times 64, sorted from low to high probability \(left to right and top to down\)\.
### D\.5Validation of the normalization constant in synthethic experiment

Analogous to the experiments of\(Guthet al\.,[2025](https://arxiv.org/html/2605.15487#bib.bib114)\), we generaten=100,000n=100\{,\}000samples from a mixture of two Gaussian distributions,12​𝒩​\(0,σ12​Id\)\+12​𝒩​\(0,σ22​Id\)\\frac\{1\}\{2\}\\mathcal\{N\}\(0,\\sigma\_\{1\}^\{2\}\\mathrm\{Id\}\)\+\\frac\{1\}\{2\}\\mathcal\{N\}\(0,\\sigma\_\{2\}^\{2\}\\mathrm\{Id\}\), withσ1=1\\sigma\_\{1\}=1andσ2=4\\sigma\_\{2\}=4, in dimensiond=1,000d=1\{,\}000\. Each sample receives a random partition of thedddimensions into two equal groupsAAandBBof sized/2d/2, with two independent noise levelsttop,tbot∼UniformLog​\(tmin,tmax\)t\_\{\\mathrm\{top\}\},t\_\{\\mathrm\{bot\}\}\\sim\\mathrm\{UniformLog\}\(t\_\{\\min\},t\_\{\\max\}\), yielding per\-dimension varianceΣt,i​i=ttop\\Sigma\_\{t,ii\}=t\_\{\\mathrm\{top\}\}ifi∈Ai\\in AandΣt,i​i=tbot\\Sigma\_\{t,ii\}=t\_\{\\mathrm\{bot\}\}ifi∈Bi\\in B\. The true \(normalized\) energy is given by

U\(y,Σt\)=−log\(∑i=12e−12​y⊤​\(σi2​Id\+Σt\)−1​y−12​log​det\(2​π​\(σi2​Id\+Σt\)\)−log⁡2\)\.U\(y,\\Sigma\_\{t\}\)=\-\\log\\\!\\mathopen\{\}\\mathclose\{\{\\left\(\\sum\_\{i=1\}^\{2\}\\mathrm\{e\}^\{\-\\frac\{1\}\{2\}y^\{\\top\}\(\\sigma\_\{i\}^\{2\}\\mathrm\{Id\}\+\\Sigma\_\{t\}\)^\{\-1\}y\\,\-\\,\\frac\{1\}\{2\}\\log\\det\\\!\\bigl\(2\\pi\(\\sigma\_\{i\}^\{2\}\\mathrm\{Id\}\+\\Sigma\_\{t\}\)\\bigr\)\\,\-\\,\\log 2\}\}\}\\right\)\.\(48\)We parameterize the energy as a mixture of quadraticsUθ\(y,ttop,tbot\)=−log\(∑i=12e−aitop​\(t\)​rtop2−aibot​\(t\)​rbot2−bi​\(t\)\),U\_\{\\theta\}\(y,t\_\{\\mathrm\{top\}\},t\_\{\\mathrm\{bot\}\}\)=\-\\log\\\!\\mathopen\{\}\\mathclose\{\{\\left\(\\sum\_\{i=1\}^\{2\}\\mathrm\{e\}^\{\-a\_\{i\}^\{\\mathrm\{top\}\}\(t\)\\,r\_\{\\mathrm\{top\}\}^\{2\}\\,\-\\,a\_\{i\}^\{\\mathrm\{bot\}\}\(t\)\\,r\_\{\\mathrm\{bot\}\}^\{2\}\\,\-\\,b\_\{i\}\(t\)\}\}\}\\right\),wherertop2=‖yA‖2r\_\{\\mathrm\{top\}\}^\{2\}=\\\|y\_\{A\}\\\|^\{2\}andrbot2=‖yB‖2r\_\{\\mathrm\{bot\}\}^\{2\}=\\\|y\_\{B\}\\\|^\{2\}are the squared norms over each group, and the functionsaitop,aibot,bia\_\{i\}^\{\\mathrm\{top\}\},a\_\{i\}^\{\\mathrm\{bot\}\},b\_\{i\}are computed by a 5\-layer MLP with a hidden dimension of 256 that takes the 2\-dimensional schedule embedding\[log⁡\(ttop\+tmin\),log⁡\(tbot\+tmin\)\]\[\\log\(t\_\{\\mathrm\{top\}\}\+t\_\{\\min\}\),\\,\\log\(t\_\{\\mathrm\{bot\}\}\+t\_\{\\min\}\)\]as input \(reducing to the isotropic case whenttop=tbott\_\{\\mathrm\{top\}\}=t\_\{\\mathrm\{bot\}\}\)\. This network is trained either with single \(space\) score matching or with dual score matching, both across noise levelsttop,tbot∈\[tmin,tmax\]t\_\{\\mathrm\{top\}\},t\_\{\\mathrm\{bot\}\}\\in\[t\_\{\\min\},t\_\{\\max\}\], where the dual loss includes one time\-derivative term per noise level\. Training is otherwise similar to the isotropic experiment, for50,00050\{,\}000training steps with a batch size of512512and an initial learning rate of0\.00010\.0001, over noise levels fromtmin=10−2t\_\{\\min\}=10^\{\-2\}totmax=102t\_\{\\max\}=10^\{2\}\. The comparison is illustrated in Fig\.[21](https://arxiv.org/html/2605.15487#A4.F21)

![Refer to caption](https://arxiv.org/html/2605.15487v1/x28.png)Figure 21:Comparison of single and dual score matching on a10001000\-dimensional Gaussian scale mixture\. Left: Dual score matching \(red dashed\) captures the true energy \(blue solid\), whereas single score matching \(green\) fails, even after global normalization\. Right: Radial score components\. Single score matching learns the score accurately within the data support \(blue bar plot\) but lacks accuracy outside of this region\.
### D\.6Reconstructed images

![Refer to caption](https://arxiv.org/html/2605.15487v1/x29.png)Figure 22:Examples of inpainting for CelebA\.![Refer to caption](https://arxiv.org/html/2605.15487v1/x30.png)Figure 23:Examples of inpainting for ImageNet\.![Refer to caption](https://arxiv.org/html/2605.15487v1/x31.png)Figure 24:Examples of deblurring for CelebA![Refer to caption](https://arxiv.org/html/2605.15487v1/x32.png)Figure 25:Examples of deblurring for ImageNet![Refer to caption](https://arxiv.org/html/2605.15487v1/x33.png)Figure 26:Examples of super\-resolution for CelebA![Refer to caption](https://arxiv.org/html/2605.15487v1/x34.png)Figure 27:Examples of inpainting for AFHQ\-cats\.![Refer to caption](https://arxiv.org/html/2605.15487v1/x35.png)Figure 28:Examples of Motion Deblurring\.![Refer to caption](https://arxiv.org/html/2605.15487v1/x36.png)Figure 29:Comparison of the different corrector samplers a function of the number of steps\.