# Visual Saliency Detection in Advanced Driver Assistance Systems

Francesco Rundo  ADG, Central R&D
STMicroelectronics
Catania, Italy
<EMAIL_ADDRESS>Michael Sebastian Rundo  Perceive Lab
University of Catania
Catania, Italy
<EMAIL_ADDRESS>Concetto Spampinato  Perceive Lab
University of Catania
Catania, Italy
<EMAIL_ADDRESS>
###### Abstract

Visual Saliency refers to the innate human mechanism of focusing on and
extracting important features from the observed environment. Recently, there
has been a notable surge of interest in the field of automotive research
regarding the estimation of visual saliency. While operating a vehicle,
drivers naturally direct their attention towards specific objects, employing
brain-driven saliency mechanisms that prioritize certain elements over others.
In this investigation, we present an intelligent system that combines a
drowsiness detection system for drivers with a scene comprehension pipeline
based on saliency. To achieve this, we have implemented a specialized 3D deep
network for semantic segmentation, which has been pretrained and tailored for
processing the frames captured by an automotive-grade external camera. The
proposed pipeline was hosted on an embedded platform utilizing the STA1295
core, featuring ARM A7 dual-cores, and embeds an hardware accelerator.
Additionally, we employ an innovative biosensor mounted on the steering wheel
to monitor the driver’s drowsiness, gathering the PhotoPlethysmoGraphy (PPG)
signal. A dedicated 1D temporal deep convolutional network has been devised to
classify the collected PPG time-series, enabling us to assess the driver’s
level of attentiveness. Ultimately, we compare the determined attention level
of the driver with the corresponding saliency-based scene classification to
evaluate the overall safety level. The efficacy of the proposed pipeline has
been validated through extensive experimental results.

###### Index Terms:

Drowsiness, Deep learning, D-CNN, Deep-LSTM, PPG (PhotoPlethySmography)

##  I Introduction

Drowsiness denotes a physiological state characterized by reduced
consciousness and challenges in maintaining wakefulness, posing a significant
risk for road accidents. Consequently, numerous studies have explored methods
to monitor the driver’s attentiveness. These investigations have identified a
correlation between human attention, often assessed through drowsiness
monitoring, and the Heart Rate Variability (HRV) index. HRV quantifies heart
activity across consecutive beats, yielding valuable insights into drowsiness
levels [1]. Specifically, HRV serves as a non-invasive means of assessing the
activity of the Autonomous Nervous System (ANS) [1]. The scientific literature
extensively recognizes the link between an individual’s attention level and
ANS activity [1]-[3]. Hence, by analyzing HRV, we can indirectly gauge the
subject’s level of attentiveness. However, it is crucial to consider the
current driving situation when evaluating the driver’s attention. For
instance, driving in low-traffic and low-speed conditions imposes less
cognitive demand compared to scenarios involving risky maneuvers such as
overtaking or lane changes. This study introduces an algorithm that integrates
saliency analysis, specifically for comprehending driving scenes, with
physiological monitoring of the car driver to assess their attention level.
Saliency analysis has found extensive utility across various domains,
including automotive applications, enabling the identification of crucial
elements within videos or images [1, 2].However, the algorithms for analyzing
driving scenes and monitoring drowsiness impose substantial computational
demands. Through the proposed saliency-based approach, it becomes possible to
characterize the driving scene by focusing solely on salient elements such as
vehicles and pedestrians. This not only reduces the computational load but
also enhances algorithmic performance. The subsequent sections of this paper
are structured as follows: Section II provides a review of previous works in
the field, exploring the utilization of physiological signal analysis for
automated drowsiness detection. Section III presents a detailed description of
the constituent blocks comprising our pipelines. Section IV outlines our
experimental setup, including information regarding the dataset and hardware
devices utilized in implementing our pipeline. Finally, Section V concludes
the paper, summarizing the findings and outlining potential directions for
future research.

Figure 1: The overall scheme of the proposed pipeline.

##  II Related Works

Numerous research endeavors have focused on developing effective methods for
simultaneously monitoring the driver’s attention level and analyzing driving
scenarios. In a notable study [3], the researchers devised a pipeline for
detecting driver drowsiness by capturing the ElectroCardioGraphy (ECG) signal
and analyzing alterations in Heart Rate Variability (HRV). Through ECG signal
stabilization and employing classical linear discriminant analysis, the
pipeline successfully discriminates between drowsy and wakeful states.
However, a key challenge in drowsiness detection lies in ECG sampling. The
minimal setup for acquiring a subject’s ECG signal necessitates the use of at
least three electrodes in contact with the skin, following the configuration
known as Einthoven’s Triangle [4]. Consequently, the driver must place their
hands on electrodes embedded in the steering wheel to capture the ECG signal,
with a third electrode typically positioned in the driver’s seat.
Unfortunately, this placement increases susceptibility to noise, thereby
affecting the reliability of the HRV signal derived from the collected ECG
signal [4]. To address these limitations, several researchers have turned to
analyzing the PhotoPlethysmoGraphic (PPG) signal, which offers a more
favorable alternative. The PPG signal requires only a single contact point for
sampling, making it more practical compared to ECG [5]. In another study [5],
he authors explored parasympathetic nervous activity by extracting HRV
features from the PPG signal, aiming to classify the driver’s fatigue level.
The obtained results demonstrated the robustness of their proposed approach.
Moreover, in a separate study [6], the authors put forth a reliable indicator
derived from the PPG signal. In previous studies [1][2], researchers have
reported encouraging outcomes by employing algorithms that leverage Pulse Rate
Variability (PRV) data processing. These investigations have demonstrated the
utility of PRV as a valuable metric for analyzing the Autonomous Nervous
System (ANS) and evaluating subject drowsiness. Another noteworthy study by
Ryu et al. [3] employed red organic light-emitting diodes (OLEDs) and organic
photodiodes (OPDs) to capture the photoplethysmographic (PPG) signal. The
evaluation of this flexible PPG sensor showcased its effectiveness in
detecting drowsiness, surpassing the performance of standard PPG probes in
terms of accuracy. In recent years, Deep Learning approaches have garnered
significant attention for drowsiness estimation, employing bio-signals and
image analysis. For instance, in studies by Hong et al. [4] and Alshaqaqi et
al. [5], methods based on image processing were implemented to track the
driver’s face, gaze, and emotions to assess their attentiveness and
drowsiness. However, it is important to note that image processing performance
can be influenced by various factors in the in-vehicle environment, such as
lighting conditions and occlusions, which can impact the accuracy of the
analysis. As highlighted, the automotive industry’s continual technological
advancements have led to the emergence of promising Drowsiness Detection
systems, as exemplified by the pipeline proposed in [6]. Furthermore, a range
of Deep Learning algorithms has been devised by researchers to assess the
level of fatigue by employing ”data fusion” techniques that combine visual and
physiological information from the car driver. Altun and Celenk [1] introduced
a vision-based driver assistance system focused on scene awareness, leveraging
visual saliency analysis. Their results demonstrated the system’s
effectiveness in enhancing the driver’s situational awareness (SA) by
facilitating adaptive road surface classification. In another study [2], Deng
et al. investigated driver drowsiness by conducting eye-tracking experiments
during driving. They recruited a group of 40 subjects, including both non-
drivers and experienced individuals. Their findings indicated that drivers
exhibit heightened attention towards the vanishing points on the road. This
observation suggests that the road’s vanishing point could serve as a valuable
cue for estimating the road area in a traffic saliency detection mode. To
incorporate drowsiness monitoring, the authors proposed an integrated solution
that extracts relevant information regarding the driving scenario [2].

##  III Methods and Materials

This section delves into the finer aspects of the proposed pipeline,
illustrated in Figure 1.

###  III-A The Attention Monitoring System

Our study involved the implementation of a system designed to capture the car
driver’s Photoplethysmographic (PPG) signal, which enables us to assess and
monitor their attention level. As mentioned earlier, the PPG signal offers a
non-invasive means of analyzing the heart’s pulse rate. Notably, it provides
valuable insights into cardiac and vascular conditions, as well as respiratory
rate [5]. The PPG waveform comprises two key components: an ’AC’ physiological
signal that correlates with cardiac-synchronous changes in blood volume, and a
’DC’ component that reflects minor variations resulting from respiration and
thermoregulation processes [5]. As the heart propels blood towards the
periphery, it generates pressure, leading to the expansion of arteries and
arterioles in the subcutaneous tissue. To capture the Photoplethysmographic
(PPG) signal, we employed a device comprising a light-emitting component and a
detector that made contact with the subject’s skin. By illuminating the skin
and gauging the amount of light scattered back, we could effectively detect
the fluctuations in volume that occur with each heartbeat. These variations
are manifested as peaks within the PPG waveform [1]. The heart’s pressure
pulses, as it pumps blood to the periphery, induce specific distension within
the subcutaneous arteries and arterioles. Consequently, our approach involved
employing a device with a light-emitting component and a detector in direct
contact with the subject’s skin, facilitating the sampling of the PPG signal.
By illuminating the skin and accurately measuring the amount of light
scattered back, we were able to identify the changes in volume associated with
each heart pressure pulse. These changes are graphically represented as peaks
within the PPG waveform. [5]. We used the PPG sampling device composed of the
Silicon Photomultiplier sensor [7, 8, 9, 10]. The proposed PPG probes show an
array detector device, called Silicon Photomultipliers (SiPMs) [8], with a
total area of $4.0\times 4.5$ $mm^{2}$ and $4871$ square microcells with $60$
µm pitch. The devices have a geometrical fill factor of $67.4\%$ and are
packaged in a surface mount housing (SMD) with about $5.1\times 5.1$ $mm^{2}$
total area [11]. We propose a Pixelteq dichroic bandpass filter with a pass-
band centered at about $540$ nm with a Full Width at Half Maximum (FWHM) of
$70$ nm and an optical transmission higher than $90-95\%$ in the pass-band
range was glued on the SMD package by using a Loctite 352TM adhesive. The SiPM
has a maximum detection efficiency of about $30\%$ at 565 nm and a PDE of
about $27.5\%$ at $540$ nm (central wavelength in the filter pass-band). As
described, the PPG detector is composed of a light emitter in combination with
a detector based on SiPM technology. The OSRAM LT M673 LEDs were used by SMD
package and InGan technology [11]. The used LEDs devices are characterized by
an area of $2.3\times 1.5$ $mm^{2}$ with a 120° angle view, a spectral
bandwidth of $33$ nm and a lower power emission (mW) in the standard range. In
order to enhance the efficiency and utility of the PPG probe, we developed a
printed circuit board (PCB) with a user-friendly interface utilizing National
Instruments (NI) instrumentation. The PCB incorporates various components such
as a 4V portable battery, power management circuits, a conditioning circuit
for output SiPMs (Silicon Photomultipliers) signals, and multiple USB
connectors to accommodate PPG probes. Additionally, the board features SMA
(SubMiniature version A) output connectors for seamless connectivity and
integration with other devices and systems.

TABLE I: HYPER LOW-PASS FILTERING SETUP (IN HZ)

F | f1 | f2 | f3 | f4 | f5 | f6 | f7 | f8 | f9 | f10 | f11
---|---|---|---|---|---|---|---|---|---|---|---
HP | 0.5 | / | / | / | / | / | / | / | / | / | /
LP | 0 | 1.4 | 2.9 | 2.5 | 3.8 | 3.9 | 4 | 4.5 | 5 | 5.3 | 6.9

TABLE II: HYPER HIGH-PASS FILTERING SETUP (IN HZ)

F | f1 | f2 | f3 | f4 | f5 | f6 | f7 | f8 | f9 | f10 | f11
---|---|---|---|---|---|---|---|---|---|---|---
HP | 0.5 | 1.2 | 2.6 | 2.7 | 3.3 | 3.5 | 4 | 4.4 | 5 | 5.7 | 6.4
LP | 7 | / | / | / | / | / | / | / | / | / | /

In our previous works [1][2][3][4], we provided further insights into the
hardware setup utilized for acquiring the Photoplethysmographic (PPG) signal.
The PPG Sensor Probe consists of the SiPM sensor along with the aforementioned
LEDs. Effective management of the SiPM device’s power consumption is
facilitated by a dedicated Power management circuit [1][2][3]. As previously
mentioned, multiple PPG probes were strategically positioned on the steering
wheel of the vehicle to capture the PPG signal. It suffices for the driver to
place only one hand over the PPG sensor to collect the physiological signal
accurately. Once the raw PPG signal is acquired, the National Instruments (NI)
device processes the collected data using an internal $24$ bit analog-to-
digital converter (ADC). Additionally, the NI device is equipped with a
Windows-based operating system, incorporating a LabView software framework
[4]. To preprocess the acquired PPG raw data, we devised a LabView algorithm
that incorporates several key steps. Firstly, we applied ad-hoc Finite Impulse
Response (FIR) filtering techniques, encompassing both low-pass and high-pass
bands, to refine the signal quality. Additionally, we performed further
processing by calculating the first and second derivatives of the PPG signal,
enabling us to evaluate the minimum and maximum values for each waveform
accurately. Lastly, we visualized the processed PPG signal on the monitor
connected to the National Instruments (NI) device, providing real-time
feedback and visualization of the results. Further details regarding the
National Instruments (NI) device and the LabView software framework can be
found in our previous works [1][2]. However, it is important to note that the
final implementation of the proposed pipeline will be transitioned to the
STA1295 embedded platform, replacing the NI device utilized solely for
development purposes [3][4]. This transition ensures a more streamlined and
optimized solution. For the purpose of this study, we specifically developed
an efficient pipeline for processing the raw Photoplethysmographic (PPG) data
using the MATLAB framework, as described in the validation session of our
research. This approach allowed us to effectively analyze and validate the
performance of our pipeline. Table I and II present comprehensive information
about the hyper-filtering configuration employed for processing the recorded
raw Photoplethysmographic (PPG) signal [1][2][3]. The hyper-filtering
technique, introduced by the authors, draws inspiration from hyper-spectral
imaging and applies it to signal processing [3]. By employing hyper-filtering,
the source PPG signal undergoes filtering and processing with various
frequency setups (as outlined in the tables) to extract discriminative
patterns associated with the drowsiness level of the subject under analysis.
Further details regarding this approach can be found in [3]. Once the hyper-
filtered PPG signal patterns have been collected [2], they are classified
using a specialized 1D Temporal Dilated Convolutional Neural Network (1D-CNN).
Our proposed Temporal Convolutional Network architecture incorporates a
Dilated Causal Convolution layer that operates on the temporal sequence of
input data. The ”causal” aspect ensures that the activations at time $t$ are
solely influenced by inputs from time $t-1$. The network architecture
comprises multiple residual blocks, each consisting of two dilated causal
convolution layers with a consistent dilation factor, followed by
normalization, ReLU activation, and spatial dropout layers. Specifically, our
implemented 1D-CNN consists of 12 blocks, culminating in a downstream softmax
layer. Each block consists of a dilated convolution layer with $3\times 3$
kernel filters, a spatial dropout layer, another dilated convolution layer, a
ReLU layer, and a final spatial dropout layer. The initial dilation size is
set to 2, incrementing it for each subsequent block [3]. The objective of this
study is to predict the drowsiness level of the car driver by processing the
hyper-filtered signal patterns derived from the acquired PPG signal using the
designed 1D-CNN architecture. The network’s output is a scalar value ranging
from 0 to 1, indicating the driver’s attention level and distinguishing
between drowsy (0) and wakeful (1) states. Our findings strongly support the
feasibility of accurately estimating driver drowsiness. To ensure effective
support during driving, we have successfully ported the proposed 1D Deep CNN
backbone to an automotive-grade STA1295 Accordo5 embedded platform [4]

Figure 2: Saliency analysis of the video representing the driving scene.

###  III-B The Video Saliency Scene Understanding Block

This section provides detailed insights into the architectural aspects of the
proposed Deep Network, specifically designed for saliency analysis in driving
scenarios.

To capture the driving scenario, we utilized an automotive-grade camera
device, which enabled the acquisition of video frames. Subsequently, we
employed a custom-designed 3D-to-2D Semantic Segmentation Fully Convolutional
Network (SS-FCN) to process the collected video frames and extract saliency
information. The architecture of our proposed segmentation model follows an
encoder/decoder structure, facilitating efficient processing of the driving
scene. The resulting saliency map highlights the most prominent object within
the captured driving scenario. Our SS-FCN architecture consists of two primary
components: an encoder and a decoder. The encoder block, known as 3D Enc Net,
focuses on extracting spatiotemporal features. It comprises five blocks, with
each block consisting of convolutional layers grouped into two sets of two
layers. The kernel size for each separable convolutional layer is 3x3x3.
Within each block, a sequence of two 3x3x3 convolutional operations is
followed by batch normalization, a Rectified Linear Unit (ReLU) layer, and a
max pooling operation with a pooling size of 1x2x2. This sequence is repeated
twice. Subsequently, another sequence of 3x3x3 convolutional layers is
performed. The remaining three blocks follow a similar pattern, with a
succession of two convolutional operations, followed by batch normalization,
another convolutional layer with a 3x3x3 kernel, batch normalization, and ReLU
activation, culminating in a downstream max-pooling layer with a size of
1x2x2. The decoder backbone, referred to as 2D Dec Net, comprises five blocks
and incorporates up-sampling layers for feature map decoding. It consists of
2D convolutional layers with a kernel size of 3x3. After each convolutional
layer, batch normalization and ReLU activation are applied. Residual
connections are introduced through convolutional blocks. The decoder section
is responsible for adjusting the size of the feature maps through up-sampling
operations. Ultimately, the output of our proposed SS-FCN architecture is the
feature map corresponding to the segmented area of the most salient object
within the driving scenario’s video frame. Semantic segmentation techniques
often highlight the most visually salient objects in an image, typically
represented by fixation points. To illustrate this, we present the saliency
maps of the driving scenarios and the corresponding training loss dynamics of
the SS-FCN model in Figures 3 and 2, respectively. For training and testing
the SS-FCN model, we utilized the DHF1K dataset [1]. Performance evaluation of
the model was conducted using several performance indexes, including the Area
Under the Curve (AUC), Similarity, Correlation Coefficient, and Normalized
Scanpath Saliency. The proposed architecture demonstrated satisfactory
performance, with the following results: AUC: 0.885, Similarity: 0.355,
Correlation Coefficient: 0.455, and Normalized Scanpath Saliency: 2.564. These
values compare favorably with similar architectures found in the literature
[1]. It is worth noting that while other architectures may exhibit superior
performance, they often require high computational resources and rely on
sophisticated hardware. In contrast, our proposed pipeline achieves notable
results while maintaining a low computational workload and without the need
for specific hardware accelerations [2].

Figure 3: An instance of saliency processing made by fully trained SS-FCN. The
network training loss dynamic is reported

Figure 4: (a) Static driving scenario requiring low attention level. The
Driver Attention Analyzer confirms an adequate level of attention (b) Dynamic
driving scenario (car overtaking) requiring High attention level. The Driver
Attention Analyzer detects an inadequate level of attention generating an
acoustic signal alert.

###  III-C The Driver Attention Analyzer

This component serves to evaluate the driver’s attention level by combining
the analysis of the PPG signal with the output of the proposed SS-FCN block,
namely the saliency map. Its main objective is to ensure that the detected
attention level aligns with the meaningful interpretation of the corresponding
salience map, which represents the dynamic nature of the driving scenario. A
static saliency map indicates a relatively low level of dynamic activity in
the driving scenario, suggesting a corresponding lower level of attention from
the driver. On the other hand, an increasing variation in the saliency maps
enable a greater demand for the driver’s attention. Hence, this block verifies
the coherence between the driver’s attention level and the saliency map,
reflecting the dynamic characteristics of the driving scenario [1]. We set ad-
hoc thresholds to verify the level of attention obtained by the Drowsiness
Monitoring systems. Both thresholds are determined during the training phase
by means of a heuristic calibration i.e. by choosing a setup that maximizes
the performance of the overall pipeline during the learning.

##  IV Experimental Results

In this section, we begin by introducing the dataset and then proceed to
discuss the training details and the key findings of our study. To evaluate
the effectiveness of our implemented pipeline, we utilized the DH1FK dataset
[1]. Additionally, we captured video sequences of various driving scenarios
using a camera with a resolution of 2.3 megapixels and a maximum frame rate of
60 frames per second. The recruitment of subjects was supervised by a group of
physiologists who ensured that the recruited individuals exhibited diverse
levels of attention while simultaneously collecting the corresponding EEG
signal, thereby corroborating their physiological responses [2]. Our
experiments involved a total of 43 subjects ranging in age from 21 to 70
years. We recorded the PPG signal from each subject at a sampling frequency of
1 KHz, with a data collection duration of 5 minutes. For training purposes, we
allocated 70% of the acquired PPG time-series and video driving scene frames,
reserving the remaining 30% for testing and validation. In summary, our
proposed pipeline aims to assess the driver’s attention level by combining the
PPG signal and the saliency map. Specifically, the Driver Attention Analyzer
block compares the level of attention required by the driving scenario with
the output of the Drowsiness Monitoring System. If the driver’s vigilance
level falls below the requirement of the current driving scenario, the system
alerts the driver through an acoustic signal. To classify the attention level,
we defined two value ranges. A value within the range of 0 to 0.6, which is
the output of the 1D-CNN-based Drowsiness Monitoring System, indicates a
medium-low attention level. Conversely, values ranging from 0.61 to 1
represent a high attention level. Finally, we set an ad-hoc normalized
threshold $\vartheta$ (0.45) to define a static scene-based saliency map (see
Eqs. (1)). The driving scene is considered as ‘dynamic’ if the values of the
normalized saliency map gradient are greater than $0.45$ (requiring a high
level of attention). For a normalized saliency map gradient lower than $0.45$
(requiring a low level of attention) the driving scene is considered ‘static’.
In Fig.4, two instances of the driving scenarios are reported.

##  V Conclusions and Discussion

The promising results obtained from our study highlight the potential for
assessing driver drowsiness and enhancing driving safety. A key advantage of
our proposed method is its independence from frequency domain analysis, unlike
other approaches that rely on HRV analysis [1]. Additionally, our method
leverages the readily available PPG signal, which can be conveniently sampled
from a range of sensors placed on the steering wheel. To evaluate the driver’s
level of attention based on PPG signal analysis, we developed a fully
convolutional deep network that generates a saliency map of the driving scene.
The results demonstrate the effectiveness of our pipeline in assessing the
required level of driver drowsiness specific to the driving scenario. The
system compares the attention level reconstructed from the driver’s PPG signal
with the calibrated attention level derived from saliency analysis. If a
discrepancy is detected between these attention levels, indicating a potential
risk, the system alerts the driver. Ongoing investigations aim to enhance the
saliency analysis by expanding the dataset to include diverse driving
scenarios from various domains. Currently, we are in the process of porting
the entire pipeline to an automotive-specific board powered by the SoC STA1295
ACCORDO 5 processor manufactured by STMicroelectronics [2]. Furthermore, our
future endeavors involve incorporating more robust domain adaptation methods,
encompassing both supervised and unsupervised approaches that have been
successfully employed in other pipelines [3][4][5][6].

## Acknowledgment

The used dataset was collected under the clinical study Ethical Committee CT1
authorization n.113 / 2018 / PO. This research was funded by the National
Funded Program 2014-2020 under grant agreement n. 1733, (ADAS + Project).

## References

  * [1] Y. Cai, Z. Liu, H. Wang, and X. Sun, “Saliency-based pedestrian detection in far infrared images,” _IEEE Access_ , vol. 5, pp. 5013–5019, 2017.
  * [2] A. Winterlich, V. Zlokolica, P. Denny, L. Kilmartin, M. Glavin, and E. Jones, “A saliency weighted no-reference perceptual blur metric for the automotive environment,” in _2013 Fifth International Workshop on Quality of Multimedia Experience (QoMEX)_. IEEE, 2013, pp. 206–211.
  * [3] J. Vicente, P. Laguna, A. Bartra, and R. Bailón, “Detection of driver’s drowsiness by means of hrv analysis,” in _2011 Computing in Cardiology_. IEEE, 2011, pp. 89–92.
  * [4] B. Abi-Saleh and B. Omar, “Einthoven’s triangle transparency: a practical method to explain limb lead configuration following single lead misplacements,” _Reviews in cardiovascular medicine_ , vol. 11, no. 1, pp. 33–38, 2019.
  * [5] F. Rundo, S. Conoci, A. Ortis, and S. Battiato, “An advanced bio-inspired photoplethysmography (ppg) and ecg pattern recognition system for medical assessment,” _Sensors_ , vol. 18, no. 2, p. 405, 2018.
  * [6] F. Rundo, S. Petralia, G. Fallica, and S. Conoci, “A nonlinear pattern recognition pipeline for ppg/ecg medical assessments,” in _Convegno Nazionale Sensori_. Springer, 2018, pp. 473–480.
  * [7] V. Vinciguerra, E. Ambra, L. Maddiona, M. Romeo, M. Mazzillo, F. Rundo, G. Fallica, F. di Pompeo, A. M. Chiarelli, F. Zappasodi _et al._ , “Ppg/ecg multisite combo system based on sipm technology,” in _Convegno Nazionale Sensori_. Springer, 2018, pp. 353–360.
  * [8] M. Mazzillo, L. Maddiona, F. Rundo, A. Sciuto, S. Libertino, S. Lombardo, and G. Fallica, “Characterization of sipms with nir long-pass interferential and plastic filters,” _IEEE Photonics Journal_ , vol. 10, no. 3, pp. 1–12, 2018\.
  * [9] F. Rundo, S. Rinella, S. Massimino, M. Coco, G. Fallica, R. Parenti, S. Conoci, and V. Perciavalle, “An innovative deep learning algorithm for drowsiness detection from eeg signal,” _Computation_ , vol. 7, no. 1, p. 13, 2019.
  * [10] F. Rundo, C. Spampinato, and S. Conoci, “Ad-hoc shallow neural network to learn hyper filtered photoplethysmographic (ppg) signal for efficient car-driver drowsiness monitoring,” _Electronics_ , vol. 8, no. 8, p. 890, 2019.
  * [11] K. Fujiwara, E. Abe, K. Kamata, C. Nakayama, Y. Suzuki, T. Yamakawa, T. Hiraoka, M. Kano, Y. Sumi, F. Masuda _et al._ , “Heart rate variability-based driver drowsiness detection and its validation with eeg,” _IEEE Transactions on Biomedical Engineering_ , vol. 66, no. 6, pp. 1769–1778, 2018.


# Equivariant plug-and-play image reconstruction

Matthieu Terris1, Thomas Moreau1, Nelly Pustelnik2, Julian Tachella2
1Université Paris-Saclay, Inria, CEA, Palaiseau, 91120, France; 2Laboratoire
de Physique, ENSL, CNRS UMR 5672, F-69342, Lyon, France.

###### Abstract

Plug-and-play algorithms constitute a popular framework for solving inverse
imaging problems that rely on the implicit definition of an image prior via a
denoiser. These algorithms can leverage powerful pre-trained denoisers to
solve a wide range of imaging tasks, circumventing the necessity to train
models on a per-task basis. Unfortunately, plug-and-play methods often show
unstable behaviors, hampering their promise of versatility and leading to
suboptimal quality of reconstructed images. In this work, we show that
enforcing equivariance to certain groups of transformations (rotations,
reflections, and/or translations) on the denoiser strongly improves the
stability of the algorithm as well as its reconstruction quality. We provide a
theoretical analysis that illustrates the role of equivariance on better
performance and stability. We present a simple algorithm that enforces
equivariance on any existing denoiser by simply applying a random
transformation to the input of the denoiser and the inverse transformation to
the output at each iteration of the algorithm. Experiments on multiple imaging
modalities and denoising networks show that the equivariant plug-and-play
algorithm improves both the reconstruction performance and the stability
compared to their non-equivariant counterparts.

##  1 Introduction

Linear inverse imaging problems are ubiquitous in imaging sciences, famous
instances of which include image restoration, magnetic resonance imaging
(MRI), computed tomography, and astronomical imaging to name a few. In this
setting, the aim is to recover an image $x\in\mathbb{R}^{n}$ from measurements
$y\in\mathbb{R}^{m}$ acquired through

$y=Ax+\epsilon,$ (1)

where $A\colon\mathbb{R}^{n}\to\mathbb{R}^{m}$ is a linear operator and
$\epsilon$ the realisation of some random noise. A myriad of methods for
solving problems of the likes of (1) have been proposed in the literature,
ranging from variational methods [12, 6] (solution of a cost function
minimisation problem) to end-to-end reconstruction with deep neural networks
[19, 22, 50] and, more recently, diffusion algorithms [9, 21, 52, 14].

In this work, we focus on approaches relying on implicit denoising priors. For
instance, plug-and-play (PnP) algorithms [43, 34, 48, 20] propose to replace
the proximity operator involved in variational methods (solved by proximal
algorithms) by a denoiser modeling an implicit image prior. Similarly,
regularization-by-denoising (RED) approaches propose to replace gradient steps
on the prior term by a denoiser [32]. While the denoiser is typically trained
as a Gaussian denoiser using grayscale or color natural images, it can be
plugged into algorithms designed to tackle a wide range of image-related
problems, without being constrained by the nature of the input images (e.g.
MRI images, CT scans, etc) [1, 26, 41, 38].

Figure 1: Instability of algorithms relying on implicit denoising priors can
be solved by incorporating equivariance. Enforcing approximate equivariance of
the denoiser at test time allows to both stabilize the algorithm and to
improve the reconstruction quality without needing to retrain the implicit
prior. Left: PnP algorithm applied to an accelerated MRI problem. Middle:
Unadjusted Langevin sampling algorithm for a motion blur problem; estimated
mean and variance of the associated Markov chain are displayed. Right: RED
algorithm on a $4\times$ super-resolution problem.

Yet, despite successful performance in various applications, these algorithms
tend to suffer from instability issues, making it difficult to transfer across
imaging tasks or to derive statistical estimates from the reconstruction. In
particular, RED and PnP algorithms often require careful finetuning [31, 27],
or departing from the original optimization algorithm for efficient
application [48, 10]. Thus, several works have been proposed to restore the
convergence of PnP and RED algorithms while establishing a clear connection
with an associated cost function [34, 40, 10, 29, 17]. However, these efforts
introduce notable constraints on the denoiser, leading to a trade-off between
stability and reconstruction performance. Typically, convergent PnP and RED
algorithms tend to transfer easily to new problems, but perform less well than
their fine-tuned and early-stopped non-convergent counterparts.

In this work, we propose to investigate theoretically and empirically the
effect of equivariance on algorithms relying on implicit denoising priors.
More precisely, we prove that enforcing equivariance on the denoiser improves
the stability of the resulting PnP algorithm and plays a symmetrizing role on
the Jacobian of the implicit prior. Our experiments show that the proposed
equivariant approach can significantly improve the quality of images
reconstructed with PnP algorithms as well as with popular algorithmic
frameworks such as RED or Langevin algorithms. We give an overview of the
possibilities offered by the proposed approach on popular algorithms relying
on implicit denoising priors in Figure 1.

##  2 Related Work

#### Stable plug-and-play algorithms

A significant amount of work has recently been proposed to provide stable PnP
algorithms. A popular line of research in that direction consists in
constraining the denoiser’s architecture to ensure its stability when plugged
into a PnP scheme. Some approaches propose to regularize the denoiser’s
training loss with a term penalizing the Lipschitz constant of the denoiser
[34, 29], which can be combined with architectural constraints [40, 17]. Other
works propose to alter the optimization algorithm itself in order to ensure
the its stability [27, 10, 48, 18]. In turn, ensuring convergence of the
algorithm allows to ensure better transferability to new imaging tasks [41],
but also to perform iteration-intensive tasks, such as sampling from the
posterior distribution [23]. However, all the aforementionned methods come at
the cost of either algorithmic modifications or strong constraints on the
design of the denoiser, as opposed to the original, desnoiser agnostic
approaches [43, 32].

#### Equivariance in imaging inverse problems

Equivariance to certain transformations, such as rotations or translation, has
often been a desired property when designing handcrafted variational priors
[36, 13, 35]. Geometric ensembling techniques have been known in the computer
vision literature, where it has been shown to improve reconstruction quality
for image super-resolution [42]. Yet, to the notable exception of [48, 41],
this strategy went unnoticed in the PnP literature, and to the best of our
knowledge, the role of equivariance has not been explored in the context of
algorithms relying on implicit denoising priors. In a different line of work,
but still in the context of imaging inverse problems, recent works have
exploited equivariance for the design of unrolled network architectures [5],
or the construction of self-supervised learning losses [7] (e.g., see the
recent review [8]).

#### Equivariant neural networks

More generally, the equivariance assumption appears as a natural assumption in
a large number of application, and there exists a vast literature on building
equivariant neural network architectures [11, 45, 4]. Typical applications
involve segmentation on spherical manifolds, robotics, point cloud analysis,
data augmentation to name a few [15, 7, 28, 33, 44, 51]. However, these
networks often fail to perform as well as other state-of-the-art
architectures. Moreover, it can be challenging to incorporate complex layers
(e.g., upsampling/downsampling, attention-like layers, etc.) without breaking
the equivariance of the resulting network. In this work, we provide a simple
method for rendering any denoiser equivariant, without any architectural
constraints.

##  3 PnP algorithms

Traditional variational approaches for solving (1) consists in reformulating
it as a minimization problem. Following a maximum-a-posteriori approach, one
can derive an estimate $\widehat{x}$ as

$\widehat{x}=\underset{x}{\operatorname{argmin}}\,\,f(x)+\lambda r(x)$ (2)

where $f$ is a data-fidelity enforcing term, $r$ is a regularization enforcing
prior knowledge about the solution, and $\lambda>0$ is a regularization
parameter.

PnP approaches propose to replace the proximity operator of $r$ [2] (implicit
gradient step) arising in algorithms for solving (2) by a denoiser
$\operatorname{D}$ [43]. For the standard case of a quadratic data-fidelity
term $f(x)=\frac{1}{2}\|Ax-y\|_{2}^{2}$, the classical PnP algorithm reads

$x_{k+1}=\operatorname{D}\left(x_{k}-\gamma A^{\top}(Ax_{k}-y)\right)$ (PnP)

where $\gamma>0$ is a stepsize.

Similarly, regularization by denoising (RED) algorithms approximate the
gradient of $r$ as $\nabla r(x)\propto x-\operatorname{D}(x)$ using Tweedie’s
formula [32]. A simple explicit gradient descent optimization based on this
definition yields

$x_{k+1}=x_{k}-\gamma
A^{\top}(Ax_{k}-y)-\gamma\lambda(x_{k}-\operatorname{D}(x_{k})).$ (RED)

We also consider the Unadjusted Langevin algorithm (ULA) [23], which aims to
obtain samples associated with the negative log posterior density $-\log
p(x|y)\propto f(x)+\lambda r(x)$, and requires adding noise to the iterates in
RED, i.e.,

$x_{k+1}=x_{k}-\gamma
A^{\top}(Ax_{k}-y)-\gamma\lambda(x_{k}-\operatorname{D}(x_{k}))+\sqrt{2\gamma}\epsilon_{k},$
(ULA)

where $\epsilon_{k}\sim\mathcal{N}(0,I)$ is a standard Gaussian vector.

While there exist many different variants of PnP, RED, and ULA, we here focus
on their most standard formulations. Interestingly, these algorithms have
shown impressive performance on a wide variety of imaging tasks while relying
on Gaussian denoisers agnostic to the imaging modality of interest. Yet, these
algorithms suffer from a lack of stability and potential divergence effects,
hurting their versatility.

###  3.1 Proposed equivariant approach

Intuitively, imaging priors should have some invariance properties with
respect to certain groups of transformations, such as rotations, translations,
and reflections. We denote these transformations associated with a group
$\mathcal{G}$, $\\{T_{g}\\}_{g\in\mathcal{G}}$ where
$T_{g}\in\mathbb{R}^{n\times n}$ is a unitary matrix111While it is possible to
define group actions with non-unitary matrices [37], here we focus on the
unitary matrices, which is the case of translations, rotations and reflections
of images.. We say that $\operatorname{D}$ is equivariant to the group action
$\\{T_{g}\\}_{g\in\mathcal{G}}$ if
$\operatorname{D}(T_{g}x)=T_{g}\operatorname{D}(x)$ for all $x$ and
$g\in\mathcal{G}$. At the algorithmic level, this requirement translates into
the equivariance of the denoiser with respect to the transforms of interest:
if $r(x)$ is a $\mathcal{G}$-invariant function, its proximal operator and
gradient (if they exist) are necessarily $\mathcal{G}$-equivariant functions
[5].

A simple way of rendering any function $\mathcal{G}$-equivariant is by
averaging over the group222This construction is known as Reynolds averaging,
see e.g. [37].. The associated averaged denoiser

$\operatorname{D}_{\mathcal{G}}(x)\overset{\text{def}}{=}\frac{1}{|\mathcal{G}|}\sum_{g\in\mathcal{G}}T_{g}^{-1}\operatorname{D}(T_{g}x).$
(3)

is equivariant by construction. For large groups, or in the case of large
denoising architectures $\operatorname{D}$, computing the averaged denoiser
might be too computationally demanding. However, in this work, we propose to
use a simple Monte Carlo approximation by sampling a single transformation at
each step of the algorithm, i.e.,

$\displaystyle g$ $\displaystyle\sim\mathcal{G}$ (4)
$\displaystyle\widetilde{\operatorname{D}}_{\mathcal{G}}(x)$
$\displaystyle=T_{g}^{-1}\operatorname{D}(T_{g}x).$

The equivariant counterpart of the PnP, RED and ULA algorithms is thus simply
obtained by replacing the denoiser $\operatorname{D}$ by a sample of the Monte
Carlo estimate (4) at each iterate of the algorithm. Explicit versions of
these algorithms can be found in the Supplementary Material (SM).

##  4 Theoretical analysis

In this section, we provide a theoretical analysis of the advantages in terms
of performance and stability of equivariant denoisers, in comparison to their
non-equivariant counterparts. In this section, we denote the Jacobian of the
denoiser as $J_{x}\overset{\text{def}}{=}\frac{\delta\operatorname{D}}{\delta
x}(x)$.

#### Optimality of equivariant denoisers

We first show that if the signal distribution is $\mathcal{G}$-invariant, then
for any denoiser $\operatorname{D}$, its averaged version 3 obtains an equal
or better denoising performance. This can be shown by computing the expected
$\ell_{2}$ error with respect to the signal and noise distributions, i.e.,

$\displaystyle\mathbb{E}\,\Big{\|}\frac{1}{|\mathcal{G}|}\sum_{g\in\mathcal{G}}T_{g}^{-1}\operatorname{D}(T_{g}(x+\varepsilon))-x\Big{\|}$
$\displaystyle\leq\frac{1}{|\mathcal{G}|}\sum_{g\in\mathcal{G}}\mathbb{E}\,\|T_{g}^{-1}\operatorname{D}(T_{g}(x+\varepsilon)-x\|$
$\displaystyle\leq\frac{1}{|\mathcal{G}|}\sum_{g\in\mathcal{G}}\mathbb{E}\,\|T_{g}^{-1}\operatorname{D}(T_{g}T_{g}^{-1}(x+\varepsilon))-T_{g}^{-1}x\|$
$\displaystyle\leq\mathbb{E}\,\|\operatorname{D}(x+\varepsilon)-x\|$

where the expectation is taken with respect to $x$ and $\varepsilon$, and
where the second line uses the triangle inequality and the third line uses
that (i) $\mathbb{E}h(x)=\mathbb{E}h(T_{g}^{-1}x)$ for any function
$h:\mathbb{R}^{n}\mapsto\mathbb{R}$ and $T_{g}$ in the group action if the
distribution of $x$ and $\varepsilon$ is $\mathcal{G}$-invariant, and that
(ii) the transformations $T_{g}$ are isometries.

#### Existence of an explicit prior

A necessary condition for a denoiser to be associated with an explicit (PnP or
RED) prior is to have a symmetric Jacobian, i.e., $J_{x}=J_{x}^{\top}$, see
[30, Theorem 1]. Unfortunately, most state-of-the-art denoisers do not exhibit
this property. Averaging a denoiser over a sufficiently large group can lead
to symmetric Jacobians. In particular, if the denoiser is linear and
$\mathcal{G}$ includes translations and reflections, then the denoiser is
assured to have a symmetric Jacobian:

######  Proposition 1.

Any linear denoiser $\operatorname{D}$ that is equivariant to the action of
2-dimensional shifts, and vertical and horizontal reflections, has a symmetric
Jacobian.

###### Proof.

Let $\operatorname{D}_{\mathcal{G}}(x)=Mx$ with Jacobian
$M\in\mathbb{R}^{n\times n}$. A matrix $M$ that is equivariant to the action
of 2-dimensional has a circulant form, i.e., $M=\text{circ}(d)$ where
$d\in\mathbb{R}^{n}$ is a filter. Thus the transposed Jacobian is also a
circulant matrix $M^{\top}=\text{circ}(d^{\prime})$ where
$d^{\prime}\in\mathbb{R}^{n}$ is the transposed filter. Since $D$ is also
equivariant to vertical and horizontal reflections, we have that $d$ is even
and real, and thus $M=M^{\top}$. ∎

While this result applies only to linear denoisers, the symmetry of non-linear
denoisers also improves when incorporating equivariance: Table 1 shows the
relative symmetry error
$\|J_{x}-J_{x}^{\top}\|_{\mathrm{F}}^{2}/\|J_{x}\|_{\mathrm{F}}^{2}$ of
popular non-linear denoisers averaged over 10 different patches of $64\times
64$ pixels. The $\mathcal{G}$-equivariant denoisers
$\operatorname{D}_{\mathcal{G}}$ have significantly smaller errors than their
non-equivariant counterparts.

|  | Lipschitz
---
DnCNN
DnCNN | SCUNet | SwinIR |  | DRUNet
---
($\sigma_{d}$=0.01)
Standard $\operatorname{D}$ | 0.014 | 0.022 | 0.954 | 0.604 | 0.030
Equivariant $\operatorname{D}_{\mathcal{G}}$ | 0.003 | 0.005 | 0.710 | 0.291 | 0.008
Table 1: Mean Jacobian Symmetry error
$\|J_{x}-J_{x}^{\top}\|_{F}^{2}/\|J_{x}\|_{F}^{2}$. Equivariant denoisers are
obtained by averaging over the group of 90-degree rotations and reflections.

#### Lipschitz constant of the denoiser

The stability of PnP algorithms depends crucially on the Lipschitz constant of
the denoiser [34, 20]. For example, if the Lipschitz constant of the denoiser
is lower than 1, both the PnP and RED iterates converge under a good choice of
step size. Since the Lipschitz constant of the sum of two mappings is smaller
or equal than the sum of their Lipschitz constants, we have that the Lipschitz
constant of the averaged equivariant denoiser is necessarily equal or lower
than the non-equivariant one. If we restrict to linear denoisers, we can show
that the equivariant denoiser will have a strictly smaller constant, as long
as the dominant singular vector is not equivariant:

######  Proposition 2.

Let $\operatorname{D}(x)=Mx$ be a linear denoiser with singular value
decomposition $M=\sum_{i=1}^{n}\lambda_{i}u_{i}v_{i}^{\top}$ and
$\lambda_{1}>\lambda_{2}\geq\dots\geq\lambda_{n}\geq 0$. If the principal
component $u_{1}v_{1}^{\top}$ is not $\mathcal{G}$-equivariant, the averaged
denoiser $\operatorname{D}_{\mathcal{G}}$ has a strictly smaller Lipschitz
constant than $\operatorname{D}$.

###### Proof.

For any $p$ matrices $A_{1},\dots,A_{p}$, we have that
$\|\frac{1}{p}\sum_{g=1}^{p}A_{g}\|=\frac{1}{p}\sum_{g=1}^{p}\|A_{g}\|$ if and
only if all matrices share the same leading left and right singular vectors.
The $\mathcal{G}$-averaged denoiser can be written as
$\operatorname{D}_{\mathcal{G}}=\frac{1}{|\mathcal{G}|}\sum_{g=1}^{|\mathcal{G}|}A_{g}$
where $A_{g}:=T_{g}MT_{g}^{-1}$. We have that $A_{g}$ has the same singular
values as $M$ since singular vectors are defined as
$u_{i}^{\prime}=T_{g}u_{i}$ and $v^{\prime}=T_{g}v_{i}$ for $i=1,\dots,n$.
Since $u_{1}v_{1}^{\top}$ is not equivariant, we have that
$T_{g}u_{1}v_{1}^{\top}T_{g}^{-1}\neq u_{1}v_{1}$ for some $g\in\mathcal{G}$.
Thus, there exist at least 2 terms in the sum $\sum_{g}T_{g}MT_{g}^{-1}$ which
do not share the same leading singular vectors, and consequently
$\|M\|>\|\frac{1}{|\mathcal{G}|}\sum_{g}T_{g}MT_{g}^{-1}\|$. ∎

In practice, we observe a significantly smaller constant for most popular non-
linear denoisers. Table 2 shows the local Lipschitz constant (i.e., the
spectral norm of the Jacobian) of various denoisers averaged over 10 different
patches of $64\times 64$ pixels. The averaged denoiser can have a
significantly smaller constant.

|  | Lipschitz
---
DnCNN
DnCNN | SCUNet | SwinIR |  | DRUNet
---
($\sigma_{d}$=0.01)
Standard $\operatorname{D}$ | 1.06 | 1.44 | 5.78 | 6.28 | 1.57
Equivariant $\operatorname{D}_{\mathcal{G}}$ | 0.92 | 1.18 | 4.19 | 4.05 | 1.26
Table 2: Local Lipschitz constant of the denoiser averaged over 16 image
patches. Equivariant denoisers are obtained by averaging over the group of 90
degree rotations and reflections.

#### Interplay between the group action and forward operator

So far we have focused on the properties associated with an equivariant
denoiser, however, the (lack of) equivariance of $A$ also plays an important
role in the convergence of PnP algorithms. The iterates in PnP can converge
even for denoisers with Lipschitz constant larger than 1, as long as the
Lipschitz constant of the composition

$\operatorname{D}\circ\,(I-\gamma A^{\top}A)$ (5)

is smaller than 1. If the spectra of the matrices $A^{\top}A$ and
$\operatorname{D}$ are incoherent, i.e., $A^{\top}A$ and the Jacobian of the
denoiser are diagonalized in different bases, the Lipschitz constant of 5 is
likely to be smaller than that of $\operatorname{D}$. A similar stability
phenomenon happens for the RED and ULA iterates due to the incoherence between
the spectrum of the forward operator and the one of the equivariant denoiser.
We formalize this intuition in the following proposition:

######  Proposition 3.

Let $x$ be a grayscale image and $\\{T_{g}\\}_{g\in\mathcal{G}}$ be a group of
transformations that includes 2-dimensional shifts and
$\operatorname{D}_{\mathcal{G}}$ be a linear $\mathcal{G}$-equivariant
denoiser. If $A^{\top}A$ is not $\mathcal{G}$-equivariant, it does not share
the same singular vectors as the Jacobian of $\operatorname{D}$.

###### Proof.

Let $\operatorname{D}(x)=Mx$ with Jacobian $M\in\mathbb{R}^{n\times n}$. A
matrix $B\in\mathbb{R}^{n\times n}$ is $\mathcal{G}$-equivariant to the action
of 2-dimensional shifts if and only if it admits a diagonal decomposition as
$B=F\text{diag}(d)F^{*}$ where $F$ is the 2-dimensional Fourier transform, see
e.g. [39, Section 4.1]. Thus, if $\operatorname{D}$ is
$\mathcal{G}$-equivariant, then $M$ is diagonal in the Fourier domain, whereas
the non-equivariant $A^{\top}A$ matrix does not admit such a diagonalization.
∎

Table 3 shows the local Lipschitz constant of the mapping in 5 for various
popular denoisers, where $A$ is a random inpainting operator, which is not
equivariant to shifts, rotations nor reflections. The constants are smaller
than those shown in Table 2 and are below 1, ensuring contraction of the PnP
iterates for that specific operator. The equivariant denoisers show smaller
constants than the non-equivariant counterparts.

|  | Lipschitz
---
DnCNN
DnCNN | SCUNet | SwinIR |  | DRUNet
---
($\sigma_{d}$=0.01)
Standard | 0.91 | 0.91 | 0.62 | 0.69 | 0.83
Equivariant | 0.79 | 0.78 | 0.52 | 0.67 | 0.70
Table 3: Local Lipschitz constant of PnP iteration $\|J_{x}(I-A^{\top}A)\|$,
with $A$ being a random inpainting operator. Equivariant denoisers are
obtained by averaging over the group of 90 rotations and reflections.

###  4.1 Non-linear example

We demonstrate some of the properties analysed in the previous subsection on a
non-linear setting of a neural network denoiser with a single hidden layer.
More precisely, we consider the case where $\operatorname{D}$ is a slight
perturbation of a proximity operator, i.e.

$\operatorname{D}(x)=B_{2}\operatorname{prox}_{\gamma\lambda\|\cdot\|_{1}}(B_{1}x)$
(6)

where $B_{1}\in\mathbb{R}^{n\times n}$ satisfies
$B_{1}B_{1}^{\top}=B_{1}^{\top}B_{1}=I$ and where $B_{2}=(B_{1}^{\top}+P)$
where $P\in\mathbb{R}^{n\times n}$ is a small random perturbation. In
particular, if $P=0$, then the denoiser is a well-defined proximal operator,
i.e.
$\operatorname{D}(x)=\operatorname{prox}_{\gamma\lambda\|B_{1}\cdot\|_{1}}(x)$,
with a well-defined prior $r(x)=\|B_{1}x\|_{1}$. We stress that, for $P\neq
0$, there exists a priori no loss function $g$ associated with the PnP
iterates with the denoiser in (6). Furthermore, we assume that $B_{1}$ is
$\mathcal{G}$-equivariant.

Figure 2: Behaviour of the PnP algorithm with an approximated proximity
operator (blue curve) and its equivariant counterpart (red curve). Contour
lines show the loss in 2 with $r(x)=\|B_{1}x\|_{1}$. Stars denote the limit
point of each sequence (when it exists) and green dots show the initialization
points.

Using the fact that $B_{1}$ and
$\operatorname{prox}_{\gamma\lambda\|\cdot\|_{1}}$ are
$\mathcal{G}$-equivariant functions, we can write the
$\mathcal{G}$-equivariant denoiser as (see a detailed derivation in the SM)

$\displaystyle\operatorname{D}_{\mathcal{G}}(x)$
$\displaystyle=(B_{1}^{\top}+P_{\mathcal{G}})\operatorname{prox}_{\gamma\lambda\|\cdot\|_{1}}(B_{1}x)$

where
$P_{\mathcal{G}}=\frac{1}{|\mathcal{G}|}\sum_{g\in\mathcal{G}}T_{g}^{-1}PT_{g}$
is the $\mathcal{G}$-averaged perturbation. If the original perturbation is
not equivariant we have $0\leq\|P_{\mathcal{G}}\|_{F}^{2}<\|P\|_{F}^{2}$, and
the equivariant denoiser will be closer to the proximal operator
$\operatorname{prox}_{\gamma\lambda\|B_{1}\cdot\|_{1}}(x)$.

Figure 2 illustrates the behaviour of the PnP sequence with and without the
group averaging in 3 for two specific choices of $A$ in a 2D toy example,
where $\mathcal{G}$ is the group of flips (see the SM for more details). In
both cases, the algorithm involving the equivariant denoiser converges to a
point close to the global minima associated with the prior $\lambda
r(x)=\lambda\|B_{1}x\|_{1}$, whereas sequence generated by the non-equivariant
algorithm diverges.

##  5 Experimental results

In this section, we evaluate the influence of the proposed equivariant
approach for different algorithms and linear inverse imaging problems.

###  5.1 Problems considered

#### Image deblurring and image super-resolution

In this setting, we set $y=h*x+\epsilon$ in (1) where $h$ is a convolutional
kernel and $*$ the circular convolution. We consider either Gaussian
deblurring, in which case $h$ is a Gaussian kernel of standard deviation 1,
and motion deblurring, in which case $h$ is the first kernel from [24]. In
both cases, $\epsilon$ is a Gaussian noise with standard deviation $0.01$. In
the image super-resolution (SR) setting, (1) writes as
$y=\left(h*x\right)_{\Downarrow S}+\epsilon$ where $h$ is a Gaussian kernel of
standard deviation 1, and $S$ denotes the undersampling factor. When $S=2$
(resp. $S=4$), $\epsilon$ is a Gaussian noise with standard deviation $0.01$
(resp. $0.05$). We test the proposed method on the Set3C dataset as well as
the BSD10 dataset, a subset of 10 images from the BSD68 dataset.

#### MRI

In this setting, we consider $y=MFx$ in (1) where $M$ is a binary mask and $F$
the 2D Fourier transform. Following [46], we consider the $\times$4 and
$\times$8 acceleration factors. In contrast to the previous problem, no noise
is added to the measurements in this case. We test the method on a subset of
10 images from the validation set of the fully acquired k-space data of [46].

#### Algorithms and backbone denoisers

We consider several backbone pretrained denoiser, namely DRUNet [48], SCUNet
[49], SwinIR [25], DiffUNet [16] and DnCNN [47] as well as its Lipschitz
version (LipDnCNN) from [29]. These architectures are representative of state-
of-the-art image reconstruction architectures, involving both convolutional
and attention layers. We stress that none of these networks is equivariant to
translations or rotations, to the exception of DnCNN that show approximate
translation equivariance. The choice of the backbone denoiser may influence
the chosen algorithm. For instance, the DnCNN and SwinIR denoisers on which
this article relies are trained for a fixed level of noise, limiting the
ability to finetune the algorithm. Moreover, DiffUNet and SCUNet architectures
can only be applied to color images, preventing their usage on the MRI
problem. Unless mentioned otherwise, each algorithm is run for $10^{4}$
iterations.

#### Baselines

In this paper, we mainly investigate the influence of equivariant priors on
the stability of PnP algorithms. We thus compare the proposed approach with
variational approaches which can be seen as a class of convergent PnP
algorithms. In particular, we use wavelet denoisers (i.e.
$\operatorname{D}(x)=\operatorname{prox}_{\lambda\|\Psi\cdot\|_{1}}(x)$ for
$\Psi$ a redundant wavelet dictionary), as well as total generalized variation
(TGV) denoisers [3]. The PnP algorithm with Lipschitz denoiser corresponds to
the method from [29] which ensures convergence of the PnP algorithm. We also
compare our approach with the state-of-the-art DPIR algorithm [48] which runs
a small number of steps of a half-quadratic splitting algorithm with fine-
tuned decaying stepsizes. As a consequence, DPIR can be seen as representative
of non-convergent, fine-tuned PnP algorithm.

| Deblur (motion) | Deblur (Gaussian) | MRI
---|---|---|---
| Set3C | BSD10 | Set3C | BSD10 |  $\times$4 |  $\times$8
LipDnCNN [29] | $31.7\pm 2.0$ | $30.9\pm 0.7$ | $32.1\pm 2.5$ | $32.6\pm 0.9$ | $30.6\pm 2.6$ | $26.4\pm 2.3$
Eq. LipDnCNN | $31.8\pm 2.0$ | $31.1\pm 0.6$ | $32.1\pm 2.5$ | $32.7\pm 0.9$ | $30.7\pm 2.7$ | $26.6\pm 2.2$
DnCNN | $29.9\pm 1.9$ | $30.4\pm 0.3$ | $22.0\pm 6.0$ | $29.6\pm 4.6$ | $28.6\pm 4.3$ | div.
Eq. DnCNN | $30.7\pm 1.6$ | $30.9\pm 0.1$ | $29.8\pm 3.1$ | $33.0\pm 0.6$ | $30.1\pm 3.9$ | $24.6\pm 2.9$
DRUNet | $10.1\pm 1.5$ | $16.5\pm 9.8$ | $14.6\pm 8.8$ | $17.5\pm 10.5$ | $27.7\pm 4.0$ | $16.9\pm 8.6$
Eq. DRUNet | $18.0\pm 6.0$ | $28.4\pm 2.3$ | $25.3\pm 10.1$ | $31.9\pm 1.1$ | $30.7\pm 3.0$ | $22.0\pm 6.0$
wavelets | $29.6\pm 2.0$ | $27.2\pm 1.0$ | $31.2\pm 2.5$ | $30.9\pm 0.8$ | $28.6\pm 2.1$ | $25.4\pm 1.8$
TGV | $29.2\pm 2.1$ | $26.7\pm 1.3$ | $30.7\pm 2.5$ | $30.5\pm 0.9$ | $28.5\pm 2.3$ | $24.8\pm 2.1$
DPIR [48] | $33.9\pm 1.9$ | $35.0\pm 0.5$ | $33.0\pm 2.5$ | $34.6\pm 0.7$ | $28.4\pm 2.2$ | $25.1\pm 1.9$
Table 4: Mean reconstruction PSNR on different problems for the PnP algorithm with different backbone denoisers. Three top rows show results with non-equivariant denoisers, while the three middle rows show the proposed equivariant modification. The last three rows provide standard reconstruction baselines. “div.” indicated divergence of the method. Groundtruth | $A^{\top}y$ | Standard | Equivariant
---|---|---|---
|  |  |
PSNR | 16.65 | 30.30 | 31.65
Figure 3: Motion deblurring on a Set3C sample with (PnP) relying on a DRUNet
backbone denoiser.

###  5.2 Stability of PnP and RED algorithms

The instability of PnP algorithms often translates into unrealistic artifacts
in the reconstructed image. Figure 3 illustrates this phenomenon and shows
that the proposed approach allows us to circumvent this drawback.

DnCNN |
---|---
DRUNet |
LipDnCNN |

---

Figure 4: Average PSNR (left) and convergence criterion
$\|x_{k+1}-x_{k}\|/\|x_{k}\|$ (right) for 3 different imaging problems as a
function of the iteration number with different backbone denoisers plugged in
the PnP algorithm. Top row: DnCNN, middle row: DRUNet, bottom row: Lipschitz
DnCNN.

Figure 4 shows the PSNR and convergence criterion
$\|x_{k+1}-x_{k}\|/\|x_{k}\|$ on different problems and backbone architecture.
For each of these problems, the denoising level $\sigma$ in (PnP) is set to
$\sigma=0.01$ (resp. $\sigma=0.015$) for DnCNN (resp. DRUNet) backbone
architectures. We notice that the (PnP) algorithm with equivariant denoiser
shows a more stable PSNR along iterations than its non-equivariant
counterpart. In particular, we observe in the right panel of Figure 4 that the
convergence rate of the (PnP) algorithm with $\mathcal{G}$-equivariant DnCNN
matches the behaviour of the algorithm with Lipschitz denoisers. Lastly, we
notice that the proposed equivariant approach also benefits convergent (PnP)
algorithms relying on Lipschitz backbone denoisers. Reconstruction metrics are
provided in Table 4.

| Groundtruth | $A^{\top}y$ | TGV (29.31) | Wavelets (29.99)
---|---|---|---|---
|  |  |  |
| DnCNN (22.34) | SCUNet (12.56) | DRUNet (8.93) | SwinIR (27.83)
Standard |  |  |  |
| DnCNN (30.52) | SCUNet (13.89) | DRUNet (29.44) | SwinIR (27.82)
Equiv. |  |  |  |
Figure 5: Gaussian deblurring on a BSD10 sample (detail) for different
denoising backbone plugged in the PnP algorithm. Middle row: standard
algorithms; bottom row: equivariant algorithms.

In the case of the (RED) algorithm, we observe similar behaviours.
Reconstructions on a $\times$2 SR problem are shown in Figure 6 with
associated convergence plots. In the case of a non-equivariant DRUNet
backbone, the reconstructed image shows important geometric artifacts that
disappear in the equivariant case.

We however stress that the proposed approach may fail to stabilize the
algorithm in certain settings. For example, in the case of a SCUNet backbone
denoiser, strong artifacts are visible from 50 iterations only in the
reconstruction in both the classical and equivariant version of the (PnP)
algorithm as seen in Figure 5. Similarly, we did not observe convergence of
the RED algorithm with the DiffUNet backbone, thus requiring early stopping of
the algorithm in order to reach good reconstruction results.

Lastly, several studies have shown that the noise level $\sigma$ plays an
important role in the stability of the PnP algorithm. Increasing $\sigma$ may
help to solve the instability issue at the cost of over-smoothed
reconstructions.

| DiffUNet | DRUNet
---|---|---
|  |
Observed | Standard (21.55) | Standard (20.75)
|  |
Groundtruth (PSNR) | Equivariant (22.26) | Equivariant (24.94)

DiffUNet | DRUNet
---|---
|

Figure 6: Left: Results of the RED algorithm on a $\times 2$ SR problem for
different backbone denoisers. Middle column: DiffUNet; right column: DRUNet.
Top row: standard algorithm; Bottom row: equivariant algorithm. Right:
Associated convergence criterion.

###  5.3 Interplay with the kernel of $A$

In the previous section, we have seen that equivariant denoisers can prevent
the emergence of artifacts along the iteration of (PnP). Despite their
unnatural aspect, these artifacts are not incompatible with good data-fidelity
measures. In fact, PnP algorithm offers no control over $\ker(A)$, which is
nontrivial by nature of the ill-posed inverse problem (1); artifacts appearing
during the reconstruction are therefore likely to belong to $\ker(A)$.

|  |  |
---|---|---|---
| Groundtruth | $A^{\top}y$ | Mask
Standard |  |  |
|  $x_{j}$, PSNR=34.10 | $|x_{j}-x_{i}|$ | $\operatorname{FFT}(x_{j}-x_{i})$
Equivariant |  |  |
|  $x_{j}$, PSNR=35.26 | $|x_{j}-x_{i}|$ | $\operatorname{FFT}(x_{j}-x_{i})$
Figure 7: Evolution of the reconstruction along the iterates of PnP with a
(non 1-Lipschitz) DnCNN backbone. The top row shows the ground truth, back-
projected data, and Fourier mask. The second (resp. third) row shows, from
left to right: reconstruction at iteration $j=10^{4}$; difference
$|x_{j}-x_{i}|$ between reconstructions at iteration $j=10^{4}$ and at
iteration $i=10^{3}$; the difference in the Fourier domain (displayed in
logarithmic scale).

This phenomenon can be illustrated in the case of MRI, where $\ker(A)$
corresponds to the non-sampled subspace of the Fourier domain. This is
illustrated in Figure 7 showing the reconstruction with and without the
proposed equivariant algorithmic update. In the standard (non-equivariant)
setting, mild artifacts appear between iteration $j=10^{4}$ and $j=10^{3}$.
These may appear unnoticed in $x_{j}$; however, they clearly appear when
plotting the difference $x_{j}-x_{i}$. Interestingly, the Fourier spectrum of
these artifacts shows significantly more energy at frequencies that were not
sampled. The same experiment with the equivariant algorithm shows no such
artifacts and a more uniform Fourier spectrum.

###  5.4 Sampling from RED prior

The gain in stability provided by the proposed approach opens the door to more
robust sampling algorithms relying on implicit denoising priors, such as ULA,
where a sufficiently large number of iterations is required in order to obtain
good estimators. We show in Figure 8 the estimated mean and variance obtained
with the ULA algorithm, with both equivariant and non-equivariant DRUNet
backbone denoisers. In the non-equivariant setting, we observe similar
artifacts to those obtained in the deterministic case on both the estimated
mean and variance; these artifacts vanish in the equivariant case.

The plots of Figure 8 further illustrates this phenomenon. We notice that
after a few hundred iterations, the reconstruction quality with the non-
equivariant ULA algorithm collapses, yielding an irrelevant MCMC chain. While
enforcing equivariance improves the situation, increasing the noise level in
the denoiser in PnP can significantly improve the situation, as shown on the
rightmost plot of Figure 8.

| Standard | Equivariant
---|---|---
|  |
Observed | Estimated mean | Estimated mean
|  |
Ground truth | Estimated variance | Estimated variance

Figure 8: Sampling with the ULA algorithm on a motion deblurring and the
DRUNet prior, in the standard and equivariant cases. Left: PSNR of a sample
from the MCMC chain for a sample of the Set3C dataset. Right: average PSNR
over the Set3C dataset as a function of $\sigma$ in the denoising prior.

##  6 Limitations

While our experiments show that the proposed method enables to always improve
both the stability of the algorithmic scheme and the reconstructed image
quality over its non-equivariant counterpart, we stress that it remains prone
to divergence and/or hallucinating artifacts. This can for example be seen in
Figure 4 in the case of the DRUNet backbone for example, which shows
decreasing mean PSNR after a certain number of iterations, or in Figure 5 for
the SCUNet backbone, which shows important reconstruction artifacts. In
particular, the proposed equivariant approach when using a SCUNet backbone
architecture does not clearly improve over its non-equivariant counterpart.

##  7 Conclusion

In this work, we have proposed a simple yet efficient method for ensuring
approximate equivariance of implicit denoising priors. In essence, the method
amounts to sampling and applying a group action at random at each step of the
algorithm. In spite of its simplicity, this method shows interesting
theoretical properties. In the case of linear denoisers, for example, it
allows us to enforce the symmetry of the Jacobian, which is a cornerstone
property when aiming to derive explicit priors from implicit denoising priors.
Furthermore, it can be shown that the Lipschitz constant of the equivariant
linear denoiser can only be lower than that of its non-equivariant
counterpart, thus improving the stability of the resulting PnP algorithm. We
showcase the symmetrization effect of equivariance on the Jacobian of several
architectures, and illustrate its stabilization effect for three families of
imaging algorithms, namely PnP, RED, and ULA. Importantly, this stabilization
procedure is not detrimental to the reconstruction quality, as often observed
in the literature. However, despite the significant improvements brought by
equivariance both in terms of stability and image reconstruction quality, the
proposed method remains prone to divergence and artifact contamination in the
reconstructions.

## References

  * Ahmad et al. [2020] Rizwan Ahmad, Charles A Bouman, Gregery T Buzzard, Stanley Chan, Sizhuo Liu, Edward T Reehorst, and Philip Schniter.  Plug-and-play methods for magnetic resonance imaging: Using denoisers for image recovery.  _IEEE signal processing magazine_ , 37(1):105–116, 2020.
  * Bauschke et al. [2017] Heinz H Bauschke, Patrick L Combettes, Heinz H Bauschke, and Patrick L Combettes.  _Correction to: Convex Analysis and Monotone Operator Theory in Hilbert Spaces_.  Springer, 2017.
  * Bredies et al. [2010] Kristian Bredies, Karl Kunisch, and Thomas Pock.  Total generalized variation.  _SIAM Journal on Imaging Sciences_ , 3(3):492–526, 2010.
  * Bronstein et al. [2021] Michael M Bronstein, Joan Bruna, Taco Cohen, and Petar Veličković.  Geometric deep learning: Grids, groups, graphs, geodesics, and gauges.  _arXiv preprint arXiv:2104.13478_ , 2021.
  * Celledoni et al. [2021] 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, 2021.
  * Chambolle and Pock [2016] Antonin Chambolle and Thomas Pock.  An introduction to continuous optimization for imaging.  _Acta Numerica_ , 25:161–319, 2016.
  * Chen et al. [2021] Dongdong Chen, Julián Tachella, and Mike E Davies.  Equivariant imaging: Learning beyond the range space.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 4379–4388, 2021.
  * Chen et al. [2023] Dongdong Chen, Mike Davies, Matthias J Ehrhardt, Carola-Bibiane Schönlieb, Ferdia Sherry, and Julián Tachella.  Imaging with equivariant deep learning: From unrolled network design to fully unsupervised learning.  _IEEE Signal Processing Magazine_ , 40(1):134–147, 2023.
  * Chung et al. [2022] Hyungjin Chung, Jeongsol Kim, Michael T Mccann, Marc L Klasky, and Jong Chul Ye.  Diffusion posterior sampling for general noisy inverse problems.  _arXiv preprint arXiv:2209.14687_ , 2022.
  * Cohen et al. [2021] Regev Cohen, Michael Elad, and Peyman Milanfar.  Regularization by denoising via fixed-point projection (red-pro).  _SIAM Journal on Imaging Sciences_ , 14(3):1374–1406, 2021.
  * Cohen and Welling [2016] Taco Cohen and Max Welling.  Group equivariant convolutional networks.  In _International conference on machine learning_ , pages 2990–2999. PMLR, 2016.
  * Combettes et al. [2014] Patrick L Combettes, Laurent Condat, J-C Pesquet, and BC Vũ.  A forward-backward view of some primal-dual optimization methods in image recovery.  In _2014 IEEE International Conference on Image Processing (ICIP)_ , pages 4141–4145. IEEE, 2014.
  * Condat [2017] Laurent Condat.  Discrete total variation: New definition and minimization.  _SIAM Journal on Imaging Sciences_ , 10(3):1258–1290, 2017.
  * De Bortoli et al. [2021] Valentin De Bortoli, James Thornton, Jeremy Heng, and Arnaud Doucet.  Diffusion schrödinger bridge with applications to score-based generative modeling.  _Advances in Neural Information Processing Systems_ , 34:17695–17709, 2021.
  * Fuchs et al. [2020] Fabian Fuchs, Daniel Worrall, Volker Fischer, and Max Welling.  Se (3)-transformers: 3d roto-translation equivariant attention networks.  _Advances in neural information processing systems_ , 33:1970–1981, 2020.
  * Ho et al. [2020] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  _Advances in neural information processing systems_ , 33:6840–6851, 2020.
  * Hurault et al. [2021] Samuel Hurault, Arthur Leclaire, and Nicolas Papadakis.  Gradient step denoiser for convergent plug-and-play.  _arXiv preprint arXiv:2110.03220_ , 2021.
  * Hurault et al. [2023] Samuel Hurault, Antonin Chambolle, Arthur Leclaire, and Nicolas Papadakis.  Convergent plug-and-play with proximal denoiser and unconstrained regularization parameter.  _arXiv preprint arXiv:2311.01216_ , 2023.
  * Jin et al. [2017] Kyong Hwan Jin, Michael T McCann, Emmanuel Froustey, and Michael Unser.  Deep convolutional neural network for inverse problems in imaging.  _IEEE transactions on image processing_ , 26(9):4509–4522, 2017.
  * Kamilov et al. [2023] Ulugbek S Kamilov, Charles A Bouman, Gregery T Buzzard, and Brendt Wohlberg.  Plug-and-play methods for integrating physical and learned models in computational imaging: Theory, algorithms, and applications.  _IEEE Signal Processing Magazine_ , 40(1):85–97, 2023.
  * Kawar et al. [2022] Bahjat Kawar, Michael Elad, Stefano Ermon, and Jiaming Song.  Denoising diffusion restoration models.  _Advances in Neural Information Processing Systems_ , 35:23593–23606, 2022.
  * Knoll et al. [2020] Florian Knoll, Jure Zbontar, Anuroop Sriram, Matthew J Muckley, Mary Bruno, Aaron Defazio, Marc Parente, Krzysztof J Geras, Joe Katsnelson, Hersh Chandarana, et al.  fastmri: A publicly available raw k-space and dicom dataset of knee images for accelerated mr image reconstruction using machine learning.  _Radiology: Artificial Intelligence_ , 2(1):e190007, 2020.
  * Laumont et al. [2022] Rémi Laumont, Valentin De Bortoli, Andrés Almansa, Julie Delon, Alain Durmus, and Marcelo Pereyra.  Bayesian imaging using plug & play priors: when langevin meets tweedie.  _SIAM Journal on Imaging Sciences_ , 15(2):701–737, 2022.
  * Levin et al. [2009] Anat Levin, Yair Weiss, Fredo Durand, and William T Freeman.  Understanding and evaluating blind deconvolution algorithms.  In _2009 IEEE conference on computer vision and pattern recognition_ , pages 1964–1971. IEEE, 2009.
  * Liang et al. [2021] Jingyun Liang, Jiezhang Cao, Guolei Sun, Kai Zhang, Luc Van Gool, and Radu Timofte.  Swinir: Image restoration using swin transformer.  In _Proceedings of the IEEE/CVF international conference on computer vision_ , pages 1833–1844, 2021.
  * Liu et al. [2022] Renhao Liu, Yu Sun, Jiabei Zhu, Lei Tian, and Ulugbek S Kamilov.  Recovery of continuous 3d refractive index maps from discrete intensity-only measurements using neural fields.  _Nature Machine Intelligence_ , 4(9):781–791, 2022.
  * Meinhardt et al. [2017] Tim Meinhardt, Michael Moller, Caner Hazirbas, and Daniel Cremers.  Learning proximal operators: Using denoising networks for regularizing inverse imaging problems.  In _Proceedings of the IEEE International Conference on Computer Vision_ , pages 1781–1790, 2017.
  * Müller et al. [2021] Philip Müller, Vladimir Golkov, Valentina Tomassini, and Daniel Cremers.  Rotation-equivariant deep learning for diffusion mri.  _arXiv preprint arXiv:2102.06942_ , 2021.
  * Pesquet et al. [2021] Jean-Christophe Pesquet, Audrey Repetti, Matthieu Terris, and Yves Wiaux.  Learning maximally monotone operators for image recovery.  _SIAM Journal on Imaging Sciences_ , 14(3):1206–1237, 2021.
  * Reehorst and Schniter [2018] Edward T Reehorst and Philip Schniter.  Regularization by denoising: Clarifications and new interpretations.  _IEEE transactions on computational imaging_ , 5(1):52–67, 2018.
  * Rick Chang et al. [2017] JH Rick Chang, Chun-Liang Li, Barnabas Poczos, BVK Vijaya Kumar, and Aswin C Sankaranarayanan.  One network to solve them all–solving linear inverse problems using deep projection models.  In _Proceedings of the IEEE International Conference on Computer Vision_ , pages 5888–5897, 2017.
  * Romano et al. [2017] Yaniv Romano, Michael Elad, and Peyman Milanfar.  The little engine that could: Regularization by denoising (red).  _SIAM Journal on Imaging Sciences_ , 10(4):1804–1844, 2017.
  * Rommel et al. [2022] Cédric Rommel, Thomas Moreau, and Alexandre Gramfort.  Deep invariant networks with differentiable augmentation layers.  _Advances in Neural Information Processing Systems_ , 35:35672–35683, 2022.
  * Ryu et al. [2019] Ernest Ryu, Jialin Liu, Sicheng Wang, Xiaohan Chen, Zhangyang Wang, and Wotao Yin.  Plug-and-play methods provably converge with properly trained denoisers.  In _International Conference on Machine Learning_ , pages 5546–5557. PMLR, 2019.
  * Saydjari and Finkbeiner [2022] Andrew K Saydjari and Douglas P Finkbeiner.  Equivariant wavelets: Fast rotation and translation invariant wavelet scattering transforms.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 45(2):1716–1731, 2022.
  * Selesnick et al. [2005] Ivan W Selesnick, Richard G Baraniuk, and Nick C Kingsbury.  The dual-tree complex wavelet transform.  _IEEE signal processing magazine_ , 22(6):123–151, 2005.
  * Serre et al. [1977] Jean-Pierre Serre et al.  _Linear representations of finite groups_.  Springer, 1977.
  * Stergiopoulou et al. [2023] Vasiliki Stergiopoulou, Subhadip Mukherjee, Luca Calatroni, and Laure Blanc-Féraud.  Fluctuation-based deconvolution in fluorescence microscopy using plug-and-play denoisers.  In _International Conference on Scale Space and Variational Methods in Computer Vision_ , pages 498–510. Springer, 2023.
  * Tachella et al. [2023] 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.
  * Terris et al. [2020] Matthieu Terris, Audrey Repetti, Jean-Christophe Pesquet, and Yves Wiaux.  Building firmly nonexpansive convolutional neural networks.  In _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 8658–8662. IEEE, 2020.
  * Terris et al. [2023] Matthieu Terris, Arwa Dabbech, Chao Tang, and Yves Wiaux.  Image reconstruction algorithms in radio interferometry: From handcrafted to learned regularization denoisers.  _Monthly Notices of the Royal Astronomical Society_ , 518(1):604–622, 2023.
  * Timofte et al. [2016] Radu Timofte, Rasmus Rothe, and Luc Van Gool.  Seven ways to improve example-based single image super resolution.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 1865–1873, 2016.
  * Venkatakrishnan et al. [2013] Singanallur V Venkatakrishnan, Charles A Bouman, and Brendt Wohlberg.  Plug-and-play priors for model based reconstruction.  In _2013 IEEE global conference on signal and information processing_ , pages 945–948. IEEE, 2013.
  * Wang et al. [2022] Dian Wang, Robin Walters, Xupeng Zhu, and Robert Platt.  Equivariant $q$ learning in spatial action spaces.  In _Conference on Robot Learning_ , pages 1713–1723. PMLR, 2022.
  * Weiler and Cesa [2019] Maurice Weiler and Gabriele Cesa.  General e (2)-equivariant steerable cnns.  _Advances in neural information processing systems_ , 32, 2019.
  * Zbontar et al. [2018] Jure Zbontar, Florian Knoll, Anuroop Sriram, Tullie Murrell, Zhengnan Huang, Matthew J Muckley, Aaron Defazio, Ruben Stern, Patricia Johnson, Mary Bruno, et al.  fastmri: An open dataset and benchmarks for accelerated mri.  _arXiv preprint arXiv:1811.08839_ , 2018.
  * Zhang et al. [2017] Kai Zhang, Wangmeng Zuo, Yunjin Chen, Deyu Meng, and Lei Zhang.  Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising.  _IEEE transactions on image processing_ , 26(7):3142–3155, 2017.
  * Zhang et al. [2021] Kai Zhang, Yawei Li, Wangmeng Zuo, Lei Zhang, Luc Van Gool, and Radu Timofte.  Plug-and-play image restoration with deep denoiser prior.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 44(10):6360–6376, 2021.
  * Zhang et al. [2022] Kai Zhang, Yawei Li, Jingyun Liang, Jiezhang Cao, Yulun Zhang, Hao Tang, Radu Timofte, and Luc Van Gool.  Practical blind denoising via swin-conv-unet and data synthesis.  _arXiv preprint arXiv:2203.13278_ , 2022.
  * Zhang et al. [2023] Yulun Zhang, Kai Zhang, Zheng Chen, Yawei Li, Radu Timofte, Junpei Zhang, Kexin Zhang, Rui Peng, Yanbiao Ma, Licheng Jia, et al.  Ntire 2023 challenge on image super-resolution (x4): Methods and results.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 1864–1883, 2023.
  * Zhu et al. [2023a] Minghan Zhu, Shizhong Han, Hong Cai, Shubhankar Borse, Maani Ghaffari, and Fatih Porikli.  4d panoptic segmentation as invariant and equivariant field prediction.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 22488–22498, 2023a.
  * Zhu et al. [2023b] Yuanzhi Zhu, Kai Zhang, Jingyun Liang, Jiezhang Cao, Bihan Wen, Radu Timofte, and Luc Van Gool.  Denoising diffusion models for plug-and-play image restoration.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 1219–1229, 2023b.

##  Appendix A Details on the non-linear example

#### Derivation of $D_{\mathcal{G}}$

We have

$\displaystyle D_{\mathcal{G}}(x)$ $\displaystyle=\frac{1}{|\mathcal{G}|}\sum
T_{g}^{-1}B_{2}\operatorname{prox}_{\gamma\lambda\|\cdot\|_{1}}(B_{1}T_{g}x)$
(7) $\displaystyle=\frac{1}{|\mathcal{G}|}\sum
T_{g}^{-1}(B_{1}+P)T_{g}\operatorname{prox}_{\gamma\lambda\|\cdot\|_{1}}(B_{1}x)$
$\displaystyle=\left(B_{1}+\frac{1}{|\mathcal{G}|}\sum
T_{g}^{-1}PT_{g}\right)\operatorname{prox}_{\gamma\lambda\|\cdot\|_{1}}(B_{1}x)$

yielding the desired result. The second step uses the fact that
$B_{2}=B_{1}+P$ and that $B_{1}$ and $\operatorname{prox}$ are
$\mathcal{G}$-equivariant functions. The third step just uses that $B_{1}$ is
a $\mathcal{G}$-equivariant function.

#### Numerical details for Figure 2

For both the leftmost and rightmost examples, we consider the group
$\mathcal{G}$ consisting of permuations of the coordinates of the vectors.
This is a group with a single element $g$, the matrix representation of its
linear application being the unitary matrix

$T_{g}=\begin{pmatrix}0&1\\\ 1&0\end{pmatrix}.$ (8)

In the leftmost example, we use $A=\operatorname{diag}(2,1)$, $B_{1}=I$
($B_{1}$ is thus $\mathcal{G}$-equivariant) and $\lambda=10$. The perturbation
and its $\mathcal{G}$-average are

$P=\begin{pmatrix}-0.228&-0.023\\\ 0.066&0.1\end{pmatrix}\quad
P_{\mathcal{G}}=\begin{pmatrix}-0.064&0.022\\\ 0.022&-0.064\end{pmatrix},$

with associated norms $\|P\|_{F}=0.26$, $\|P_{\mathcal{G}}\|_{F}=0.10$. The
PnP algorithm is ran with $\gamma=5e-2$.

In the rightmost example, we use $A=\operatorname{diag}(2,5e-4)$, $B_{1}=I$
and $\lambda=2$. The perturbation and its $\mathcal{G}$-average are

$P=\begin{pmatrix}0.0275&0.0244\\\
0.0112&-0.1842\end{pmatrix},P_{\mathcal{G}}=\begin{pmatrix}-0.0783&0.0178\\\
0.0178&-0.0783\end{pmatrix},$

with associated norms $\|P\|_{F}=0.0469$, $\|P_{\mathcal{G}}\|_{F}=0.0366$.
The PnP algorithm is ran with $\gamma=0.2$.

##  Appendix B Equivariant algorithms

The equivariant counterpart of (PnP) is

$\displaystyle\text{Sample }\,g_{k}\sim\mathcal{G}$ (eq. PnP)
$\displaystyle\text{Set
}\,\widetilde{\operatorname{D}}_{\mathcal{G},k}(x)=T_{g_{k}}^{-1}\operatorname{D}(T_{g_{k}}x)$
$\displaystyle
x_{k+1}=\widetilde{\operatorname{D}}_{\mathcal{G},k}\left(x_{k}-\gamma
A^{\top}(Ax_{k}-y)\right).$

The equivariant counterpart of (RED) is

$\displaystyle\text{Sample }\,g_{k}\sim\mathcal{G}$ (eq. RED)
$\displaystyle\text{Set
}\,\widetilde{\operatorname{D}}_{\mathcal{G},k}(x)=T_{g_{k}}^{-1}\operatorname{D}(T_{g_{k}}x)$
$\displaystyle x_{k+1}=x_{k}-\gamma A^{\top}(Ax_{k}-y)$
$\displaystyle\qquad\qquad\qquad\qquad-\gamma\lambda(x_{k}-\widetilde{\operatorname{D}}_{\mathcal{G},k}(x_{k})).$

The equivariant counterpart of (ULA) is

$\displaystyle\text{Sample }\,g_{k}\sim\mathcal{G}$ (eq. ULA)
$\displaystyle\text{Set
}\,\widetilde{\operatorname{D}}_{\mathcal{G},k}(x)=T_{g_{k}}^{-1}\operatorname{D}(T_{g_{k}}x)$
$\displaystyle x_{k+1}=x_{k}-\gamma A^{\top}(Ax_{k}-y)$
$\displaystyle\qquad\qquad-\gamma\lambda(x_{k}-\widetilde{\operatorname{D}}_{\mathcal{G},k}(x_{k}))+\sqrt{2\gamma}\epsilon_{k}.$