# An analytic theory of convolutional neural network inverse problems solvers

Minh Hai Nguyen, Quoc Bao Do, Edouard Pauwels, Pierre Weiss

<figure><img src="image_1.png" /></figure>

## To cite this version:

Minh Hai Nguyen, Quoc Bao Do, Edouard Pauwels, Pierre Weiss. An analytic theory of convolutional neural network inverse problems solvers. 2026. hal-05459284

HAL Id: hal-05459284

https://hal.science/hal-05459284v1

Preprint submitted on 15 Jan 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
Minh-Hai Nguyen<fnref n="1" /> Quoc-Bao Do² Edouard Pauwels³ Pierre Weiss<fnref n="1" />

Abstract

Supervised convolutional neural networks (CNNs) are widely used to solve imaging inverse problems, achieving state-of-the-art performance in numerous applications. However, despite their empirical success, these methods are poorly understood from a theoretical perspective and often treated as black boxes. To bridge this gap, we analyze trained neural networks through the lens of the Minimum Mean Square Error (MMSE) estimator, incorporating functional constraints that capture two fundamental inductive biases of CNNs: translation equivariance and locality via finite receptive fields. Under the empirical training distribution, we derive an analytic, interpretable, and tractable formula for this constrained variant, termed Local-Equivariant MMSE (LE-MMSE). Through extensive numerical experiments across various inverse problems (denoising, inpainting, deconvolution), datasets (FFHQ, CIFAR-10, FashionMNIST), and architectures (U-Net, ResNet, PatchMLP), we demonstrate that our theory matches the neural networks outputs (PSNR ≳ 25dB). Furthermore, we provide insights into the differences between physics-aware and physics-agnostic estimators, the impact of high-density regions in the training (patch) distribution, and the influence of other factors (dataset size, patch size, etc.).

1. Introduction

Let $A: \mathbb{R}^N \to \mathbb{R}^M$ be a linear or nonlinear mapping, $x$ a random vector in $\mathbb{R}^N$. We consider the forward model with additive white Gaussian noise

$$y = Ax + e \quad \text{where} \quad e \sim \mathcal{N}(0, \sigma^2 I_M) \quad (1)$$

<figure><img src="image_3.png" /><figcaption>Figure 1: Our analytic theory accurately predicts neural network outputs across settings. We consider three inverse problems: denoising, inpainting, and deconvolution (left to right), on FFHQ, CIFAR10, and FashionMNIST datasets (top to bottom) with varying noise levels $\sigma$ (columns). For each setting, we show the measurements (top row), our analytic LE-MMSE estimator (second row), and outputs of trained UNet, ResNet, and PatchMLP models (last three rows). Theory closely matches network outputs.</figcaption></figure>

<footnote n="1">¹IRIT, Toulouse University, France ²INSA Toulouse, France ³Toulouse School of Economics, France. Correspondence to: Minh-Hai Nguyen <minh-hai.nguyen@irit.fr>.</footnote>
and aim to construct an estimator $\hat{x}(y)$ of $x$ given $y$. A variety of neural network architectures have been proposed for this task. Some are physics-agnostic, i.e. do not depend on $A$ explicitly and often rely on fully connected layers (Zhu et al., 2018). Others are physics-aware, and rely on an approximate inverse of $A$. In this work, we focus on linear operators $A$ in the numerical experiments for didactic reasons, but our theory applies to nonlinear operators as well. We consider estimators of the following form

$$ \hat{x}(y) = N_w(By), \quad (2) $$

where $N_w$ is a neural network with weights $w \in \Theta$ and $B \in \mathbb{R}^{N \times M}$ is an operator such as the identity when $M = N$ (physics-agnostic), or an approximate inverse of $A$ (physics-aware) such as the transpose $A^\top$, the pseudo-inverse $A^+$, or a Tikhonov-regularized inverse. Such designs were popularized in (Kim et al., 2016; Jin et al., 2017; Schlemper et al., 2017; Würfl et al., 2016) and are still widely used in practice (McCann et al., 2017; Zhou et al., 2022; Wang et al., 2020). They can be considered as a robust baseline, and be interpreted as a single-step variant of unrolled neural networks (Adler & Öktem, 2018b; Celle-doni et al., 2021), which often achieve top performance in empirical benchmarks (Muckley et al., 2021).

## Empirical MMSE estimators

**Definition 1.1** (Constrained empirical MMSE estimators). Let $\mathcal{D}$ be a dataset of clean images. We define the empirical MMSE estimators as $\hat{x}_{\mathcal{M}} \stackrel{\text{def}}{=} \phi_{\mathcal{M}}^* \circ B$, where

$$ \phi_{\mathcal{M}}^* \stackrel{\text{def}}{=} \arg \min_{\phi \in \mathcal{M}} \frac{1}{2} \mathbb{E} \left[ \|\phi(By) - x\|^2 \right]. \quad (3) $$

Here, $x$ follows the empirical distribution of the dataset $p_D = \frac{1}{|\mathcal{D}|} \sum_{x \in \mathcal{D}} \delta_x$ and $y$ defined by (1). The set $\mathcal{M}$ encodes the range of functions reachable by the neural network, $\mathcal{M} = \{N_w : w \in \Theta\}$, with $\Theta$ a set of admissible weights.

A precise characterization of the set $\mathcal{M}$ is not available for a given neural network architecture $N_w$. However, we can relate the trained neural network with statistical estimators constrained to function classes encoding architectural constraints of the neural networks as follows.

**Definition 1.2** (MMSE, E-MMSE, LE-MMSE). Table 1 defines the constraints $\mathcal{M}$ and estimators studied in this work.

<table><caption>Table 1: Different constraint sets correspond to different variants of the MMSE estimator.</caption><thead><tr><th>CONSTRAINT <i>M</i></th><th>ESTIMATOR</th><th>SEE</th></tr></thead><tbody><tr><td>(1) ANY MEASURABLE FUNC.</td><td>MMSE − <i>x̂</i><sub>MMSE</sub></td><td>(2.1)</td></tr><tr><td>(2) (1) + TRANSLATION EQUIV.</td><td>E-MMSE − <i>x̂</i><sub>τ</sub></td><td>(3.5)</td></tr><tr><td>(2) + LOCALITY</td><td>LE-MMSE − <i>x̂</i><sub>τ,LOC</sub></td><td>(3.8)</td></tr></tbody></table>

These function classes capture different architectural constraints of the neural network $N_w$. It is well known that sufficiently deep and wide multi-layer perceptrons (MLPs) can approximate any measurable function (Hornik et al., 1989), hence the MMSE $\hat{x}_{\text{MMSE}}$ is a proxy for unconstrained MLPs. The E-MMSE $\hat{x}_{\tau}$ approximates CNNs without constraints on the receptive field, while the LE-MMSE $\hat{x}_{\tau,\text{loc}}$ approximates the set of functions reachable by CNNs with finite receptive fields.

In what follows, for results involving neural networks, we let $\mathcal{M}$ denote the function class reachable by the considered architecture, and for results involving analytical formulas, we let $\mathcal{M}$ denote the abstract function class in Table 1.

**Related works** In the context of generative diffusion models, (Kamb & Ganguli, 2025; Scarvelis et al., 2025) derive closed-form expressions for the MMSE under additive Gaussian noise (denoising) with architectural constraints such as equivariance and locality, and use them to study memorization and generalization. Locality has also been examined in (Kadkhodaie et al., 2023) for modeling patch distributions. For imaging inverse problems, equivariance has been widely explored in self-supervised learning (Chen et al., 2023; Terris et al., 2024). In the supervised setting for linear inverse problems, it is standard to assume that a trained neural network approximates the MMSE (Adler & Öktem, 2018a). To the best of our knowledge, for general inverse problems, our work is the first to introduce additional functional constraints to more accurately characterize the action of neural networks and to derive closed-form expressions for the resulting constrained MMSE estimators.

We provide an analytic LE-MMSE formula and show that it precisely predicts trained CNN outputs across tasks, datasets, and architectures, see Figure 1 for an overview, Figures 5, 11 and 12 for quantitative comparisons and Appendix C for more results.

**Contributions** In the context of inverse problems outlined above, our contributions and outline are as follows:

1. We review the MMSE estimator (Section 2).

2. We define constrained MMSE variants that capture inductive biases of CNNs and derive their closed-form expressions (Section 3). Remarkably, we find tractable and interpretable formulas that generalize those of (Kamb & Ganguli, 2025) to arbitrary linear inverse problems.

3. We analyze some theoretical properties of these estimators (Section 3). While the MMSE and E-MMSE memorize the training set, the LE-MMSE gives patchworks of training patches and generalizes better. We also discuss the influence of the pre-inverse $B$ on the estimators.

4. We validate our theory on extensive numerical exper-
iments (Section 4) with various architectures, inverse problems and dataset. We find strong agreement between network outputs and the MMSE estimator (PSNR $\gtrsim 25$ dB for all tasks). See Figure 1.

5. We show empirically that different choices of architectures converge to the same solution, which is shown to be very close to our analytical formula (Section 4.2).

6. We characterize where the theory better matches the practice – high-density regions of the *patch distribution*, and explain how CNNs generalize. The influence of additional factors (training set size, noise level, receptive field size) is also investigated (Section 4.3).

In what follows, we let $\bar{x}$ denote a fixed but unknown signal to recover, and $y = A\bar{x} + e$ be a realization of the measurement $y$ defined in (1), with noise $e \sim \mathcal{N}(0, \sigma^2 I_M)$.

## 2. Unconstrained MMSE estimator

It is well known that the MMSE estimator coincides with the conditional mean almost surely (Kay, 1993):

$$ \hat{x}_{\text{MMSE}}(y) = \mathbb{E}[\boldsymbol{x} | \boldsymbol{B} \boldsymbol{y} = \boldsymbol{B} y] = \int x \cdot p(x | \boldsymbol{B} y) \, dx. $$

When the data distribution $p_x$ is replaced by the empirical measure $p_D = \frac{1}{|D|} \sum_{x \in D} \delta_x$, with a finite dataset $D$, the empirical MMSE estimator becomes a kernel regressor (Nadaraya, 1964; Watson, 1964):

**Proposition 2.1 (Closed-form of the MMSE).** *The MMSE estimator in Definition 1.2 writes, for any $y \in \mathbb{R}^M$,*

$$ \hat{x}_{\text{MMSE}}(y) = \sum_{x \in D} x \cdot w(x|y), \quad (4) $$

$$ \text{where } w(x|y) = \frac{\mathcal{N}(By; BAx, \sigma^2 BB^\top)}{\sum_{x' \in D} \mathcal{N}(By; BAx', \sigma^2 BB^\top)}. $$

A proof is provided in B.1. $\mathcal{N}(\cdot; \mu, \Sigma)$ refers to the density of a (possibly degenerate) Gaussian distribution with mean $\mu$ and covariance $\Sigma$, see Definition A.1. The weights $w(x|y)$ is the posterior probability that training sample $x$ explains $By$. The denominator term ensures that they sum to 1.

**Remark 2.2.** When $B = I$, $\hat{x}_{\text{MMSE}}$ is exactly the *posterior mean*: $\hat{x}_{\text{MMSE}}(y) = \mathbb{E}[\boldsymbol{x}|y]$.

**Physics-aware estimator** To elucidate the role of $B$, we remark that the weights have the following explicit form:

$$ w(x|y) \propto \exp \left( -\frac{1}{2\sigma^2} \| \Pi_{\text{Im}(B^\top)} (A(\bar{x} - x) + e) \|^2 \right) $$

Taking $B = A^+$ or $A^\top$, $\text{Im}(B^T) = \text{Im}(A)$ and the inner term is $A(\bar{x} - x) + \Pi_{\text{Im}(A)}e$. The noise $e$ is projected onto the image of $A$, reducing its magnitude. Therefore, for the

<figure><img src="image_4.png" /><figcaption>Figure 2: The MMSE and E-MMSE estimators memorize (yield the nearest neighbor), while the LE-MMSE estimator recombines training patches to give good reconstruction.</figcaption></figure>

MMSE, physics-aware estimators is preferable to physics-agnostic ones, which is consistent with a popular belief and the usual practice. However, beyond the projection effect, the choice of $B$ only plays a little role, for example, the estimator is the same for any $B$ such that $\text{Im}(B^T) = \text{Im}(A)$. We discuss this further in Section 3.5.

**Memorization of the empirical MMSE** The weights in (4) satisfy $w(x|y) \ge 0$ and $\sum_{x \in D} w(x|y) = 1$, hence $\hat{x}_{\text{MMSE}}(y) \in \text{conv}(D)$ — the convex envelop of the dataset $D$, for every $y$. In particular, as $\sigma \to 0$, $\hat{x}_{\text{MMSE}}(y)$ converges to $x \in D$ such that $BAx$ is the closest to $By$ (1-nearest neighbor); ties being averaged. Therefore, the empirical MMSE “memorizes”: it reconstructs by re-weighting training examples and never extrapolates outside $\text{conv}(D)$, see Figure 2. This aligns with recent findings related to memorization in generative models (Kamb & Ganguli, 2025; Bertrand et al., 2025; Scarvelis et al., 2025).

## 3. Constrained MMSE estimators

We first provide details on the constraint classes and then describe our main theoretical results: analytical formulas and properties of constrained MMSE estimators.

### 3.1. Functional constraints: equivariance and locality

Motivated by the success of CNNs in inverse problems, we study function classes $\mathcal{M}$ that encode two core inductive biases – equivariance to geometric transformations ( Cohen & Welling, 2016; Veeling et al., 2018) and locality via bounded receptive field (Kadkhodaie et al., 2023). We start with formal definition of these concepts (see Appendix A.1 for precise definitions). Let $x \in \mathbb{R}^N$ represent an image defined on a discrete 2D grid of size $H \times W = N$.

**Definition 3.1 (Translation equivariant).** Let $\mathcal{T} = \mathbb{Z}_H \times \mathbb{Z}_W$ be the group of 2D cyclic translations, where $\mathbb{Z}_H = \{0, 1, \dots, H-1\}$ and $\mathbb{Z}_W = \{0, 1, \dots, W-1\}$. For each $g = (g_h, g_w) \in \mathcal{T}$, its action is represented by a permutation matrix $T_g \in \mathbb{R}^{N \times N}$, where $T_g x$ shifts the image $x \in \mathbb{R}^N$ by $g_h$ pixels vertically and $g_w$ pixels horizontally with periodic boundary conditions. A measurable map $\phi: \mathbb{R}^N \to \mathbb{R}^N$ is
said to be translation equivariant if it satisfies $\phi(T_g x) = T_g \phi(x)$ for all $x \in \mathbb{R}^N$ and $g \in \mathcal{T}$. We denote by $\mathcal{M}_\mathcal{T}$ the set of all such maps.

We now add locality constraint. For each pixel $n$, consider the patch extractor $\Pi_n : x \in \mathbb{R}^N \mapsto \Pi_n x = x[\omega_n] \in \mathbb{R}^P$ that extracts a square patch $x[\omega_n]$ of size $\sqrt{P} \times \sqrt{P}$ centered at pixel $n$, and $\omega_n$ is the set of pixel indices in the patch with circular boundary conditions.

**Definition 3.2 (Local and translation equivariant).** A measurable map $\phi: \mathbb{R}^N \to \mathbb{R}^N$ is local and translation equivariant if there exists a measurable map $f: \mathbb{R}^P \to \mathbb{R}$ such that, for all $x \in \mathbb{R}^N$ and any pixel $n$, the output of $\phi$ at pixel $n$ denoted by $\phi(x)[n]$, depends only on the local patch $x[\omega_n]$, that is $\phi(x)[n] = f(\Pi_n x)$. We denote by $\mathcal{M}_{\mathcal{T},\text{loc}}$ the set of all such functions.

The class $\mathcal{M}_{\mathcal{T},\text{loc}}$ captures standard CNNs with small kernels and weight sharing or patch-based models, e.g. MLPs acting on patches (Khorashadizadeh et al., 2025).

## 3.2. Preliminary remarks

**Constraints and projection** The constrained MMSE has a simple geometric interpretation, see proof in A.6.

**Proposition 3.3.** For a closed set $\mathcal{M}$, the $\mathcal{M}$-constrained MMSE estimator in Definition 1.1 is the orthogonal projection in $\mathbb{L}^2$ of the posterior mean $\mathbb{E}[\mathbf{x}|\mathbf{y}]$ onto the subspace of random vectors of the form $\mathcal{X} = \{\phi(\mathbf{By}) : \phi \in \mathcal{M}\}$, that is $\hat{x}_{\mathcal{M}}(\mathbf{y}) = \Pi_{\mathcal{X}}(\mathbb{E}[\mathbf{x}|\mathbf{y}])$.

**Architecture equivariance and data augmentation** A standard approach to obtain equivariant estimators is by using data augmentation: the set $\mathcal{D}$ is replaced by $\mathcal{T}(\mathcal{D}) = \{T_g x : x \in \mathcal{D}, g \in \mathcal{T}\}$ at training time. This promotes reconstruction equivariance (Chen et al., 2021) of $\hat{x}_{\mathcal{M}}$:

$$ \hat{x}_{\mathcal{M}}(AT_g \bar{x} + e) = T_g \hat{x}_{\mathcal{M}}(A\bar{x} + T_g^{-1}e). \quad (5) $$

The data-augmented MMSE estimators, are defined below.

**Definition 3.4 (Data-augmented MMSE estimator).** The data-augmented MMSE estimator $\hat{x}_{\mathcal{M}}^{\text{aug}}$ is defined as the MMSE estimator $\hat{x}_{\mathcal{M}}$ with respect to the empirical measure on $\mathcal{T}(\mathcal{D})$.

Reconstruction equivariance is often a desired property, but it differs from the structural equivariance of $\phi$ (Chen et al., 2020; Nordenfors & Flinth, 2025), as illustrated in Figure 3.

## 3.3. Translation Equivariant MMSE estimator

We start by the closed-form of the E-MMSE estimator $\hat{x}_{\mathcal{T}}$ and its properties. See Appendix B.2 for proofs.

**Theorem 3.5 (Closed-form of E-MMSE).** The E-MMSE

estimator writes, for any $y \in \mathbb{R}^M$

$$ \hat{x}_{\mathcal{T}}(y) = \sum_{x \in \mathcal{D}, g \in \mathcal{T}} T_g x \cdot w_g(x|y), \quad (6) $$

where $w_g(x|y) \propto \mathcal{N}(T_g^{-1} By; BAx, \sigma^2 BB^\top)$.

**Remark 3.6.** This result holds true for any group $\mathcal{G}$ of orthogonal transformations, not just translations.

The normalization constant of the weights $w_g(x|y)$ ensures that they sum to 1 over all $x \in \mathcal{D}$ and $g \in \mathcal{T}$, we ignore it here for brevity. The E-MMSE estimator $\hat{x}_{\mathcal{T}}$ is a weighted average of the augmented dataset $\mathcal{T}(\mathcal{D})$, which is similar to $\hat{x}_{\text{MMSE}}^{\text{aug}}$ in Definition 3.4. In particular, for denoising ($A = B = I$), they coincide. However, depending on the forward operator $A$ and the pre-inverse $B$, it can differ significantly from $\hat{x}_{\text{MMSE}}^{\text{aug}}$. Some properties of the E-MMSE $\hat{x}_{\mathcal{T}}$ are presented in Corollary 3.7 and illustrated in Figure 3.

**Corollary 3.7.** The E-MMSE estimator $\hat{x}_{\mathcal{T}}$ in Theorem 3.5 satisfies the following properties:

* The E-MMSE estimator $\hat{x}_{\mathcal{T}}$ memorizes the augmented dataset: reconstructed images live in $\text{conv}(\mathcal{T}(\mathcal{D}))$.

* Data-augmentation and architecture equivariance are identical ($\hat{x}_{\mathcal{T}} = \hat{x}_{\text{MMSE}}^{\text{aug}}$) if $A$ and $B$ are circular convolution, with $B$ invertible. In this case, the E-MMSE is reconstruction equivariant (satisfies (5)), physics-agnostic and physics-aware solvers are identical.

* This is a near equivalence: $w_g(x|y) = w(T_g x|y)$ only holds if $A$ and $B$ are invertible circular convolutions.

<figure><img src="image_5.png" /><figcaption>Figure 3: Architectural equivariance does not always guarantee reconstruction equivariance (5). Input $y$ in 2nd row is shifted. The E-MMSE estimator $\hat{x}_{\mathcal{T}}$ is reconstruction equivariant for deconvolution but not inpainting. The LE-MMSE estimator $\hat{x}_{\mathcal{T},\text{loc}}$ is reconstruction equivariant for deconvolution and shows reduced sensitivity to shifts for inpainting.</figcaption></figure>

## 3.4. Locality and Translation Equivariance

We now analyze the effect of local receptive fields. See Appendix B.3 for proofs.

**Theorem 3.8 (Closed-form of LE-MMSE).** Suppose that the $N$ matrices $Q_n \stackrel{\text{def}}{=} \Pi_n B \in \mathbb{R}^{P \times M}$ have the same rank $r > 0$ for any $n^1$. The LE-MMSE estimator $\hat{x}_{\mathcal{T},\text{loc}}$ admits,

<footnote n="1">¹This assumption can be relaxed: a general formula by rank stratification is provided in Proposition B.5 for completeness.</footnote>
for any $y \in \mathbb{R}^M$ the following expression, defined pixel-wise for each pixel $n'$:

$$ \hat{x}_{\mathcal{T}, \text{loc}}(y)[n'] = \sum_{x \in \mathcal{D}} \sum_{n=1}^{N} x[n] \cdot w_{n',n}(x|y), \quad (7) $$

$$ \text{where } w_{n',n}(x|y) \propto \mathcal{N}(Q_{n'}y; Q_n Ax, \sigma^2 Q_n Q_n^\top). $$

The value at pixel $n'$ of the LE-MMSE estimator is a weighted average of all the pixels in the training images. We show that the resulting estimator can be interpreted as a patchwork of training patches in Figures 23 and 24. This recombination can produce images outside $\text{conv}(\mathcal{D})$ or $\text{conv}(\mathcal{T}(\mathcal{D}))$ (contrary to the MMSE or E-MMSE); see Figure 2 and Section 4. The weights $w_{n',n}(x|y)$ can be seen as the posterior probability of the pixel $x[n]$ given the patch $(By)[\omega_{n'}] = Q_{n'}y$.

Recall that $y = A\bar{x} + e$, with $\bar{x}$ the signal to recover and let $\Delta_{n',n}(\bar{x}, x) \stackrel{\text{def}}{=} (BA\bar{x})[\omega_{n'}] - (BAx)[\omega_n]$, the weights can be rewritten as $w_{n',n}(x|y) \propto \exp(-\frac{\eta^2}{2\sigma^2})$ with:

$$ \eta \stackrel{\text{def}}{=} \|Q_n^+ \Delta_{n',n}(\bar{x}, x) + Q_n^+ Q_{n'} e\|. \quad (8) $$

Hence, the weights measure the similarity between local patches of measured-then-reconstructed images with the metric $Q_n^+$. It is perturbed by the noise term $Q_n^+ Q_{n'} e$, which is responsible for the estimator's variance. Some properties of the LE-MMSE $\hat{x}_{\mathcal{T}, \text{loc}}$ are given in Corollary 3.9. See Appendix B.3 for proofs and Appendix D.4 for further discussions and illustrations.

**Corollary 3.9.** The LE-MMSE estimator $\hat{x}_{\mathcal{T}, \text{loc}}$ in Theorem 3.8 satisfies the following properties:

* The LE-MMSE estimator $\hat{x}_{\mathcal{T}, \text{loc}}$ is a patchwork (see Appendix D.5) of training patches. As $\sigma \to 0$, each pixel is set to the central pixel of the best matching patch from $\mathcal{D}$.

* It is not a posterior mean: if $A$ is the identity (denoising), for a generic database $\mathcal{D}$, there exists no prior distribution of $x$ such that $\hat{x}_{\mathcal{T}, \text{loc}}(y) = E[x|y=y]$.

### 3.5. Physics-agnostic and physics-aware estimators

The expectation of $\eta^2$ in (8) is given by

$$ \mathbb{E}_e[\eta^2] = \|Q_n^+ \Delta_{n',n}(\bar{x}, x)\|^2 + \sigma^2 \cdot \text{Tr}(\text{Cov}_{n',n}) $$

with $\text{Cov}_{n',n} = Q_n^+ Q_{n'} Q_n^\top Q_n^+$. This expression reveals that the pre-inverse $B$ plays two distinct roles:

* Signal discrimination: the term $Q_n^+ \Delta_{n',n}(\bar{x}, x)$ measures the amplification/reduction of difference between dissimilar/similar patches for the pre-inverse $B$. A good choice of $B$ should improve this discrimination.

* Noise robustness: the second term $\text{Tr}(\text{Cov}_{n',n})$ measures noise amplification/reduction for $B$. A good choice should reduce this amplification.

Ideally, one would like to choose $B$ to optimize both criteria, but unfortunately they can be conflicting: preserving the signal $BAx \approx x$ may amplify the noise, in relation to the classical bias-variance decomposition. Most importantly this effect is operator dependent, as we now discuss.

For inpainting, choosing a physics-aware estimator with $B = A^+$ is an effective way to reduce the noise sensitivity while keeping signal's information. Indeed, with this choice, the noise within the masked region is eliminated, while the signal is preserved in the unmasked region. For deconvolution, however, the situation is different. Choosing $B = A^+$ can result in significant noise amplification (i.e. $\text{Tr}(\text{Cov}_{n',n})$ is large), which is not counter-balanced by a better signal's discrimination. In this setting, a regularized inverse, balancing both aspects should likely be preferred. This effect is illustrated in Figure 4. When designing reconstruction algorithms, this trade-off between data-consistency and noise amplification must be carefully considered, in relation to the inverse problem at hand.

<figure><img src="image_6.png" /><figcaption>Figure 4: Physics-aware ($B=A^+$, bottom) has lower variance for inpainting (left), while physics-agnostic ($B=I$, top) has lower variance for deconvolution (right). Mean and pixel-wise variance are computed w.r.t 50 noise realizations.</figcaption></figure>

## 4. Numerical Experiments

While the MMSE and E-MMSE estimators are interesting from a theoretical perspective, the LE-MMSE estimator is more relevant to practical neural network architectures used in imaging inverse problems. We validate our theoretical findings about the LE-MMSE on three representative inverse problems: denoising, inpainting, and deconvolution.

### 4.1. Experimental setup

We implement three different local and translation equivariant neural network architectures: UNet2D (Ronneberger et al., 2015), ResNet (He et al., 2016) and PatchMLP (an MLP acting on patches).

Implementing the analytical formula in Theorem 3.8 requires computing a huge amount of pairwise distances, which is computationally demanding. We therefore restrict our experiments to 32 × 32 color images and datasets of
$10^4$ images. To accumulate Gaussian weights robustly and stably across batches, we use a careful online log-sum-exp implementation of the weighted averages.

For inpainting we use a square mask of size $15 \times 15$ at the center. For deconvolution, we use an isotropic Gaussian kernel of standard deviation of 1.0. The noise level $\sigma$ is varying uniformly between 0 and 1.0 during training. Full details are described in Appendix C and the implementation is available at https://github.com/mh-nguyen712/analytical_mmse.

## 4.2. Analytical formula and neural networks

**Case-by-case neural outputs prediction** We verify that trained networks approximate the LE-MMSE estimator, by comparing their outputs to the formula in Theorem 3.8. The PSNR between the trained UNet2D, the formula and the ground truth are reported in Figure 5. The reconstruction quality (orange and blue curves) degrades as noise level increases, which is expected as noise makes the problem harder. However, the PSNR between neural networks and the analytical formula (green curves) remains high ($\gtrsim$ 25 dB) across all noise levels, indicating a strong alignment between the two. Similar behavior is observed in Table 2 for other architectures and datasets. We provide extensive

<figure><img src="image_7.png" /><figcaption>Figure 5: The green curves reveal a strong agreement (PSNR $\gtrsim$ 25 dB) between the trained UNet2D and the analytical formula for all inverse problems. Median and interquartile range (IQR) using 50 images per $\sigma$, $P = 5 \times 5$ and $B=1$. See Figures 11 and 12 for other architectures.</figcaption></figure>

numerical results on various settings (architectures, datasets, tasks, physics-agnostic/aware, patch sizes) in Appendix D.1. These results confirm that trained neural networks closely approximate the analytical LE-MMSE estimator. It is remarkable that different architectures consistently yield a similar output, as shown in Figures 1 and 13 to 16.

**Low-density regions** A noticeable train-test gap remains in Figure 5 and Table 2, most pronounced at low noise level. We attribute this to the measurement distribution $p(y)$ induced by the empirical distribution. It is a Gaussian mixture, centered at the training measurements $Ax: p(y) = \frac{1}{|D|} \sum_{x \in D} N(y; Ax, \sigma^2 I_M)$. Its density is high near the

centers $Ax, x \in D$, or for large noise levels $\sigma$, creating overlap between Gaussians. It becomes negligible far from the centers in low noise regimes.

<figure><img src="image_8.png" /><figcaption>Figure 6: Higher density yield better alignment. We select test images from FFHQ-32, compute $-\log p(y)$ and plot it against the PSNR between the outputs of the LE-MMSE formula and a trained UNet2D, $\sigma = 0.05$ and $P = 5 \times 5$.</figcaption></figure>

Both LE-MMSE and the networks are optimized with respect to the measurement distribution $p(y)$. The networks and the theoretical formula best align in high-density regions of $p(y)$. This phenomenon is illustrated in Figure 6, where we show that the alignment degrades as $-\log p(y)$ increases (or equivalently as $p(y)$ decreases). Notice that for large $\sigma$, the density is higher as the Gaussians overlap more, resulting in a better agreement.

**Out-of-distribution: how do networks generalize?** Neural networks often generalize well to out-of-distribution (OOD) data, but understanding this behavior remains an open question. Interestingly, our analytical formulas provide insights to explain this phenomenon. We consider here images $\bar{x}$ that lie in another dataset $D'$ disjoint from the training set $D$. In Figure 7, we compare the output of UNet2D

<figure><img src="image_9.png" /><figcaption>Figure 7: Our theory can predict the neural network outputs for out-of-distribution data: both UNet2D and the LE-MMSE formula are trained (or computed) on $D = \text{FFHQ-32}$. They are tested on $D' = \text{CIFAR10}$.</figcaption></figure>

and the LE-MMSE formula. For large $\sigma$, we observe very close outputs between the neural network and the analytical formula. For small $\sigma$, as the OOD images lie in low-density regions of $p(y)$, the alignment degrades as expected from the previous discussion. However, even in this regime, the PSNR between the two outputs remains reasonably high (see Appendix D.6 for quantitative results). This suggests that the network approximates the LE-MMSE estimator
<figure><img src="image_10.png" /><figcaption>Table 2: Theorem 3.8 is verified across architectures, tasks, and datasets. We show the median of the PSNR between neural networks and the analytical formula with $P = 5 \times 5$ and $B=I$. The PSNR on FashionMNIST is consistently higher ($2 \sim 3$ dB) than on FFHQ and CIFAR10, which we attribute to the lower complexity of FashionMNIST images. A lower PSNR on the test set in low-noise regimes is explained by the low-density regions.</figcaption></figure>

even on OOD data: generalization of the neural networks can be understood through the lens of the LE-MMSE estimator – it leverages local patches from the training set to reconstruct unseen images.

### 4.3. Hyperparameter influence

**Patch size** The patch size $P$ in the LE-MMSE estimator plays a crucial role and should be chosen carefully depending on the inverse problem and noise level $\sigma$. For large noise levels, the inverse problem requires a higher degree of regularization, which can be achieved by using larger patches. On the contrary, small patches are preferred for low noise levels to capture fine details. This effect is illustrated in Figure 19, where $5 \times 5$ patches should be preferred for the test set at low noise levels, while $11 \times 11$ patches perform better at high noise levels.

On the other hand, the match between the analytical formula and the neural networks output is higher for small patches, as they depend on the density of the measurement patches distribution $p(y[\omega])$, see Figure 18. There are $N \cdot |\mathcal{D}|$ patches, but the space dimension is $P$. Consequently, for a fixed dataset size $|\mathcal{D}|$, increasing $P$ lowers the density of $p(y[\omega])$ significantly. A similar conclusion was already reached in denoising diffusion models, where (Kamb & Ganguli, 2025) proposed to vary the patch sizes from large to small across noise levels to guarantee a better match.

**Dataset size** We further investigate the influence of the dataset size $|\mathcal{D}|$ on the alignment between neural networks and the LE-MMSE formula. Specifically, we vary $|\mathcal{D}|$ from $10^3$ to $5 \times 10^4$ images on FFHQ-32 and trained UNet2D with $P = 5 \times 5$ and $B = I$. As shown in Figure 26, our theory holds for all dataset sizes, with small variations in PSNR between the neural network and the formula.

### 4.4. Smoothed LE-MMSE and spectral bias

Deep neural networks also exhibit a spectral bias (Xu et al., 2019; Rahaman et al., 2019), i.e. a preference for low-frequency functions, which yields smoother reconstructions in practice. Describing this spectral bias with functional constraints is missing from our derivation. We model this effect by considering a randomized smoothing (Duchi et al., 2012) variant: for a smoothing parameter $\epsilon > 0$ and $z \sim \mathcal{N}(0, I_M)$, define $\hat{x}_{T,\text{loc}}^{\text{smooth}}(y) \stackrel{\text{def}}{=} \mathbb{E}[\hat{x}_T(y + \epsilon \cdot z)]$. This averages the estimates over random perturbations, acting as a nonlinear low-pass filter that attenuates high-frequency variability in the output. It also mirrors stan-

<figure><img src="image_11.png" /><figcaption>Figure 8: Smoothed LE-MMSE ($\epsilon = 0.05$) better matches neural network and improves reconstruction quality.</figcaption></figure>

dard neural network training, where different noise realizations are seen across epochs. In Figure 8, we compare the smoothed estimator $\hat{x}_{T,\text{loc}}^{\text{smooth}}$ with the original $\hat{x}_{T,\text{loc}}$, both computed with $P = 5 \times 5$, $\sigma = 0.05$ and $B = I$ on FFHQ-32 images. We observe that the smoothed LE-MMSE better matches the neural network output ($1 \sim 3$ dB) and improves reconstruction quality ($\sim 1$dB). This suggests that we can better capture neural network behavior by incorporating carefully designed smoothness constraints.

### 4.5. Validation on higher resolution images

To further validate our theoretical findings, we conduct experiments on FFHQ images of size $64 \times 64$ with $10^4$ samples. We train a UNet2D architecture with patch size $P = 11 \times 11$
<figure><img src="image_12.png" /><figcaption>Figure 9: Our theory is validated at higher resolution: comparison of UNet2D output and the analytical formula on FFHQ-64.</figcaption></figure>

for denoising, inpainting, and deconvolution tasks. We observe that the trained neural network closely approximates the analytical LE-MMSE estimator, further strengthening our conclusions (see Figures 9 and 27). Details about architecture, training and additional results are provided in Appendix D.8. Note that computing the analytical formula at this resolution is computationally intensive, limiting the extent of experiments compared to the 32 × 32 case.

# 5. Discussions and conclusion

## 5.1. Neural networks: compress and rearticulate data

Although our analysis focuses on the local and translation equivariant MMSE estimator, the underlying conclusions extend to far more general learning settings. Neural networks compress and rearticulate the training data in an efficient manner. This viewpoint is consistent with recent findings in deep networks (Arpit et al., 2017), including large language models (Biderman et al., 2023) and image generation models (Somepalli et al., 2023), where new samples are synthesized by recombining and transforming elements implicitly stored from the training distribution.

In our experiments, neural networks with roughly $4 \times 10^6$ parameters are trained on datasets containing $10^4$ images of resolution 3×32×32, with approximately $3.1 \times 10^7$ pixels. The networks compress this information into a comparatively small number of parameters (roughly 13% of the number of pixels) while still retaining the ability to reconstruct images by rearticulating training set patches. Furthermore, neural network inference is substantially faster than evaluating the analytical estimator: ~4ms versus ~0.7s for denoising, resulting in ~175× speedup. This gain is even more pronounced for physics-aware estimator or other operators: the formula inverses the operators $Q_n$ for each patch while the neural networks learn to approximate it from training data. This demonstrates the practical advantage of neural networks in approximating complex estimators such as the LE-MMSE: they not only compress the data but also enable efficient, near-instantaneous evaluation. Originally considered

as black boxes, our analysis provides a theoretical foundation for understanding how neural networks achieve this remarkable feature.

## 5.2. Perspective and future works

**Reconstruction guarantee and stability** Neural networks often lack theoretical guarantees regarding reconstruction accuracy and stability. By modeling CNNs as LE-MMSE approximations, our closed-form expressions enable the theoretical analysis of these properties. This framework moves beyond purely empirical validation, offering a rigorous basis for understanding the strengths and limitations of neural inverse problem solvers.

**Extensions beyond current assumptions** While we assumed in this work additive Gaussian noise and local, translation equivariant functional classes, the framework can be broadened. For example, vision transformers (Dosovitskiy, 2020) with self-attention capture long-range dependencies, but it is permutation equivariant (Xu et al., 2024). The Gaussian assumption can be relaxed (e.g., exponential-family) by using the appropriate likelihood $p_{y|x}(y|x)$. Finally, beyond MSE, practical training often uses $\ell_1$ (yielding posterior median), perceptual, or mixed losses. Analyzing the MMSE analogs under these alternatives is an interesting future venue.

## 5.3. Conclusion

We presented an analytic framework for CNN-based inverse problem solvers by considering MMSE with functional constraints. We derived closed-form expressions that make explicit the impact the forward operator $A$ and a pre-inverse $B$ on the reconstruction. The theory separates memorization and generalization behaviors: while the MMSE and E-MMSE memorize, the LE-MMSE forms patchwork reconstructions and closely matches trained CNN outputs across settings. This perspective turns black-box inverse problem CNNs into predictable estimators, enabling future
theoretical analyzes of generalization, stability, and princi-
pled choices of pre-inverses.

Acknowledgement

The authors acknowledge a support from the ANR Micro-
Blind (ANR-620 21-CE48-0008) and from the ANR
CLEAR-Microscopy (ANR-25-CE45-3780). This work was
performed using HPC resources from GENCI-IDRIS (Grant
AD011012210)

**References**

Adler, J. and Öktem, O. Deep bayesian inversion. *arXiv preprint arXiv:1811.05910*, 2018a.

Adler, J. and Öktem, O. Learned primal-dual reconstruction. *IEEE transactions on medical imaging*, 37(6):1322–1332, 2018b.

Arpit, D., Jastrzębski, S., Ballas, N., Krueger, D., Bengio,
E., Kanwal, M. S., Maharaj, T., Fischer, A., Courville,
A., Bengio, Y., and Lacoste-Julien, S. A closer look
at memorization in deep networks. In *Proceedings of
the 34th International Conference on Machine Learning*,
volume 70 of *Proceedings of Machine Learning Research*,
pp. 233–242. PMLR, 06–11 Aug 2017.

Bertrand, Q., Gagneux, A., Massias, M., and Emonet, R. On
the closed-form of flow matching: Generalization does
not arise from target stochasticity. *Advances in neural
information processing systems*, 2025.

Biderman, S., Prashanth, S., Sutawika, L., Schoelkopf, H., Anthony, Q., Purohit, S., and Raff, E. Emergent and predictable memorization in large language models. In Oh, A., Naumann, T., Globerson, A., Saenko, K., Hardt, M., and Levine, S. (eds.), *Advances in Neural Information Processing Systems*, volume 36, pp. 28072–28090. Curran Associates, Inc., 2023. URL https://openreview.net/pdf?id=Iq0DvhB4Kf.

Celledoni, E., Ehrhardt, M. J., Etmann, C., Owren, B., Schönlieb, C.-B., and Sherry, F. Equivariant neural networks for inverse problems. *Inverse Problems*, 37(8): 085006, July 2021. ISSN 1361-6420. doi: 10.1088/1361-6420/ac104f. URL http://dx.doi.org/10.1088/1361-6420/ac104f.

Chen, D., Tachella, J., and Davies, M. E. Equivariant imaging: Learning beyond the range space. In *Proceedings of the IEEE/CVF International Conference on Computer Vision*, pp. 4379–4388, 2021.

Chen, D., Davies, M., Ehrhardt, M. J., Schönlieb, C.-B., Sherry, F., and Tachella, J. Imaging with equivariant deep learning: From unrolled network design to fully unsupervised learning. *IEEE Signal Processing Magazine*, 40(1): 134–147, 2023.

Chen, S., Dobriban, E., and Lee, J. H. A group-theoretic framework for data augmentation. *Journal of Machine Learning Research*, 21(245):1–71, 2020.

Cohen, T. and Welling, M. Group equivariant convolutional networks. In *International conference on machine learning*, pp. 2990–2999. PMLR, 2016.
Dosovitskiy, A. An image is worth 16x16 words: Transformers for image recognition at scale. *arXiv preprint arXiv:2010.11929*, 2020.

Duchi, J. C., Bartlett, P. L., and Wainwright, M. J. Randomized smoothing for stochastic optimization. *SIAM Journal on Optimization*, 22(2):674–701, 2012.

Federer, H. *Geometric Measure Theory*. Classics in Mathematics. Springer-Verlag Berlin Heidelberg, Berlin, Heidelberg, 1 edition, 1996. ISBN 978-3-540-60656-7. doi: 10.1007/978-3-642-62010-2. Originally published as volume 153 in the series: Grundlehren der mathematischen Wissenschaften.

He, K., Zhang, X., Ren, S., and Sun, J. Deep residual learning for image recognition. In *Proceedings of the IEEE conference on computer vision and pattern recognition*, pp. 770–778, 2016.

Hornik, K., Stinchcombe, M., and White, H. Multilayer feedforward networks are universal approximators. *Neural networks*, 2(5):359–366, 1989.

Jin, K. H., McCann, M. T., Froustey, E., and Unser, M. Deep convolutional neural network for inverse problems in imaging. *IEEE transactions on image processing*, 26 (9):4509–4522, 2017.

Kadkhodaie, Z., Guth, F., Mallat, S., and Simoncelli, E. P. Learning multi-scale local conditional probability models of images. In *The Eleventh International Conference on Learning Representations*, 2023. URL https://openreview.net/forum?id=VZX2I_VVJKH.

Kamb, M. and Ganguli, S. An analytic theory of creativity in convolutional diffusion models. In *Forty-second International Conference on Machine Learning*, 2025. URL https://openreview.net/forum?id=ilpL2qAC1a.

Karras, T., Laine, S., and Aila, T. A style-based generator architecture for generative adversarial networks. In *Proceedings of the IEEE/CVF conference on computer vision and pattern recognition*, pp. 4401–4410, 2019.

Kay, S. M. *Fundamentals of statistical signal processing: estimation theory*. Prentice-Hall, Inc., 1993.

Khorashadizadeh, A., Debarnot, V., Liu, T., and Dokmanić, I. Glimpse: Generalized locality for scalable and robust ct. *IEEE Transactions on Medical Imaging*, 2025. URL https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=11018464.

Kim, J., Lee, J. K., and Lee, K. M. Accurate image super-resolution using very deep convolutional networks. In *Proceedings of the IEEE conference on computer vision and pattern recognition*, pp. 1646–1654, 2016.

Kingma, D. P. and Ba, J. Adam: A method for stochastic optimization, 2014. URL https://arxiv.org/abs/1412.6980.

Krizhevsky, A., Hinton, G., et al. Learning multiple layers of features from tiny images. (2009), 2009.

McCann, M. T., Jin, K. H., and Unser, M. Convolutional neural networks for inverse problems in imaging: A review. *IEEE Signal Processing Magazine*, 34(6):85–95, 2017.

Mityagin, B. The zero set of a real analytic function. *Mathematical Notes*, 107(3):529–530, 2020.

Muckley, M. J., Riemenschneider, B., Radmanesh, A., Kim, S., Jeong, G., Ko, J., Jun, Y., Shin, H., Hwang, D., Mostapha, M., et al. Results of the 2020 fastmri challenge for machine learning mr image reconstruction. *IEEE transactions on medical imaging*, 40(9):2306–2317, 2021.

Nadaraya, E. A. On estimating regression. *Theory of Probability & Its Applications*, 9(1):141–142, 1964.

Nordenfors, O. and Flinth, A. Data augmentation and regularization for learning group equivariance. In *15th International Conference on Sampling Theory and Applications*, 2025.

Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al. Pytorch: An imperative style, high-performance deep learning library. *Advances in neural information processing systems*, 32, 2019.

Rahaman, N., Baratin, A., Arpit, D., Draxler, F., Lin, M., Hamprecht, F., Bengio, Y., and Courville, A. On the spectral bias of neural networks. In *International conference on machine learning*, pp. 5301–5310. PMLR, 2019.

Rao, C. R. *Linear statistical inference and its applications*, volume 2. Wiley New York, 1973.

Ronneberger, O., Fischer, P., and Brox, T. U-net: Convolutional networks for biomedical image segmentation. In *International Conference on Medical image computing and computer-assisted intervention*, pp. 234–241. Springer, 2015.

Scarvelis, C., de Ocáriz Borde, H. S., and Solomon, J. Closed-form diffusion models. *Transactions on Machine Learning Research*, 2025.

Schlemper, J., Caballero, J., Hajnal, J. V., Price, A. N., and Rueckert, D. A deep cascade of convolutional neural networks for dynamic mr image reconstruction. *IEEE transactions on Medical Imaging*, 37(2):491–503, 2017.
Somepalli, G., Singla, V., Goldblum, M., Geiping, J., and Goldstein, T. Diffusion art or digital forgery? investigating data replication in diffusion models. In *Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)*, pp. 6048–6058, June 2023.

Terris, M., Moreau, T., Pustelnik, N., and Tachella, J. Equivariant plug-and-play image reconstruction. In *Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition*, pp. 25255–25264, 2024.

Veeling, B. S., Linmans, J., Winkens, J., Cohen, T., and Welling, M. Rotation equivariant cnns for digital pathology. In *International Conference on Medical image computing and computer-assisted intervention*, pp. 210–218. Springer, 2018.

Wang, G., Ye, J. C., and De Man, B. Deep learning for tomographic image reconstruction. *Nature machine intelligence*, 2(12):737–748, 2020.

Watson, G. S. Smooth regression analysis. *Sankhyā: The Indian Journal of Statistics, Series A*, pp. 359–372, 1964.

Würfl, T., Ghesu, F. C., Christlein, V., and Maier, A. Deep learning computed tomography. In *International conference on medical image computing and computer-assisted intervention*, pp. 432–440. Springer, 2016.

Xiao, H., Rasul, K., and Vollgraf, R. Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms. *arXiv preprint arXiv:1708.07747*, 2017.

Xu, H., Xiang, L., Ye, H., Yao, D., Chu, P., and Li, B. Permutation equivariance of transformers and its applications. In *Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition*, pp. 5987–5996, 2024.

Xu, Z.-Q. J., Zhang, Y., Luo, T., Xiao, Y., and Ma, Z. Frequency principle: Fourier analysis sheds light on deep neural networks. *arXiv preprint arXiv:1901.06523*, 2019.

Zhou, B., Chen, X., Zhou, S. K., Duncan, J. S., and Liu, C. Dudodr-net: Dual-domain data consistent recurrent network for simultaneous sparse view and metal artifact reduction in computed tomography. *Medical Image Analysis*, 75:102289, 2022.

Zhu, B., Liu, J. Z., Cauley, S. F., Rosen, B. R., and Rosen, M. S. Image reconstruction by domain-transform manifold learning. *Nature*, 555(7697):487–492, 2018.
# A. Mathematical formalism

## A.1. Preliminaries

**Notation conventions** In what follows, we use the following notation:

* Sets are indicated by calligraphic letters, e.g. $D, M, T, G$.

* Random vectors are denoted by bold lowercase letters, e.g. $x, y, z$.

* $A \in \mathbb{R}^{M \times N}$ denotes a linear operator from $\mathbb{R}^N$ to $\mathbb{R}^M$ (e.g. convolution, inpainting).

* $\mathcal{N}(x; \mu, \Sigma)$ denotes the probability density function (PDF) of a Gaussian distribution with mean $\mu$ and covariance $\Sigma$ evaluated at $x$. The covariance matrix $\Sigma$ can be non-singular or singular (see Definition A.1).

* The n-th coordinate of a vector $x \in \mathbb{R}^N$ is denoted either $x_n$ or $x[n]$. Similarly, the coordinates of a multivalued function $\phi: \mathbb{R}^N \to \mathbb{R}^N$ can be denoted either $\phi_n(x)$ or $\phi(x)[n]$.

* The empirical data distribution $p_D$ is defined by

$$p_D = \frac{1}{|D|} \sum_{x \in D} \delta_x, \qquad (9)$$

where $D$ is the training dataset of finite size $|D|$.

**Degenerate Gaussian distribution** We can define the multivariate Gaussian distribution for positive semi-definite covariance matrices (Rao, 1973). Below, we recall a simple definition of the positive semi-definite covariance multivariate Gaussian, sometimes called the degenerate or singular multivariate Gaussian.

**Definition A.1 (Degenerate Gaussian distribution).** A random vector $z \in \mathbb{R}^N$ has a Gaussian distribution with mean $\mu \in \mathbb{R}^N$ and covariance matrix $\Sigma \in \mathbb{R}^{N \times N}$, $\Sigma \succeq 0$ if its probability density function (PDF) is given by:

$$\mathcal{N}(z; \mu, \Sigma) = \frac{1}{(2\pi)^{r/2} \sqrt{|\Sigma|_+}} \exp \left( -\frac{1}{2} (z - \mu)^{\top} \Sigma^{+} (z - \mu) \right) \mathbb{1}_{\text{supp}(z)}(z), \qquad (10)$$

on the support $z \in \text{supp}(z) = \mu + \text{Im}(\Sigma)$ and zero elsewhere. In this equation, $r = \text{rank}(\Sigma)$, $\Sigma^+$ denotes the Moore-Penrose pseudo-inverse of $\Sigma$ and $|\Sigma|_+$ is the pseudo-determinant of $\Sigma$ (product of non-zero eigenvalues).

The support supp($z$) of the distribution is the $r$-dimensional affine subspace of $\mathbb{R}^N$, $r$ is the rank of $\Sigma$. The density is with respect to the Lebesgue measure restricted to supp($z$), constructed via the pushforward measure of the standard Lebesgue measure on $\mathbb{R}^r$ by the affine map $v \mapsto \mu + \Sigma_r v$, where $\Sigma = \Sigma_r \Sigma_r^\top$. This measure also coincides (up to a normalization constant, equal to the pseudo-determinant of $\Sigma$) with the $r$-dimensional Hausdorff measure. For completeness, we provide a derivation of the density of $z$ with respect to the $r$-dimensional Hausdorff measure $\mathcal{H}^r$ on the affine support supp($z$).

Suppose $\Sigma_r \in \mathbb{R}^{n \times r}$ satisfies $\Sigma = \Sigma_r \Sigma_r^\top$ and $\text{rank}(\Sigma) = r$. Define the random vector $z = \mu + \Sigma_r w$, where $w \sim \mathcal{N}(0, I_r)$. Consider the affine map $\Phi: \mathbb{R}^r \to \mathbb{R}^N$ defined by $\Phi(w) = \mu + \Sigma_r w$. This map is a smooth bijection from $\mathbb{R}^r$ onto the support supp($z$) = $\mu + \text{Im}(\Sigma)$.

The probability measure of $z$, denoted $\mathbb{P}_z$, is the pushforward of the standard Gaussian measure $\gamma^r$ via $\Phi$. The standard Gaussian density on $\mathbb{R}^r$ is $g(w) = (2\pi)^{-r/2}e^{-\frac{1}{2}\|w\|^2}$. For any measurable set $A \subset \text{supp}(z)$, we have:

$$\mathbb{P}_z(A) = \gamma^r(\Phi^{-1}(A)) = \int_{\Phi^{-1}(A)} g(w) d\lambda^r(w). \qquad (11)$$

To find the density with respect to the Hausdorff measure $\mathcal{H}^r$ (the standard volume measure on the support), we apply the Area Formula (Federer, 1996) (change of variables). The Jacobian determinant of $\Phi$ is:

$$J\Phi = \sqrt{\det(\Sigma_r^\top \Sigma_r)} = \sqrt{|\Sigma|_+}.$$
Using the change of variables $z = \Phi(w)$, the volume elements relate via $d\mathcal{H}^r(z) = J\Phi d\lambda^r(w)$. Therefore:

$$d\lambda^r(w) = \frac{1}{\sqrt{|\Sigma|_+}} d\mathcal{H}^r(z).$$

Substituting this into the integral and using $w = \Phi^{-1}(z) = \Sigma_r^+(z - \mu)$:

$$ \begin{aligned} \mathbb{P}_z(A) &= \int_A g(\Phi^{-1}(z)) \frac{1}{\sqrt{|\Sigma|_+}} d\mathcal{H}^r(z) \\ &= \int_A \frac{1}{(2\pi)^{r/2} \sqrt{|\Sigma|_+}} \exp \left( -\frac{1}{2} \|\Sigma_r^+(z - \mu)\|^2 \right) d\mathcal{H}^r(z). \end{aligned} $$

Noting that $\|\Sigma_r^+(z - \mu)\|^2 = (z - \mu)^\top (\Sigma_r^+)^\top \Sigma_r^+(z - \mu) = (z - \mu)^\top \Sigma^+(z - \mu)$, we recover the density given in Definition A.1.

When the covariance matrix is non-singular, we recover the standard PDF of a Gaussian distribution.

**Degenerate Multivariate Gaussian and Mahalanobis Distance** The Mahalanobis distance quantifies the distance from $z$ to $\mu$ relative to the covariance structure and is defined as:

$$D_{\Sigma}(z, \mu) = \sqrt{(z - \mu)^{\top} \Sigma^{+} (z - \mu)}. \qquad (12)$$

When $z \notin \text{supp}(z)$, the density is zero, corresponding to an infinite effective distance outside the support.

*   **Eigenvalue Decomposition.** Consider the spectral decomposition of the covariance matrix:

$$\Sigma = U\Lambda U^{\top},$$

where:

- $U$ is an orthogonal matrix whose columns are the eigenvectors of $\Sigma$.

- $\Lambda = \text{diag}(\lambda_1, \lambda_2, \dots, \lambda_d)$ is a diagonal matrix with eigenvalues $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_r > 0 = \lambda_{r+1} = \dots = \lambda_d$, sorted in descending order, with zeros corresponding to the degenerate directions.

The pseudo-inverse is:

$$\Sigma^{+} = U\Lambda^{+}U^{\top}, \qquad (13)$$

where $\Lambda^{+} = \text{diag}(1/\lambda_1, \dots, 1/\lambda_r, 0, \dots, 0)$.

Transform the coordinates to the eigen-basis by defining $y = U^{\top}(x - \mu)$. The Mahalanobis distance simplifies to:

$$D_{\Sigma}(x, \mu) = \sqrt{y^{\top} \Lambda^{+} y} = \sqrt{\sum_{i=1}^{r} \frac{y_i^2}{\lambda_i}}, \qquad (14)$$

since $y_i = 0$ for $i > r$ on the support $\text{supp}(z)$. The exponent in the density becomes:

$$-\frac{1}{2} \sum_{i=1}^{r} \frac{y_i^2}{\lambda_i}. \qquad (15)$$

*   **Distance in Different Directions**

The effect of deviations from the mean $\mu$ along the directions of the eigenvectors (principal axes) depends on the corresponding eigenvalues:

- **Directions with large eigenvalues ($\lambda_i \gg 0$):** The term $\frac{y_i^2}{\lambda_i}$ grows slowly as $|y_i|$ increases. Deviations along these high-variance directions contribute less to reducing the density, effectively scaling the distance by $\sqrt{\lambda_i}$. This allows larger Euclidean deviations in these directions.
- **Directions with small positive eigenvalues** ($0 < \lambda_i \ll 1$): The term $\frac{y_i^2}{\lambda_i}$ grows rapidly even for small $|y_i|$. Deviations are heavily penalized, scaling the distance by $1/\sqrt{\lambda_i}$, making small deviations appear “far” in the Mahalanobis sense.

- **Directions with zero eigenvalues** ($\lambda_i = 0$): These correspond to the null space of $\Sigma$. Here, $y_i$ must be exactly zero for the density to be non-zero, enforced by the Dirac delta measure. Any non-zero deviation results in zero density, equivalent to an infinite Mahalanobis distance, indicating no variability in these directions.

**Integration on linear subspaces** The following lemma is useful for change of variable with linear mapping (Federer, 1996).

**Lemma A.2 (Integration on linear subspaces).** Let $\mathcal{N}(y; \mu, \Sigma)$ be the PDF of a Gaussian distribution in $\mathbb{R}^N$, with mean $\mu \in \mathbb{R}^N$ and covariance matrix $\Sigma \in \mathbb{R}^{N \times N}$ of rank $r_0 \le N$. For any linear application $B: \mathbb{R}^N \to \mathbb{R}^P$, and any measurable $f$, we have:

$$ \int_{\mathbb{R}^N} f(By)\mathcal{N}(y; \mu, \Sigma) d\mathcal{H}^{r_0}(y) = \int_{\mathbb{R}^P} f(v)\mathcal{N}(v; B\mu, B\Sigma B^\top) d\mathcal{H}^r(v), \quad (16) $$

where $\mathcal{N}(v; B\mu, B\Sigma B^\top)$ is the PDF of a (possibly degenerate) Gaussian distribution with respect to the $r$-dimensional Hausdorff measure $\mathcal{H}^r$, with $r$ being the rank of $B\Sigma B^\top$. This PDF is supported on the $r$-dimensional subspace $B\mu + \text{Im}(B\Sigma B^\top) \subset \text{Im}(B) \subset \mathbb{R}^P$.

*Proof.* The proof of this result is relatively straightforward by using the density of a degenerate Gaussian distribution Definition A.1. Let $\mathbf{y} \sim \mathcal{N}(\mu, \Sigma)$ and $\mathbf{z} = By$, then $\mathbf{z} \sim \mathcal{N}(B\mu, B\Sigma B^\top)$. We have

$$ \int_{\mathbb{R}^N} f(By)\mathcal{N}(y; \mu, \Sigma) d\mathcal{H}^{r_0}(y) = \mathbb{E}[f(By)] \quad (17) $$

$$ = \mathbb{E}[f(\mathbf{z})] \quad (18) $$

$$ = \int_{\mathbb{R}^P} f(v)\mathcal{N}(v; B\mu, B\Sigma B^\top) d\mathcal{H}^r(v). \quad (19) $$

□

*Remark A.3.* When $\Sigma$ is full rank, the Hausdorff measure $\mathcal{H}^{r_0}$ coincides with the Lebesgue measure on $\mathbb{R}^N$.

## A.2. Minimum Mean Square Error (MMSE) estimator

The MMSE estimator is the optimal estimator in the following sense

**Definition A.4 (MMSE estimator).** Given two random vectors $\mathbf{x} \in \mathbb{R}^N$, $\mathbf{y} \in \mathbb{R}^M$ and a linear operator $B: \mathbb{R}^M \to \mathbb{R}^N$. The MMSE estimator of $\mathbf{x}$ given $By$ is the best approximation random variable $\phi^*(By)$ to $\mathbf{x}$, in the least-square sense:

$$ \hat{x}_{\text{MMSE}} = \phi^* \circ B \quad \text{where} \quad \phi^* = \arg\min_{\phi: \mathbb{R}^N \to \mathbb{R}^M} \mathbb{E} \left[ \| \phi(By) - \mathbf{x} \|^2 \right] \quad (20) $$

It is well-known that the MMSE estimator coincides with the conditional expectation. That is, for any $y$ such that $p(y) > 0$, we have:

$$ \hat{x}_{\text{MMSE}}(y) = \mathbb{E}[\mathbf{x}|By = By] = \int x \cdot p(x|By)dx \quad (21) $$

Composing $\hat{x}_{\text{MMSE}}$ with the random variable $\mathbf{y}$, we get

$$ \hat{x}_{\text{MMSE}}(\mathbf{y}) = \mathbb{E}[\mathbf{x}|By]. $$

In particular, if $B=I$ (in this case $M=N$), the MMSE estimator reduces to the classical posterior mean: $\hat{x}_{\text{MMSE}}(y) = \mathbb{E}[\mathbf{x}|y=y]$. Note that both $\mathbb{E}[\mathbf{x}|y=y]$ and $\mathbb{E}[\mathbf{x}|y]$ are often called condition expectation, but these are different objects. In particular, $\mathbb{E}[\mathbf{x}|y=\cdot] = \hat{x}_{\text{MMSE}}(\cdot)$ is a function $\mathbb{R}^M \to \mathbb{R}^N$ while $\mathbb{E}[\mathbf{x}|y]$ is a random variable assuming values in $\mathbb{R}^N$. However, finding the MMSE estimator amounts to finding the optimal function $\phi^*$.
### A.3. Constrained Minimum Mean Square Error (MMSE) estimator

The classical MMSE estimator is defined as the best approximation function over the space of measurable function from $\mathbb{R}^M$ to $\mathbb{R}^N$. When adding functional constraints to the estimator, we would like to find the best approximation function over a subspace $\mathcal{M}$.

$$ \min_{\phi \in \mathcal{M}} \mathbb{E} \left[ \| \phi(B\mathbf{y}) - \mathbf{x} \|^2 \right] \qquad (22) $$

For example, the subspace $\mathcal{M}$ could be the set of measurable functions from $\mathbb{R}^N \to \mathbb{R}^N$ and translation equivariant.

### A.4. Optimality condition

We first state a simple first-order sufficient and necessary optimality condition for solving the Constrained MMSE.

**Proposition A.5 (Optimality condition).** Let $\mathbf{x} \in \mathbb{R}^N$, $\mathbf{y} \in \mathbb{R}^M$ be two random variables and $\mathcal{M}$ be a vector space of measurable functions from $\mathbb{R}^N$ to $\mathbb{R}^N$, $B$ be a linear operator from $\mathbb{R}^M$ to $\mathbb{R}^N$. Then $\phi^*$ is a minimizer of

$$ \min_{\phi \in \mathcal{M}} \mathbb{E} \left[ \| \phi(B\mathbf{y}) - \mathbf{x} \|^2 \right] $$

if and only if

$$ \mathbb{E} [\langle \phi(B\mathbf{y}), \phi^*(B\mathbf{y}) - \mathbf{x} \rangle] = 0 \quad \text{for all } \phi \in \mathcal{M}. \qquad (23) $$

*Proof.* Let $J(\phi) = \mathbb{E} \left[ \|\phi(B\mathbf{y}) - \mathbf{x}\|^2 \right]$. For all $t \in \mathbb{R}$ and for all $\phi \in \mathcal{M}$, we have

$$ \begin{aligned} J(\phi^* + t\varphi) &= \mathbb{E} \left[ \| (\phi^* + t\varphi)(B\mathbf{y}) - \mathbf{x} \|^2 \right] \\ &= J(\phi^*) + 2t \underbrace{\mathbb{E} [\langle \phi(B\mathbf{y}), \phi^*(B\mathbf{y}) - \mathbf{x} \rangle]}_{a} + t^2 \underbrace{\mathbb{E} [\|\phi(B\mathbf{y})\|^2]}_{b} \\ &= J(\phi^*) + 2at + bt^2 \end{aligned} $$

Therefore, $\phi^*$ is a minimizer of $J$ on $\mathcal{M}$ if and only if $J(\phi^* + t\varphi) \ge J(\phi^*)$ for all $t \in \mathbb{R}$ and all $\varphi \in \mathcal{M}$. This is equivalent to the condition that $2at + bt^2 \ge 0$, for all $t \in \mathbb{R}$ and all $\varphi \in \mathcal{M}$. Since this difference term is a quadratic function in $t$ and $b \ge 0$, it is non-negative for all $t \in \mathbb{R}$ if and only if $a = 0$. Therefore, we have the sufficient and necessary condition that $\mathbb{E} [\langle \phi(B\mathbf{y}), \phi^*(B\mathbf{y}) - \mathbf{x} \rangle] = 0$ for all $\varphi \in \mathcal{M}$. $\square$

The optimality condition Proposition A.5 states that the residual $\phi^*(By) - x$ is orthogonal, in the $L^2$ sense, to every perturbation in the feasible set $\mathcal{M}$. Equivalently, $\phi^*(By)$ is the orthogonal projection of $x$ onto $\mathcal{M}$ in the $L^2$ sense. When $\mathcal{M}$ is the space of all square-integrable functions of $y$ (i.e. no constraints) and $B = I$, this projection yields the classical MMSE estimator (posterior mean), $\mathbb{E}[x|y]$. In the constrained case, $\phi^*(By)$ can also be viewed as the orthogonal projection of the conditional expectation $\mathbb{E}[x|y]$ onto the subspace $\{\phi(By) : \phi \in \mathcal{M}\}$.

**Proposition A.6.** [Proposition 3.3 in the main paper] Given a closed set $\mathcal{M}$. The $\mathcal{M}$-constrained MMSE estimator in Definition 1.1 is the orthogonal projection (in $L^2$ sense) of the posterior mean $\mathbb{E}[x|y]$ onto the subspace of $y$-measurable random vectors of form $\chi = \{\phi(By) : \phi \in \mathcal{M}\}$. That is,

$$ \hat{x}_{\mathcal{M}}(\mathbf{y}) = \prod_{\chi} \mathbb{E}[x|y]. $$

*Proof.* The classes $\mathcal{M}_T$ and $\mathcal{M}_{T,\text{loc}}$ are linear subspaces of the space of measurable functions from $\mathbb{R}^N$ to $\mathbb{R}^N$: they are closed under addition and scalar multiplication, hence the projection is well-defined. The MMSE estimator in Proposition 2.1 is the posterior mean $\mathbb{E}[x|y]$, which is the orthogonal projection of $x$ onto the space of $y$-measurable random vectors. Similarly, the $\mathcal{M}$-constrained MMSE estimator is the projection of $x$ onto the subspace $\{\phi(By) : \phi \in \mathcal{M}\}$. Using the Pythagorean decomposition, for any $\phi \in \mathcal{M}$, we have

$$ \mathbb{E} \left[ \| \phi(B\mathbf{y}) - \mathbf{x} \|^2 \right] = \mathbb{E} \left[ \| \phi(B\mathbf{y}) - \mathbb{E}[x|y] \|^2 \right] + \mathbb{E} \left[ \| \mathbb{E}[x|y] - \mathbf{x} \|^2 \right] \qquad (24) $$
Therefore

$$ \arg\min_{\phi \in \mathcal{M}} \mathbb{E} \left[ \|\phi(B\mathbf{y}) - \mathbf{x}\|^2 \right] = \arg\min_{\phi \in \mathcal{M}} \mathbb{E} \left[ \|\phi(B\mathbf{y}) - \mathbb{E}[\mathbf{x}|\mathbf{y}]\|^2 \right] \quad (25) $$

and the $\mathcal{M}$-constrained MMSE in Definition 1.1 is the projection of the posterior mean $\mathbb{E}[\mathbf{x}|\mathbf{y}]$ onto the subspace $\{\phi(B\mathbf{y}) : \phi \in \mathcal{M}\}$.

Moreover, the decomposition in Equation (24) has interesting interpretation: it isolates the sources of error in the estimation. The first term $\mathbb{E}[\|\phi(B\mathbf{y}) - \mathbb{E}[\mathbf{x}|\mathbf{y}]\|^2]$ is the approximation error due to the restriction to the subspace $\mathcal{M}$, while the second term $\mathbb{E}[\|\mathbb{E}[\mathbf{x}|\mathbf{y}] - \mathbf{x}\|^2]$ is the irreducible Bayes error. $\square$

### A.5. Structural constraint: equivariant functions

The below definitions are taken from (Celledoni et al., 2021).

**Definition A.7 (Group).** A group, to be denoted $\mathcal{G}$, is a set equipped with an associative operator $\cdot : \mathcal{G} \times \mathcal{G} \to \mathcal{G}$, which satisfies the following conditions:

1. If $g_1, g_2 \in \mathcal{G}$ then $g_2 \cdot g_1 \in \mathcal{G}$

2. If $g_1, g_2, g_3 \in \mathcal{G}$ then $(g_1 \cdot g_2) \cdot g_3 = g_1 \cdot (g_2 \cdot g_3)$

3. There exists $\iota \in \mathcal{G}$ such that $e \cdot g = g \cdot \iota = g$ for all $g \in \mathcal{G}$.

4. If $g \in \mathcal{G}$ there exists $g^{-1} \in \mathcal{G}$ such that $g^{-1} \cdot g = g \cdot g^{-1} = \iota$.

**Definition A.8 (Group action).** Given a group $\mathcal{G}$ and a set $\mathbb{R}^N \subset \mathbb{R}^N$, we say that $\mathcal{G}$ acts on $\mathbb{R}^N$ if there exists a function $T : \mathcal{G} \times \mathbb{R}^N \to \mathbb{R}^N$ (we denote by $T_g(x)$ for $g \in \mathcal{G}$ and $x \in \mathbb{R}^N$) that satisfies:

$$ T_{g_1} \circ T_{g_2} = T_{g_1 \cdot g_2} \quad \text{and} \quad T_{\iota} = \text{id} $$

Given a general group $\mathcal{G}$. A function $\phi : \mathbb{R}^N \to \mathbb{R}^N$ and group action $T$ of $\mathcal{G}$ on $\mathbb{R}^N$. A function $\phi$ is called $\mathcal{G}$-equivariant if it satisfies

$$ \phi(T_g(y)) = T_g \phi(y) \quad \text{for all } x \in \mathbb{R}^N \text{ and for all } g \in \mathcal{G}. $$

**Proposition A.9.** If the group $\mathcal{G}$ is finite, the following properties hold true:

1. **Invariance:** for any function $\phi : \mathbb{R}^N \to \mathbb{R}$, the function $\bar{\phi} = \sum_g \phi \circ T_g$ is invariant. (And similarly for function defined in $\mathbb{R}^M$).

2. **Equivariance:** for any function $\phi : \mathbb{R}^N \to \mathbb{R}^N$, the function $\bar{\phi} = \sum_g T_g^{-1} \circ \phi \circ T_g$ is equivariant. This is called **Reynolds averaging**.

For image data $x \in \mathbb{R}^N$, let $H, W \in \mathbb{N}$ be the dimensions of a discrete grid $\Omega = \mathbb{Z}_H \times \mathbb{Z}_W$ with $H \times W = N$.

**Definition A.10 (Translation equivariant functions).** Let $\mathcal{T} = \mathbb{Z}_H \times \mathbb{Z}_W$ be the group of 2D cyclic translations. For every group element $g = (g_h, g_v) \in \mathcal{T}$, we define the translation operator $T_g : \mathbb{R}^N \to \mathbb{R}^N$ as the permutation matrix that acts on an image $x \in \mathbb{R}^N$ by shifting its indices:

$$ (T_g x)[i, j] = x[(i - g_h) \mod H, (j - g_v) \mod W] $$

for all $(i, j) \in \Omega$. A measurable map $\phi : \mathbb{R}^N \to \mathbb{R}^N$ is said to be translation equivariant if it commutes with the translation operator for all $g \in \mathcal{T}$:

$$ \phi(T_g x) = T_g \phi(x) \quad \text{for all } x \in \mathbb{R}^N. $$
### A.6. Structural constraints: local and translation equivariant functions

We first define precisely the patch extractor. Let $n \in \Omega$ be a pixel coordinate on the grid. We define the patch extractor $\Pi_n: x \in \mathbb{R}^N \to \Pi_n x = x[\omega_n] \in \mathbb{R}^P$ which extracts a square patch $x[\omega_n]$ of size $\sqrt{P} \times \sqrt{P}$ centered at $n$. The extraction uses circular boundary conditions, such that the patch is given by the grid values at indices:

$$ \omega_n = \{n + \delta \ (\text{mod}(H, W)) \mid \delta \in \Delta\} $$

where $\Delta$ is the set of offsets defining the square neighborhood centered at zero.

**Definition A.11 (Local and translation-equivariant functions).** A measurable map $\phi: \mathbb{R}^N \to \mathbb{R}^N$ is said to be local if it can be represented as a sliding window operation. Specifically, if there exist a measurable function $f: \mathbb{R}^P \to \mathbb{R}$ such that for all $x \in \mathbb{R}^N$ and all $n \in \Omega$:

$$ \phi(x)[n] = f(\Pi_n x). $$

By construction, any such function $\phi$ is also translation equivariant due to the circular boundary conditions of the patch extractor $\Pi_n$. We denote by $M_{\mathcal{T},\text{loc}}$ the set of all such maps.

This class captures standard CNNs with finite kernels and weight sharing or patch-based methods, e.g. MLPs acting on patches (Khorashadizadeh et al., 2025): the output of $\phi$ at pixel $n$, denoted $\phi(x)[n]$ or $\phi_n(x)$, depends only on the local patch (receptive field) $x[\omega_n]$. Note that by construction, functions in Definition A.11 are translation equivariant.

### B. Analytical solution for inverse problems

In the following, we will derive the analytical solution for the MMSE estimator under these constraints.

Consider a random variable $x \sim p_x$, $x \in \mathbb{R}^N$ and the forward (measurement) model

$$ y = Ax + e \quad \text{where} \quad e \sim \mathcal{N}(0, \sigma^2 I). \qquad (26) $$

We would like to find the MMSE estimator of $x$ given $y$, with or without constraints. In this section, we will derive the analytical solution for the MMSE estimator under various constraints, when the true underlying distribution of $x$ is replaced by the empirical distribution $p_D$. Under Gaussian noise, the likelihood of the measurement $y$ given $x$ is given by

$$ p(y|x) = \mathcal{N}(y; Ax, \sigma^2 I) \propto \exp\left(-\frac{\|y - Ax\|^2}{2\sigma^2}\right). \qquad (27) $$

Therefore, for any function $\phi^*$ and $\phi$, we have

$$ \begin{aligned} \mathbb{E}[\langle\phi(B\mathbf{y}), \phi^*(B\mathbf{y}) - \mathbf{x}\rangle] &= \mathbb{E}_{\mathbf{x}\sim p_D}\mathbb{E}_{\mathbf{y}|\mathbf{x}}[\langle\phi(B\mathbf{y}), \phi^*(B\mathbf{y}) - \mathbf{x}\rangle] \\ &= \frac{1}{|\mathcal{D}|}\sum_{\mathbf{x} \in \mathcal{D}} \int_{\mathbb{R}^M} \langle\phi(By), \phi^*(By) - \mathbf{x}\rangle p(y|x)dy \\ &= \frac{1}{|\mathcal{D}|}\sum_{\mathbf{x} \in \mathcal{D}} \int_{\mathbb{R}^M} \langle\phi(By), \phi^*(By) - \mathbf{x}\rangle \mathcal{N}(y; Ax, \sigma^2 I) dy \end{aligned} \qquad (28) $$

This simplified expression (28) is useful for deriving the analytical solution of the MMSE estimator under various constraints.

#### B.1. Unconstrained MMSE estimator

We start with the unconstrained MMSE estimator, which is the optimal estimator in the least-square sense. Then, we will derive the equivariant MMSE estimator, which is the optimal estimator under the constraint of equivariance to translation. Then, we will derive the local MMSE estimator, which is the optimal estimator under the constraint of locality. Finally, we will also show that combining both constraints leads to a local and equivariant MMSE estimator.

**Proposition B.1 (Unconstrained MMSE estimator).** When the data distribution is replaced by the empirical distribution $p_D$, the unconstrained MMSE estimator of $x$ given $y$ is given by

$$ \hat{x}_{\text{MMSE}}(y) = \phi^*(By) \quad \text{where} \quad \phi^*(v) = \frac{\sum_{x \in \mathcal{D}} x \mathcal{N}(v; BAx, \sigma^2 BB^\top)}{\sum_{x \in \mathcal{D}} \mathcal{N}(v; BAx, \sigma^2 BB^\top)} \quad \text{for any } v \in \text{Im}(B). \qquad (29) $$

The estimator is well-defined for all $y \in \mathbb{R}^M$.
*Proof.* We will verify that $\phi^*$ satisfies the optimality condition from Proposition A.5. Using Equation (28), for any $\varphi$, we have:

$$
\begin{align}
\mathbb{E} [\langle \varphi(B\boldsymbol{y}), \phi^*(B\boldsymbol{y}) - \boldsymbol{x} \rangle] &= \frac{1}{|\mathcal{D}|} \sum_{x \in \mathcal{D}} \int_{\mathbb{R}^M} \langle \varphi(B\boldsymbol{y}), \phi^*(B\boldsymbol{y}) - x \rangle \mathcal{N}(y, A x, \sigma^2 I_M) dy \nonumber \\
&= \frac{1}{|\mathcal{D}|} \sum_{x \in \mathcal{D}} \int_{\mathbb{R}^N} \langle \varphi(v), \phi^*(v) - x \rangle \mathcal{N}(v, B A x, \sigma^2 B B^\top) d\mathcal{H}^r(v) \tag{30} \\
&= \frac{1}{|\mathcal{D}|} \int_{\mathbb{R}^N} \left\langle \varphi(v), \phi^*(v) \sum_{x \in \mathcal{D}} \mathcal{N}(v, B A x, \sigma^2 B B^\top) - x \sum_{x \in \mathcal{D}} \mathcal{N}(v, B A x, \sigma^2 B B^\top) \right\rangle d\mathcal{H}^r(v) \nonumber \\
&= 0 \nonumber
\end{align}
$$

where Equation (30) comes from the change of variable $v = By$ and Lemma A.2. The last equality holds by definition of $\phi^*$.
□

## B.2. Equivariant MMSE estimator

**Proposition B.2 (Translation equivariant MMSE estimator).** The translation equivariant MMSE estimator is as follows, for any $y \in \mathbb{R}^M$:

$$ \hat{x}_{\mathcal{T}}(y) = \phi^*(By) \quad \text{where} \quad \phi^*(v) = \frac{\sum_{x \in \mathcal{D}, g \in \mathcal{T}} T_g x \cdot \mathcal{N}(T_g^{-1}v; B A x, \sigma^2 B B^\top)}{\sum_{x \in \mathcal{D}, g \in \mathcal{T}} \mathcal{N}(T_g^{-1}v; B A x, \sigma^2 B B^\top)} \quad (31) $$

for any $v \in \bigcup_g T_g \mathrm{Im}(B) \subset \mathbb{R}^N$.

*Proof.* We will verify that $\phi^*$ is admissible and satisfies the optimality condition Proposition A.5.

**Admissibility.** Let $q(v, x) = \mathcal{N}(v; B A x, \sigma^2 B B^\top)$ denote the PDF of a Gaussian distribution with mean $B A x$ and covariance $\sigma^2 B B^\top$. The optimal estimator can be written as:

$$ \phi^*(v) = \frac{\sum_{x \in \mathcal{D}, g \in \mathcal{T}} T_g x q(T_g^{-1} v, x)}{\sum_{x \in \mathcal{D}, g \in \mathcal{T}} q(T_g^{-1} v, x)} \quad (32) $$

The denominator $\sum_{x \in \mathcal{D}, g \in \mathcal{T}} q(T_g^{-1}v, x) = \sum_{g \in \mathcal{T}} \bar{q}_1(T_g^{-1}v)$ is $\mathcal{T}$-invariant by Proposition A.9, where $\bar{q}_1(v) = \sum_{x \in \mathcal{D}} q(v, x)$.

The nominator $\sum_{x \in \mathcal{D}, g \in \mathcal{T}} T_g x \cdot q(T_g^{-1}v, x) = \sum_{g \in \mathcal{T}} T_g \bar{q}_2(T_g^{-1}v)$ is $\mathcal{T}$-equivariant, where $\bar{q}_2(v) = \sum_{x \in \mathcal{D}} x \cdot q(v, x)$. Therefore, $\phi^*$ is admissible, that is $\phi^* \in \mathcal{M}_{\mathcal{T}}$.

**Optimality condition.** By using Equation (28), for any $\varphi \in \mathcal{M}_{\mathcal{T}}$, we have:

$$
\begin{align}
\mathbb{E} [\langle \varphi(B\boldsymbol{y}), \phi^*(B\boldsymbol{y}) - \boldsymbol{x} \rangle] \\
&= \frac{1}{|\mathcal{D}|} \sum_{x \in \mathcal{D}} \int_{\mathbb{R}^M} \langle \varphi(B\boldsymbol{y}), \phi^*(B\boldsymbol{y}) - x \rangle \mathcal{N}(y; A x, \sigma^2 I_M) dy \nonumber \\
&= \frac{1}{|\mathcal{D}|} \sum_{x \in \mathcal{D}} \int_{\mathbb{R}^N} \langle \varphi(v), \phi^*(v) - x \rangle \mathcal{N}(v; B A x, \sigma^2 B B^\top) d\mathcal{H}^r(v) \tag{33}
\end{align}
$$

$$ = \frac{1}{|\mathcal{D}||\mathcal{T}|} \sum_{g \in \mathcal{T}, x \in \mathcal{D}} \int_{\mathbb{R}^N} \langle T_g^{-1} \varphi(T_g v), T_g^{-1} \phi^*(T_g v) - x \rangle \mathcal{N}(v; B A x, \sigma^2 B B^\top) d\mathcal{H}^r(v) \quad (34) $$

$$ = \frac{1}{|\mathcal{D}||\mathcal{T}|} \sum_{g \in \mathcal{T}, x \in \mathcal{D}} \int_{\mathbb{R}^N} \langle \varphi(T_g v), \phi^*(T_g v) - T_g x \rangle \mathcal{N}(v; B A x, \sigma^2 B B^\top) d\mathcal{H}^r(v) \quad (35) $$

$$ = \frac{1}{|\mathcal{D}||\mathcal{T}|} \sum_{g \in \mathcal{T}, x \in \mathcal{D}} \int_{\mathbb{R}^N} \langle \varphi(v), \phi^*(v) - T_g x \rangle \mathcal{N}(T_g^{-1}v; B A x, \sigma^2 B B^\top) d\mathcal{H}^r(v) \quad (35) $$
$$ 
\begin{aligned}
&= \frac{1}{|\mathcal{D}||\mathcal{T}|} \int_{\mathbb{R}^N} \left\langle \varphi(v), \phi^*(v) \sum_{g \in \mathcal{T}, x \in \mathcal{D}} q(T_g^{-1}v, n, x) - \sum_{g \in \mathcal{T}, x \in \mathcal{D}} T_g x \cdot q(T_g^{-1}v, n, x) \right\rangle d\mathcal{H}^r(v) \\
&= 0
\end{aligned}
$$

where $q(T_g^{-1}v, n, x) \stackrel{\text{def}}{=} \mathcal{N}(T_g^{-1}v; BAx, \sigma^2 BB^\top)$ and

* In Equation (33), we apply Lemma A.2 for $B$ in the change of variables.

* In Equation (34), we use the fact that $\varphi = \varphi \circ T_g \circ T_g^{-1} = T_g \circ \varphi \circ T_g^{-1}$ for all $g \in \mathcal{T}$ and $\varphi \in \mathcal{M}_G$. Moreover, $T_g^{-1} = T_g^\top$ for all $g \in \mathcal{T}$.

* In Equation (35), we use the change of variable $T_g v \to v$ and the fact that $T_g$ is an isometry.

□

We next state and prove that the augmented MMSE estimator is reconstruction equivariant.

**Lemma B.3 (Reconstruction equivariance of $\hat{x}_{\text{MMSE}}^{\text{aug}}$).** Let $A : \mathbb{R}^N \to \mathbb{R}^M$ be a circular convolution operator, i.e., $AT_g = T_g A$ for all $g \in \mathcal{T}$. Then, the data-augmented MMSE estimator $\hat{x}_{\text{MMSE}}^{\text{aug}}$ defined in Definition 3.4 satisfies the reconstruction equivariance property in Equation (5):

$$ \hat{x}_{\text{MMSE}}^{\text{aug}}(AT_g \bar{x} + e) = T_g \hat{x}_{\text{MMSE}}^{\text{aug}}(A\bar{x} + T_g^{-1}e) $$

for all $\bar{x} \in \mathbb{R}^N$, all $e \in \mathbb{R}^M$, and all $g \in \mathcal{T}$.

*Proof.* Let $y = AT_g \bar{x} + e$ be the input observation. Recall the definition of the data-augmented MMSE estimator in Definition 3.4:

$$ \hat{x}_{\text{MMSE}}^{\text{aug}}(y) = \frac{\sum_{x \in \mathcal{T}(D)} x \cdot \exp\left(-\frac{1}{2\sigma^2} \|y - Ax\|^2\right)}{\sum_{x \in \mathcal{T}(D)} \exp\left(-\frac{1}{2\sigma^2} \|y - Ax\|^2\right)}. $$

Substituting the input $y = AT_g \bar{x} + e$ into the estimator:

$$ \hat{x}_{\text{MMSE}}^{\text{aug}}(y) = \frac{\sum_{x \in \mathcal{T}(D)} x \cdot \exp\left(-\frac{1}{2\sigma^2} \|AT_g \bar{x} + e - Ax\|^2\right)}{\sum_{x \in \mathcal{T}(D)} \exp\left(-\frac{1}{2\sigma^2} \|AT_g \bar{x} + e - Ax\|^2\right)}. $$

Since $\mathcal{T}(D)$ is invariant under the group action, we can perform a change of variable $z = T_g^{-1}x$. As $x$ traverses $\mathcal{T}(D)$, $z$ also traverses $\mathcal{T}(D)$. Note that $x = T_g z$. Substituting $x$ with $T_g z$ in the numerator and denominator:

$$ \hat{x}_{\text{MMSE}}^{\text{aug}}(y) = \frac{\sum_{z \in \mathcal{T}(D)} T_g z \cdot \exp\left(-\frac{1}{2\sigma^2} \|AT_g \bar{x} + e - AT_g z\|^2\right)}{\sum_{z \in \mathcal{T}(D)} \exp\left(-\frac{1}{2\sigma^2} \|AT_g \bar{x} + e - AT_g z\|^2\right)}. $$

Using the assumption that $A$ is a circular convolution ($AT_g = T_g A$) and $T_g$ is linear, the term inside the norm becomes:

$$ 
\begin{aligned}
AT_g \bar{x} + e - AT_g z &= T_g A \bar{x} + T_g (T_g^{-1} e) - T_g Az \\
&= T_g (A \bar{x} + T_g^{-1} e - Az).
\end{aligned}
$$

Since $T_g$ is unitary, it preserves the norm, i.e., $\|T_g u\| = \|u\|$. Therefore:

$$ \|AT_g \bar{x} + e - AT_g z\|^2 = \|A \bar{x} + T_g^{-1} e - Az\|^2. $$
Substituting this back into the estimator expression and factoring the linear operator $T_g$ out of the sum in the numerator:

$$ \begin{aligned} \hat{x}_{\text{MMSE}}^{\text{aug}}(y) &= \frac{T_g \sum_{z \in \mathcal{T}(D)} z \cdot \exp\left(-\frac{1}{2\sigma^2} \|(A\bar{x} + T_g^{-1}e) - Az\|_2^2\right)}{\sum_{z \in \mathcal{T}(D)} \exp\left(-\frac{1}{2\sigma^2} \|(A\bar{x} + T_g^{-1}e) - Az\|_2^2\right)} \\ &= T_g \left( \frac{\sum_{z \in \mathcal{T}(D)} z \cdot \mathcal{N}(A\bar{x} + T_g^{-1}e; Az, \sigma^2\mathbf{I})}{\sum_{z' \in \mathcal{T}(D)} \mathcal{N}(A\bar{x} + T_g^{-1}e; Az', \sigma^2\mathbf{I})} \right). \end{aligned} $$

The term in the parentheses is exactly the estimator evaluated at the input $A\bar{x} + T_g^{-1}e$. Thus:

$$ \hat{x}_{\text{MMSE}}^{\text{aug}}(AT_g\bar{x} + e) = T_g \hat{x}_{\text{MMSE}}^{\text{aug}}(A\bar{x} + T_g^{-1}e). $$

<figure><img src="image_13.png" /></figure>

Now we state and prove the properties of the E-MMSE estimator presented in Corollary 3.7.

**Proof of Corollary 3.7.** The first point is a direct consequence of Theorem 3.5. We prove the equivalence in the second point by establishing both directions. We first express the weights of both estimators. The two estimators admit the same form:

$$ \hat{x}_{\mathcal{T}}(y) = \sum_{x \in \mathcal{D}, g \in \mathcal{T}} T_g x \cdot w_g(x|y) \quad \text{and} \quad \hat{x}_{\text{MMSE}}^{\text{aug}}(y) = \sum_{x \in \mathcal{D}, g \in \mathcal{T}} T_g x \cdot w_g^{\text{aug}}(x|y). $$

The weights of the E-MMSE estimator $\hat{x}_{\mathcal{T}}$ corresponding to the component $(x, g)$ are given by:

$$ \begin{aligned} w_g(x|y) &= \frac{\mathcal{N}(T_g^{-1}By; BAx, \sigma^2BB^\top)}{\sum_{x' \in \mathcal{D}, g' \in \mathcal{T}} \mathcal{N}(T_{g'}^{-1}By; BAx', \sigma^2BB^\top)} \\ &= \frac{\exp\left(-\|B^{-1}(T_g^{-1}By - BAx)\|^2/(2\sigma^2)\right)}{\sum_{x' \in \mathcal{D}, g' \in \mathcal{T}} \exp\left(-\|B^{-1}(T_{g'}^{-1}By - BAx')\|^2/(2\sigma^2)\right)} \\ &= \frac{\exp\left(-\|B^{-1}T_g^{-1}By - Ax\|^2/(2\sigma^2)\right)}{\sum_{x' \in \mathcal{D}, g' \in \mathcal{T}} \exp\left(-\|B^{-1}T_{g'}^{-1}By - Ax'\|^2/(2\sigma^2)\right)} \end{aligned} $$

where we used the invertibility of $B$ to write $B^+ = B^{-1}$ and simplified the Mahalanobis distance. The weights of the data-augmented estimator $\hat{x}_{\text{MMSE}}^{\text{aug}}$ are:

$$ \begin{aligned} w_g^{\text{aug}}(x|y) &= \frac{\mathcal{N}(By; BAT_g x, \sigma^2 BB^\top)}{\sum_{x' \in \mathcal{D}, g' \in \mathcal{T}} \mathcal{N}(By; BAT_{g'} x', \sigma^2 BB^\top)} \\ &= \frac{\exp\left(-\|B^{-1}(By - BAT_g x)\|^2/(2\sigma^2)\right)}{\sum_{x' \in \mathcal{D}, g' \in \mathcal{T}} \exp\left(-\|B^{-1}(By - BAT_{g'} x')\|^2/(2\sigma^2)\right)} \\ &= \frac{\exp\left(-\|y - AT_g x\|^2/(2\sigma^2)\right)}{\sum_{x' \in \mathcal{D}, g' \in \mathcal{T}} \exp\left(-\|y - AT_{g'} x'\|^2/(2\sigma^2)\right)}. \end{aligned} $$

**Sufficiency** (=$\Leftarrow$) Assume $A$ and $B$ are circular convolutions (they commute with $T_g$) and $B$ is invertible.

Commutativity implies $B^{-1}T_g^{-1} = T_g^{-1}B^{-1}$. We can rearrange the term in the exponent of the above weights as follows:

$$ \|B^{-1}T_g^{-1}By - Ax\|^2 = \|B^{-1}BT_g^{-1}y - B^{-1}BAx\|^2 = \|T_g^{-1}y - Ax\|^2. $$
Since $T_g$ is unitary and $A$ commutes with $T_g$, we have:

$$ \|T_g^{-1}y - Ax\|^2 = \|T_g(T_g^{-1}y - Ax)\|^2 = \|y - T_gAx\|^2 = \|y - AT_gx\|^2. $$

Thus $\hat{x}_{\mathcal{T}} = \hat{x}_{\text{MMSE}}^{\text{aug}}$. Furthermore, since the weights are independent of $B$, the physics-agnostic and physics-aware solvers are identical. The reconstruction equivariance follows immediately from Lemma B.3.

**Necessity (⇒)** Assume that $w_g^{\text{aug}}(x|y) = w_g(x|y)$ for all $y \in \mathbb{R}^M$ and all $x \in \mathbb{R}^N$. This implies that:

$$ \|B^{-1}T_g^{-1}By - Ax\|^2 = \|y - AT_gx\|^2 + c(y) \quad (36) $$

for some function $c(y)$ independent of $x$ and $g$. We expand the squared norms and the terms depending solely on $y$ can be absorbed into $c(y)$.

$$ \text{LHS} = \underbrace{\|B^{-1}T_g^{-1}By\|^2}_{\text{depends only on } y} - 2\langle B^{-1}T_g^{-1}By, Ax \rangle + \|Ax\|^2, $$

$$ \text{RHS} = \underbrace{\|y\|^2}_{\text{depends only on } y} - 2\langle y, AT_gx \rangle + \|AT_gx\|^2 + c(y). $$

For the equality in Equation (36) to hold for all $x \in \mathbb{R}^N$ and $y \in \mathbb{R}^M$, the cross-terms (the bilinear forms coupling $y$ and $x$) must be identical. Moreover, since the estimator $\hat{x}_{\text{MMSE}}^{\text{aug}}$ is independent of $B$, the Equation (36) must hold for all invertible $B$. In particular, taking $B = I$, we deduce that

$$ \langle T_g^{-1}y, Ax \rangle = \langle y, AT_gx \rangle \quad \text{for all } x \in \mathbb{R}^N, g \in \mathcal{T}, y \in \mathbb{R}^M. $$

This equality yields $T_g A = AT_g$, or that $A$ commutes with $T_g$. Plugging this result into the cross-term equality for a general invertible $B$, we get:

$$ \langle B^{-1}T_g^{-1}By, Ax \rangle = \langle y, AT_gx \rangle = \langle T_g^{-1}y, Ax \rangle \quad \text{for all } x \in \mathbb{R}^N, g \in \mathcal{T}, y \in \mathbb{R}^M. $$

This implies that $A^\top B^{-1}T_g^{-1}B = A^\top T_g^{-1}$ for all $g$. For invertible $A$, this implies that $B^{-1}T_g^{-1}B = T_g^{-1}$ for all $g \in \mathcal{T}$. Multiplying by $B$ on each side of the equality implies that $B$ commutes with $T_g$. $\square$

### B.3. Local and Translation Equivariant MMSE estimator

**Proposition B.4 (Local and translation equivariant MMSE estimator).** Suppose that the $N$ matrices $Q_n = \Pi_n B \in \mathbb{R}^{P \times M}$ have constant rank $r > 0$. The local and equivariant MMSE is defined for any $y \in \mathbb{R}^M$ by:

$$ \hat{x}_{\mathcal{T},\text{loc}}(y) = \phi^*(By), \quad (37) $$

where

$$ \phi_n^* = f_{\phi^*} \circ \Pi_n \quad \text{with} \quad f_{\phi^*}(v) = \frac{\sum_{x \in \mathcal{D}} \sum_{n=1}^{N} x_n \mathcal{N}(v; Q_n Ax, \sigma^2 Q_n Q_n^\top)}{\sum_{x \in \mathcal{D}} \sum_{n=1}^{N} \mathcal{N}(v; Q_n Ax, \sigma^2 Q_n Q_n^\top)}. \quad (38) $$

*Proof.* We will verify that $\phi^*$ is admissible and satisfies the optimality condition Proposition A.5.

**Admissibility.** Firstly, we show that $\phi^* \in \mathcal{M}_{\mathcal{T},\text{loc}}$. By construction, we have $\phi^* = (\phi_1^*, \dots, \phi_N^*) = (f_{\phi^*} \circ \Pi_1, \dots, f_{\phi^*} \circ \Pi_N)$, therefore $\phi^*$ is local. For any $g \in \{1, \dots, N\}$, the group transformation (translation) $T_g$ is defined as $T_g : (x_1, \dots, x_N) \in \mathbb{R}^N \mapsto (x_{1-g}, \dots, x_{N-g}) \in \mathbb{R}^N$, where $x_{n-g} = x_{n-g} \equiv_N$ (circular boundary). For any $g$, we have

$$ \begin{aligned} \phi^* \circ T_g &= (\phi_1^* \circ T_g, \dots, \phi_N^* \circ T_g) \\ &= (f_{\phi^*} \circ \Pi_1 \circ T_g, \dots, f_{\phi^*} \circ \Pi_N \circ T_g) \\ &= (f_{\phi^*} \circ \Pi_{1-g}, \dots, f_{\phi^*} \circ \Pi_{N-g}) \quad \text{where } \Pi_{n-g} = \Pi_{n-g} \equiv_N \text{ (circular boundary)} \\ &= T_g \circ \phi^* \end{aligned} $$
Therefore, $\phi^*$ is translation equivariant. We deduce that $\phi^* \in \mathcal{M}_{\mathcal{T},\text{loc}}$.

**Optimality condition.** We will verify that this estimator satisfies the optimality condition Proposition A.5. By using Equation (28), for any $\varphi \in \mathcal{M}_{\mathcal{T},\text{loc}}$, we have:

$$ 
\begin{aligned}
& \mathbb{E} [\langle \varphi(B\mathbf{y}), \phi^*(B\mathbf{y}) - \mathbf{x} \rangle] \\
&= \frac{1}{|\mathcal{D}|} \sum_{x \in \mathcal{D}} \int_{\mathbb{R}^M} \langle \varphi(By), \phi^*(By) - x \rangle \mathcal{N}(y; Ax, \sigma^2 I_N) dy \\
&= \frac{1}{|\mathcal{D}|} \sum_{x \in \mathcal{D}} \sum_{n=1}^N \int_{\mathbb{R}^M} f_{\varphi}(\Pi_n By) (f_{\phi^*}(\Pi_n By) - x_n) \mathcal{N}(y; Ax, \sigma^2 I_N) dy \\
&= \frac{1}{|\mathcal{D}|} \sum_{x \in \mathcal{D}} \sum_{n=1}^N \int_{\mathbb{R}^P} f_{\varphi}(v) (f_{\phi^*}(v) - x_n) \underbrace{\mathcal{N}(v; Q_n Ax, \sigma^2 Q_n Q_n^\top)}_{q(v,n,x)} d\mathcal{H}^r(v) \\
&= \frac{1}{|\mathcal{D}|} \int_{\mathbb{R}^P} f_{\varphi}(v) \left( f_{\phi^*}(v) \sum_{x \in \mathcal{D}} \sum_{n=1}^N q(v, n, x) - \sum_{x \in \mathcal{D}} \sum_{n=1}^N x_n q(v, n, x) \right) d\mathcal{H}^r(v) \\
&= 0
\end{aligned}
$$

□

The constant rank assumption in Proposition B.4 can be relaxed by stratifying the image space according to the rank of the matrices $Q_n$ as follows.

**Proposition B.5 (Local and translation equivariant MMSE estimator – rank stratification).** Let $Q_n = \Pi_n B \in \mathbb{R}^{P \times M}$. For each rank $r \in \{0, 1, \dots, P\}$, define the index set and the union of subspaces:

$$ I_r = \{n \in [1:N] : \text{rank}(Q_n) = r\}, \quad E_r = \bigcup_{n \in I_r} \text{Im}(Q_n) \subset \mathbb{R}^P. $$

We define the disjoint strata recursively: $\bar{E}_r = E_r \setminus \bigcup_{r' < r} E_{r'}$. Then $\bigcup_{r=0}^P \bar{E}_r = \bigcup_{n=1}^N \text{Im}(Q_n)$ is a disjoint union and, for any $r' < r$, we have $\mathcal{H}^r(E_{r'}) = 0$.

Define the weighted density sums for $v \in \mathbb{R}^P$:

$$ 
\begin{aligned}
a_r(v) &= \sum_{x \in \mathcal{D}} \sum_{n \in I_r} x_n \mathcal{N}(v; Q_n Ax, \sigma^2 Q_n Q_n^\top), \\
b_r(v) &= \sum_{x \in \mathcal{D}} \sum_{n \in I_r} \mathcal{N}(v; Q_n Ax, \sigma^2 Q_n Q_n^\top).
\end{aligned}
$$

The local and translation equivariant MMSE estimator $\phi^*$ is given component-wise by $\phi_n^* = f_{\phi^*} \circ \Pi_n$, where:

$$ f_{\phi^*}(v) = \begin{cases} \sum_{r=0}^{P} \frac{a_r(v)}{b_r(v)} \mathbb{1}_{\bar{E}_r}(v) & \text{if } v \in \bigcup \text{Im}(Q)_n \text{ and } b_r(v) > 0, \\ 0 & \text{otherwise.} \end{cases} \quad (39) $$

**Proof.** Admissibility (locality and translation equivariance).

* Locality. By definition, the estimator is defined component-wise by $\phi_n^* = f_{\phi^*} \circ \Pi_n$, so it is local.

* Translation equivariance on $\text{Im}(B)$. Let $T_g$ denote the circular translation on $\mathbb{R}^N$, the selection operators commute with $T_g$ by a simple index check

$$ \Pi_n \circ T_g = \Pi_{n-g}, \quad \forall n, g \in [1, N]. $$

Then, for any $v \in \text{Im}(B)$ and any $n$,

$$ [\phi^*(T_g v)]_n = f_{\phi^*}(\Pi_n T_g v) = f_{\phi^*}(\Pi_{n-g} v) = [\phi^*(v)]_{n-g} = [T_g \phi^*(v)]_n. \quad (40) $$

Hence $\phi^* \circ T_g = T_g \circ \phi^*$ on $\text{Im}(B)$, i.e., it is translation equivariant. Combining with locality gives $\phi^* \in \mathcal{M}_{\mathcal{T},\text{loc}}$.
• Well-definedness. On each $\bar{E}_r$, $f_{\phi^*}$ is defined by the ratio $a_r/b_r$. If $b_r(v) = 0$ on a negligible set (w.r.t. $\mathcal{H}^r$), assign any fixed value (e.g. 0); this does not affect admissibility nor optimality.

**Optimality.** For any $\varphi \in \mathcal{M}_{\mathcal{T},\text{loc}}$, we will verify that the optimality condition Proposition A.5 holds. By using Equation (28), we have:

$$
\begin{align*}
& \mathbb{E} [\langle \varphi(B\boldsymbol{y}), \phi^*(B\boldsymbol{y}) - \boldsymbol{x} \rangle] \\
&= \frac{1}{|\mathcal{D}|} \sum_{x \in \mathcal{D}} \int_{\mathbb{R}^M} \langle \varphi(By), \phi^*(By) - x \rangle \mathcal{N}(y; Ax, \sigma^2 I_N) dy \\
&= \frac{1}{|\mathcal{D}|} \sum_{x \in \mathcal{D}} \sum_{n=1}^N \int_{\mathbb{R}^M} f_{\varphi}(\Pi_n By)(f_{\phi^*}(\Pi_n By) - x_n) \mathcal{N}(y; Ax, \sigma^2 I) dy \\
&\stackrel{(\star)}{=} \frac{1}{|\mathcal{D}|} \sum_{x \in \mathcal{D}} \sum_{n=1}^N \int_{\mathbb{R}^P} f_{\varphi}(v)(f_{\phi^*}(v) - x_n) \mathcal{N}(v; Q_n Ax, \sigma^2 Q_n Q_n^\top) d\mathcal{H}^{r_n}(v),
\end{align*}
$$

where $(\star)$ uses Lemma A.2 and $r_n = \text{rank}(Q_n)$. We decompose the summation over $n$ by rank. Note that for $n \in I_r$, the domain is $\text{Im}(Q)_n$. We observe that $\text{Im}(Q)_n \setminus \bar{E}_r \subseteq \bigcup_{r' < r} E_{r'}$. Since the union of lower-dimensional subspaces has $\mathcal{H}^r$-measure zero, we can restrict the integration domain from $\text{Im}(Q)_n$ to $\text{Im}(Q)_n \cap \bar{E}_r$ without changing the value of the integral. Thus, the previous expression becomes:

$$
\sum_{r=0}^{P} \int_{\bar{E}_r} f_{\varphi}(v) \underbrace{\left( f_{\phi^*}(v) \sum_{x \in \mathcal{D}} \sum_{n \in I_r} \mathcal{N}(v; Q_n Ax, \sigma^2 Q_n Q_n^\top) \right)}_{b_r(v)} \quad (41)
$$

$$
- \underbrace{\sum_{x \in \mathcal{D}} \sum_{n \in I_r} x_n \mathcal{N}(v; Q_n Ax, \sigma^2 Q_n Q_n^\top)}_{a_r(v)} d\mathcal{H}^r(v). \quad (42)
$$

Choosing $f_{\phi^*}(v) = a_r(v)/b_r(v)$ on $\bar{E}_r$ cancels each integrand, hence the optimality condition Proposition A.5 holds. $\square$

**Remark B.6.** If all $Q_n$ have the same rank $r$, then $\bar{E}_r = \bigcup_{n=1}^N \text{Im}(Q_n)$ and the formula reduces to the constant-rank case in Proposition B.4, with a single ratio over $n = 1, \dots, N$.

**Remark B.7 (Singular limit and rank stratification).** Consider a regularization of $B$ as $B^{(\epsilon)} = U \Sigma_{\epsilon} V^{\top}$ and $Q_n^{(\epsilon)} = \Pi_n B^{(\epsilon)}$, where $B = U \Sigma V^{\top}$ is the SVD of $B$ and the diagonal matrix $\Sigma_{\epsilon}$ is constructed by adding a small positive number $\epsilon$ to the singular values in $\Sigma$. For $\epsilon > 0$ (and $M \ge P$), each $Q_n^{(\epsilon)}$ has full row rank $P$. For each $(n, x)$, define the probability measure $\mu_{n,x}^{(\epsilon)}$ on $\mathbb{R}^P$ with density (Radon–Nikodým derivative) with respect to Lebesgue measure $\lambda^P$ given by

$$
\frac{d\mu_{n,x}^{(\epsilon)}}{d\lambda^P}(v) = \mathcal{N}\left(v; Q_n^{(\epsilon)}Ax, \sigma^2 Q_n^{(\epsilon)}Q_n^{(\epsilon)T}\right).
$$

As $\epsilon \to 0$, $Q_n^{(\epsilon)} \to Q_n$ and some ranks $r_n = \text{rank}(Q)_n$ may drop. Then $\mu_{n,x}^{(\epsilon)}$ converges weakly to a probability measure $\mu_{n,x}$ supported on the linear subspace $\text{Im}(Q)_n$, which is absolutely continuous with respect to the Hausdorff measure $\mathcal{H}^{r_n}$ on $\text{Im}(Q)_n$, with density

$$
\frac{d\mu_{n,x}}{d\mathcal{H}^{r_n}}(v) = \mathcal{N}(v; Q_n Ax, \sigma^2 Q_n Q_n^\top), \quad v \in \text{Im}(Q_n),
$$

where, $\mathcal{N}(\cdot; \mu, \Sigma)$ denotes the degenerate Gaussian density Definition A.1 with respect to the appropriate Hausdorff measure on its support (and vanishes off that support).

For $\epsilon > 0$, the estimator reads (constant-rank case, similar to Proposition B.4)

$$
f_{\phi^*}^{(\epsilon)}(v) = \frac{\sum_{x \in \mathcal{D}} \sum_{n=1}^{N} x_n \mathcal{N}(v; Q_n Ax, \sigma^2 Q_n Q_n^\top)}{\sum_{x \in \mathcal{D}} \sum_{n=1}^{N} \mathcal{N}(v; Q_n Ax, \sigma^2 Q_n Q_n^\top)}.
$$
In the singular limit, for each $r$ and $\mathcal{H}^r$-a.e. $v \in \bar{E}_r$ with $b_r(v) > 0$,

$$ \lim_{\epsilon \to 0} f_{\phi^*}^{(\epsilon)}(v) = \frac{a_r(v)}{b_r(v)}, $$

i.e. the $\epsilon$-regularized estimator converges (stratum-wise, $\mathcal{H}^r$-a.e.) to the rank-stratified formula

$$ f_{\phi^*}(v) = \sum_{r=0}^{P} \frac{a_r(v)}{b_r(v)} \mathbb{1}_{\bar{E}_r}(v), $$

interpreted up to $\mathcal{H}^r$-null sets on each $\bar{E}_r$. If all $Q_n$ have the same rank $r$, only the stratum $\bar{E}_r$ is nonempty, and the above reduces to the constant-rank expression in Proposition B.4.

*Proof of i) of Corollary 3.9—Physics-agnostic LE-MMSE estimator.* When $B = I$, the weights of the LE-MMSE estimator in Theorem 3.8 simplifies to

$$ w_{n',n}(x|y) \propto \mathcal{N}(\Pi_{n'}y; \Pi_{n'}Ax, \sigma^2\Pi_n\Pi_n^\top) \propto \exp\left(-\frac{1}{2\sigma^2}\|y[\omega_{n'}] - (Ax)[\omega_n]\|^2\right) $$

Therefore, the LE-MMSE estimator at pixel $n'$ reads:

$$ \hat{x}_{\mathcal{T},\text{loc}}(y)[n'] = \frac{\sum_{x \in \mathcal{D}} \sum_{n=1}^{N} x_n \exp\left(-\frac{1}{2\sigma^2} \|y[\omega_{n'}] - (Ax)[\omega_n]\|^2\right)}{\sum_{x \in \mathcal{D}} \sum_{n=1}^{N} \exp\left(-\frac{1}{2\sigma^2} \|y[\omega_{n'}] - (Ax)[\omega_n]\|^2\right)}. $$

For small noise level $\sigma \to 0$, it returns the central pixel value of the patches in the dataset $\mathcal{D}$ whose degraded version $(Ax)[\omega_n]$ is closest to the observed patch $y[\omega_{n'}]$. Hence, the LE-MMSE estimator is a patch-work of training patches. $\square$

*Proof of ii) of Corollary 3.9—LE-MMSE is not a posterior mean.* We focus on the simplest case of denoising, that is $y = x + e$. Recall that the posterior mean satisfies the so-called Tweedie formula:

$$ \mathbb{E}[x|y] = y + \sigma^2 \nabla \log p_y(y). $$

Taking the gradient of both sides w.r.t $y$ yields:

$$ \nabla_y \mathbb{E}[x|y] = I + \sigma^2 \nabla^2 \log p_y(y). $$

Therefore, the Jacobian of any pure MMSE estimator (posterior mean) must be symmetric since the Hessian of $\log p_y$ is symmetric. To simplify the notation, we use $f$ instead of $\hat{x}_{\mathcal{T},\text{loc}}$ in what follows. Using the fact that the scaling by $\sigma^{-2}$ does not affect the symmetry (note that $p_y$ is the convolution of $p_x$ and a Gaussian so it's $\mathcal{C}^\infty$), if $\hat{x}_{\mathcal{T},\text{loc}}$ were a posterior mean, we must have:

$$ \frac{\partial}{\partial y_m} f_n(y) = \frac{\partial}{\partial y_n} f_m(y), \quad \forall n, m \text{ and } \forall y. \qquad (43) $$

For notation simplicity, we let $e_{n,l}(x, y) = \exp\left(-\frac{\|y[\omega_n] - x[\omega_l]\|^2}{2\sigma^2}\right)$. We have the weights of the LE-MMSE estimator in Theorem 3.8 becomes:

$$ w_{n,l}(x|y) = \exp\left(-\frac{\|y[\omega_n] - x[\omega_l]\|}{2\sigma^2}\right) / Z_n(y) = \frac{e_{n,l}(x, y)}{Z_n(y)}. $$

where $Z_n(y) = \sum_{x' \in \mathcal{D}} \sum_{l'=1}^{N} e_{n,l'}(x', y)$ is the normalization constant. Therefore, the LE-MMSE estimator at pixel $n$ reads:

$$ f_n(y) = \frac{1}{Z_n(y)} \sum_{x \in \mathcal{D}} \sum_{l=1}^{N} x_l \cdot e_{n,l}(x, y). $$

That is:

$$ f_n(y) = \frac{S_n(y)}{Z_n(y)} \quad \text{where} \quad S_n(y) = \sum_{x \in \mathcal{D}} \sum_{l=1}^{N} x_l \cdot e_{n,l}(x, y). $$
A key observation is that this estimator is real-analytic for $\sigma > 0$. To prove it, we remind that the product of polynomials and exponentials are real-analytic functions, and that the sum and quotient (with non-vanishing denominator) of real-analytic functions are also real-analytic. The denominator $Z_n(y)$ does not vanish for $\sigma > 0$, which concludes the proof of real-analyticity.

We will use the following classical result (see e.g. (Federer, 1996, Section 3.1.24) and (Mityagin, 2020) for an elementary proof):

**Lemma B.8 (Zeros of real-analytic functions).** Let $f: \mathbb{R}^p \to \mathbb{R}$ be a real-analytic function for some $p \in \mathbb{N}^*$. If $f$ is not identically zero, then the set of zeros of $f$ has Lebesgue measure zero.

Applying this lemma, it therefore suffices to find a dataset $\mathcal{D}$ and a point $y$ such that $\frac{\partial}{\partial y_m} f_n(y) - \frac{\partial}{\partial y_n} f_m(y) \neq 0$. To this end, we consider the case of overlapping patches $\omega_n \cap \omega_m \neq \emptyset$ with $n \neq m$, which is always possible for patch sizes $P \geq 2$. Now consider the case where $n \in \omega_m$ and $m \in \omega_n$ (they are equivalent).

The chain rule gives

$$ \frac{\partial e_{n,l}}{\partial y_m}(x, y) = -\sigma^{-2}(y_m - x_{l+m-n}) \cdot e_{n,l}(x, y) $$

and the quotient rule gives:

$$
\begin{align*}
\frac{\partial w_{n,l}}{\partial y_m}(x, y) &= \frac{\partial e_{n,l}}{\partial y_m}(x, y) \cdot \frac{1}{Z_n(y)} - e_{n,l}(x, y) \cdot \frac{\partial Z_n}{\partial y_m}(y) \cdot \frac{1}{Z_n(y)^2} \\
&= -\sigma^{-2}(y_m - x_{l+m-n}) \cdot \frac{e_{n,l}(x, y)}{Z_n(y)} - \frac{e_{n,l}(x, y)}{Z_n(y)^2} \cdot \frac{\partial Z_n}{\partial y_m}(y) \\
&= -\sigma^{-2}(y_m - x_{l+m-n}) \cdot w_{n,l}(x, y) - w_{n,l}(x, y) \cdot \frac{1}{Z_n(y)} \cdot \left( \sum_{x' \in \mathcal{D}} \sum_{l'=1}^{N} \frac{\partial e_{n,l'}}{\partial y_m}(x', y) \right) \\
&= -\sigma^{-2} \left( (y_m - x_{l+m-n}) \cdot w_{n,l}(x, y) + w_{n,l}(x, y) \cdot \frac{1}{Z_n(y)} \cdot \left( \sum_{x' \in \mathcal{D}} \sum_{l'=1}^{N} (y_m - x_{l'+m-n}) \cdot e_{n,l'}(x', y) \right) \right) \\
&= -\sigma^{-2} \left( (y_m - x_{l+m-n}) \cdot w_{n,l}(x, y) + w_{n,l}(x, y) \cdot \frac{1}{Z_n(y)} \cdot \left( y_m \cdot Z_n(y) - \sum_{x' \in \mathcal{D}} \sum_{l'=1}^{N} x_{l'+m-n} \cdot e_{n,l'}(x', y) \right) \right) \\
&= \sigma^{-2} \cdot \left( x_{l+m-n} \cdot w_{n,l}(x, y) - w_{n,l}(x, y) \cdot \frac{1}{Z_n(y)} \cdot \sum_{x' \in \mathcal{D}} \sum_{l'=1}^{N} x_{l'+m-n} \cdot e_{n,l'}(x', y) \right) \\
&= \sigma^{-2} \cdot \left( x_{l+m-n} \cdot w_{n,l}(x, y) - w_{n,l}(x, y) \cdot \sum_{x' \in \mathcal{D}} \sum_{l'=1}^{N} x_{l'+m-n} \cdot w_{n,l'}(x', y) \right) \\
&= \sigma^{-2} w_{n,l}(x, y) \left( x_{l+m-n} - \bar{x}_n^{(m-n)}(y) \right)
\end{align*}
$$

where we define $\bar{x}_n^{(k)}(y) = \sum_{x' \in \mathcal{D}} \sum_{l'=1}^{N} x_{l'+k} \cdot w_{n,l'}(x', y)$ the local posterior mean of pixel $n$ with offset $k$. Finally, differentiating $f_n(y)$ gives:

$$
\begin{align*}
\frac{\partial f_n}{\partial y_m}(y) &= \sum_{x \in \mathcal{D}} \sum_{l=1}^{N} x_l \cdot \frac{\partial w_{n,l}}{\partial y_m}(x, y) \\
&= \sigma^{-2} \sum_{x \in \mathcal{D}} \sum_{l=1}^{N} x_l \cdot w_{n,l}(x, y) \left( x_{l+m-n} - \bar{x}_n^{(m-n)}(y) \right)
\end{align*}
$$

Set $n = 2$ and $m = 3$. Set $\mathcal{D} = \{x, 0, \dots, 0\}$ to be dataset containing only one non-zero image $x$ and the rest are zeros. We
choose $x$ such that $x_2 > 0$ and $x_i = 0$ for $i \neq 2$. The derivatives simplify to:

$$ \frac{\partial f_2}{\partial y_3}(y) = \sigma^{-2}x_2w_{2,2}(x,y)(x_3 - \bar{x}_2^{(1)}(y)) = -\sigma^{-2}x_2^2w_{2,2}(x,y)w_{2,1}(x,y) $$

$$ \frac{\partial f_3}{\partial y_2}(y) = \sigma^{-2}x_2w_{3,2}(x,y)(x_1 - \bar{x}_3^{(-1)}(y)) = -\sigma^{-2}x_2^2w_{3,2}(x,y)w_{3,3}(x,y) $$

Fix an index $1 \le i \le N$ such that $y_i$ is a variable appearing in the vector $y[\omega_3]$ but not in the vector $y[\omega_2]$ (note that $i \neq 3$, otherwise $y_i = y_3$ is the center of the patch $y[\omega_3]$ and it's also in the patch $y[\omega_2]$ since the patch size $P > 1$). We have

$$ \frac{\partial w_{2,2}(x,y)}{\partial y_i} = \frac{\partial w_{2,1}(x,y)}{\partial y_i} = 0 \quad \text{and} \quad \frac{\partial^2 f_2}{\partial y_3 \partial y_i}(x,y) = 0. $$

On the other hand, since $i \neq 3$, we have

$$ \frac{\partial w_{3,2}(x,y)}{\partial y_i} = \sigma^{-2}w_{3,2}(x,y)(x_{2+i-3} - \bar{x}_3^{(i-3)}(y)) = -\sigma^{-2}x_2w_{3,2}(x,y)w_{3,2-(i-3)}(x,y) $$

$$ \frac{\partial w_{3,3}(x,y)}{\partial y_i} = \sigma^{-2}w_{3,3}(x,y)(x_i - \bar{x}_3^{(i-3)}(y)) = -\sigma^{-2}x_2w_{3,3}(x,y)w_{3,2-(i-3)}(x,y). $$

We obtain

$$ \frac{\partial^2 f_3(x,y)}{\partial y_2 \partial y_i} = \sigma^{-4}x_2^3 w_{3,2}(x,y)w_{3,2-(i-3)}w_{3,3}(x,y) + \sigma^{-4}x_2^3 w_{3,3}(x,y)w_{3,2-(i-3)}(x,y)w_{3,2}(x,y) > 0. $$

Therefore, for this value of $x$, $\frac{\partial f_3}{\partial y_2}$ and $\frac{\partial f_2}{\partial y_3}$ are independent, when seen as functions of $y$. We deduce that the equation $\frac{\partial f_2}{\partial y_3} = \frac{\partial f_3}{\partial y_2}$ is a nondegenerate analytic equation in variables $D, y$. Using Lemma B.8 and by Fubini's theorem, the set $\{D, \forall y, \frac{\partial f_2}{\partial y_3}(y) = \frac{\partial f_3}{\partial y_2}(y)\}$ also has zero measure.

In other words, for almost every dataset $D$, there exists $y$ such that the symmetry condition Equation (43) does not hold, and the LE-MMSE estimator is not a posterior mean. $\square$
<table><caption>Table 3: Details of neural network architectures with various receptive fields used in our experiments for images at 32 × 32 resolution.</caption><thead><tr><th></th><th>RECEPTIVE FIELD (PATCH SIZE)</th><th colspan="2">ARCHI. HYPER-PARAMETERS</th><th>NUM. PARAMETERS</th></tr></thead><tbody><tr><td rowspan="4">UNET2D</td><td></td><td>block_out_channels</td><td>kernel_size</td><td></td></tr><tr><td>5</td><td>(96, 224, 480)</td><td>(3, 1, 1, 1)-1</td><td>4.6M</td></tr><tr><td>7</td><td>(64, 192, 448)</td><td>(3, 1, 1, 1)-3</td><td>5.1M</td></tr><tr><td>9</td><td>(96, 192, 288)</td><td>(3, 1, 1, 1)-5</td><td>4.3M</td></tr><tr><td></td><td>11</td><td>(64, 96, 160, 288)</td><td>(3, 1, 1, 1, 3)-1</td><td>4.3M</td></tr><tr><td rowspan="4">RESNET</td><td></td><td colspan="2">num_res_blocks</td><td></td></tr><tr><td>5</td><td>1</td><td>640</td><td>4.5M</td></tr><tr><td>7</td><td>2</td><td>448</td><td>4.2M</td></tr><tr><td>9</td><td>3</td><td>384</td><td>4.5M</td></tr><tr><td></td><td>11</td><td>4</td><td>328</td><td>4.4M</td></tr><tr><td rowspan="4">PATCHMLP</td><td></td><td colspan="2">hidden_dim</td><td></td></tr><tr><td>5</td><td>7168</td><td>5</td><td>3.5M</td></tr><tr><td>7</td><td>6144</td><td>5</td><td>3.7M</td></tr><tr><td>9</td><td>5120</td><td>5</td><td>4.3M</td></tr><tr><td></td><td>11</td><td>3072</td><td>5</td><td>4.0M</td></tr></tbody></table>

# C. Numerical experiments

## C.1. Datasets

We use the following datasets in our experiments: FFHQ (Karras et al., 2019) downscaled to 32 × 32 or 64 × 64, CIFAR-10 (Krizhevsky et al., 2009) and Fashion-MNIST (Xiao et al., 2017). For each dataset, we randomly select 10,000 images for training, which is denoted by $\mathcal{D}$ in the main paper.

## C.2. Network architectures

For the local and translation equivariant estimator, we examine 3 different architectures with noise level $\sigma$ conditioning:

* UNet2D: we use a state-of-the-art UNet2D (Ronneberger et al., 2015) from the diffusers<fnref n="2" /> library. The model has several downsampling and upsampling layers and skip connections, with varying channel dimensions (defined by `block_out bidirectional` and `kernel size` for down-blocks and mid-blocks). The architecture is modified slightly (circular boundary conditions and kernel sizes of convolutional layers) to have a desired receptive field and to ensure translation equivariance. The noise level $\sigma$ is conditioned using a time embedding module with linear layers.

* ResNet: we use a minimal ResNet with residual block (He et al., 2016), with a 1 × 1 convolutional layer at the beginning and at the end, similar to (Kamb & Ganguli, 2025). Each residual block contains a convolutional layer with 3 × 3 kernels at a channel dimension defined by `num bidirectional channels`, followed by a batch normalization and a ReLU nonlinearity. The noise level $\sigma$ is conditioned using sine-cosine positional embedding.

* PatchMLP: a fully local MLP acting on patches. The MLP has 5 residual blocks, each containing two linear layers with hidden dimension of `hidden dim`, followed by a layer normalization and GELU activation. The noise level $\sigma$ is conditioned using sine-cosine positional embedding.

All convolutional layers use circular padding to ensure translation equivariance and they have approximately 4 million parameters. Details of the architectures with various receptive fields are provided in Table 3.

## C.3. Training procedure

All models are trained on a single NVIDIA A100 GPU with the following settings:

<footnote n="2">²https://huggingface.co/docs/diffusers/en/api/models/unet2d</footnote>
* Optimizer: Adam optimizer (Kingma & Ba, 2014)

* Learning rate: starting at $10^{-4}$ with cosine decay schedule and minimum learning rate at $10^{-6}$.

* Number of epochs: 600 for images at 32 × 32 resolution and 900 for images at 64 × 64 resolution.

* Batch size: 256

* Exponential moving average (EMA) with decay rate 0.99 from the 1000-th training step and updates every 5 steps. The EMA weights are used for evaluation.

### C.4. Forward operators

We consider the 3 representative inverse problems as forward operators A:

* Denoising: the forward operator is simply the identity $A = I$ and only the physics-agnostic estimator is applicable in this case.

* Inpainting: we consider the inpainting operator with a center square mask of size 15 × 15. The forward operator A is therefore a diagonal matrix with 0 on the masked pixels and 1 elsewhere. In this case, we define the physics-aware estimator with $B = A^+ = A$: it simply removes noise inside the masked region and keeps the observed pixels unchanged.

* Deconvolution: we consider an isotropic Gaussian blur kernel with standard deviation 1.0 with circular boundary conditions. The same kernel is used for all color channels. In this case, we build the full matrix A as a block-circulant matrix representing the convolution operation and B is its pseudo-inverse, which is computed once, in double precision.

### C.5. Analytical formula implementation

Implementing the analytical formulas of the E-MMSE Theorem 3.5 and the LE-MMSE Theorem 3.8 estimators requires computing many distance terms between images or patches in the dataset $\mathcal{D}$ and the observed measurement $y$. We process by batch to avoid memory overflow. For numerical stability and avoidance of overflow/underflow, the exponential terms are accumulated using an online log-sum-exp trick. All theoretical estimators are computed in an *exact* manner without any approximation, modular finite precision arithmetic. We use PyTorch (Paszke et al., 2019) for implementation and all computations are performed in single precision (FP32) on a single A-100 GPU, unless otherwise specified. Full implementation details are available at https://github.com/mh-nguyen712/analytical_mmse.
# D. Additional numerical results

## D.1. Comparison between trained neural networks and the analytical LE-MMSE estimator

We show in Figures 10 to 12 additional PSNR results between trained neural networks and the analytical LE-MMSE estimator for different architectures (UNet2D, ResNet, PatchMLP) on both training and test sets of various datasets. We recover a consistent conclusion across different settings: the trained neural networks closely approximate the analytical LE-MMSE estimator, with PSNR values exceeding 20 in most cases and often reaching above 30 dB.

<figure><img src="image_14.png" /><figcaption>Figure 10: Additional PSNR between trained UNet2D and the analytical formula of the LE-MMSE for different inverse problems on both training and test sets of various datasets. The patch size is $P = 11 \times 11$. Left: physics-agnostic estimator with $B=I$. Right: physics-aware estimator with $B=A^+$. We recover the same conclusions as in Figure 5.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 11: PSNR between trained ResNet and the analytical formula of the LE-MMSE for different inverse problems on both training and test sets of various datasets. The patch size is $P = 11 \times 11$. Left: physics-agnostic estimator with $B=1$. Right: physics-aware estimator with $B=A^+$. We recover the same conclusions across architectures and settings.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 12: PSNR between trained PatchMLP neural network and the analytical formula of the LE-MMSE for different inverse problems on both training and test sets of various datasets. The patch size is $P = 11 \times 11$. We recover the same conclusions as in Figures 5, 10 and 11. An exception is observed for the deconvolution task with physics-aware models (rightmost column). Here, the noise is amplified by the inversion of the blurring operator and the local PatchMLP architecture struggles to accurately reconstruct fine details, leading to a lower PSNR comparing to CNNs. We hypothesize that CNNs have other inductive biases that are more suited to handle such challenging tasks.</figcaption></figure>
## D.2. Addition qualitative comparison

We provide additional qualitative comparisons between the analytical LE-MMSE estimator and trained neural networks (UNet2D, ResNet, PatchMLP) for different patch sizes in Figures 13 to 16 on the FFHQ, CIFAR10 and FashionMNIST datasets. Overall, the neural networks closely match the analytical solution on both training and test sets across all datasets and inverse problems, confirming our theoretical findings.

<figure><img src="image_17.png" /><figcaption>Figure 13: Qualitative comparison. The patch size is $P = 5 \times 5$ and $B=1$. The noise level is $\sigma = 0.05, 0.2, 0.8$ from left to right for each column. The neural networks (UNet2D, ResNet, and PatchMLP) closely match the analytical LE-MMSE on both training and test sets across all datasets (FFHQ, CIFAR10, and FashionMNIST) and inverse problems, confirming our theoretical findings.</figcaption></figure>
An analytic theory of convolutional neural network inverse problems solvers

<figure><img src="image_18.png" /><figcaption>Figure 14: Qualitative comparison between the analytical formula LE-MMSE, UNet2D, ResNet, and PatchMLP, with $B=I$ on FFHQ, CIFAR10 and FashionMNIST. The patch size is $P = 7 \times 7$. The noise level is $\sigma = 0.05, 0.2, 0.8$ from left to right for each column. The neural networks closely match the analytical solution on both training and test sets across all datasets and inverse problems, confirming our theoretical findings. Yet, some discrepancies can be observed, especially at low noise levels on the test set, which we attribute to generalization issues discussed in Section 4.2.</figcaption></figure>
An analytic theory of convolutional neural network inverse problems solvers

<figure><img src="image_19.png" /><figcaption>Figure 15: Qualitative comparison between the analytical formula LE-MMSE, UNet2D, ResNet, and PatchMLP, with $B=I$ on FFHQ, CIFAR10 and FashionMNIST. The patch size is $P=9$. The noise level is $\sigma = 0.05, 0.2, 0.8$ from left to right for each column. The neural networks closely match the analytical solution on both training and test sets across all datasets and inverse problems, confirming our theoretical findings. Yet, some discrepancies can be observed, especially at low noise levels on the test set, which we attribute to generalization issues discussed in Section 4.2.</figcaption></figure>
An analytic theory of convolutional neural network inverse problems solvers

<figure><img src="image_20.png" /><figcaption>Figure 16: Qualitative comparison between the analytical formula LE-MMSE, UNet2D, ResNet, and PatchMLP, with $B=I$ on FFHQ, CIFAR10 and FashionMNIST. The patch size is $P=11$. The noise level is $\sigma = 0.05, 0.2, 0.8$ from left to right for each column. The neural networks closely match the analytical solution on both training and test sets across all datasets and inverse problems, confirming our theoretical findings. Yet, some discrepancies can be observed, especially at low noise levels on the test set, which we attribute to generalization issues discussed in Section 4.2.</figcaption></figure>

### D.3. Influence of the patch size (receptive field)

**Density of the patch distribution** We analyze here the influence of the patch size on the density of the patch distribution in FFHQ-32 dataset. We compute the negative-log-density of patches as a function of the patch size, for points on either the training set or the test set. The patch density is exactly the denominator term in the analytical LE-MMSE formula Equation (7). The results are shown in Figure 17.
<figure><img src="image_21.png" /><figcaption>Figure 17: The negative-log-density of patches in FFHQ-32 training set (top) and test set (bottom) as a function of the patch size. As the patch size increases, the density of patches decreases significantly, indicating a sparser coverage of the patch space by the dataset. In the training set, the negative-log-density is relatively lower (around $10^2$) than in the test set (around $10^4$ for small noise levels), indicating a better coverage of the patch space by the training set. This observation helps explain the drop in PSNR between trained neural networks and the analytical LE-MMSE estimator in the test set for low noise levels.</figcaption></figure>

**Patch size influence** In Figure 18, we show the PSNR between the trained UNet2D and the analytical LE-MMSE estimator for different patch sizes on both training and test sets of FFHQ-32 dataset. For low noise levels and in the test set, the PSNR decreases as the patch size increases, which we attribute to the sparser coverage of the patch space by the dataset for larger patches, as shown in Figure 17.

<figure><img src="image_22.png" /><figcaption>Figure 18: PSNR between trained UNet2D and the analytical formula of the LE-MMSE versus the patch size, on the FFHQ-32 dataset. Top: on the training set. Bottom: on the test set.</figcaption></figure>

In Figure 19, we show the PSNR between the analytical LE-MMSE estimator and the ground truth as a function of the patch size on both training and test sets of FFHQ-32 dataset. We observe that – on the test set – smaller patch size are preferable for low noise levels, while larger patch sizes yield better performance for high noise levels. The behavior on the training set is different: for low noise levels, the analytical formula copy-pastes exactly the right patches, explaining the blow up of PSNR at the origin.
<figure><img src="image_23.png" /><figcaption>Figure 19: PSNR between the analytical formula of the LE-MMSE and the ground truth versus the patch size, on the FFHQ-32 dataset. Top: on the training set. Bottom: on the test set.</figcaption></figure>

### D.4. Mass concentration

The LE-MMSE estimator Theorem 3.8 is a weighted average of the central pixel values of patches in the dataset. To understand the behavior of the estimator, we analyze how many patches contribute significantly to the estimate at each pixel location. In Figure 20, we show the number of patches contributing to 99% of the mass of the LE-MMSE estimator as a function of the noise level $\sigma$. We observe that for low noise levels, the estimator concentrates its mass on fewer patches (even nearest neighbors). As the noise level increases, the number of contributing patches increases significantly, with a critical value of $\sigma$ where the number of patches starts to increase rapidly. This behavior shows that although the MMSE estimator is typically associated with averaging, it can behave like a nearest-neighbor estimator when the noise is small, due to the strong concentration of mass on a very limited subset of patches. We also observe a clear difference between physics-agnostic and physics-informed settings for deconvolution, where the latter has a significantly higher number of contributing patches due to the amplification of noise by the pseudo-inverse of the blurring operator. We also visualize

<figure><img src="image_24.png" /><figcaption>Figure 20: For low noise levels, the LE-MMSE estimator concentrates its mass on fewer patches (even nearest neighbors). Here we show the median and IQR (over all pixels of 50 samples on the test set, for each $\sigma$) of the number of patches contributing to 99% of the mass of the LE-MMSE estimator. There is a significant increase in number of contributed patches as the noise level $\sigma$ increases, and a critical value of $\sigma$ where the number of patches starts to increase rapidly. The patch size is $P = 5 \times 5$, on FFHQ-32. Note that due to computational constraints, we only keep the top $10^4$ nearest patches (over $32 \times 32 \times 10^4 \approx 10^7$ patches) when computing mass concentration.</figcaption></figure>

in Figures 21 and 22 the $\log_{10}$ of number of patches contributing to 99% of the mass of the LE-MMSE estimator at each pixel location for different inverse problems on FFHQ-32 dataset. We observe that for low noise levels, the nearest patch (or a very small number of patches) contributes 99% of the mass at each pixel location.
<figure><img src="image_25.png" /><figcaption>Figure 21: Visualization of the number of patches (in log<sub>10</sub> scale) contributing to 99% of the mass of the physic-agnostic LE-MMSE estimator at each pixel location for different inverse problems on FFHQ-32 dataset.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure 22: Visualization of the number of patches (in $\log_{10}$ scale) contributing to 99% of the mass of the physic-informed LE-MMSE estimator at each pixel location for different inverse problems on FFHQ-32 dataset. For deconvolution, the number of contributing patches is significantly higher than for the physic-informed case, due to the amplification of noise by the pseudo-inverse of the blurring operator.</figcaption></figure>

### D.5. LE-MMSE is a patchwork

From the LE-MMSE formula Theorem 3.8, we see that the estimate at each pixel location is a weighted average of the central pixel values of patches in the dataset. In fact, the estimator often concentrates its mass on very few patches, as shown in Appendix D.4. Moreover, contiguous pixels may come from the central pixels of patches of the same image in the dataset, leading to a patchwork behavior.

We visualize in Figures 23 and 24 this behavior for different inverse problems on FFHQ-32 and FashionMNIST datasets. More precisely, at each pixel location, the source index is the index of the image in the dataset that contains the patch whose central pixel contributes at least 50% of the mass of the LE-MMSE estimator at that pixel location. We observe that, large contiguous regions of the estimate come from the same source image in the dataset, which we can interpret as the estimator to be a patchwork of training patches. As the noise level increases, more pixels become white, meaning that the corresponding pixel value is not mostly due to a single patch: the “patchwork effect” diminishes, as more patches contribute to the estimate at each pixel location.

Note that we use a patch size of $P = 11 \times 11$ to better visualize the patchwork behavior. This choice results in slightly weaker reconstruction performance, as discussed in Appendix D.3.
<figure><img src="image_27.png" /><figcaption>Figure 23: Illustration of the patchwork behavior of the LE-MMSE estimator (with $B=1$) on FFHQ-32 dataset for different inverse problems and noise levels. We use a patch size of $P = 11 \times 11$. We observe contiguous regions coming from the same source image in the dataset for low noise levels. White pixels correspond to locations where no patch contributes more than 50% of the mass.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 24: Similar illustration as in Figure 23 but on FashionMNIST dataset.</figcaption></figure>

## D.6. Out-of-distribution data

We provide additional results on out-of-distribution (OOD) data using UNet2D architecture with receptive field of size 5 × 5. The models and the formula are trained / evaluated on $\mathcal{D}$ = FFHQ-32. They are then evaluated on OOD dataset $\mathcal{D}'$ = CIFAR10. We show in Figure 25 the PSNR between trained UNet2D and the analytical LE-MMSE estimator. For large noise levels, as the Gaussians overlap significantly even on OOD data, the value of $-\log p(y)$ is low and the trained neural network closely matches the analytical LE-MMSE estimator with PSNR values above 25 dB. For low noise levels, the PSNR is about 3 dB lower than in the in-distribution case shown in Figure 5, which we attribute to the low-density of the measurement density $p(y)$, as discussed in Section 4.2 and can be seen in Figure 25.
<figure><img src="image_29.png" /></figure>

Figure 25: Comparison of UNet2D and the analytical LE-MMSE estimator on OOD dataset CIFAR10 when both are trained on FFHQ-32. Median and IQR using 50 images per $\sigma$, $P = 5 \times 5$ and $B=I$. On top: PSNR between UNet2D and the analytical LE-MMSE estimator. On bottom: negative-log-density of measurements $y$ in CIFAR10 under the measurement distribution induced by FFHQ-32.

## D.7. Dataset size influence

We analyze here the influence of the dataset size on the alignment between trained neural networks and the analytical LE-MMSE estimator. We train UNet2D models with receptive field of size $P = 5 \times 5$ and $B=I$ on various dataset sizes from $10^3$ to $5 \times 10^4$ images from FFHQ-32. The PSNR between trained UNet2D and the analytical LE-MMSE estimator is reported in Figure 26. We observe that the dataset size has limited influence on the alignment between trained neural networks and the analytical LE-MMSE formula, with a slight improvement when increasing the dataset size.

<figure><img src="image_30.png" /></figure>

Figure 26: Dataset size has limited influence on the alignment between neural networks and the LE-MMSE formula. Median and IQR using 50 images per $\sigma$, $P = 5 \times 5$ and $B=I$.
### D.8. Results on 3 × 64 × 64 images

We provide additional results on images at 3 × 64 × 64 resolution using UNet2D architecture with receptive field of size 11 × 11. The models are trained on 10⁴ images from FFHQ downscaled to 64 × 64, with the same training procedure as in Appendix C.3.

<table><caption>Table 4: Details of neural network architectures with various receptive fields used in our experiments for images at 3 × 64 × 64 resolution.</caption><thead><tr><th>RECEPTIVE FIELD (PATCH SIZE)</th><th colspan="2">ARCHI. HYPER-PARAMETERS</th><th>NUM. PARAMETERS</th></tr></thead><tbody><tr><td>UNET2D</td><td colspan="2">block_out bid channels kernel_size</td><td>13.3M</td></tr><tr><td></td><td colspan="2">(64, 128, 256, 512)</td><td>(3, 1, 1, 1, 3)-1</td></tr></tbody></table>

<figure><img src="image_31.png" /><figcaption>Figure 27: Additional qualitative comparison between UNet2D and the analytical LE-MMSE estimator on FFHQ-64 across tasks.</figcaption></figure>