{'fraction_non_alphanumeric': 0.040708651234967023, 'fraction_numerical': 0.010474589421957843, 'mean_word_length': 4.6372649074209065, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Deep Equilibrium Architectures for Inverse Problems in Imaging

Recent efforts on solving inverse problems in imaging via deep neural networks use architectures inspired by a fixed number of iterations of an optimization method. The number of iterations is typically quite small due to difficulties in training networks corresponding to more iterations; the resulting solvers cannot be run for more iterations at test time without incurring significant errors. This paper describes an alternative approach corresponding to an infinite number of iterations, yielding a consistent improvement in reconstruction accuracy above state-of-the-art alternatives and where the computational budget can be selected at test time to optimize context-dependent trade-offs between accuracy and computation. The proposed approach leverages ideas from Deep Equilibrium Models, where the fixed-point iteration is constructed to incorporate a known forward model and insights from classical optimization-based reconstruction methods.


Introduction
A collection of recent efforts surveyed in [1] consider the problem of using training data to solve inverse problems in imaging. Specifically, imagine we observe a corrupted set of measurements y of an image x under a linear measurement operator A with some noise ε according to Our task is to compute an estimate of x given measurements y and knowledge of A. This task is particularly challenging when the inverse problem is ill-posed, i.e., when the system is underdetermined or ill-conditioned, in which case simple methods such as least squares estimation may not have a unique solution or may produce estimates that are highly sensitive to noise. Decades of research has explored geometric models of image structure that can be used to regularize solutions to this inverse problem, including [2,3,4] and many others. More recent efforts have focused instead on using large collections of training images, {x * i } n i=1 , to learn effective regularizers. One particularly popular and effective approach involves augmenting standard iterative inverse problem solvers with learned deep networks. This approach, which we refer to as deep unrolling (DU), is reviewed in §2.1. The basic idea is to build an architecture that mimics a small number of iterations of an iterative algorithm. In practice, the number of iterations is quite small (typically 5-10) because of issues stability, memory, and numerical issues arising in backpropagation. This paper sidesteps this key limitation of deep unrolling methods with a novel approach based on deep equilibrium models (DEMs) [5], which are designed for training arbitrarily deep networks. The result is a novel approach to training networks to solve inverse problems in imaging that yields a consistent improvement in performance above state-ofthe-art alternatives and where the computational budget can be selected at test time to optimize context-dependent tradeoffs between accuracy and computation. The key empirical findings, which are detailed in §6.4, are illustrated in Fig. 2(a).

Contributions
This paper presents a novel approach to machine learning-based methods for solving linear inverse problems in imaging. Unlike most state-of-the-art methods, which are based on unrolling a small number of iterations of an iterative reconstruction scheme ("deep unrolling"), our method is based on deep equilibrium models that correspond to a potentially infinite number of iterations. This framework yields more accurate reconstructions that the current state-of-the-art across a range of inverse problems and gives users the ability to navigate a tradeoff between reconstruction computation time and accuracy during inference. Furthermore, because our formulation is based on finding a fixed points of a operator, we can use standard acceleration techniques to speed inference computations -something that is not possible with deep unrolling methods. In addition, our approach inherits provable convergence guarantees depending on the "base" algorithm used to select a fixed point equation for the deep equilibrium framework. Experimental results also show that our proposed initialization for Deep Equilibrium Models based on pre-training is superior to random initialization, and the proposed approach is more robust to noise than past methods. Overall, the proposed DEM approach is a unique bridge between conventional fixed-point methods in numerical analysis and learning-based techniques for inverse problems in imaging.

Review of Deep Unrolling Methods
Deep unrolling methods for solving inverse problems in imaging consist of a fixed number of architecturally identical "blocks," which are often inspired by particular optimization algorithm. These methods represent the current state-of-the-art in MRI reconstruction, with most top submissions to the fastMRI challenge [6] being some sort of unrolled net. Deep unrolling architectures have also been successfully applied to other inverse problems in imaging, such as low-dose CT [7], light-field photography [8], and emission tomography [9].
We describe here a specific deep unrolling method based on the gradient descent algorithm, although many other variants exist based on alternative optimization or fixed point iteration schemes [1]. Suppose we have a known regularization function r that could be applied to an image x; e.g., in Tikhonov regularization,   Figure 2: (a) PSNR of reconstructed images for an MRI reconstruction problem as a function of iterations used to compute the reconstruction. Unrolled methods are optimized for a fixed computational budget during training, and running additional steps at inference yields a significant drop in performance. Our deep equilibrium methods can achieve the same PSNR for the optimal computational budget of an unrolled method, but can trade slightly more computation time for a significant increase in the PSNR, allowing a user to choose the desired computational budget and reconstruction quality. (b) Standard unrolled deep optimization networks typically require choosing some fixed number of iterates during training. Deviating from this fixed number at inference incurs a significant penalty in PSNR. The forward model here is 8x accelerated single-coil MRI reconstruction, and the unrolled algorithm is unrolled proximal gradient descent with K iterates, labeled PROX-K (Fig. 4). For further experimental details see §6.4.
r(x) = λ 2 x 2 2 for some scalar λ > 0. Then we could compute an image estimate x by solving the optimization problem If r is differentiable, this can be accomplished via gradient descent. That is, we start with an initial estimate x (0) such as x (0) = A y and choose a step size η > 0, such that for iteration k = 1, 2, 3, . . ., we set where ∇r is the gradient of the regularizer. The basic idea behind deep unrolled methods is to fix some number of iterations K (typically K ranges from 5 to 10), declare that x (K) will be our estimate x, and model ∇r with a neural network, denoted R θ (x), whose weights θ can be learned from training data. For example, we may define the unrolled gradient descent estimate to be x (K) (y; θ) := x (K) where x (0) = A y and for k = 0, . . . , K − 1 we have the recursive update Training attempts to minimize the cost function with respect to the network weights θ. This form of training is often called "end-to-end"; that is, we do not train the network R θ that replaces ∇r in isolation, but rather on the quality of the resulting estimate x (K) , which depends on the forward model A.
Above we assume that all instances of R θ have identical weights θ, although other works explore variants where the R θ has iteration dependent weights [10].
The number of iterations in deep unrolling methods is kept small for two reasons. First, at deployment, these systems are optimized to compute image estimates quickly -a desirable property we wish to retain in developing new methods. Second, it is challenging to train deep unrolled networks for many iterations due to memory limitations of GPUs because the memory required to calculate the backpropagation updates scales linearly with the number of unrolled iterations.
As one potential workaround, suppose we train a deep unrolled method for small number of iterations K (e.g., K = 5), then extract the learned regularizer gradient R θ and at inference time run the iterative scheme (3) until convergence (i.e., for more iterations K than used in training). Our numerical results highlight how poorly this approach performs in practice ( §6.4). Choosing a sufficiently large number of iterations K (and hence the computational budget for inference) at training time is essential. As we illustrate in Fig. 2(b), one cannot deviate from the choice of K used in training and expect good performance.

Review of Deep Equilibrium Models
In [5], the authors propose a method for training arbitrarily-deep networks given by the repeated application of a single layer. More precisely, consider an L-layer network with input y and weights θ. Letting x (k) denote the output of the k th hidden layer, we may write where k is the layer index and f (k) θ is a nonlinear transformation such as inner products followed by the application of a nonlinear activation function. Recent prior work explored forcing this transformation at each layer to be the same (i.e. weight tying), so that f (k) θ = f θ for all k and showed that such networks still yield competitive performance [11,12]. Under weight tying, we have the recursion The limit of x (K) as K → ∞, provided it exists, is a fixed point of the operator f θ (·, y). In [5] the authors show that the network weights θ can be learned with constant memory using implicit differentiation, bypassing computation and numerical stability issues associated with related techniques on large-scale problems [13,14]. This past work focused on sequence models and time-series tasks, assuming that each f θ was a single layer of a neural network, and did not explore the image reconstruction task that is the focus of this paper. Following the posting of a preprint of this paper, [15] propose "fixed-point networks" using a strategy similar to ours, independently verifying the potential of this framework in image reconstruction.

Plug-and-Play and Regularization by Denoising Methods
Initiated by [16], a collection of methods based on the plug-and-play (PnP) framework have been proposed, allowing denoising algorithms to be used as priors for model-based image reconstruction. The starting point of PnP is to write reconstructed image as the minimizer of a cost function given by a sum of a data-fit term and a regularizer as in (2). Applying alternating directions method of multipliers (ADMM, [17,18]) to this minimization problem gives a collection of update equations, one of which has the form arg min where r(x) is the regularizer and σ > 0 is a parameter; this update can be considered as a "denoising" of the image z. PnP methods replace this explicit optimization step with a "plugged-in" denoising method.
Notably, some state-of-the-art denoisers (e.g., BM3D [4] and U-nets [19]) do not have an explicit r associated with them, but nevertheless empirically work well within the PnP framework. A related framework called Regularization by Denoising (RED) [20] is based on a similar philosophy as PnP, but instead considers an explicit regularizer of the form where ρ(x) corresponds to an image denoising function. Recent PnP and RED efforts focuses on using training data to learn denoisers [21,22,23,24,25]. In contrast to the unrolling methods described in §2.1, these methods are not trained end-to-end; rather, the denoising module is trained independent of the inverse problem at hand (i.e., independent of the forward model A). As described by [1], decoupling the training of the learned component from A results in a reconstruction system that is flexible and does not need to be re-trained for each new A, but can require substantially more training samples to achieve the reconstruction accuracy of a method trained end-to-end for a specific A.

Proposed Approach
Our approach is to design an iteration map f θ (· ; y) so that a fixed-point x (∞) satisfying is a good estimate of the image x given its measurements y.
Here we describe choices of f θ (and hence of the implicit infinite-depth neural network architecture) that explicitly account for the forward model A and generally for the inverse problem at hand. Specifically, we propose choosing f θ based on different optimization algorithms applied to regularized least squares problem (2). This approach is similar to a DU approach (see 2.1), but where the number of iterations is effectively infinite -a paradigm that has been beyond the reach of all previous deep unrolling architectures for solving inverse problems in imaging. Below we consider three specific choices of f θ , but we note that many other options are possible.

Deep Equilibrium Gradient Descent (DE-Grad)
Connecting the unrolled gradient descent iterations in (3) with the deep equilibrium model in (4), we let Recall that in this setting R θ is a trainable network that replaces the gradient of the regularizer. See Figure  3 for a block diagram illustrating this choice of f θ .

Deep Equilibrium Proximal Gradient Descent (DE-Prox)
Proximal gradient methods [26] use a proximal operator associated with a function h: Specifically, the proximal gradient descent algorithm applied to the optimization problem in (2) yields the iterates where η > 0 is a step size. Similar to the deep unrolling approach of [27], we consider replacing prox ηr with a trainable network R θ : R n → R n , which gives the iteration map See Figure 4 for a block diagram illustrating this choice of f θ . Figure

Deep Equilibrium Alternating Directions Method of Multipliers (DE-ADMM)
The Alternating Directions Method of Multipliers (ADMM, [17]) is an efficient first-order algorithm for large-scale constrained optimization problems. ADMM can be applied to the uncontrained optimization problem (2) by rewriting it as the equivalent constrained problem The augmented Lagrangian (in its "scaled form" -see [17]) associated with this problem is given by where u is an additional auxiliary variable and α > 0 is a user-defined parameter. The ADMM iterates are then Here the z-and x-updates simplify as As in the DE-Prox approach, prox αr (·) can be replaced with a learned network, denoted R θ . Making this replacement, and substituting z (k+1) directly into the expressions for x (k+1) and u (k+1) gives: Note that the updates for x (k+1) and u (k+1) depend only on the previous iterates x (k) and u (k) . Therefore, the above updates can be interpreted as fixed-point iterations on the joint variable q = (x, u), where the iteration map f θ (q; y) is implicitly defined as the map that satisfies Here we take the estimated image to be ) is a fixed-point of f θ (·; y). See Figure 5 for a block diagram illustrating this choice of f θ .

Calculating forward passes and gradient updates
Given a choice of iteration map f θ (·; y) defining a DEM, we confronted the following obstacles. (1) Forward calculation: given an observation y and network weights θ, we need to be able to compute a fixed point of f θ (·; y) efficiently.
(2) Training: given a collection of training samples {x i } n i=1 , we need to find the optimal network weights θ.

Calculating Fixed-Points
Both training and inference in a DEM require calculating a fixed point of the iteration map f θ (·; y) given some initial point y. The most straightforward approach is to use fixed-point iterations given in (4). Convergence of this scheme for specific f θ designs is discussed in Section 5.
However, fixed-point iterations may not converge quickly. By viewing unrolled deep networks as fixedpoint iterations, we inherit the ability to accelerate inference with standard fixed-point accelerators. To our knowledge, this work is the first time iterative inversion methods incorporating deep networks have been accelerated using fixed-point accelerators.
Anderson Acceleration Anderson acceleration [28] 1 utilizes past iterates to identify promising directions to move during the iterations. This takes the form of identifying a vector α (k) ∈ R m and setting, for β > 0 The vector α (k) is the solution to the optimization problem: where G is a matrix whose i th column is the (vectorized) residual The optimization problem in (12) admits a least-squares solution, adding negligible computational overhead when m is small (e.g., m = 5).
In Section 6.4 we compare the performance and time characteristics of solving using Anderson acceleration with standard fixed-point iterations, as well as a technique which uses Broyden's method (a quasi-Newton algorithm) to find fixed points, as proposed in [5].
An important practical consideration is that accelerating fixed-point iterations arising from optimization algorithms with auxiliary variables (like ADMM) is non-trivial. Our implementation of DE-ADMM accelerates ADMM using the results of [30]. However, in general acceleration is not required to learn to solve inverse problems, and for other algorithms or settings standard fixed-point iterations may be attractive for their simplicity of implementation.

Gradient Calculation
In this section, we provide a brief overview of the training procedure used to train all networks in Section 6.4. We use stochastic gradient descent to find network parameters θ that (locally) minimize a cost function of the form 1 is a given loss function, x * i is the ith training image with paired measurements y i , and x (∞) (y i ; θ) denotes the reconstructed image given as the fixed-point of f θ (· ; y i ). For our image reconstruction experiments, we use the mean-squared error (MSE) loss: To simplify the calculations below, we consider gradients of the cost function with respect to a single training measurement/image pair, which we denote (y, x * ). Following [5], we leverage the fact that x (∞) := x (∞) (y; θ) is a fixed-point of f θ (·; y) to find the gradient of the loss with respect to the network parameters θ without backpropagating through an arbitrarily-large number of fixed-point iterations. We summarize this approach below.
First, abbreviating (x (∞) , x ) by , then by the chain rule the gradient of with respect to the network parameters is given by where ∂x (∞) ∂θ is the Jacobian of x (∞) with respect to θ, and ∂ ∂x (∞) is the gradient of with respect to its first argument evaluated at x (∞) . Since we assume is the MSE loss, the gradient ∂ ∂x (∞) is simply the residual between x and the equilibrium point: ∂ ∂x (∞) = x (∞) − x . Now, in order to compute the Jacobian ∂x (∞) ∂θ we start with the fixed point equation: Differentiating both sides of this equation, and solving for ∂x (∞) Plugging this expression into (14) gives This converts the memory-intensive task of backpropagating through many iterations of f θ (· ; y) to the problem of calculating an inverse Jacobian-vector product. To approximate the inverse Jacobian-vector product, first we define the vector β (∞) by Following [29], we note that β = β (∞) is a fixed point of the equation and the same machinery used to calculate the fixed point x (∞) may be used to calculate β (∞) . For analysis purposes, we note that the limit of fixed-point iterations for solving (16) with initial iterate β (0) = 0 is equivalent to the Neumann series: Convergence of the above Neumann series is discussed in Section 5. Conventional autodifferentiation tools permit quickly computing the vector-Jacobian products in (16) and (17). Once an accurate approximation to β (∞) is calculated, the gradient in (14) is given by The gradient calculation process is summarized in the following steps, assuming a fixed point x (∞) of f θ (· ; y) is known: 1. Compute the residual r = x ∞ − x * .

Convergence Theory
Here we study convergence of the proposed deep equilibrium models to a fixed-point at inference time, i.e., given the iteration map f θ (· ; y) : R d → R d we give conditions that guarantee the convergence of the iterates x (k+1) = f θ (x (k) ; y) to a fixed-point x (∞) as k → ∞. Classical fixed-point theory ensures that the iterates converge to a unique fixed-point if the iteration map Below we give conditions on the regularization network R θ : R d → R d (replacing the gradient or proximal mapping of a regularizer) used in the DE-Grad, DE-Prox and DE-ADMM models that that ensure the resulting iteration map is contractive and thus the fixed-point iterations for these models converge.
In particular, following [22], we assume that the regularization network R θ satisfies the following condition: there exists an > 0 such that for all x, x ∈ R d we have where In other words, we assume the map R θ − I is -Lipschitz.
If we interpret R θ as a denoising or de-artifacting network, then R θ − I is the map that outputs the noise or artifacts present in a degraded image. In practice, often R θ is implemented with a residual "skip-connection", such that R θ = I + N θ , where N θ is, e.g., a deep U-net. Therefore, in this case, (19) is equivalent to assuming the trained network N θ is -Lipschitz.
First, we have the following convergence result for DE-Grad: Theorem 1 (Convergence of DE-Grad). Assume that R θ − I is -Lipschitz (19), and let L = λ max (A A) and µ = λ min (A A), where λ max (·) and λ min (·) denote the maximum and minimum eigenvalue, respectively. If the step-size parameter η > 0 is such that η < 1/(L + 1), then the DE-Grad iteration map f θ (·; y) defined in (6) satisfies for all x, x ∈ R d . The coefficient γ is less than 1 if < 1 + µ, in which case the the iterates of DE-Grad converge.
Proof. Let f θ (x; y) be the iteration map for DE-Grad. The Jacobian of f θ (x; y) with respect to x ∈ R d , denoted by ∂ x f θ (x; y), is given by To prove f θ (· ; y) is contractive it suffices to show ∂ x f θ (x; y) < 1 for all x ∈ R d where · denotes the spectral norm. Towards this end, we have where λ i denotes the ith eigenvalue of A A, and in the final inequality (20) we used our assumption that the map (R θ − I)(x) := R θ (x) − x is -Lipschitz, and therefore the spectral norm of its Jacobian ∂ x R θ (x) − I is bounded by . Finally, by our assumption η < 1 1+L where L := max i λ 1 , we have η < 1 1+λi for all i, which implies (1 − η) − ηλ i > 0 for all i. Therefore, the maximum in (20) is obtained at µ := min i λ i , which gives This shows f θ is γ-Lipschitz with γ = 1 − η(1 + µ) + η , proving the claim.
Convergence of PnP approaches PnP-Prox and PnP-ADMM is studied in [22]. At inference time, the proposed DE-Prox and DE-ADMM methods are equivalent to the corresponding PnP method but with a retrained denoising network R θ . Therefore, the convergence results in [22] apply directly to DE-Prox and DE-ADMM. To keep the paper self-contained, we restate these results below, specialized to the case of the quadratic data-fidelity term assumed in (2).
Theorem 3 (Convergence of DE-ADMM). Assume that R θ − I is -Lipschitz (19), and let L = λ max (A A) and µ = λ min (A A) > 0, where λ max (·) and λ min (·) denote the maximum and minimum eigenvalue, respectively. Then the iteration map f θ (·; y) for DE-ADMM defined in (11) is contractive if the ADMM step-size parameter α parameter satisfies See Corollary 1 of [22]. Unlike the convergence result for DE-Grad given in Theorem 1, the convergence results for DE-Prox and DE-ADMM in Theorem 2 and Theorem 3 make the assumption that λ min (A A) > 0, i.e., A has a trivial nullspace. This is condition is satisfied for certain inverse problems, such as denoising or deblurring, but violated in many others, including compressed sensing and undersampled MRI. However, in practice we observe that the iterates of DE-Prox and DE-ADMM still appear to converge even in situations where A has a nontrivial nullspace, indicating this assumption may be stronger than necessary.
Finally, an important practical concern when training deep equilibrium models is whether the fixed-point iterates used to compute gradients (as detailed in Section 4.2) will converge. Specifically, the gradient of the loss at the training pair (y, x * ) involves computing the truncated Neumann series in (17). This series converges if the Jacobian ∂ x f θ (x; y) has spectral norm strictly less than 1 when evaluated at any x ∈ R d , which is true if and only if the iteration map f θ (· ; y) is contractive. Therefore, the same conditions in Theorems 1-3 that ensure the iteration map is contractive also ensure that the Neumann series in (17) used to compute gradients converges.

Comparison Methods and Inverse Problems
Our numerical experiments include comparisons with a variety of models and methods. Total-variation Regularized Least Squares (TV) is an important baseline that does not use any training data but rather leverages geometric models of image structure [3,31,32]. The PnP and RED methods are described in §2.3; we consider both the original ADMM variant of [16] PnP-ADMM and a proximal gradient PnP-Prox method as described in [22]. We utilize the ADMM formulation of RED. Deep Unrolled methods (DU) are described in §2.1; we consider DU using gradient descent, proximal gradient, and ADMM. The preconditioned Neumann network [33] does not have simple Deep Equilibrium or Plug-and-Play analogues, and is included as an alternative deep unrolled method. All deep unrolled methods have tied weights, i.e. the network used at each iteration is the same. This is done to ensure that the number of parameters in the deep unrolled, deep equilibrium, and Plug-and-Play/RED methods are the same. Moreover, on modestly-sized datasets recent work has shown that tied weights result in better reconstructions [34].
We compare the above approaches across three inverse problems: Image deblurring (Deblur), compressed sensing (CS), and accelerated MRI reconstruction (MRI). In the Deblur setting we simulate blurry images  (2) refers the setting of additive white Gaussian noise with standard deviation σ = 0.0001 (low noise). In the CS setting, we take linear measurements of the image by forming inner products with random Gaussian vectors whose entries i.i.d. standard normals, and use an undersampling factor of 4×, i.e., the corresponding forward model A has 4× is a random Gaussian matrix fewer rows than columns. The measurements are then corrupted with additive white Gaussian noise with standard deviation σ = 0.01. In the MRI setting, we investigate recovery at 4× and 8× acceleration, where "ρ× acceleration" indicates an undersampling factor of ρ in k-space (not to be confused with the acceleration techniques used in finding fixed points). Our MRI experiments focus on the case of (virtual) single-coil MR data acquired on a Cartesian grid in k-space, where corresponding forward model A is a subsampling of rows of the discrete Fourier transform matrix. We additionally add complex white Gaussian noise to the undersampled k-space measurements with standard deviation σ = 0.01. For the Deblur and CS problems, we utilize a subset of the Celebrity Faces with Attributes (CelebA) dataset [35], which consists of centered human faces. We train on a subset of 10000 of the training images. All images are resized to 128×128. For the MRI problem, we use a random subset of size 2000 of the fastMRI single-coil knee dataset [36] for training. We trim the ground truth MR images to a 320×320 pixel region-of-interest.

Architecture Specifics
For our learned network, we utilize a DnCNN architecture as in [22]. (We also experimented with U-Nets, but found that DnCNN yielded superior performance for both our proposed deep eq methods and the comparison methods.) For both the CelebA and fastMRI datasets, we train six DnCNN denoisers with noise variances σ 2 = 0.1, 0.05, 0.02, 0.01, 0.005, 0.001 on the training split. Training follows the methodology of [22]. Specifically, to ensure contractivity of the learned component, we add spectral normalization to all layers, ensuring that each layer has a Lipschitz constant bounded above by 1. This normalization is enforced during pretraining as well as during the Deep Equilibrium training phrase.
During training, we utilize Anderson acceleration for both the forward and backward pass fixed-point iterations. In backward passes, the number of fixed-point iterations was limited to 50 due to memory constraints, but fixed-point iterations in forward passes were run until convergence was observed (defined to be when the relative norm difference between iterations is less than 10 −3 ). Test time results were produced using the same method as the forward pass during training, but over the test set. We compare different fixed-point calculation methods for the forward and backward passes in Section 6.5.
Further details on settings, parameter choices, and data may be found in the supplementary materials and in our publicly-available code. 2

Parameter Tuning and Pretraining
Each of the iterative optimization algorithms we test has its own set of hyperparameters to choose, e.g., the step size η in DE-Grad, plus any parameters used to calculate the initial estimate x (0) . Tuning hyperparameters requires choosing a particular regularization network R θ during tuning. We choose from a collection of R θ that have been pretrained for Gaussian denoising at different noise levels. Pretraining can be done on the training dataset (e.g., training on MRI images directly) or using an independent dataset (e.g., the BSD500 image dataset [37]). We use the former approach in our experiments.
To tune hyperparameters, we first choose parameters to optimize the performance of PnP on a validation set via a grid search over pretrained R θ as well as algorithm-specific hyperparameters (such as the η stepsize parameter in gradient descent approaches). Then, we use the PnP hyperparameter settings as initial hyperparameter values when training Deep Equilibrium or Deep Unrolling methods.

Main Results
We present the main reconstruction accuracy comparison in Table 1. Each entry for Deep Equilibrium (DE), Regularization by Denoising (RED), and Plug-and-Play (PnP) approaches is the result of running fixed-point iterations until the relative change between iterations is less than 10 −3 . During training, all DE models were limited to a maximum 100 forward updates, but terminate iterations on the relative norm difference between iterations falling below 10 −3 . The DU models are tested at the number of iterations for which they were trained and all parameters for TV reconstructions (including number of TV iterations) are cross-validated to maximize PSNR. Performance as a function of iteration is shown in Figs. 2(a), 6(a), and 6(b), with example reconstructions in Fig. 7. Further example reconstructions are available for qualitative evaluation in the supplementary materials.
We observe our Deep Equilibrium-based approaches consistently outperform Deep Unrolled approaches across different choices of base algorithm (i.e., Grad,Prox,ADMM). Among choices of iterative reconstruction architectures for f θ , there does not appear to be an obvious winner, suggesting the optimal choice may be problem-or setting-dependent.
As seen in Figs. 2(a), 6(a), and 6(b), the Deep Equilibrium approach generally outperforms Deep Unrolled solvers. Figs. 2(a), 6(a), and 6(b) show that our approach requires no more computation than Deep Unrolled networks to achieve the same performance level and has an increasing advantage with further computation.

Effect of Acceleration
Here we demonstrate the effect of using different fixed-point solvers during both the training and inference procedures. Leveraging acceleration can decrease computational costs during both training and inference and result in better empirical performance at inference. Table 6.5 compares Anderson accelerated Deep Equilibrium approaches with Deep Unrolling, Plug and Play, Deep Equilibrium utilizing Broyden's Method (as was used in [5]), and non-accelerated Deep Equilibrium. All results were determined using PyTorch utilizing an NVidia RTX 2080 Ti GPU. As mentioned previously, non-accelerated Deep Equilibrium at inference time has identical per-iteration cost and memory requirements as Plug and Play and Deep Unrolling.
We observe that while finding the fixed-points using Broyden's method and Anderson acceleration requires more time per iterate, convergence occurs faster than a standard fixed-point iteration so the net time spent at inference is less. Since the Broyden solution was slightly worse in terms of PSNR, Anderson acceleration was used for all other experiments.
For practical matters, the memory cost of each of the compared methods may also be an important factor to consider. At train time, Deep Unrolled methods require memory scaling linearly with the number of iterations used, while Deep Equilibrium methods require only the memory necessary to compute the gradient in (16), which is what permits training at convergence. Plug and Play methods are the least memory-intensive of all to train. At inference time, each method only needs to store at most a constant number of iterations, so all methods are cheap in terms of memory to evaluate.

Effect of Pre-Training
Here we compare the effect of initializing the learned component R θ in our deep equilibrium models with a pretrained denoiser versus initializing with random weights. We use the same hyperparameter tuning scheme described in Section 6.3. We present our results on Deep Equilibrium Proximal Gradient Descent (DE-Prox) in Figure 8(a). We observe an improvement in reconstruction quality when utilizing our pretraining method compared to a random initialization. We also note that pretraining enables a simple choice of algorithm-specific hyperparameters, such as the initial internal step size for DE-Prox.

Noise Sensitivity
We observe empirically that the Deep Equilibrium approach to training achieves competitive reconstruction quality and increased flexibility with respect to allocating computation budget at inference time. Recent work in deep inversion has questioned these methods' robustness to noise and unexpected inputs [38,39,40]. To examine whether the Deep Equilibrium approach is brittle to simple changes in the noise distribution, we varied the level of Gaussian noise added to the observations at test time and observed the effect on reconstruction quality in a setting where DE-Prox and DU-Prox perform similarly. Fig. 8(b) demonstrates that the Deep Equilibrium model DE-Prox is more robust to variation in the noise level than the analogous Deep Unrolled approach DU-Prox. The forward model used in Fig. 8(b) is 8× MRI reconstruction.

Conclusions
This paper illustrates non-trivial quantitative benefits to using implicitly-defined infinite-depth networks for solving linear inverse problems in imaging. These empirical benefits complement convergence guarantees that are unavailable to widely-used deep unrolling methods. Other recent work has focused on such implicit networks akin to the deep equilibrium models considered here (e.g. [41]). Whether these models could lead to additional advances in image reconstruction remains an open question for future work. Furthermore, while the exposition in this work focused on linear inverse problems, nonlinear inverse problems may be solved with iterative approaches just as well. The conditions under which deep equilibrium methods proposed here may be used on such iterative approaches are an active area of investigation.

Further Qualitative Results
In this section, we provide further visualizations of the reconstructions produced by Deep Equilibrium models and the corresponding Deep Unrolled approaches, beyond those shown in the main body. Figures 11, 10, and 9 are best viewed electronically, and contain the ground-truth images, the measurements (projected back to image space in the case of MRI and compressed sensing), and reconstructions by DU-Prox and DE-Prox.
We also visualize the intermediate iterates in the fixed-point iterations, to further demonstrate the convergence properties of DEMs for image reconstruction. We find that DEMs converge quickly to reasonable reconstructions, and maintain high-quality reconstructions after more than one hundred iterations.

Visualizing Iterates
In Figures 13 and 12 we visualize the outputs of the K'th iteration of the mapping f θ in DE-Prox. We observe that across forward problems, the reconstructions converge to good reconstructions.
We illustrate 90 iterations for compressed sensing and 31 for MRI reconstructions (as iterations terminate at 31 iterations).

Further Experimental Details
In this section we provide further details related to the experimental setup. The input to the deblurring algorithms is the preconditioned measurement (A A + λI) −1 A y, where λ is set to be equal to the noise level σ. For MRI reconstruction and compressed sensing experiments, the input is instead simply A y. The masks used in the MRI reconstruction experiments are based on a Cartesian sampling pattern, as in the standard fastMRI setting. For both 4× and 8×, the center 4% of frequencies are fully sampled, and further frequencies are sampled according to a Gaussian distribution centered at 0 frequency with σ = 1.
The compressed sensing design matrices have entries sampled and scaled so that each entry is drawn from a Gaussian distribution with variance 1/m, where A ∈ R m×n . The same design matrix is used for all learned methods.
Optimization algorithm parameters for RED, Plug-and-Play, and all Deep Equilibrium approaches are all chosen via a logarithmic grid search from 10 −4 to 10 1 with 20 elements in each dimension of the grid. All DU methods were trained for 10 iterations. All testing was done on an NVidia RTX 2080 Ti. All networks were trained on a cluster with a variety of computing resources 3 . Every experiment was run utilizing a single GPU-single CPU setup with less than 12 GB of GPU memory.

K=0
K=10 K=20 K=30 K=31 Figure 12: Sample images and reconstructions for MRI reconstruction with acceleration 4× and additive Gaussian noise with σ = 0.01. Each image represents the output of iterate number K. Below each image is the residual between iterate K and the previously-visualized iterate, or in the case of K = 0, between the input to the network and the output of the initial iterate. In this case, the algorithm stops iterating (as the relative norm between iterations drops below 10 −3 ) at iteration 31. The ground truth may be viewed in the initial column of Figure 9. Best viewed electronically.