# Equivariant Splitting: Self-supervised learning from
incomplete data

Victor Sechaud, Jérémy Scanvic, Quentin Barthélemy, Patrice Abry, Julián
Tachella

► To cite this version:

Victor Sechaud, Jérémy Scanvic, Quentin Barthélemy, Patrice Abry, Julián Tachella. Equivariant Splitting: Self-supervised learning from incomplete data. 2026. hal-05301498v2

HAL Id: hal-05301498

https://hal.science/hal-05301498v2

Preprint submitted on 7 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# EQUIVARIANT SPLITTING: SELF-SUPERVISED LEARNING FROM INCOMPLETE DATA

Victor Sechaud*<fnref n="1" />, Jérémy Scanvic*<fnref n="1" />,², Quentin Barthélemy², Patrice Abry<fnref n="1" />, and Julián Tachella<fnref n="1" />

¹LPENSL, CNRS, ENS de Lyon, France

²Prysm, Lyon, France

## ABSTRACT

Self-supervised learning for inverse problems allows to train a reconstruction network from noise and/or incomplete data alone. These methods have the potential of enabling learning-based solutions when obtaining ground-truth references for training is expensive or even impossible. In this paper, we propose a new self-supervised learning strategy devised for the challenging setting where measurements are observed via a single incomplete observation model. We introduce a new definition of equivariance in the context of reconstruction networks, and show that the combination of self-supervised splitting losses and equivariant reconstruction networks results in unbiased estimates of the supervised loss. Through a series of experiments on image inpainting, accelerated magnetic resonance imaging, sparse-view computed tomography, and compressive sensing, we demonstrate that the proposed loss achieves state-of-the-art performance in settings with highly rank-deficient forward models<fnref n="1" />.

## 1 INTRODUCTION

Inverse problems are ubiquitous in many sensing and imaging applications. They are written as

$$y = Ax + \varepsilon \quad (1)$$

where $A \in \mathbb{R}^{m \times n}$ is the known forward matrix, $x \in \mathbb{R}^n$ is the ground truth image to be estimated, $y \in \mathbb{R}^m$ is the observed measurement vector, and $\varepsilon \in \mathbb{R}^m$ is the unknown noise, generally assumed to follow a Gaussian distribution. This model is suitable for many imaging modalities, including magnetic resonance imaging (MRI) for medical imaging (Zbontar et al., 2019), computed tomography (CT) (Withers et al., 2021), microscopy (Ragone et al., 2023), remote sensing (Fassnacht et al., 2024), and astronomical imaging (Vojtekova et al., 2021).

The number of linearly independent measurements is often smaller than the number of pixels in the target images due to physical and practical constraints. In this case, the forward matrix is rank-deficient, the forward process discards some of the information present in the target image. Further information about the target images is needed to solve the problem and different methods make different assumptions about the signal distribution.

Modern learning-based solvers generally obtain state-of-the-art performance by training on supervised pairs of ground-truth images and measurements. For certain applications, it is expensive or even impossible to obtain enough ground-truth data for supervised training (Belthangady & Royer, 2019). This is notably the case in astronomical imaging, microscopy and medical imaging.

Recent self-supervised methods overcome this limitation by learning to reconstruct without ground-truth data, relying only on a dataset of measurements, and they generally differ in the assumption they make on the forward model. For denoising problems where the forward matrix is the identity mapping, certain methods rely on knowing the exact noise distribution (Eldar, 2009; Pang et al., 2021; Monroy et al., 2025), others only assume it is entry-wise independent (Krull et al., 2019), while some make intermediate assumptions (Tachella et al., 2025a).

<footnote>*The first two authors have contributed equally to this paper.</footnote>

<footnote n="1">¹The code is available at https://github.com/vsechaud/Equivariant-Splitting</footnote>
In settings where measurements are observed via multiple incomplete forward operators, such as accelerated MRI with masks varying across acquisitions or inpainting problems with missing pixels varying across images, the main approaches are splitting (Yaman et al., 2020; Millard & Chiew, 2023) and consistency across operators (Tachella et al., 2022). Splitting losses divide the measurements into input and target components, and are unbiased estimators of the supervised loss if there is enough diversity of operators in the dataset (Daras et al., 2023; Millard & Chiew, 2023).

In the more challenging case of measurement data obtained via a single incomplete forward operator, the main self-supervised approach is equivariant imaging (Chen et al., 2021; 2022) which makes the assumption that the target distribution is invariant to a certain group of transformations including geometric transformations (Wang & Davies, 2024; Scanvic et al., 2025) and range transformations such as intensity scalings (Sechaud et al., 2024). Experimental results show that it can obtain competitive performances to supervised methods even though it does not require ground-truth references for training (Chen et al., 2022). However, training with EI is typically slower than supervised learning, as it requires two to three evaluations of the model for every iteration, and can obtain subpar performances if the operator is highly incomplete.

In this work, we propose equivariant splitting (ES), a new self-supervised method for learning from measurements obtained via a single forward operator that combines the invariance to transformations assumption of equivariant imaging and the simplicity and computational efficiency of splitting methods. The key idea behind our method is the use of recent developments in equivariant architectures (Chaman & Dokmanić, 2021a; Puny et al., 2022) to design a training loss that performs implicit ground truth data augmentation without any transformation overhead. Our theoretical results show that our method yields (in expectation) the minimum mean squared error (MMSE) estimator as long as the model is expressive enough, which is not guaranteed for equivariant imaging and previous splitting methods. Additionally, our experiments demonstrate state-of-the-art performance on a wide range of self-supervised imaging problems including compressed sensing, image inpainting, accelerated MRI and sparse-view CT. This work is additional evidence that architectural constraints built upon equivariance are a powerful tool to solve ill-posed imaging inverse problems.

Our contributions are the following:

1. We propose a new definition for equivariance in inverse problems, and propose architectures that satisfy this property, including unrolled architectures.

2. We propose a new self-supervised loss that leverages equivariant networks (according to our definition above) whose global minimizer is the gold standard MMSE estimator under the assumption of an invariant signal distribution.

3. We demonstrate the performance of our method on a wide range of inverse problems including compressed sensing, inpainting, accelerated MRI and sparse-view computed tomography.

## 2 RELATED WORK

**Measurement splitting** Various self-supervised losses consist of dividing measurement vectors into two components, one used as input and the other as target. It has been used to solve oversampled single-operator inverse problems including full-view CT (Hendriksen et al., 2020), and also undersampled multi-operator inverse problems with theoretical guarantees of producing estimates equivalent to supervised estimates in expectation, notably accelerated MRI and image inpainting with varying masks (Daras et al., 2023; Millard & Chiew, 2023). To the best of our knowledge, this work is the first to extend these methods to the more challenging single-operator undersampled setting, which notably includes sparse-view CT and accelerated MRI (with fixed mask).

**Equivariant imaging** It is possible to learn from incomplete data associated with a single degradation operator, as long as the underlying signal distribution is invariant to a group of transformations (Tachella et al., 2023). Equivariant imaging (Chen et al., 2021; 2022) assumes that the distribution of clean images remains unchanged under certain transformations, including translations, rotations and flips, and introduces a training loss that enforces the equivariance to these transformations of the entire measurement-reconstruction process, thereby constraining the set of learnable models. It has been used effectively to solve various inverse problems using adequate groups of transformations (Wang & Davies, 2024; Sechaud et al., 2024), but it is computationally expensive
due to the two to three network evaluations. Moreover, the EI loss is not necessarily an unbiased estimator of the supervised loss, and it is unclear whether this approach recovers the optimal MMSE estimator in expectation.

**Equivariant neural networks** The design of equivariant networks is an active research topic and many different approaches exist. Some rely on data augmentation to make neural networks more equivariant to rotations and flips, while other rely on averaging on the group of transformations at test time (Rivera et al., 2021; Puny et al., 2022; Kaba et al., 2023; Sannai et al., 2024). A third line of work focuses on architectures that are equivariant by design: Cohen & Welling (2016; 2017) propose convolutional layers that are equivariant to rotations and flips, and other works also rely on the design of more equivariant layers to improve translation-equivariance. Zhang (2019); Chaman & Dokmanić (2021b) propose equivariant pooling and downsampling/upsampling layers, and Karras et al. (2021); Michaeli et al. (2023) equivariant non-linearities and activation layers. Closer to our work, in the specific setting of inverse problems, Celledoni et al. (2021) propose the use of equivariant denoising blocks within unrolled architectures, and Terris et al. (2024) similarly propose to render plug-and-play denoisers more equivariant by averaging over transformations at test time. We take these analyses further by providing a clear definition of equivariance in inverse problems and showing which popular posterior estimators and related architectures verify these properties.

# 3 BACKGROUND

Solving the inverse problem in eq. (1) amounts to designing a reconstruction function $f(y, A) \approx x$ estimating a ground truth signal $x \in \mathbb{R}^n$ from its measurement vector $y \in \mathbb{R}^m$ and forward matrix $A \in \mathbb{R}^{m \times n}$. In practice, it is often implemented as a neural network parametrized by a set of weights. Supervised methods assume the existence of a finite dataset containing pairs of ground truth signals and measurements $\{(x_i, y_i)\}_{i \in \mathcal{I}}$ that are used to learn the reconstruction function $f(y, A)$. The main approach is to minimize a training loss equal to the mean squared error (Ongie et al., 2019)

$$ \min_{f} \left\{ \frac{1}{|\mathcal{I}|} \sum_{i \in \mathcal{I}} \mathcal{L}_{\text{SUP}}(x_i, y_i, A, f) \right\}, \quad \mathcal{L}_{\text{SUP}}(x, y, A, f) = \|f(y, A) - x\|^2. \quad (2) $$

While this approach obtains state-of-the-art performance, it cannot be used in the absence of ground-truth data. Self-supervised methods overcome this limitation using a finite dataset containing only measurements $\{y_i\}_{i \in \mathcal{I}}$, and a training loss $\mathcal{L}(y, A, f)$ that need no ground truth data to be evaluated, and which is designed to approximate well the supervised objective in eq. (2).

In denoising problems ($A = I$) with Gaussian noise of known variance, Recorrupted2Recorrupted (R2R) (Pang et al., 2021) and SURE (Metzler et al., 2020) provide unbiased estimators of the supervised loss. The R2R loss is computed by adding synthetic Gaussian noise $\omega \sim N(0, \sigma^2 I)$ to the measurements, creating input-target pairs as $(y + \alpha\omega, y - \frac{\omega}{\alpha})$ for some $\alpha \in (0, +\infty)$. However, if measurements are observed by an incomplete operator $A$ with a non-trivial nullspace, these losses fail to learn in the nullspace, approximating $f(y, A) = A^\dagger A \mathbb{E}_{x|y,A}\{x\} + v(y, A)$, with $v$ being any function taking values in the nullspace of $A$ (Chen et al., 2021).

In the rest of this section, we present the two main self-supervised losses that can learn beyond the nullspace, measurement splitting in Section 3.1 and equivariant imaging in Section 3.2. By combining their main ideas, we obtain our new self-supervised loss introduced in Section 4.

## 3.1 MEASUREMENT SPLITTING

One strategy to address the limitations imposed by the nullspace of a single operator is to employ multiple operators (Millard & Chiew, 2023). The key idea is that different operators generally do not share the same nullspace; thus, observing measurements through the image spaces of multiple operators allows access to the whole space $\mathbb{R}^n$. Formally, they assume that measurements are obtained according to $y \sim p(y|Ax)$ where the measurement operator $A$ is itself drawn from a distribution $p(A)$, and differ for each acquisition.

Splitting losses divide the measurements into two components $y = [y_1^\top, y_2^\top]^\top$ with corresponding operators $A = [A_1^\top, A_2^\top]^\top$, where $A_1 \in \mathbb{R}^{m_1 \times n}$ and $A_2 \in \mathbb{R}^{m_2 \times n}$ with $m_1 + m_2 = m$.
The network is then trained to predict the entire measurements vector $\mathbf{y}$ from only one of its two components $\mathbf{y}_1$, using the training loss<fnref n="2" />

$$ \mathcal{L}_{\text{SPLIT}}(\mathbf{y}, \mathbf{A}, f) = \mathbb{E}_{\mathbf{y}_1, \mathbf{A}_1 | \mathbf{y}, \mathbf{A}} \{ \| \mathbf{A} f(\mathbf{y}_1, \mathbf{A}_1) - \mathbf{y} \|^2 \}, \quad (3) $$

where $p(\mathbf{y}_1, \mathbf{A}_1 | \mathbf{y}, \mathbf{A})$ is a random splitting distribution, which is chosen on a per-problem basis. The loss encourages the model to learn in the nullspace of each operator by predicting the unobserved part. In practice, the expectation in eq. (3) is estimated using a single split $\mathbf{y}_1$ for each training batch.

## 3.2 EQUIVARIANT IMAGING

Equivariant imaging (EI) relies on the assumption that the distribution of images is invariant under a group of transformations $T_g \in \mathbb{R}^{n \times n}$, $g \in \mathcal{G}$, to learn beyond the nullspace of measurements obtained via a single operator (Chen et al., 2021; 2022). In this setting, the reconstruction model is expected to be able to estimate ground truth images $\mathbf{x}$ as well as their transformations $T_g \mathbf{x}$ in a coherent manner, i.e., such that the entire measurement-reconstruction pipeline $f(\mathbf{A} \mathbf{x})$ is equivariant with respect to the transformations $f(\mathbf{A} T_g \mathbf{x}, \mathbf{A}) = T_g f(\mathbf{A} \mathbf{x}, \mathbf{A})$. In order to achieve this, they propose a self-supervised loss which consists in a traditional measurement consistency term (replaced by SURE in the presence of noise), along with an equivariance-promoting term

$$ \mathcal{L}_{\text{EI}}(\mathbf{y}, \mathbf{A}, f) = \| \mathbf{A} f(\mathbf{y}, \mathbf{A}) - \mathbf{y} \|^2 + \lambda \mathbb{E}_g \{ \| T_g f(\mathbf{y}, \mathbf{A}) - f(\mathbf{A} T_g f(\mathbf{y}, \mathbf{A}), \mathbf{A}) \|^2 \}, \quad (4) $$

where $\lambda > 0$ is a trade-off coefficient. Even though it has been shown to be particularly effective on a wide variety of problems (Wang & Davies, 2024; Sechaud et al., 2024), it typically requires from two to three evaluations of the neural network which makes it very computationally expensive (Xu et al., 2025). Moreover, the equivariant loss is only effective at enforcing equivariance when the learned estimator achieves almost perfect reconstructions, $f(\mathbf{y}, \mathbf{A}) \approx \mathbf{x}$ (Chen et al., 2021), which is not the case for very ill-conditioned problems and which leads to the method having multiple possible solutions in general.

## 4 METHOD

In this section, we present our method that combines measurement splitting and EI introduced in Section 3. In order to learn from incomplete measurements obtained via a single operator, we rely on the same assumption of EI:

**Assumption 1.** The distribution of ground truth images $p(\mathbf{x})$ is invariant to the transformations $\{T_g\}_{g \in \mathcal{G}}$

$$ p(T_g \mathbf{x}) = p(\mathbf{x}), \quad \forall g \in \mathcal{G}, \forall \mathbf{x} \in \mathbb{R}^n. \quad (5) $$

The set of transformations is a design choice of the method to be chosen on a per-problem basis. Corollary 1 helps to choose them for specific problems, and we use it in the experiments. This assumption applies in many different settings, natural image distributions are generally invariant to rigid transformations (translations, rotations, flips, scalings), especially microscopic, aerial or remote sensing images which have no privileged orientation. Moreover, our experiments show that our method performs well even in the presence of approximate invariance, e.g., for medical images. Refer to Appendix A.2 for more details about how to choose the transformations for a given application.

In Section 4.1, we present our proposed loss and state its optimality under mild assumptions in Theorem 1 and Proposition 1. In Section 4.2, we present a new definition of equivariance for reconstruction functions and state sufficient conditions for common architectures to be equivariant in Theorem 2. In Section 4.3, we finally present a computational synergy between our loss and these equivariant architectures stated in Theorem 3. See Appendix A.3 for a description of the end-to-end algorithm, and Appendix C for the detailed proofs.

<footnote n="2">²We use the notation $\mathbb{E}_{a|b}\{g(a)\}$ for $\int_a g(a)p(a|b)da$.</footnote>
## 4.1 PROPOSED LOSS

Under Assumption 1, the measurements can be understood in a different way than they are traditionally. Indeed, measurements $y$ are generally thought of as being associated to the ground truth image $x$ and the forward matrix $A$, but they can equally be understood as being associated to the virtual ground truth image $x_g = T_g^{-1}x$ and the virtual forward matrix $A_g = AT_g$, with

$$y = Ax + \varepsilon = AT_g T_g^{-1}x + \varepsilon = A_g x_g + \varepsilon. \quad (6)$$

The implicit multi-operator structure (with operators {$AT_g$}$_{g \in \mathcal{G}}$) of the problem hints that we should be able to leverage the splitting approaches presented in Section 3. Moreover, since we use the same invariance assumption as EI, we can combine the two approaches to obtain equivariant splitting (ES), a new self-supervised loss $\mathcal{L}_{\text{ES}}$ that has the advantages of both methods.

**Noiseless measurements** The ES self-supervised loss is expressed as

$$\mathcal{L}_{\text{ES}}(y, A, f) \triangleq \mathbb{E}_g \{ \mathcal{L}_{\text{SPLIT}}(y, AT_g, f) \} \quad (7)$$

$$= \mathbb{E}_g \{ \mathbb{E}_{y_1, A_1 | y, AT_g} \{ \| AT_g f(y_1, A_1) - y \|^2 \} \} \quad (8)$$

where $A_1 \sim p(A_1|AT_g) \triangleq p(A_1|g)$ is a random splitting of $AT_g$.

**Theorem 1.** In the case of noiseless measurements with $p(x) \mathcal{G}$-invariant (Assumption 1), if the matrix $Q_{A_1} \triangleq \mathbb{E}_{g|A_1} \{(AT_g)^\top AT_g\} \in \mathbb{R}^{n \times n}$ has full rank for some split $A_1$, then the splitting method yields the same MMSE-optimal reconstructions as the supervised method, i.e.,

$$f^*(y_1, A_1) = \mathbb{E}_{x|y_1, A_1} \{x\}. \quad (9)$$

While it is sufficient for the matrix $Q_{A_1}$ to be invertible for the reconstructions to be almost or exactly optimal, in practice the spectrum (number of its non-negligible eigenvalues) of the matrix determines how close to optimal they are. In the experiments, we use a single (random) split per batch element and we average the reconstructions corresponding to 10 splits at inference.

**Proposition 1.** If the matrix $\bar{Q}_A \triangleq \mathbb{E}_{A_1|A} \{Q_{A_1}\} \in \mathbb{R}^{n \times n}$ is invertible and $f$ minimizes $\mathbb{E}_y \{\mathcal{L}_{\text{ES}}(y, A, f)\}$. Then the reconstruction function

$$\bar{f}(y, A) \triangleq \mathbb{E}_{y_1, A_1 | y, A} \{ \bar{Q}_A^{-1} Q_{A_1} f(y_1, A_1) \} \quad (10)$$

satisfies

$$\bar{f}(y, A) = \mathbb{E}_{y_1, A_1 | y, A} \{ \bar{Q}_A^{-1} Q_{A_1} \mathbb{E}_{x|y_1, A_1} \{x\} \}. \quad (11)$$

where eq. (11) is a convex combination of MMSE estimators for different splittings.

In practice, often neither $\bar{Q}_A$ nor $Q_{A_1}$ can be computed in closed-form, and we use a non-weighted average over random splittings

$$\bar{f}(y, A) := \frac{1}{J} \sum_{j=1}^{J} f(y_1^{(j)}, A_1^{(j)}) \quad \text{with} \quad (y_1^{(j)}, A_1^{(j)}) \sim p(y_1 A_1 | y, AT_g) \quad (12)$$

where $g$ is chosen randomly over the group of transformations for each split.

As with EI, the forward operator should not be equivariant with respect to the choice of transformations, in order to learn beyond the nullspace of the operator:

**Corollary 1.** In order for the matrices $Q_{A_1}$ or $\bar{Q}_A$ to have full rank, it is necessary that $A$ is not equivariant:

$$\exists g \in \mathcal{G}, AT_g \neq T_g A. \quad (13)$$

**Noisy measurements** The ES loss can be split into two separate terms, one enforcing measurement consistency, and the other prediction accuracy:

$$\mathcal{L}_{\text{ES}}(y, A, f) = \mathbb{E}_g \{ \mathbb{E}_{y_1, A_1 | y, AT_g} \{ \| A_1 f(y_1, A_1) - y_1 \|^2 + \| A_2 f(y_1, A_1) - y_2 \|^2 \} \},$$
where $A_1$ and $A_2$ are a splitting of $AT_g$. If the measurements are noisy, the first term can be replaced by a self-supervised denoising loss. In particular, if measurements are corrupted by Gaussian noise of standard deviation $\sigma$, we replace the first term by the R2R loss, yielding:

$$ \mathcal{L}_{\text{G-ES}}(\mathbf{y}, \mathbf{A}, f) = \mathbb{E}_{g, y_1, A_1, \omega | \mathbf{y}, \mathbf{A} T_g} \left\{ \left\| A_1 f(\mathbf{y}_1 + \alpha \omega, A_1) - \left( \mathbf{y}_1 - \frac{\omega}{\alpha} \right) \right\|^2 \right. \\ \left. + \left\| A_2 f(\mathbf{y}_1 + \alpha \omega, A_1) - \mathbf{y}_2 \right\|^2 \right\} $$

with $\omega \sim \mathcal{N}(0, \sigma^2 I)$ and a hyper-parameter $\alpha \in (0, +\infty)$. Since R2R provides an unbiased estimate of the clean measurement consistency term (Pang et al., 2021), we can apply Theorem 1 to show that minimizing this loss (in expectation) also results in MMSE estimators (if the conditions on $Q_{A_1}$ or $Q_A$ are verified). In the case of non-Gaussian noise, the R2R loss can be replaced by its non-Gaussian extension (Monroy et al., 2025). As with random splits, the expectation over $\omega$ is computed using a random realization per batch. At test time, we modify eq. (12) to average over both splits and synthetic noise additions.

## 4.2 EQUIVARIANT RECONSTRUCTORS

The ES loss requires a model evaluation for every mask and transformation. We show that, instead of sampling a random transformation each evaluation, imposing architectural equivariance constraints removes the need to explicitly compute the transforms.

Image-to-image functions $\phi(\mathbf{x})$ are equivariant if they satisfy (Cohen & Welling, 2016)

$$ \phi(T_g \mathbf{x}) = T_g \phi(\mathbf{x}), \quad \forall \mathbf{x} \in \mathbb{R}^n, \forall g \in \mathcal{G}. \qquad (14) $$

In this work, we introduce an extension of this definition to reconstruction functions $f(\mathbf{y}, \mathbf{A})$. To the best of our knowledge, this is the first work that introduces this definition.

**Definition 1.** We say that the reconstruction function $f(\mathbf{y}, \mathbf{A})$ is an equivariant reconstructor if

$$ f(\mathbf{y}, \mathbf{A} T_g) = T_g^{-1} f(\mathbf{y}, \mathbf{A}), \quad \forall \mathbf{y} \in \mathbb{R}^m, \forall g \in \mathcal{G}, \forall \mathbf{A} \in \mathbb{R}^{m \times n}. \qquad (15) $$

This property is very general and the class of classical reconstruction functions that satisfy it is large.

**Theorem 2.** The reconstruction functions defined in points below are all equivariant as in eq. (15).

1. **Artifact removal network.** For a denoiser $\phi(\mathbf{x})$ equivariant in the sense of eq. (14),

$$ f(\mathbf{y}, \mathbf{A}) = \phi(\mathbf{A}^{\top} \mathbf{y}), \quad \text{or} \quad f(\mathbf{y}, \mathbf{A}) = \phi(\mathbf{A}^{\dagger} \mathbf{y}). \qquad (16) $$

2. **Unrolled network.** For $\phi(\mathbf{x})$ equivariant, any $\gamma \in \mathbb{R}$ and data fidelity $d(\mathbf{A}\mathbf{x}, \mathbf{y})$, with

$$ \mathbf{x}_0 = \mathbf{0}, \quad \mathbf{x}_{k+1} = \phi(\mathbf{x}_k - \gamma \nabla_{\mathbf{x}_k} d(\mathbf{A}\mathbf{x}_k, \mathbf{y})) \qquad (17) $$

for $k = 0, \dots, L-1$ and $f(\mathbf{y}, \mathbf{A}) = \mathbf{x}_L$.

3. **Reynolds averaging.** For a possibly non-equivariant reconstructor $r(\mathbf{y}, \mathbf{A})$, with

$$ f(\mathbf{y}, \mathbf{A}) = \frac{1}{|\mathcal{G}|} \sum_{g \in \mathcal{G}} T_g r(\mathbf{y}, \mathbf{A} T_g). \qquad (18) $$

4. **Maximum a posteriori (MAP).** For a distribution $p(\mathbf{x})$ invariant as in eq. (5), with

$$ f(\mathbf{y}, \mathbf{A}) = \underset{\mathbf{x} \in \mathbb{R}^n}{\operatorname{argmax}} \left\{ p(\mathbf{x} | \mathbf{y}, \mathbf{A}) \right\}. \qquad (19) $$

5. **Minimum mean squared error (MMSE).** For a distribution $p(\mathbf{x})$ invariant as in eq. (5),

$$ f(\mathbf{y}, \mathbf{A}) = \mathbb{E}_{\mathbf{x} | \mathbf{y}, \mathbf{A}} \{ \mathbf{x} \}. \qquad (20) $$

For additional motivation and details about these reconstructor architectures, see Appendix A.
## 4.3 EFFICIENT LOSS EVALUATION WITH EQUIVARIANT RECONSTRUCTORS

For equivariant reconstructors, the ES loss in eq. (7) reduces to the splitting loss in eq. (3).

**Theorem 3.** If $f(\mathbf{A}, \mathbf{x})$ is an equivariant reconstructor, then ES is equivalent to the splitting loss

$$ \mathcal{L}_{\text{ES}}(\mathbf{y}, \mathbf{A}, f) = \mathcal{L}_{\text{SPLIT}}(\mathbf{y}, \mathbf{A}, f). \qquad (21) $$

We emphasize that the condition for a reconstructor to be equivariant is different from the condition enforced by the EI loss, i.e., $f(\mathbf{A}T_g\mathbf{x}, \mathbf{A}) = T_g f(\mathbf{A}\mathbf{x}, \mathbf{A})$. They are only equivalent if the reconstruction function is a perfect one-to-one mapping over all possible images.

In our experiments, we build equivariant reconstructors using i) artifact removal networks with a translation equivariant UNet denoiser Chaman & Dokmanić (2021b) and ii) unrolled networks with a denoiser architecture equivariant to rotations and flips via averaging (Sannai et al., 2021).

# 5 EXPERIMENTS

We assess the effectiveness of the proposed self-supervised loss using experiments conducted on different inverse problems. For each experiment, we train a model corresponding to our method as well as baseline methods and compare their performance. The inverse problems we consider are 1) inpainting, 2) compressive sensing, 3) accelerated MRI and 4) sparse-view CT. We also validate our theoretical predictions by testing the effect of using an equivariant architecture. In Section 5.1 we detail our experiments on compressive sensing, in Section 5.2 on image inpainting, in Section 5.3 on accelerated MRI and in Section 5.4 on sparse-view CT, and in Section 5.5 we present an ablation study on the effect of equivariant architectures. For additional details about the experiments, see Appendix B.

In each experiment, we compare against multiple baselines: a supervised baseline using the supervised loss described in eq. (2), the EI (Chen et al., 2021) baseline described in eq. (4), a measurement consistency baseline (MC and SURE) Eldar (2009) and a learning-free baseline which are either the measurements directly, or their image under the adjoint or the pseudo-inverse of the forward operator. We use the same architecture and the same training procedure for every method to ensure a fair comparison. We report the peak signal-to-noise ratio (PSNR) and the structural similarity index measure (SSIM) (Wang et al., 2004) of the final reconstructions for the different methods. They are distortion metrics indicating how close the reconstructions are to the ground truth images. We do not include perception metrics which are known to be at odds with them (Blau & Michaeli, 2018). In each case, we also compute equivariance metrics (EQUIV) for translations or rotations and flips.

We design a model architecture from the principles introduced in Section 4, with a variant equivariant to shifts, one equivariant to rotations and flips, and one without equivariance. It uses an existing unrolled architecture (Aggarwal et al., 2019) with a prior step implemented as a standard UNet (Ronneberger et al., 2015), or that of Chaman & Dokmanić (2021b) to enforce the equivariance to shifts. It also optionally uses Reynolds’ averaging to enforce the equivariance to rotations and flips. We use a single equivariant variant dictated by Corollary 1 for each problem, the shift-equivariant one for compressive sensing and inpainting, and that equivariant to rotations and flips for MRI and CT.

## 5.1 COMPRESSIVE SENSING

The 28 × 28 ground truth images are obtained from the MNIST dataset and are measured without additional noise through compressive matrices $\mathbf{A} \in \mathbb{R}^{m \times n}$, with $m < n$ varying from training to training to assess the impact of the compression rate. These matrices are obtained by sampling $A_{i,j} \sim \mathcal{N}(0, 1/m)$ for $i = 1, \dots, m$ and $j = 1, \dots, m$, where $n = 28 \times 28$. Figure 1 shows the performance of the different methods for the different compression rates. Our method performs almost as well as the supervised baseline, while the equivariant imaging baseline performs close to the supervised baseline only for higher compression rates.

## 5.2 IMAGE INPAINTING

The dataset consists of 128 × 128 images from DIV2K (Agustsson & Timofte, 2017) measured without additional noise through a single subsampling matrix $\mathbf{A} \in \mathbb{R}^{m \times n}$ selecting about 30% of
the pixels. This matrix is obtained by sampling $a_1, ..., a_n \sim B(0.3)$ where $n = 3 \times 128 \times 128$ and letting $A_{i,j} = \delta_{j_i, j}$ for $i = 1, ..., m$ and $j = 1, ..., n$, where $m \approx 0.3n$ is the number of nonzero values in $a$, and $j_i$ is the $i$-th index in $a$ corresponding to a nonzero value. Among the 900 images in the dataset, 800 are used for training while the remaining 100 are used for testing. For the supervised method, we use different crops at each evaluation. Table 2 and Figure 2 show that ES performs almost as well as the supervised baseline, and better than EI.

<table><caption>Table 1: **Medical imaging results.** ES (ours) performs better than EI, SURE and MC (baselines), while performing almost as well as the supervised baseline in reconstruction quality (PSNR, SSIM) and measured equivariance (EQUIV). In **bold**, the best self-supervised metrics (avg ± st.d.).</caption><thead><tr><th rowspan="2">Method</th><th colspan="3">MRI (x8 Accel., 40 dB SNR)</th></tr><tr><th>PSNR ↑</th><th>SSIM ↑</th><th>EQUIV ↑</th></tr></thead><tbody><tr><td>Supervised</td><td>28.74 ± 2.81</td><td>0.6445 ± 0.1094</td><td>31.71 ± 2.83</td></tr><tr><td>ES (Ours)</td><td><strong>28.54 ± 2.75</strong></td><td><strong>0.6195 ± 0.1188</strong></td><td><strong>31.53 ± 2.74</strong></td></tr><tr><td>EI</td><td>27.88 ± 2.64</td><td>0.5731 ± 0.1299</td><td>30.79 ± 2.64</td></tr><tr><td>SURE</td><td>24.45 ± 1.86</td><td>0.5479 ± 0.0740</td><td>27.35 ± 1.90</td></tr><tr><td>IDFT</td><td>23.62 ± 1.90</td><td>0.5052 ± 0.0900</td><td>25.99 ± 1.94</td></tr></tbody><thead><tr><th rowspan="2">Method</th><th colspan="3">Real MRI measurements (x8 Accel.)</th></tr><tr><th>PSNR ↑</th><th>SSIM ↑</th><th>EQUIV ↑</th></tr></thead><tbody><tr><td>Supervised</td><td>28.81 ± 2.85</td><td>0.6480 ± 0.1103</td><td>31.81 ± 2.84</td></tr><tr><td>ES (Ours)</td><td><strong>28.30 ± 2.64</strong></td><td><strong>0.6151 ± 0.1179</strong></td><td><strong>31.29 ± 2.62</strong></td></tr><tr><td>EI</td><td>27.88 ± 2.61</td><td>0.5740 ± 0.1290</td><td>30.80 ± 2.61</td></tr><tr><td>MC</td><td>23.63 ± 1.90</td><td>0.5061 ± 0.0904</td><td>26.00 ± 1.94</td></tr><tr><td>IDFT</td><td>23.63 ± 1.90</td><td>0.5060 ± 0.0904</td><td>26.00 ± 1.94</td></tr></tbody><thead><tr><th rowspan="2">Method</th><th colspan="3">CT (50 views, 50 dB SNR)</th></tr><tr><th>PSNR ↑</th><th>SSIM ↑</th><th>EQUIV ↑</th></tr></thead><tbody><tr><td>Supervised</td><td>33.99 ± 2.48</td><td>0.8819 ± 0.0585</td><td>34.00 ± 2.49</td></tr><tr><td>ES (Ours)</td><td><strong>32.62 ± 2.16</strong></td><td><strong>0.8570 ± 0.0596</strong></td><td><strong>32.60 ± 2.17</strong></td></tr><tr><td>EI</td><td>28.61 ± 1.28</td><td>0.7400 ± 0.0466</td><td>28.61 ± 1.29</td></tr><tr><td>IDFT</td><td>25.59 ± 0.69</td><td>0.4805 ± 0.0363</td><td>25.59 ± 0.70</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>Figure 1: **Compressive sensing results.** ES (ours) performs similarly as the supervised baseline, unlike EI (baseline) whose performance gap widens with higher compression levels.</figcaption></figure>

## 5.3 MAGNETIC RESONANCE IMAGING

The dataset contains 320 ×320 images from FastMRI (Zbontar et al., 2019) subsampled in the Fourier domain by a single binary mask corresponding to an acceleration of 8, as well as Gaussian noise with a standard deviation of 0.005 corresponding to a signal-to-noise ratio (SNR) of 40 dB. Mathematically, the forward operator $A \in \mathbb{R}^{m \times n}$ is expressed as $A = MF$ where $F \in \mathbb{R}^{n \times n}$ denotes the $n \times n$ discrete Fourier transform matrix and where $M \in \mathbb{R}^{m \times n}$ is the subsampling mask defined as $M_{i,j} = \delta_{j_i, j}$ for $i = 1, ..., m$ and $j = 1, ..., n$, where $j_i$ denotes the $i$-th component in $\mathbb{R}^n$ corresponding to a pixel in one of the subsampled vertical lines in a random Gaussian mask. Out
<table><caption>Table 2: **Inpainting results.** ES (ours) performs better than EI (baseline), both in terms of reconstruction quality (PSNR, SSIM) and measured equivariance (EQUIV), while performing competitively against the supervised baseline. In **bold**, the best self-supervised metrics (avg ± st.d.).</caption><thead><tr><th>Method</th><th>PSNR ↑</th><th>SSIM ↑</th><th>EQUIV ↑</th></tr></thead><tbody><tr><td>Supervised</td><td>28.46 ± 2.97</td><td>0.8982 ± 0.0411</td><td>28.46 ± 2.97</td></tr><tr><td>ES (Ours)</td><td><b>27.45 ± 2.86</b></td><td><b>0.8737 ± 0.0461</b></td><td><b>27.46 ± 2.85</b></td></tr><tr><td>EI</td><td>25.89 ± 2.65</td><td>0.8332 ± 0.0521</td><td>25.89 ± 2.65</td></tr><tr><td>MC</td><td>8.22 ± 2.47</td><td>0.0983 ± 0.0551</td><td>8.22 ± 2.47</td></tr><tr><td>Incomplete image</td><td>8.22 ± 2.47</td><td>0.0973 ± 0.0542</td><td>N/A</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2: **Sample reconstructions for image inpainting.** ES (ours) produces images perceptually closer to the supervised baseline than EI (baseline) which appears blurry.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: **Sample reconstructions for MRI (x8 Accel., 40 dB SNR).** Unlike EI (baseline) which suffers from dot-shaped artifacts, ES (ours) is perceptually closer to the supervised baseline. In line with the theoretical predictions, SURE and IDFT (baselines) fail to recover information beyond the observed frequencies.</figcaption></figure>
of the 973 images in the full dataset, 900 are selected for training and the remaining 73 are used for testing. We also test our method on a noise dominated setting using a different mask corresponding to an acceleration of 6, with a higher noise level of 0.1 corresponding to an SNR of only 10 dB, see these results in Appendix B.4. The variant of our proposed loss we use in this experiment is the R2R one introduced in Section 4 with $\alpha = 0.5$. Table 1 shows the performance of the different methods on the test set. Table 3 shows the synergy of our method with the equivariant architecture. Figure 3 shows sample reconstructions from the trained models.

We also evaluate our method on real MRI measurements from FastMRI instead of synthetic measurements. We normalize the k-space data that are sampled on different grids for different scans by resampling them using aliasing-free sinc interpolation on the 320 × 320 sampling grid of the ground truth scans. The resulting k-space data are subsampled using the same ×8 acceleration mask used in the setting with synthetic measurements. Table 1 shows that our method performs competitively.

## 5.4 COMPUTED TOMOGRAPHY

The dataset consists in pairs of ground truth CT scans and corresponding sinograms. In order to obtain the CT scans, we resize to 256 × 256 pixels and clip between -1,000 and 1,000 HUs the ground truth scans from the LIDC-IDRI dataset (Armato III et al., 2011). For each scan, we compute the corresponding sinogram using a discrete Radon transform with 50 views and additive white Gaussian noise with a standard deviation of 0.001 corresponding to a SNR of about 50 dB. The resulting 1,010 pairs are further split into 900 training pairs and 110 test pairs. Table 1 shows that ES performs almost as well as the supervised baseline, and better than EI.

## 5.5 ABLATION STUDY

Table 3 shows that architectures designed to be equivariant are measurably more equivariant across imaging modalities and training losses. Moreover, it shows that networks trained using the splitting loss perform better for equivariant architectures than for non-equivariant architectures, with an even greater gap than for supervised inpainting baselines, which confirms the theoretical analysis made in Section 4. Finally, we observe that non-equivariant architectures still lead to fairly measurably equivariant models. While surprising, this phenomenon has already been witnessed and is usually referred to as learned equivariance, whereby the training data and inductive biases lead to fairly equivariant learned models (Gruver et al., 2024). We believe that this learned equivariance is responsible for the high performance of splitting methods even when using non-equivariant architectures. For additional results on the impact of equivariant architectures, see Appendix B.4.

<table><caption>Table 3: **Impact of using equivariant architectures.** In accordance with the theoretical results described in Section 4, there is a synergy between the splitting loss and equivariant architectures resulting in higher performance. Non-equivariant models have surprisingly high equivariance measures (EQUIV) which might explain their high performance when using the splitting loss. Eq. arch. denotes whether the architecture is equivariant. In **bold**, the best self-supervised metrics (avg ± st.d.).</caption><thead><tr><th rowspan="2">Training loss</th><th rowspan="2">Eq. arch.</th><th colspan="3">Image inpainting</th></tr><tr><th>PSNR ↑</th><th>SSIM ↑</th><th>EQUIV ↑</th></tr></thead><tbody><tr><td rowspan="2">Supervised</td><td>✓</td><td>28.46 ± 2.97</td><td>0.8982 ± 0.0411</td><td>28.46 ± 2.97</td></tr><tr><td>×</td><td>28.62 ± 3.03</td><td>0.9002 ± 0.0414</td><td>27.85 ± 2.71</td></tr><tr><td rowspan="2">Splitting (Ours)</td><td>✓</td><td><b>27.45 ± 2.86</b></td><td><b>0.8737 ± 0.0461</b></td><td><b>27.46 ± 2.85</b></td></tr><tr><td>×</td><td>27.20 ± 2.83</td><td>0.8652 ± 0.0461</td><td>26.52 ± 2.60</td></tr></tbody></table><table><thead><tr><th rowspan="2">Training loss</th><th rowspan="2">Eq. arch.</th><th colspan="3">MRI (×8 Accel., 40 dB SNR)</th></tr><tr><th>PSNR ↑</th><th>SSIM ↑</th><th>EQUIV ↑</th></tr></thead><tbody><tr><td rowspan="2">Supervised</td><td>✓</td><td>28.74 ± 2.81</td><td>0.6445 ± 0.1094</td><td>31.71 ± 2.83</td></tr><tr><td>×</td><td>28.48 ± 2.68</td><td>0.6381 ± 0.1082</td><td>28.78 ± 1.95</td></tr><tr><td rowspan="2">Splitting (Ours)</td><td>✓</td><td><b>28.54 ± 2.75</b></td><td><b>0.6195 ± 0.1188</b></td><td><b>31.53 ± 2.74</b></td></tr><tr><td>×</td><td>28.18 ± 2.58</td><td>0.6104 ± 0.1176</td><td>27.28 ± 2.10</td></tr></tbody></table>
# 6 CONCLUSION

In this work, we propose a new self-supervised loss for solving inverse problems which bridges the gap between existing equivariance and splitting-based self-supervised losses. We motivate the design of our loss by showing that minimizing the expected loss results in MMSE estimators. We further validate our method using numerical simulations on different image distributions and imaging modalities, including inpainting of natural images, MRI and CT. These results suggest that the proposed method compares favorably to the equivariant imaging baseline and is close to supervised methods. To the best of our knowledge, this work is the first to leverage equivariant networks to learn from incomplete data alone, going beyond the usual goal of improving the generalization of the networks to unseen transformations at test time. Our method provides a new way to evaluate the benefits of using different equivariant architectures, and can benefit from future advances made in this field. More broadly, our work is further evidence that invariance is a promising prior for learning from incomplete data.

## REPRODUCIBILITY STATEMENT

We share the implementation of our method and experiments to make our work easier to reproduce.

## REFERENCES

- Hemant Kumar Aggarwal, Merry P. Mani, and Mathews Jacob. MoDL: Model Based Deep Learning Architecture for Inverse Problems. *IEEE Trans. Med. Imaging*, 38(2):394–405, February 2019.

- Eirikur Agustsson and Radu Timofte. NTIRE 2017 Challenge on Single Image Super-Resolution: Dataset and Study. In *CVPRW*, pp. 1122–1131, July 2017.

- Samuel G. Armato III, Geoffrey McLennan, Luc Bidaut, Michael F. McNitt-Gray, Charles R. Meyer, Anthony P. Reeves, Binsheng Zhao, Denise R. Aberle, Claudia I. Henschke, Eric A. Hoffman, Ella A. Kazerooni, Heber MacMahon, Edwin J. R. van Beek, David Yankelevitz, Alberto M. Biancardi, Peyton H. Bland, Matthew S. Brown, Roger M. Engelmann, Gary E. Laderach, Daniel Max, Richard C. Pais, David P.-Y. Qing, Rachael Y. Roberts, Amanda R. Smith, Adam Starkey, Poonam Batra, Philip Caligiuri, Ali Farooqi, Gregory W. Gladish, C. Matilda Jude, Reginald F. Munden, Iva Petkovska, Leslie E. Quint, Lawrence H. Schwartz, Baskaran Sundaram, Lori E. Dodd, Charles Fenimore, David Gur, Nicholas Petrick, John Freymann, Justin Kirby, Brian Hughes, Alessi Vande Casteele, Sangeeta Gupte, Maha Sallam, Michael D. Heath, Michael H. Kuhn, Ekta Dharaiya, Richard Burns, David S. Fryd, Marcos Salganicoff, Vikram Anand, Uri Shreter, Stephen Vastagh, Barbara Y. Croft, and Laurence P. Clarke. The Lung Image Database Consortium (LIDC) and Image Database Resource Initiative (IDRI): A Completed Reference Database of Lung Nodules on CT Scans. *Medical Physics*, 38(2):915–931, 2011. ISSN 2473-4209. doi: 10.1118/1.3528204.

- Chinmay Belthangady and Loic A Royer. Applications, promises, and pitfalls of deep learning for fluorescence image reconstruction. *Nature methods*, 16(12):1215–1225, 2019.

- Yochai Blau and Tomer Michaeli. The Perception-Distortion Tradeoff. In *CVPR*, pp. 6228–6237, 2018.

- Elena Celledoni, Matthias J. Ehrhardt, Christian Etmann, Brynjulf Owren, Carola-Bibiane Schönlieb, and Ferdia Sherry. Equivariant neural networks for inverse problems. *Inverse Problems*, 37(8):085006, July 2021.

- Anadi Chaman and Ivan Dokmanić. Truly shift-invariant convolutional neural networks. In *CVPR*, pp. 3773–3783, 2021a.

- Anadi Chaman and Ivan Dokmanić. Truly shift-equivariant convolutional neural networks with adaptive polyphase upsampling. In *Asilomar Conf Signals Syst Comput*, pp. 1113–1120, 2021b.

- Dongdong Chen, Julián Tachella, and Mike E. Davies. Equivariant Imaging: Learning Beyond the Range Space. In *ICCV*, pp. 4379–4388, 2021.
Dongdong Chen, Julián Tachella, and Mike E. Davies. Robust Equivariant Imaging: A fully unsupervised framework for learning to image from noisy and partial measurements. In CVPR, pp. 5647-5656, 2022.

Taco S. Cohen and Max Welling. Group Equivariant Convolutional Networks. In ICML, pp. 2990-2999, 2016.

Taco S. Cohen and Max Welling. Steerable CNNs. In ICLR, 2017.

Giannis Daras, Kulin Shah, Yuval Dagan, Aravind Gollakota, Alexandros G. Dimakis, and Adam Klivans. Ambient Diffusion: Learning Clean Distributions from Corrupted Data. In NeurIPS, pp. 288-313, 2023.

Leo Davy, Luis M. Briceno-Arias, and N. Pustelnik. Restarted contractive operators to learn at equilibrium, June 2025.

Yonina C. Eldar. Generalized SURE for Exponential Families: Applications to Regularization. IEEE Trans Signal Process, 57(2):471-481, 2009.

Fabian Ewald Fassnacht, Joanne C White, Michael A Wulder, and Erik Næsset. Remote sensing in forestry: current challenges, considerations and directions. Forestry: An International Journal of Forest Research, 97(1):11-37, 2024.

Nate Gruver, Marc Finzi, Micah Goldblum, and Andrew Gordon Wilson. The Lie Derivative for Measuring Learned Equivariance, 2024.

Allard A. Hendriksen, Daniel M. Pelt, and K. Joost Batenburg. Noise2Inverse: Self-supervised deep convolutional denoising for tomography. IEEE Trans. Comput. Imaging, 6:1320-1335, 2020.

Kyong Hwan Jin, Michael T McCann, Emmanuel Froustey, and Michael Unser. Deep convolutional neural network for inverse problems in imaging. IEEE Trans Image Process, 26(9):4509-4522, 2017.

Sékou-Oumar Kaba, Arnab Kumar Mondal, Yan Zhang, Yoshua Bengio, and Siamak Ravanbakhsh. Equivariance with Learned Canonicalization Functions. In ICML, pp. 15546-15566, 2023.

Tero Karras, Miika Aittala, Samuli Laine, Erik Härkönen, Janne Hellsten, Jaakko Lehtinen, and Timo Aila. Alias-free generative adversarial networks. In NeurIPS, volume 34, pp. 852-863, 2021.

Achim Klenke. Probability theory: a comprehensive course. Springer, 2008.

Alexander Krull, Tim-Oliver Buchholz, and Florian Jug. Noise2Void - Learning Denoising From Single Noisy Images. In CVPR, pp. 2129-2137, 2019.

Ilya Loshchilov and Frank Hutter. Decoupled Weight Decay Regularization, 2019.

Christopher A. Metzler, Ali Mousavi, Reinhard Heckel, and Richard G. Baraniuk. Unsupervised Learning with Stein's Unbiased Risk Estimator, July 2020.

Hagay Michaeli, Tomer Michaeli, and Daniel Soudry. Alias-Free Convnets: Fractional Shift Invariance via Polynomial Activations. In CVPR, pp. 16333-16342, 2023.

Charles Millard and Mark Chiew. A theoretical framework for self-supervised MR image reconstruction using sub-sampling via variable density Noisier2Noise. IEEE Trans Comput Imaging, 9:707-720, 2023.

Brayan Monroy, Jorge Bacca, and Julián Tachella. Generalized Recorrupted-to-Recorrupted: Self-Supervised Learning Beyond Gaussian Noise. In CVPR, pp. 28155-28164, 2025.

Greg Ongie, Rebecca Willett, Daniel Soudry, and Nathan Srebro. A Function Space View of Bounded Norm Infinite Width ReLU Nets: The Multivariate Case, October 2019.

Tongyao Pang, Huan Zheng, Yuhui Quan, and Hui Ji. Recorrupted-to-Recorrupted: Unsupervised Deep Learning for Image Denoising. In CVPR, pp. 2043-2052, June 2021.
Javier Portilla and Eero P Simoncelli. A parametric texture model based on joint statistics of com-
plex wavelet coefficients. *International journal of computer vision*, 40(1):49–70, 2000.

Omri Puny, Matan Atzmon, Heli Ben-Hamu, Ishan Misra, Aditya Grover, Edward J. Smith, and Yaron Lipman. Frame Averaging for Invariant and Equivariant Network Design, March 2022.

Marco Ragone, Reza Shahabazian-Yassar, Farzad Mashayek, and Vitaliy Yurkiv. Deep learning modeling in microscopy imaging: A review of materials science applications. *Progress in Materials Science*, 138:101165, 2023.

Sebastián Rivera, Iván Ortíz, Tatiana Gelvez, Fernando Rojas, and Henry Arguello. Rotation invariant deep learning approach for image inpainting. In *STSIVA*, pp. 1–5, 2021.

Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-Net: Convolutional Networks for Biomedical Image Segmentation. In *MICCAI*, pp. 234–241, 2015.

Leonid I. Rudin, Stanley Osher, and Emad Fatemi. Nonlinear total variation based noise removal algorithms. *Physica D: Nonlinear Phenomena*, 60(1):259–268, November 1992.

Akiyoshi Sannai, Makoto Kawano, and Wataru Kumagai. Equivariant and Invariant Reynolds Networks, October 2021.

Akiyoshi Sannai, Makoto Kawano, and Wataru Kumagai. Invariant and equivariant reynolds networks. *Journal of Machine Learning Research*, 25(42):1–36, 2024.

Jérémy Scanvic, Mike Davies, Patrice Abry, and Julián Tachella. Scale-Equivariant Imaging: Self-Supervised Learning for Image Super-Resolution and Deblurring, March 2025.

Victor Sechaud, Laurent Jacques, Patrice Abry, and Julián Tachella. Equivariance-based self-supervised learning for audio signal recovery from clipped measurements. In *EUSIPCO*, August 2024.

Julián Tachella, Dongdong Chen, and Mike Davies. Unsupervised learning from incomplete measurements for inverse problems. *NeurIPS*, 35:4983–4995, 2022.

Julián Tachella, Dongdong Chen, and Mike Davies. Sensing Theorems for Unsupervised Learning in Linear Inverse Problems. *Journal of Machine Learning Research*, 24(39):1–45, 2023. ISSN 1533-7928.

Julián Tachella, Mike Davies, and Laurent Jacques. UNSURE: Self-supervised learning with Un-known Noise level and Stein's Unbiased Risk Estimate. In *ICLR*, 2025a.

Julián Tachella, Matthieu Terris, Samuel Hurault, Andrew Wang, Dongdong Chen, Minh-Hai Nguyen, Maxime Song, Thomas Davies, Leo Davy, Jonathan Dong, Paul Escande, Johannes Hertrich, Zhiyuan Hu, Tobías I. Liaudat, Nils Laurent, Brett Levac, Mathurin Massias, Thomas Moreau, Thibaut Modrzyk, Brayan Monroy, Sebastian Neumayer, Jérémy Scanvic, Florian Sar-ron, Victor Sechaud, Georg Schramm, Chao Tang, Romain Vo, and Pierre Weiss. DeepInverse: A Python package for solving imaging inverse problems with deep learning, May 2025b.

Matthieu Terris, Thomas Moreau, Nelly Pustelnik, and Julian Tachella. Equivariant plug-and-play image reconstruction. In *CVPR*, pp. 25255–25264, 2024.

Singanallur V Venkatakrishnan, Charles A Bouman, and Brendt Wohlberg. Plug-and-play priors for model based reconstruction. In *IEEE Glob Conf Signal Inf Process*, pp. 945–948, 2013.

Antonia Vojtekova, Maggie Lieu, Ivan Valtchanov, Bruno Altieri, Lyndsay Old, Qifeng Chen, and Filip Hroch. Learning to denoise astronomical images with u-nets. *Monthly Notices of the Royal Astronomical Society*, 503(3):3204–3215, 2021.

Andrew Wang and Mike Davies. Perspective-equivariance for unsupervised imaging with camera geometry. In Alessio Del Bue, Cristian Canton, Jordi Pont-Tuset, and Tatiana Tommasi (eds.), *ECCV Workshops*, pp. 116–134, 2024.

Andrew Wang, Steven McDonagh, and Mike Davies. Benchmarking Self-Supervised Learning Methods for Accelerated MRI Reconstruction, May 2025.
Zhou Wang, A.C. Bovik, H.R. Sheikh, and E.P. Simoncelli. Image quality assessment: From error visibility to structural similarity. *IEEE Trans Image Process*, 13(4):600–612, April 2004.

Philip J Withers, Charles Bouman, Simone Carmignato, Veerle Cnudde, David Grimaldi, Charlotte K Hagen, Eric Maire, Marena Manley, Anton Du Plessis, and Stuart R Stock. X-ray computed tomography. *Nature Reviews Methods Primers*, 1(1):18, 2021.

Guixian Xu, Jinglai Li, and Junqi Tang. Fast Equivariant Imaging: Acceleration for Unsupervised Learning via Augmented Lagrangian and Auxiliary PnP Denoisers, July 2025.

Xuemiao Xu, Huaidong Zhang, Guoqiang Han, Kin Chung Kwan, Wai-Man Pang, Jiaming Fang, and Gansen Zhao. A two-phase space resection model for accurate topographic reconstruction from lunar imagery with pushbroomscanners. *Sensors*, 16(4):507, 2016.

Burhaneddin Yaman, Seyed Amir Hossein Hosseini, Steen Moeller, Jutta Ellermann, Kâmil Uğurbil, and Mehmet Akçakaya. Self-Supervised Learning of Physics-Guided Reconstruction Neural Networks without Fully-Sampled Reference Data. *Magnetic Resonance in Med*, 84(6):3172–3191, 2020.

Jure Zbontar, Florian Knoll, Anuroop Sriram, Tullie Murrell, Zhengnan Huang, Matthew J. Muckley, Aaron Defazio, Ruben Stern, Patricia Johnson, Mary Bruno, Marc Parente, Krzysztof J. Geras, Joe Katsnelson, Hersh Chandarana, Zizhao Zhang, Michal Drozdzal, Adriana Romero, Michael Rabbat, Pascal Vincent, Nafissa Yakubova, James Pinkerton, Duo Wang, Erich Owens, C. Lawrence Zitnick, Michael P. Recht, Daniel K. Sodickson, and Yvonne W. Lui. fastMRI: An Open Dataset and Benchmarks for Accelerated MRI, December 2019.

Richard Zhang. Making Convolutional Networks Shift-Invariant Again. In *ICML*, pp. 7324–7334, 2019.

# A DETAILS ABOUT THE METHOD

## A.1 EQUIVARIANT RECONSTRUCTOR

<figure><img src="image_5.png" /><figcaption>Figure 4: **Compressive sensing results.** Adds the pseudo-inverse reconstruction (pseudo-inv.) as a baseline to Figure 1</figcaption></figure>

In Theorem 2, we consider reconstruction functions of 4 different structures. In this section, we give additional details about them.

The artifact removal reconstructor architecture (Jin et al., 2017) in eq. (16) consists in a projection step that maps the measurements back into the image space using the adjoint or the pseudo-inverse of the forward matrix, which is immediately followed by a very general trainable network, often a UNet or another encoder-decoder type of network, or sometimes simply a fully convolutional network. Theorem 2 states that this reconstructor architecture produces equivariant reconstructors as long as the denoiser architecture is, itself, equivariant in the sense of eq. (14).

Reynolds averaging for possibly non-equivariant reconstruction functions $r(y, A)$ in eq. (18)

$$f(y, A) = \frac{1}{|G|} \sum_{g \in G} T_g r(y, AT_g). \qquad (18)$$
$$ \text{Equivariant reconstructor} \\ y \longrightarrow \frac{1}{|\mathcal{G}|} \sum_{g \in \mathcal{G}} T_g^r(y, AT_g) \longrightarrow f(y, A) $$

<figure><img src="image_6.png" /><figcaption>Figure 5: **Reconstructer equivariant to rotations and flips.** It is the Reynolds averaging of 90° rotations and horizontal and vertical flips as in eq. (18) of a non-equivariant reconstructor of the MAP type, as in eq. (19), implemented using a MoDL unrolled algorithm (Aggarwal et al., 2019) with 3 iterations, shared weights, and using a non-equivariant residual UNet (Ronneberger et al., 2015) as the denoiser architecture.</figcaption></figure>

has, as far as we know, not been defined in previous works. It is a natural extension of Reynolds averaging for possibly non-equivariant image-to-image functions or denoisers $\psi(x)$ (Sannai et al., 2024; Terris et al., 2024),

$$ \phi(x) = \frac{1}{|\mathcal{G}|} \sum_{g \in \mathcal{G}} T_g^{-1} \psi(T_g x) \qquad (22) $$

which makes them equivariant in the sense of eq. (14), to reconstruction functions in order to make them equivariant in the sense of eq. (15). It is a fairly simple way to make a reconstructor equivariant and it is relatively inexpensive for small groups of transformations such as the group of 90° rotations and horizontal and vertical flips (Cohen & Welling, 2016). It is however too expensive to be used in practice when the group is relatively large, like the group of shifts, since it would require to evaluate the neural network for as many times as there are pixels in the input image, for each image.

The maximum a posteriori (MAP) reconstructor architecture defined in eq. (19)

$$ f(y, A) = \underset{x \in \mathbb{R}^n}{\operatorname{argmax}} \{p(x | y, A)\}. \qquad (19) $$

is a classical reconstructor architecture that has been used in different ways, including iterative algorithms with a hand-crafted prior (Rudin et al., 1992; Davy et al., 2025), plug-and-play architectures using an iterative approach but with the hand-crafted prior replaced with a pre-trained denoiser (Venkatakrishnan et al., 2013), and unrolled architectures where the optimization problem is done with a fixed number of steps and where parts of the algorithm are replaced with trainable modules (Aggarwal et al., 2019). It is often interpreted as a Bayesian maximum a posteriori estimator, but also commonly outside of a Bayesian framework as a variational approach with a data-fidelity term and a regularization term. It is one of the reconstructor architectures that we use for our experiments in Section 5. Theorem 2 states that these reconstructors are equivariant as long as the prior distribution is invariant in the sense of eq. (5), or equivalently, if the associated regularization function or negative log-prior is invariant.

The minimum mean squared error (MMSE) estimator in eq. (20)

$$ f(y, A) = \mathbb{E}_{x \in y, A} \{x\} \qquad (20) $$

is generally the target theoretical reconstructor as it achieves the highest theoretical PSNR (Tachella et al., 2023). It is generally estimated by reconstructors trained with a mean squared error
loss (Chen et al., 2021), which is notably how we train the supervised baselines in the experiments in Section 5. Theorem 2 states that as long as the prior distribution is invariant in the sense of eq. (5), the MMSE reconstructor is equivariant in the sense of eq. (15).

<table><caption>Table 4: **MRI results in another setting.** Supplementary results for a different MRI problem (×6 Accel., 10 dB SNR) than in Table 1 In **bold**, the best self-supervised metrics. Values: avg ± st.d.</caption><thead><tr><th rowspan="2">Method</th><th colspan="3">MRI (×6 Accel., 10 dB SNR)</th></tr><tr><th>PSNR ↑</th><th>SSIM ↑</th><th>EQUIV ↑</th></tr></thead><tbody><tr><td>Supervised</td><td>27.39 ± 2.44</td><td>0.5243 ± 0.1373</td><td>30.38 ± 2.43</td></tr><tr><td>ES (Ours)</td><td><b>27.33 ± 2.45</b></td><td><b>0.5126 ± 0.1444</b></td><td><b>30.32 ± 2.44</b></td></tr><tr><td>EI</td><td>27.23 ± 2.41</td><td>0.5110 ± 0.1421</td><td>30.21 ± 2.40</td></tr><tr><td>SURE</td><td>27.08 ± 2.29</td><td>0.5097 ± 0.1372</td><td>30.06 ± 2.28</td></tr><tr><td>IDFT</td><td>23.85 ± 1.05</td><td>0.3878 ± 0.0272</td><td>25.14 ± 0.79</td></tr></tbody></table>

## A.2 HOW TO CHOOSE TRANSFORMATIONS FOR A GIVEN INVERSE PROBLEM

<table><caption>Table 5: **Decision table for the transformations.** Corollary 1 shows that not all transformations are well-suited for all problems. Namely, transformations for which the operator is equivariant introduce no additional information and should not be used. This table specifies which transformations are well-suited for which operator.</caption><thead><tr><th>Operator</th><th>Translation</th><th>Rotation</th><th>Permutation</th><th>Amplitude</th></tr></thead><tbody><tr><td>Isotropic blur</td><td>×</td><td>×</td><td>✓</td><td>×</td></tr><tr><td>Image inpainting</td><td>✓</td><td>✓</td><td>✓</td><td>×</td></tr><tr><td>Sparse-view CT</td><td>×</td><td>✓</td><td>✓</td><td>×</td></tr><tr><td>Accelerated MRI</td><td>×</td><td>✓</td><td>✓</td><td>×</td></tr><tr><td>Compressive sensing</td><td>✓</td><td>✓</td><td>✓</td><td>×</td></tr></tbody></table>

There are two major criteria for choosing the transformations for a given application. First, the image distribution of interest should be invariant to the chosen transformations. Aerial, remote sensing and microscopic images are invariant to translations and rotations as the scenes and subjects they measure exhibit no privileged position and orientation with respect to the image plane. Natural image distributions and texture distributions (Portilla & Simoncelli, 2000) are also generally invariant to translations but they are less invariant to rotations as natural images are typically oriented upward and texture distributions might be anisotropic.

Second, the transformations should also be chosen in accordance with the measurement operator. Corollary 1 shows that transformations for which the operator is equivariant do not improve the reconstruction process. It is a well-known criterion introduced in the original work on equivariant imaging (Chen et al., 2021) which remains correct in our setting where measurement splitting is added to the theoretical analysis. Table 5 lists correct choices of transformations for common measurement operators.

## A.3 END-TO-END ALGORITHM

In this section, we present the end-to-end ES algorithm. It consists in a training step where an equivariant reconstructor is trained using backpropagation against a training dataset of measurements only, after which it is applied to obtain the reconstructions associated to the test measurements. The ES loss used in the training step is computed using the expression in Equation (3). The expectations are estimated using Monte Carlo sampling where a single sample is used at training time and $T = 10$ samples are used at inference. In the experiments, we use the splitting ratio $s = 0.8$ corresponding to $m_1 = 0.8m$ and $m_2 = 0.2m$. Algorithms 1 and 2 show the detailed algorithms in pseudo-code.
**Algorithm 1:** Equivariant Splitting (Training procedure)

**Input:** Dataset $D = \{y_i\}_{i \in I}$, forward operator $A \in \mathbb{R}^{m \times n}$, split ratio $s$, learning rate $\eta$, number of epochs $E$, equivariant reconstructor $f_\theta$

**Output:** Trained model $f_\theta$

for epoch = 1 to E do

Shuffle dataset $D$;

foreach mini-batch $B \subset D$ do

init: $L = 0$

foreach $y$ in $B$ do

init: $M \in \mathbb{R}^{m_1 \times m}$ a random split s.t $m_1 = s \times m$.

Split measurement and operator: $y_1 = My, A_1 = MA$;

Compute predictions: $\hat{x} = f_\theta(y_1, A_1)$;

Compute sample loss $\ell = \|A\hat{x} - y\|_2^2$;

Aggregate loss: $L \leftarrow L + \ell$;

end

Mean: $L \leftarrow \frac{1}{|B|}L$;

Backpropagation: compute gradients $\nabla_\theta L$;

Update parameters: $\theta \leftarrow \text{ADAM}(\theta, \nabla_\theta L, \eta)$;

end

end

**Algorithm 2:** Equivariant Splitting (Inference)

**Input:** Measurement $y \in \mathbb{R}^{m \times n}$, forward operator $A \in \mathbb{R}^{m \times n}$, split ratio $s$, number of samples $T$, trained reconstructor $f_\theta$

**Output:** Reconstructed images $\hat{x} = \bar{f}_\theta(y, A)$

init: $\hat{x} = 0$

for $i = 1$ to $T$ do

init: $M \in \mathbb{R}^{m_1 \times m}$ a random split s.t $m_1 = s \times m$.

Split measurement and operator: $y_1 = My, A_1 = MA$;

Compute one predictions: $\hat{x}_i = f_\theta(y_1, A_1)$;

Aggregate: $\hat{x} \leftarrow \hat{x} + \hat{x}_i$;

end

Mean: $\hat{x} \leftarrow \frac{1}{T}\hat{x}$;

## B DETAILS ABOUT THE EXPERIMENTS

We use the optimizer AdamW (Loshchilov & Hutter, 2019) for every training with different learning rates for the different inverse problems, a weight decay of $10^{-8}$, beta coefficients equal to 0.9 and 0.999 and without the AMSGrad option. For longer trainings, we use step schedulers that divide the learning rate by a factor ranging from 2 to 10 at specific epochs, up to 3 or 4 times.

In our experiments, we make extensive use of the DeepInverse library (Tachella et al., 2025b) that provides an implementation of the various forward operators and training losses that we use. Every model is trained for up to 50 hours on a single GPU, either an NVIDIA H100, GH200 or RTX 4090. See Appendix B.4 for more details about the training durations.

### B.1 DATASETS

**Image inpainting** The ground truth images are images obtained from the dataset DIV2K (Agustsson & Timofte, 2017) which contains pictures of natural scenes (landscapes, animals) by first resizing them to a resolution of 256 × 256 pixels before extracting a central 128 × 128 pixels crop. We synthesize the measurements by corrupting the ground truth images using a single binary mask sampled from a pixel-wise Bernoulli distribution with a 30% chance of keeping each pixel value. In this setting, we do not corrupt the measurements further with additional noise.
<figure><img src="image_7.png" /><figcaption>Figure 6: Examples of k-space measurements in MRI. (left) a less noisy problem with less measurements, (right) a noisier problem with more measurements.</figcaption></figure>

**MRI** The dataset consists in 973 pairs of ground truth images from the FastMRI dataset (Zbontar et al., 2019) and associated k-space measurements synthesized using a single coil sensitivity map. The ground truth images have a size of 320 × 320 pixels and each corresponds to the middle slice of a different 3d knee acquisition. The k-spaces are synthesized as discrete Fourier transforms subsampled on a single non-regular grid and further corrupted with additive white Gaussian noise with a standard deviation of 0.005 corresponding to an average SNR of about 40 dB. The subsampling grid models the coil sensitivity map, is sampled from a Gaussian distribution and corresponds to an acceleration of 8. The entire dataset is finally split into a train/val split containing 900 images and a test split containing the remaining 73 images. Our implementation uses the work from Wang et al. (2025).

We also consider an additional dataset obtained in the same way except using a different mask corresponding to an acceleration of 6 and with a higher noise level corresponding to a standard deviation of 0.1 or an average SNR of 10 dB in the k-space domain.

## B.2 NETWORK ARCHITECTURES

The unrolled architecture uses 3 iterations and the weights are shared across different iterations. Every UNet is residual, has 4 scales and has no normalization layer as we find them to be detrimental to the performance.

The transforms we use in the model architecture and in the metrics are grid-preserving: grid-aligned shifts, 90° rotations, and vertical and horizontal flips. The transforms we use in the EI are grid-aligned shifts and 1° rotations following the prior art Chen et al. (2021). Reynolds' averaging is implemented using an unbiased Monte-Carlo estimator whereby a single random transform is sampled at every evaluation to save on computational cost.

## B.3 METRICS

In the experiments, we use three different performance metrics including two standard distortion metrics (PSNR, SSIM) and a new equivariance metric for reconstructors (EQUIV). The peak signal-to-noise ratio is defined, for images with a dynamic range normalized to [0, 1], as the mean squared error expressed in decibels (dB)

$$ \mathrm{PSNR} = \mathbb{E}_{x,y} \left\{ -10 \log_{10} \left( \|f(\mathbf{y}, \mathbf{A}) - \mathbf{x}\|^2 \right) \right\}. \quad (23) $$

The structural similarity index measure (Wang et al., 2004) is a more perceptual metric which is a combination of the empirical means, standard deviations and correlation coefficient associated to the reference and compared images. The complete definition of the metric is too long to be included in this work and we refer the reader to the original publication for more details. In addition to these two standard distortion metrics, we use a new equivariance metric similar to that used by Chaman & Dokmanić (2021b) adapted for our proposed definition of equivariant reconstructors. It is the average mean squared error associated with eq. (15) and expressed in dB for readability

$$ \mathrm{EQUIV} = -10 \log_{10} \left( \mathbb{E}_{y,g} \left\{ \| f(\mathbf{y}, \mathbf{A} \mathbf{T}_g) - \mathbf{T}_g^{-1} f(\mathbf{y}, \mathbf{A}) \|^2 \right\} \right). \quad (24) $$
It can be roughly understood as a PSNR for equivariance. In every experiment, we use the same group of transformations for EQUIV as we use for the equivariant reconstructor architectures.

## B.4 RESULTS

In this section, we provide more details about the main experiments and present additional experiments.

For each training done in the main experiments, we report the average epoch duration and the number of epochs until the model has finished training. For the sake of comparability, we made sure to conduct the different trainings for the same imaging modality on the same GPU. Namely, we used a single NVIDIA RTX 3090 Ti GPU for every inpainting experiment and a single NVIDIA H100 GPU for every MRI experiment. Table 10 shows that self-supervised methods including ours have generally longer epochs and require more epochs than the fully supervised gold standard. We believe that this is due to a fundamental trade-off whereby the use of ground truth data accelerates the learning procedure while also enabling the use of simpler training algorithms. Moreover, unlike standard implementations of EI and SURE, ES requires only a single network pass per iteration resulting in significantly faster epochs. Overall, ES is computationally more efficient than EI, in addition to being more performant in terms of reconstruction quality.

Figure 4 shows additional compressive sensing results, Table 4 and Figure 7 show results on a MRI experiment with settings different from the main one, Table 6 shows extended results for the ablation study on equivariant architectures, and Figure 6 shows sample k-spaces from the MRI experiments. Figure 5 shows the equivariant reconstructor architecture that we adopt in the MRI experiments.

<table><caption>Table 6: Extended results on the impact of equivariant architectures. Adds to Table 3 the results for EI with a non-equivariant architecture for the inpainting task, results for the noise-dominated MRI task. In **bold**, the best self-supervised metrics. Values: avg ± st.d.</caption><thead><tr><th rowspan="2">Training loss</th><th rowspan="2">Eq. arch.</th><th colspan="3">Image inpainting</th></tr><tr><th>PSNR ↑</th><th>SSIM ↑</th><th>EQUIV ↑</th></tr></thead><tbody><tr><td rowspan="2">Supervised</td><td>✓</td><td>28.46 ± 2.97</td><td>0.8982 ± 0.0411</td><td>28.46 ± 2.97</td></tr><tr><td>×</td><td>28.62 ± 3.03</td><td>0.9001 ± 0.0415</td><td>27.85 ± 2.71</td></tr><tr><td rowspan="2">Splitting (Ours)</td><td>✓</td><td><b>27.45 ± 2.86</b></td><td>0.8737 ± 0.0461</td><td>27.46 ± 2.85</td></tr><tr><td>×</td><td>27.20 ± 2.83</td><td>0.8651 ± 0.0463</td><td>26.52 ± 2.60</td></tr><tr><td rowspan="2">EI loss</td><td>✓</td><td>25.89 ± 2.65</td><td>0.8332 ± 0.0521</td><td>25.89 ± 2.65</td></tr><tr><td>×</td><td>26.33 ± 2.81</td><td>0.8451 ± 0.0536</td><td>25.58 ± 2.52</td></tr><tr><td rowspan="2">MC loss</td><td>✓</td><td>8.22 ± 2.47</td><td>0.098 ± 0.055</td><td>8.22 ± 2.47</td></tr><tr><td>×</td><td>8.24 ± 2.48</td><td>0.100 ± 0.056</td><td>8.24 ± 2.48</td></tr></tbody></table><table><thead><tr><th rowspan="2">Training loss</th><th rowspan="2">Eq. arch.</th><th colspan="3">MRI (×8 Accel., 40 dB SNR)</th></tr><tr><th>PSNR ↑</th><th>SSIM ↑</th><th>EQUIV ↑</th></tr></thead><tbody><tr><td rowspan="2">Supervised</td><td>✓</td><td>28.74 ± 2.81</td><td>0.6445 ± 0.1094</td><td>31.71 ± 2.83</td></tr><tr><td>×</td><td>28.48 ± 2.68</td><td>0.6381 ± 0.1082</td><td>28.78 ± 1.95</td></tr><tr><td rowspan="2">Splitting (Ours)</td><td>✓</td><td><b>28.54 ± 2.75</b></td><td><b>0.6195 ± 0.1188</b></td><td><b>31.53 ± 2.74</b></td></tr><tr><td>×</td><td>28.18 ± 2.58</td><td>0.6104 ± 0.1176</td><td>27.28 ± 2.10</td></tr></tbody></table><table><thead><tr><th rowspan="2">Training loss</th><th rowspan="2">Eq. arch.</th><th colspan="3">MRI (×6 Accel., 10 dB SNR)</th></tr><tr><th>PSNR ↑</th><th>SSIM ↑</th><th>EQUIV ↑</th></tr></thead><tbody><tr><td rowspan="2">Supervised</td><td>✓</td><td>27.39 ± 2.44</td><td>0.5243 ± 0.1373</td><td>30.38 ± 2.43</td></tr><tr><td>×</td><td>27.33 ± 2.42</td><td>0.5174 ± 0.1410</td><td>29.73 ± 2.20</td></tr><tr><td rowspan="2">Splitting (Ours)</td><td>✓</td><td><b>27.33 ± 2.45</b></td><td><b>0.5126 ± 0.1444</b></td><td><b>30.32 ± 2.44</b></td></tr><tr><td>×</td><td>27.20 ± 2.38</td><td>0.5095 ± 0.1430</td><td>28.66 ± 1.76</td></tr></tbody></table>

We present an additional inpainting experiment, designed to simulate a more realistic acquisition scenario and study the consistency of our method. Specifically, we consider a noisy inpainting set-
<figure><img src="image_8.png" /><figcaption>Figure 7: Sample reconstructions for noise dominated MRI (×6 Accel., 10 dB SNR). In the noise dominated setting, the different models perform more similarly than in the less noisy setting shown in Figure 3.</figcaption></figure>

<table><caption>Table 7: **Inpainting with push-broom masks.** We consider an unevenly-spaced mask sampled randomly to test our method in a realistic inpainting setting, and an evenly-spaced mask to verify the claim in Corollary 1 empirically in a case of almost-equivariance.</caption><thead><tr><th rowspan="2">Method</th><th colspan="2">Image inpainting (Randomly-spaced push-broom mask)</th></tr><tr><th>PSNR ↑</th><th>SSIM ↑</th></tr></thead><tbody><tr><td>Supervised</td><td>23.72 ± 2.10</td><td>0.743 ± 0.051</td></tr><tr><td>ES (Ours)</td><td>23.04 ± 1.81</td><td>0.734 ± 0.050</td></tr><tr><td>EI</td><td>23.05 ± 2.13</td><td>0.707 ± 0.066</td></tr><tr><td>Incomplete image</td><td>9.44 ± 2.31</td><td>0.141 ± 0.048</td></tr></tbody></table><table><thead><tr><th rowspan="2">Method</th><th colspan="2">Image inpainting (Evenly-spaced push-broom mask)</th></tr><tr><th>PSNR ↑</th><th>SSIM ↑</th></tr></thead><tbody><tr><td>Supervised</td><td>28.37 ± 2.15</td><td>0.873 ± 0.035</td></tr><tr><td>ES (Ours)</td><td>21.94 ± 2.12</td><td>0.617 ± 0.091</td></tr><tr><td>Incomplete image</td><td>9.61 ± 2.42</td><td>0.152 ± 0.061</td></tr></tbody></table>
ting in which entire image columns are randomly removed. Such sampling patterns naturally arise in satellite imaging systems based on push-broom scanners (PBS) (Xu et al., 2016). As reported in Table 7, the proposed method achieves performance comparable to that of the EI baseline, confirming its consistency under this more practical acquisition model.

To clarify the behavior of ES when the forward operator $\mathbf{A}$ is equivariant or nearly equivariant, we conducted the following experiment: we performed an PBS inpainting operation in which every other column of the image is removed. In this setting, the operator $\mathbf{A}$ becomes “almost” equivariant in the sense that for any even horizontal shift $g$, we have $AT_g = T_g\mathbf{A}$, and the same property holds for all vertical shifts. Table 7 shows the reconstruction performance deteriorates significantly.

<figure><img src="image_9.png" /><figcaption>Figure 8: Performance evolution during training for inpainting. Splitting methods perform better than EI independent of the network architecture.</figcaption></figure>

### B.4.1 TESTING THE EFFECT OF UNKNOWN NOISE DISTRIBUTIONS

In addition to the main experiments where we use the knowledge of the noise distribution in the reconstruction algorithm, we test a more realistic scenario where the noise distribution is unknown and is estimated to be zero for a lack of a better estimate. To do so, we train models using variants of ES and EI where the measurements are assumed to be noiseless even though the training and testing data are nonetheless corrupted with noise. Table 8 shows that ES also performs better than EI when the noise distribution is unknown. Moreover, it shows that the performance of ES tends to be lower when the noise distribution is unknown, but it does not drop exceedingly which demonstrates the stability of ES to unexpected noise.

### B.4.2 EMPIRICAL VERIFICATION OF THE EQUIVARIANCE OF MAP RECONSTRUCTORS

We verify empirically that MAP reconstructors are equivariant as long as the prior is itself equivariant, i.e., the claim made in Theorem 2. Since they cannot be computed exactly in general, we consider a specific scenario where they can. We assume that 1) $x \sim \mathcal{N}(0, \tau^2\mathbf{I}_n)$, 2) $A \in \mathbb{R}^{m \times n}$ is the two-dimensional decimation operator with decimation rate 2, 3) $y | Ax \sim \mathcal{N}(0, \sigma^2\mathbf{I}_m)$, and 4) that $T_g$ denotes the rotation by angle $g \in \{0^\circ, 90^\circ, 180^\circ, 270^\circ\}$. Under these assumptions, the prior distribution is equivariant and the MAP estimator in eq. (19) can be expressed in closed-form as

$$f(\mathbf{y}, \mathbf{A}T_g) = \frac{\tau^2}{\tau^2 + \sigma^2} T_g^{-1} \mathbf{A}^\top \mathbf{y}, \qquad (25)$$

with $f(\mathbf{y}, \mathbf{A})$ being the special case where $T_g = \mathbf{I}_n$. In the experiment, we set $n = 128 \times 128$ for a grayscale image with 128 rows and 128 columns and we compute the equivariance metric in eq. (24) (EQUIV) for the MAP reconstructor using 256 i.i.d. samples from the joint distribution. Table 9 shows the results for every angle and for the average over all angles. As predicted theoretically, perfect equivariance is achieved.
<table><caption>Table 8: Performance when the noise distribution is unknown. In the unknown noise scenarios, we use the variants of ES and EI corresponding to assuming that the measurements are noiseless.</caption><thead><tr><th rowspan="2">Method</th><th rowspan="2">Known</th><th colspan="2">Image inpainting</th></tr><tr><th>PSNR ↑</th><th>SSIM ↑</th></tr></thead><tbody><tr><td>Supervised</td><td></td><td>23.72 ± 2.10</td><td>0.743 ± 0.051</td></tr><tr><td>ES (Ours)</td><td>✓</td><td>23.04 ± 1.81</td><td>0.734 ± 0.050</td></tr><tr><td>ES</td><td>×</td><td>22.05 ± 1.36</td><td>0.59 ± 0.061</td></tr><tr><td>EI</td><td>✓</td><td>23.05 ± 2.13</td><td>0.707 ± 0.066</td></tr><tr><td>EI</td><td>×</td><td>22.01 ± 1.71</td><td>0.59 ± 0.051</td></tr><tr><td>Incomplete image</td><td></td><td>9.44 ± 2.31</td><td>0.141 ± 0.048</td></tr><tr><th rowspan="2">Method</th><th rowspan="2">Known</th><th colspan="2">MRI (×8 Accel., 40 dB SNR)</th></tr><tr><th>PSNR ↑</th><th>SSIM ↑</th></tr><tr><td>Supervised</td><td></td><td>28.74 ± 2.81</td><td>0.6445 ± 0.1094</td></tr><tr><td>ES (Ours)</td><td>✓</td><td>28.54 ± 2.75</td><td>0.6195 ± 0.1188</td></tr><tr><td>ES</td><td>×</td><td>28.52 ± 2.75</td><td>0.6194 ± 0.1191</td></tr><tr><td>EI</td><td>✓</td><td>27.88 ± 2.64</td><td>0.5731 ± 0.1299</td></tr><tr><td>EI</td><td>×</td><td>27.89 ± 2.59</td><td>0.5755 ± 0.1285</td></tr><tr><td>IDFT</td><td></td><td>23.62 ± 1.90</td><td>0.5052 ± 0.0900</td></tr><tr><th rowspan="2">Method</th><th rowspan="2">Known</th><th colspan="2">MRI (×6 Accel., 10 dB SNR)</th></tr><tr><th>PSNR ↑</th><th>SSIM ↑</th></tr><tr><td>Supervised</td><td></td><td>27.39 ± 2.44</td><td>0.5243 ± 0.1373</td></tr><tr><td>ES (Ours)</td><td>✓</td><td>27.33 ± 2.45</td><td>0.5126 ± 0.1444</td></tr><tr><td>ES</td><td>×</td><td>25.73 ± 1.49</td><td>0.4566 ± 0.0622</td></tr><tr><td>EI</td><td>✓</td><td>27.23 ± 2.41</td><td>0.5110 ± 0.1421</td></tr><tr><td>EI</td><td>×</td><td>26.02 ± 1.65</td><td>0.4706 ± 0.0873</td></tr><tr><td>IDFT</td><td></td><td>23.85 ± 1.05</td><td>0.3878 ± 0.0272</td></tr></tbody></table>
<table><caption>Table 9: Empirical validation of the equivariance of MAP estimators.</caption><thead><tr><th></th><th>0°</th><th>90°</th><th>180°</th><th>270°</th><th>Average</th></tr></thead><tbody><tr><th>EQUIV</th><td>∞</td><td>∞</td><td>∞</td><td>∞</td><td>∞</td></tr></tbody></table>

<table><caption>Table 10: **Training durations.** For each training, we report the average epoch duration and the number of epochs until the model is trained. Inpainting trainings are conducted on a single NVIDIA RTX 3090 Ti GPU and MRI trainings on a NVIDIA H100 GPU.</caption><tbody><tr><th colspan="3">Image inpainting</th></tr><tr><th>Method</th><th>Epoch duration (s)</th><th>Epochs</th></tr><tr><td>Supervised</td><td>12</td><td>200</td></tr><tr><td>ES (Ours)</td><td>12</td><td>1000</td></tr><tr><td>EI</td><td>14</td><td>1000</td></tr><tr><th colspan="3">MRI (×8 Accel., 40 dB SNR)</th></tr><tr><th>Method</th><th>Epoch duration (s)</th><th>Epochs</th></tr><tr><td>Supervised</td><td>29</td><td>200</td></tr><tr><td>ES (Ours)</td><td>24</td><td>13800</td></tr><tr><td>EI</td><td>53</td><td>7800</td></tr><tr><td>SURE</td><td>36</td><td>5700</td></tr><tr><th colspan="3">MRI (×6 Accel., 10 dB SNR)</th></tr><tr><th>Method</th><th>Epoch duration (s)</th><th>Epochs</th></tr><tr><td>Supervised</td><td>19</td><td>70</td></tr><tr><td>ES (Ours)</td><td>19</td><td>3100</td></tr><tr><td>EI</td><td>75</td><td>2400</td></tr><tr><td>SURE</td><td>35</td><td>1200</td></tr></tbody></table>

## C PROOFS

For the sake of clarity, we state the propositions and theorems a second time before their proofs. We also state and prove the additional Lemma 1 which helps prove Theorem 1.

**Lemma 1.** *The minimization problem*

$$ \min_{f} \mathbb{E}_{x,y} \{ \| \mathbf{A} f(\mathbf{y}) - \mathbf{A} \mathbf{x} \|^2 \} \qquad (26) $$

admits as solutions the functions of the form

$$ f(\mathbf{y}) = \mathbf{A}^{\dagger} \mathbf{A} \mathbb{E}_{x|\mathbf{y}} \{\mathbf{x}\} + (\mathbf{I} - \mathbf{A}^{\dagger} \mathbf{A}) v(\mathbf{y}) \qquad (27) $$

where $v(\mathbf{y})$ is any function.

*Proof.* Let's start by stating that $f(\mathbf{y}) = \mathbb{E}_{x|\mathbf{y}} \{\mathbf{x}\}$ is the only solution of (Klenke, 2008)

$$ \min_{f} \mathbb{E}_{x,y} \{ \| f(\mathbf{y}) - \mathbf{x} \|^2 \}. \qquad (28) $$

For $f$ any solution of eq. (26), applying it to $\tilde{f}(\mathbf{y}) = \mathbf{A}f(\mathbf{y})$ and $\tilde{x} = \mathbf{A}\mathbf{x}$ gives

$$ \mathbf{A} f(\mathbf{y}) = \mathbf{A} \mathbb{E}_{x|\mathbf{y}} \{\mathbf{x}\}, \qquad (29) $$

and applying $f(\mathbf{y}) = \mathbf{A}^{\dagger} \mathbf{A} f(\mathbf{y}) + (\mathbf{I} - \mathbf{A}^{\dagger} \mathbf{A}) f(\mathbf{y})$ with $v(\mathbf{y}) := f(\mathbf{y})$ yields eq. (27). Conversely, the objective in eq. (26) has the same value no matter the $f$ satisfying eq. (27) and since at least one of them is solution of eq. (26), they all are. $\square$

**Theorem 1.** In the case of noiseless measurements with $p(\mathbf{x}) \mathcal{G}$-invariant (Assumption 1), if the matrix $Q_{\mathbf{A}_1} \triangleq \mathbb{E}_{g|\mathbf{A}_1} \{(\mathbf{A}T_g)^{\top} \mathbf{A}T_g\} \in \mathbb{R}^{n \times n}$ has full rank for some split $\mathbf{A}_1$, then the splitting method yields the same MMSE-optimal reconstructions as the supervised method, i.e.,

$$ f^{*}(\mathbf{y}_1, \mathbf{A}_1) = \mathbb{E}_{x|\mathbf{y}_1, \mathbf{A}_1} \{\mathbf{x}\}. \qquad (9) $$
Proof.

$$
\begin{align*}
\mathbb{E}_y\{\mathcal{L}_{\text{ES}}(\mathbf{y}, \mathbf{A}, f)\} \\
&= \mathbb{E}_y\{\mathbb{E}_g\{\mathbb{E}_{y_1, \mathbf{A}_1|\mathbf{y}, \mathbf{A}_{\text{T}_g}}\{\|\mathbf{A}_{\text{T}_g}f(\mathbf{y}_1, \mathbf{A}_1) - \mathbf{y}\|^2\}\}\} \\
&= \mathbb{E}_{y,g}\{\mathbb{E}_{y_1, \mathbf{A}_1|\mathbf{y}, g}\{\|\mathbf{A}_{\text{T}_g}f(\mathbf{y}_1, \mathbf{A}_1) - \mathbf{y}\|^2\}\} \\
&= \mathbb{E}_{y, y_1, \mathbf{A}_1, g}\{\|\mathbf{A}_{\text{T}_g}f(\mathbf{y}_1, \mathbf{A}_1) - \mathbf{y}\|^2\} \\
&= \mathbb{E}_{x, y_1, \mathbf{A}_1, g}\{\|\mathbf{A}_{\text{T}_g}f(\mathbf{y}_1, \mathbf{A}_1) - \mathbf{A}_{\text{T}_g}x\|^2\} \\
&= \mathbb{E}_{y_1, \mathbf{A}_1, x}\mathbb{E}_{g|\mathbf{y}_1, \mathbf{A}_1}\{\|\mathbf{A}_{\text{T}_g}(f(\mathbf{y}_1, \mathbf{A}_1) - x)\|^2\} \\
&= \mathbb{E}_{y_1, \mathbf{A}_1, x}\{(f(\mathbf{y}_1, \mathbf{A}_1) - x)^{\top}\mathbf{Q}_{\mathbf{A}_1}(f(\mathbf{y}_1, \mathbf{A}_1) - x)\},
\end{align*}
$$

where the third line use that $p(\mathbf{y}_1, \mathbf{A}_1 | \mathbf{y}, \mathbf{A}_{\text{T}_g}) = p(\mathbf{y}_1, \mathbf{A}_1 | \mathbf{y}, g)$ as $\mathbf{A}$ is fixed. The fifth line use the noiseless measurements assumption and the invariance of the distribution $p(\mathbf{x})$. The last line uses definition of $\mathbf{Q}_{\mathbf{A}_1}$. By applying Lemma 1, the global minimizer of the expected loss is given by:

$$f^*(\mathbf{y}_1, \mathbf{A}_1) = \mathbf{Q}_{\mathbf{A}_1}^{\dagger} \mathbf{Q}_{\mathbf{A}_1} \mathbb{E}_{x|\mathbf{y}_1, \mathbf{A}_1} \{\mathbf{x}\} + (\mathbf{I} - \mathbf{Q}_{\mathbf{A}_1}^{\dagger} \mathbf{Q}_{\mathbf{A}_1}) v(\mathbf{y}_1) \quad (30)$$

where $v: \mathbb{R}^n \to \mathbb{R}^n$ is any function. Moreover, since $\mathbf{Q}$ has full-rank, $\mathbf{Q}_{\mathbf{A}_1}^{\dagger} = \mathbf{Q}_{\mathbf{A}_1}^{-1}$ and then

$$f^*(\mathbf{y}_1, \mathbf{A}_1) = \mathbb{E}_{x|\mathbf{y}_1, \mathbf{A}_1} \{\mathbf{x}\}. \quad (31)$$

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

**Proposition 1.** If the matrix $\bar{\mathbf{Q}}_\mathbf{A} \triangleq \mathbb{E}_{\mathbf{A}_1|\mathbf{A}}\{\mathbf{Q}_{\mathbf{A}_1}\} \in \mathbb{R}^{n \times n}$ is invertible and $f$ minimizes $\mathbb{E}_y\{\mathcal{L}_{\text{ES}}(\mathbf{y}, \mathbf{A}, f)\}$. Then the reconstruction function

$$\bar{f}(\mathbf{y}, \mathbf{A}) \triangleq \mathbb{E}_{\mathbf{y}_1, \mathbf{A}_1|\mathbf{y}, \mathbf{A}} \{\bar{\mathbf{Q}}_{\mathbf{A}}^{-1} \mathbf{Q}_{\mathbf{A}_1} f(\mathbf{y}_1, \mathbf{A}_1)\} \quad (10)$$

satisfies

$$\bar{f}(\mathbf{y}, \mathbf{A}) = \mathbb{E}_{\mathbf{y}_1, \mathbf{A}_1|\mathbf{y}, \mathbf{A}} \{\bar{\mathbf{Q}}_{\mathbf{A}}^{-1} \mathbf{Q}_{\mathbf{A}_1} \mathbb{E}_{x|\mathbf{y}_1, \mathbf{A}_1} \{\mathbf{x}\}\}. \quad (11)$$

where eq. (11) is a convex combination of MMSE estimators for different splittings.

Proof. By applying eq. (30) to $f$ in the definition of $\bar{f}$ we obtain:

$$
\begin{align*}
\bar{f}(\mathbf{y}, \mathbf{A}) &= \mathbb{E}_{\mathbf{y}_1, \mathbf{A}_1|\mathbf{y}, \mathbf{A}} \left\{ \bar{\mathbf{Q}}_{\mathbf{A}}^{-1} \mathbf{Q}_{\mathbf{A}_1} \left( \mathbf{Q}_{\mathbf{A}_1}^{\dagger} \mathbf{Q}_{\mathbf{A}_1} \mathbb{E}_{x|\mathbf{y}_1, \mathbf{A}_1} \{\mathbf{x}\} + (\mathbf{I} - \mathbf{Q}_{\mathbf{A}_1}^{\dagger} \mathbf{Q}_{\mathbf{A}_1}) v(\mathbf{y}_1) \right) \right\} \\
&= \mathbb{E}_{\mathbf{y}_1, \mathbf{A}_1|\mathbf{y}, \mathbf{A}} \left\{ \bar{\mathbf{Q}}_{\mathbf{A}}^{-1} \mathbf{Q}_{\mathbf{A}_1} \mathbb{E}_{x|\mathbf{y}_1, \mathbf{A}_1} \{\mathbf{x}\} \right\} + \mathbb{E}_{\mathbf{y}_1, \mathbf{A}_1|\mathbf{y}, \mathbf{A}} \left\{ \bar{\mathbf{Q}}_{\mathbf{A}}^{-1} \mathbf{Q}_{\mathbf{A}_1} \left( (\mathbf{I} - \mathbf{Q}_{\mathbf{A}_1}^{\dagger} \mathbf{Q}_{\mathbf{A}_1}) v(\mathbf{y}_1) \right) \right\} \\
&= \mathbb{E}_{\mathbf{y}_1, \mathbf{A}_1|\mathbf{y}, \mathbf{A}} \left\{ \bar{\mathbf{Q}}_{\mathbf{A}}^{-1} \mathbf{Q}_{\mathbf{A}_1} \mathbb{E}_{x|\mathbf{y}_1, \mathbf{A}_1} \{\mathbf{x}\} \right\}.
\end{align*}
$$

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

**Corollary 1.** In order for the matrices $\mathbf{Q}_{\mathbf{A}_1}$ or $\bar{\mathbf{Q}}_\mathbf{A}$ to have full rank, it is necessary that $\mathbf{A}$ is not equivariant:

$$\exists g \in \mathcal{G}, \mathbf{A}_{\text{T}_g} \neq \mathbf{T}_g \mathbf{A}. \quad (13)$$

Proof. Let's assume by contradiction that $\mathbf{A}$ is equivariant with respect to $\mathbf{T}_g$

$$\mathbf{A}_{\text{T}_g} = \mathbf{T}_g \mathbf{A}, \quad (32)$$

and let $\mathbf{x} \in \ker(\mathbf{A})$.

$$
\begin{align*}
\mathbf{Q}_{\mathbf{A}_1} \mathbf{x} &= \left( \mathbb{E}_{g|\mathbf{A}_1} \{ (\mathbf{A}_{\text{T}_g})^{\top} \mathbf{A}_{\text{T}_g} \} \right) \mathbf{x} \\
&= \mathbb{E}_{g|\mathbf{A}_1} \{ (\mathbf{A}_{\text{T}_g})^{\top} \mathbf{A}_{\text{T}_g} \mathbf{x} \} \\
&= \mathbb{E}_{g|\mathbf{A}_1} \{ (\mathbf{A}_{\text{T}_g})^{\top} \mathbf{T}_g \mathbf{A} \mathbf{x} \} \\
&= \mathbb{E}_{g|\mathbf{A}_1} \{ (\mathbf{A}_{\text{T}_g})^{\top} \mathbf{T}_g \mathbf{0} \} \\
&= \mathbf{0}
\end{align*}
$$

Therefore,

$$\ker(\mathbf{Q}_{\mathbf{A}_1}) \supseteq \ker(\mathbf{A}) \supseteq \{\mathbf{0}\}. \quad (33)$$

The matrix $\mathbf{Q}_{\mathbf{A}_1}$ has a non-trivial nullspace and thus cannot have full rank. Moreover, since this non-trivial nullspace is the same for all virtual operators $\mathbf{A}_{\text{T}_g}$, then $\bar{\mathbf{Q}}_\mathbf{A}$ shares the same non-trivial nullspace. $\square$
**Theorem 2.** The reconstruction functions defined in points below are all equivariant as in eq. (15).

1. **Artifact removal network.** For a denoiser $\phi(\mathbf{x})$ equivariant in the sense of eq. (14),
$$f(\mathbf{y}, \mathbf{A}) = \phi(\mathbf{A}^{\top}\mathbf{y}), \text{ or } f(\mathbf{y}, \mathbf{A}) = \phi(\mathbf{A}^{\dagger}\mathbf{y}). \qquad (16)$$

2. **Unrolled network.** For $\phi(\mathbf{x})$ equivariant, any $\gamma \in \mathbb{R}$ and data fidelity $d(\mathbf{A}\mathbf{x}, \mathbf{y})$, with
$$\mathbf{x}_0 = \mathbf{0}, \quad \mathbf{x}_{k+1} = \phi(\mathbf{x}_k - \gamma \nabla_{\mathbf{x}_k} d(\mathbf{A}\mathbf{x}_k, \mathbf{y})) \qquad (17)$$
for $k = 0, \dots, L-1$ and $f(\mathbf{y}, \mathbf{A}) = \mathbf{x}_L$.

3. **Reynolds averaging.** For a possibly non-equivariant reconstructor $r(\mathbf{y}, \mathbf{A})$, with
$$f(\mathbf{y}, \mathbf{A}) = \frac{1}{|\mathcal{G}|} \sum_{g \in \mathcal{G}} T_g r(\mathbf{y}, \mathbf{A}T_g). \qquad (18)$$

4. **Maximum a posteriori (MAP).** For a distribution $p(\mathbf{x})$ invariant as in eq. (5), with
$$f(\mathbf{y}, \mathbf{A}) = \underset{\mathbf{x} \in \mathbb{R}^n}{\operatorname{argmax}} \{p(\mathbf{x} \mid \mathbf{y}, \mathbf{A})\}. \qquad (19)$$

5. **Minimum mean squared error (MMSE).** For a distribution $p(\mathbf{x})$ invariant as in eq. (5),
$$f(\mathbf{y}, \mathbf{A}) = \mathbb{E}_{\mathbf{x}|\mathbf{y}, \mathbf{A}} \{\mathbf{x}\}. \qquad (20)$$

*Proof.* We prove each case separately.

1. Denoting $\mathbf{A}^\times := \mathbf{A}^\top$ or $\mathbf{A}^\times := \mathbf{A}^\dagger$, eq. (16) gives, as $(\mathbf{A}T_g)^\times = T_g^{-1}\mathbf{A}^\times$,
$$f(\mathbf{y}, \mathbf{A}T_g) = \phi(T_g^{-1}\mathbf{A}^\times\mathbf{y}), \qquad (34)$$
and since $\phi(\mathbf{x})$ is equivariant, i.e., eq. (14) holds, it simplifies to eq. (15).

2. We start by making the notation show the explicit dependency on $\mathbf{y}$ and $\mathbf{A}$:
$$\mathbf{x}_0(\mathbf{y}, \mathbf{A}) = \mathbf{0}, \quad \mathbf{x}_{k+1}(\mathbf{y}, \mathbf{A}) = \phi(\mathbf{x}_k(\mathbf{y}, \mathbf{A}) - \gamma \nabla_{\mathbf{x}_k(\mathbf{y}, \mathbf{A})} d(\mathbf{A}\mathbf{x}_k(\mathbf{y}, \mathbf{A}), \mathbf{y})). \qquad (35)$$
and proceed to show that for $k = 0, \dots, L$ it holds that
$$\mathbf{x}_k(\mathbf{y}, \mathbf{A}T_g) = T_g^{-1}\mathbf{x}_k(\mathbf{y}, \mathbf{A}). \qquad (36)$$
For $k = 0$, it holds as
$$\mathbf{x}_0(\mathbf{y}, \mathbf{A}T_g) = \mathbf{0} = T_g^{-1}\mathbf{x}_0(\mathbf{y}, \mathbf{A}). \qquad (37)$$

Let's assume that eq. (36) holds for $k < L$. Applying it and the chain rule in eq. (35) yields
$$\mathbf{x}_{k+1}(\mathbf{y}, \mathbf{A}T_g) = \phi(T_g^{-1}(\mathbf{x}_k(\mathbf{y}, \mathbf{A}) - \gamma \nabla_{\mathbf{x}_k(\mathbf{y}, \mathbf{A})} d(\mathbf{A}\mathbf{x}_k(\mathbf{y}, \mathbf{A}), \mathbf{y}))). \qquad (38)$$

Finally, applying eq. (14) in this equation gives
$$\mathbf{x}_{k+1}(\mathbf{y}, \mathbf{A}T_g) = T_g^{-1}\mathbf{x}_{k+1}(\mathbf{y}, \mathbf{A}), \qquad (39)$$
and by induction, as $f(\mathbf{y}, \mathbf{A}) = \mathbf{x}_L(\mathbf{y}, \mathbf{A})$, eq. (15) holds.

3. From eq. (18), it holds that
$$f(\mathbf{y}, \mathbf{A}T_g) = \frac{1}{|\mathcal{G}|} \sum_{h \in \mathcal{G}} T_h r(\mathbf{y}, \mathbf{A}T_g T_h), \qquad (40)$$

which, as the group action property holds $T_g T_h = T_{gh}$, rewrites as

$$f(\mathbf{y}, \mathbf{A}T_g) = \frac{1}{|\mathcal{G}|} \sum_{h \in \mathcal{G}} T_h r(\mathbf{y}, \mathbf{A}T_{gh}), \qquad (41)$$

Applying the change of variable $h' = gh$ in this equation gives

$$f(\mathbf{y}, \mathbf{A}T_g) = \frac{1}{|\mathcal{G}|} \sum_{h \in \mathcal{G}} T_{g^{-1}h} r(\mathbf{y}, \mathbf{A}T_h), \qquad (42)$$

which finally, using group action property again $T_{g^{-1}h} = T_g^{-1}T_g$, gives eq. (15).
4. Taking the negative natural logarithm in eq. (19) and using Bayes' theorem gives

$$f(\mathbf{y}, \mathbf{A}) = \underset{\mathbf{x} \in \mathbb{R}^n}{\operatorname{argmin}} \left\{ d(\mathbf{A}\mathbf{x}, \mathbf{y}) + \rho(\mathbf{x}) \right\}, \quad (43)$$

where $d(\mathbf{A}\mathbf{x}, \mathbf{y}) = -\log p(\mathbf{y} | \mathbf{A}\mathbf{x})$ and $\rho(\mathbf{x}) = -\log p(\mathbf{x})$. Applying $\mathbf{x}' = \mathbf{T}_g \mathbf{x}$,

$$f(\mathbf{y}, \mathbf{A}\mathbf{T}_g) = \mathbf{T}_g^{-1} \underset{\mathbf{x} \in \mathbb{R}^n}{\operatorname{argmin}} \left\{ d(\mathbf{A}\mathbf{x}, \mathbf{y}) + \rho(\mathbf{T}_g^{-1}\mathbf{x}) \right\}, \quad (44)$$

and eq. (5) makes $\rho(\mathbf{x})$ invariant as well $\rho(\mathbf{T}_g^{-1}\mathbf{x}) = \rho(\mathbf{x})$. Therefore, eq. (15) holds.

5. Let's assume that $p(\mathbf{x})$ and $p(\mathbf{A})$ are invariant in the sense of eq. (5). We first prove that

$$p(\mathbf{y}, \mathbf{A}\mathbf{T}_g) = p(\mathbf{y}, \mathbf{A}). \quad (45)$$

Using eq. (5), the invariance of $p(\mathbf{A})$ and the independence of $\mathbf{x}$ and $\mathbf{A}$, we compute

$$\begin{aligned}
p(\mathbf{y}, \mathbf{A}\mathbf{T}_g) &= \mathbb{E}_{\mathbf{x}} \{p(\mathbf{y}, \mathbf{A}\mathbf{T}_g | \mathbf{x})\} = \mathbb{E}_{\mathbf{x}} \{p(\mathbf{y}, | \mathbf{A}\mathbf{T}_g, \mathbf{x})p(\mathbf{A}\mathbf{T}_g | \mathbf{x})\} \\
&= \mathbb{E}_{\mathbf{x}} \{p(\mathbf{y} | \mathbf{A}\mathbf{T}_g, \mathbf{x})p(\mathbf{A}\mathbf{T}_g)\} = \mathbb{E}_{\mathbf{x}} \{p(\mathbf{y} | \mathbf{A}\mathbf{T}_g\mathbf{x})p(\mathbf{A})\} \\
&= \mathbb{E}_{\mathbf{x}} \{p(\mathbf{y} | \mathbf{A}\mathbf{x})p(\mathbf{A})\} = \mathbb{E}_{\mathbf{x}} \{p(\mathbf{y} | \mathbf{A}, \mathbf{x})p(\mathbf{A} | \mathbf{x})\} \\
&= \mathbb{E}_{\mathbf{x}} \{p(\mathbf{y}, \mathbf{A} | \mathbf{x})\} = p(\mathbf{y}, \mathbf{A}).
\end{aligned}$$

Next we start from

$$f(\mathbf{y}, \mathbf{A}\mathbf{T}_g) = \mathbb{E}_{\mathbf{x}|\mathbf{y}, \mathbf{A}\mathbf{T}_g} \{\mathbf{x}\}. \quad (46)$$

and applying the integral formula for expectations gives

$$f(\mathbf{y}, \mathbf{A}\mathbf{T}_g) = \int \mathbf{x} p(\mathbf{x} | \mathbf{y}, \mathbf{A}\mathbf{T}_g) \, \mathrm{d}\mathbf{x}. \quad (47)$$

Using Bayes' formula and $p(\mathbf{y} | \mathbf{A}, \mathbf{x}) = p(\mathbf{y} | \mathbf{A}\mathbf{x})$, it becomes

$$f(\mathbf{y}, \mathbf{A}\mathbf{T}_g) = \int \mathbf{x} \frac{p(\mathbf{y} | \mathbf{A}\mathbf{T}_g\mathbf{x})}{p(\mathbf{A}\mathbf{T}_g, \mathbf{y})} p(\mathbf{A}\mathbf{T}_g, \mathbf{x}) \, \mathrm{d}\mathbf{x}. \quad (48)$$

By using eq. (45), the invariance of $p(\mathbf{A})$ and the independence of $\mathbf{A}$ with $\mathbf{x}$, we obtain

$$f(\mathbf{y}, \mathbf{A}\mathbf{T}_g) = \int \mathbf{x} \frac{p(\mathbf{y} | \mathbf{A}\mathbf{T}_g\mathbf{x})}{p(\mathbf{A}, \mathbf{y})} p(\mathbf{A})p(\mathbf{x}) \, \mathrm{d}\mathbf{x}. \quad (49)$$

With the change of variable $\mathbf{x}' = \mathbf{T}_g\mathbf{x}$, and since $\mathbf{T}_g$ is unitary, we arrive at

$$f(\mathbf{y}, \mathbf{A}\mathbf{T}_g) = \int \mathbf{T}_g^{-1}\mathbf{x} \frac{p(\mathbf{y} | \mathbf{A}\mathbf{x})}{p(\mathbf{A}, \mathbf{y})} p(\mathbf{A})p(\mathbf{T}_g^{-1}\mathbf{x}) \, \mathrm{d}\mathbf{x}. \quad (50)$$

Finally, applying eq. (5) in this equation yields eq. (15).

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

**Theorem 3.** If $f(\mathbf{A}, \mathbf{x})$ is an equivariant reconstructor, then ES is equivalent to the splitting loss

$$\mathcal{L}_{\mathrm{ES}}(\mathbf{y}, \mathbf{A}, f) = \mathcal{L}_{\mathrm{SPLIT}}(\mathbf{y}, \mathbf{A}, f). \quad (21)$$

*Proof.* We start from eq. (7)

$$\mathcal{L}_{\mathrm{ES}}(\mathbf{y}, \mathbf{A}, f) \triangleq \mathbb{E}_g \{ \mathcal{L}_{\mathrm{SPLIT}}(\mathbf{y}, \mathbf{A}\mathbf{T}_g, f) \} \quad (51)$$

$$= \mathbb{E}_g \{ \mathbb{E}_{\mathbf{y}_1, \mathbf{A}_1 | \mathbf{y}, \mathbf{A}\mathbf{T}_g} \{ \| \mathbf{A}\mathbf{T}_g f(\mathbf{y}_1, \mathbf{A}_1) - \mathbf{y} \|^2 \} \} \quad (52)$$

As $\mathbf{A}_1$ is a splitting of $\mathbf{A}\mathbf{T}_g$, we can write $\mathbf{A}_1 = \mathbf{M}\mathbf{A}\mathbf{T}_g$ for $\mathbf{M}$ a splitting matrix. We obtain

$$\mathcal{L}_{\mathrm{ES}}(\mathbf{y}, \mathbf{A}, f) = \mathbb{E}_g \{ \mathbb{E}_{\mathbf{M}|\mathbf{y},g} \{ \| \mathbf{A}\mathbf{T}_g f(\mathbf{M}\mathbf{y}, \mathbf{M}\mathbf{A}\mathbf{T}_g) - \mathbf{y} \|^2 \} \} \quad (53)$$

Applying eq. (15) and cancelling out $\mathbf{T}_g$ with $\mathbf{T}_g^{-1}$ yields

$$\mathcal{L}_{\mathrm{ES}}(\mathbf{y}, \mathbf{A}, f) = \mathbb{E}_g \{ \mathbb{E}_{\mathbf{M}|\mathbf{y},g} \{ \| \mathbf{A} f(\mathbf{M}\mathbf{y}, \mathbf{M}\mathbf{A}) - \mathbf{y} \|^2 \} \}. \quad (54)$$

By dropping the expectation in $g$ and rewriting $(\mathbf{M}\mathbf{y}, \mathbf{M}\mathbf{A})$ as $(\mathbf{y}_1, \mathbf{A}_1)$, this yields in eq. (21).

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