Thomas Eboli
email: thomas.eboli@inria.fr

Jian Sun
email: jiansun@xjtu.edu.cn

Jean Ponce
email: jean.ponce@inria.fr

End-to-end Interpretable Learning of Non-blind Image Deblurring

Keywords: Non-blind deblurring, preconditioned fixed-point method, end-toend learning

Non-blind image deblurring is typically formulated as a linear leastsquares problem regularized by natural priors on the corresponding sharp picture's gradients, which can be solved, for example, using a half-quadratic splitting method with Richardson fixed-point iterations for its least-squares updates and a proximal operator for the auxiliary variable updates. We propose to precondition the Richardson solver using approximate inverse filters of the (known) blur and natural image prior kernels. Using convolutions instead of a generic linear preconditioner allows extremely efficient parameter sharing across the image, and leads to significant gains in accuracy and/or speed compared to classical FFT and conjugate-gradient methods. More importantly, the proposed architecture is easily adapted to learning both the preconditioner and the proximal operator using CNN embeddings. This yields a simple and efficient algorithm for non-blind image deblurring which is fully interpretable, can be learned end to end, and whose accuracy matches or exceeds the state of the art, quite significantly, in the non-uniform case.

Introduction

This presentation addresses the problem of non-blind image deblurring-that is, the recovery of a sharp image given its blurry version and the corresponding uniform or non-uniform motion blur kernel. Applications range from photography [START_REF] Hu | Image deblurring using smartphone inertial sensors[END_REF] to astronomy [START_REF] Starck | Astronomical Image and Data Analysis, Second Edition[END_REF] and microscopy [START_REF] Goodman | Introduction to Fourier optics[END_REF]. Classical approaches to this problem include leastsquares and Bayesian models, leading to Wiener [START_REF] Wiener | The Extrapolation, Interpolation, and Smoothing of Stationary Time Series[END_REF] and Lucy-Richardson [START_REF] Richardson | Bayesian-based iterative method of image restoration[END_REF] deconvolution techniques for example. Since many sharp images can lead to the same blurry one, blur removal is an ill-posed problem. To tackle this issue, variational methods [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] inject a priori knowledge over the set of solutions using penalized least-squares. Geman and Yang [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF] introduce an auxiliary variable to solve this problem by iteratively evaluating a proximal operator [START_REF] Parikh | Proximal algorithms[END_REF] and solving a least-squares problem. The rest of this presentation builds on this half-quadratic splitting approach. Its proximal part has received a lot of attention through the design of complex model-based [START_REF] Krishnan | Fast image deconvolution using hyper-laplacian priors[END_REF][START_REF] Sun | Image super-resolution using gradient profile prior[END_REF][START_REF] Xu | Unnatural L0 sparse representation for natural image deblurring[END_REF][START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] or learning-based priors [START_REF] Roth | Fields of experts[END_REF]. Far less attention had been paid to the solution of the companion least-squares problem, typically relying on techniques such as conjugate gradient (CG) descent [START_REF] Boyd | Convex Optimization[END_REF] or fast Fourier transform (FFT) [START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF][START_REF] Xu | Inverse kernels for fast spatial deconvolution[END_REF]. CG is relatively slow in this context, and it does not exploit the fact that the linear operator corresponds to a convolution. FFT exploits this property but is only truly valid under periodic conditions at the boundaries, which are never respected by real images.

We propose instead to use Richardson fixed-point iterations [START_REF] Kelley | Iterative Methods for Linear and Nonlinear Equations[END_REF] to solve the leastsquares problem, using approximate inverse filters of the (known) blur and natural image prior kernels as preconditioners. Using convolutions instead of a traditional linear preconditioner allows efficient parameter sharing across the image, which leads to significant gains in accuracy and/or speed over FFT and conjugate-gradient methods. To further improve performance and leverage recent progress in deep learning, several recent approaches to denoising and deblurring unroll a finite number of proximal updates and least-squares minimization steps [START_REF] Aljadaany | Douglas-rachford networks: Learning both the image prior and data fidelity terms for blind image deconvolution[END_REF][START_REF] Chen | Trainable nonlinear reaction diffusion: A flexible framework for fast and effective image restoration[END_REF][START_REF] Kobler | Variational networks: Connecting variational methods and deep learning[END_REF][START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF][START_REF] Schmidt | Shrinkage fields for effective image restoration[END_REF]. Compared to traditional convolutional neural networks (CNNs), these algorithms use interpretable components and produce intermediate feature maps that can be directly supervised during training [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF][START_REF] Schmidt | Shrinkage fields for effective image restoration[END_REF].

We propose a solver for non-blind deblurring, also based on the splitting scheme of [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF] but, in addition to learning the proximal operator as in [START_REF] Zhang | Learning fully convolutional networks for iterative non-blind deconvolution[END_REF], we also learn parameters in the fixed-point algorithm by embedding the preconditioner into a CNN whose bottom layer's kernels are the approximate filters discussed above. Unlike the algorithm of [START_REF] Zhang | Learning fully convolutional networks for iterative non-blind deconvolution[END_REF], our algorithm is trainable end to end, and achieves accuracy that matches or exceeds the state of the art. Furthermore, in contrast to other state-of-the-art CNN-based methods [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF][START_REF] Zhang | Learning fully convolutional networks for iterative non-blind deconvolution[END_REF] relying on FFT, it operates in the pixel domain and thus easily extends to non-uniform blurs scenarios.

Related work

Uniform image deblurring. Classical priors for natural images minimize the magnitude of gradients using the 2 [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], 1 [START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF] and hyper-Laplacian [START_REF] Krishnan | Fast image deconvolution using hyper-laplacian priors[END_REF] (semi) norms or parametric potentials [START_REF] Sun | Image super-resolution using gradient profile prior[END_REF]. Instead of restoring the whole image at once, some works focus on patches by learning their probability distribution [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] or exploiting local properties in blurry images [START_REF] Michaeli | Blind deblurring using internal patch recurrence[END_REF]. Handcrafted priors are designed so that the optimization is feasible and easy to carry out but they may ignore the characteristics of the images available. Data-driven priors, on the other hand, can be learned from a training dataset (e.g., new regularizers). Roth and Black [START_REF] Roth | Fields of experts[END_REF] introduce a learnable total variation (TV) prior whose parameters are potential/filter pairs. This idea has been extended in shallow neural networks in [START_REF] Schmidt | Shrinkage fields for effective image restoration[END_REF][START_REF] Schmidt | Discriminative non-blind deblurring[END_REF] based on the splitting scheme of [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF]. Deeper models based on Roth and Black's learnable prior had been proposed ever since [START_REF] Chen | Trainable nonlinear reaction diffusion: A flexible framework for fast and effective image restoration[END_REF][START_REF] Kobler | Variational networks: Connecting variational methods and deep learning[END_REF][START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF]. The proximal operator can also be replaced by a CNN-based denoiser [START_REF] Meinhardt | Learning proximal operators: Using denoising networks for regularizing inverse imaging problems[END_REF][START_REF] Zhang | Learning deep CNN denoiser prior for image restoration[END_REF] or a CNN specifically trained to mimic a proximal operator [START_REF] Zhang | Learning fully convolutional networks for iterative non-blind deconvolution[END_REF] or the gradient of the regularized [START_REF] Gong | Learning deep gradient descent optimization for image deconvolution[END_REF]. More generally, CNNs are now used in various image restoration tasks, including non-blind deblurring by refining a low-rank decomposition of a deconvolution filter kernel [START_REF] Xu | Deep convolutional neural network for image deconvolution[END_REF], using a FFT-based solver [START_REF] Schuler | Learning to deblur[END_REF], or to improve the accuracy of splitting techniques [START_REF] Aljadaany | Douglas-rachford networks: Learning both the image prior and data fidelity terms for blind image deconvolution[END_REF].

Non-uniform image deblurring. Non-uniform deblurring is more challenging than its uniform counterpart [START_REF] Chakrabarti | Analyzing spatially-varying blur[END_REF][START_REF] Cho | Removing non-uniform motion blur from images[END_REF]. Hirsch et al. [START_REF] Hirsch | Efficient filter flow for space-variant multiframe blind deconvolution[END_REF] consider large overlapping patches and suppose uniform motion on their supports before removing the local blur with an FFTbased uniform deconvolution method. Other works consider pixelwise locally linear motions [START_REF] Brooks | Learning to synthesize motion blur[END_REF][START_REF] Couzinie-Devy | Learning to estimate and remove nonuniform image blur[END_REF][START_REF] Kim | Segmentation-free dynamic scene deblurring[END_REF][START_REF] Sun | Learning a convolutional neural network for non-uniform motion blur removal[END_REF] as simple elements representing complex global motions and solve penalized least-squares problems to restore the image. Finally, geometric non-uniform blur can be used in the case of camera shake to predict motion paths [START_REF] Tai | Richardson-lucy deblurring for scenes under a projective motion path[END_REF][START_REF] Whyte | Non-uniform deblurring for shaken images[END_REF].

Main contributions

Our contributions can be summarized as follows.

-We introduce a convolutional preconditioner for fixed-point iterations that efficiently solves the least-squares problem arising in splitting algorithms to minimize penalized energies. It is faster and/or more accurate than FFT and CG for this task, with theoretical convergence guarantees. -We propose a new end-to-end trainable algorithm that implements a finite number of stages of half-quadratic splitting [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF] and is fully interpretable. It alternates between proximal updates and preconditioned fixed-point iterations. The proximal operator and linear preconditioner are parameterized by CNNs in order to learn these functions from a training set of clean and blurry images. -We evaluate our approach on several benchmarks with both uniform and nonuniform blur kernels. We demonstrate its robustness to significant levels of noise, and obtain results that are competitive with the state of the art for uniform blur and significantly outperforms it in the non-uniform case.

Proposed method

Let y and k respectively denote a blurry image and a known blur kernel. The deconvolution (or non-blind deblurring) problem can be formulated as

min x 1 2 ||y -k x|| 2 F + λΩ( n i=1 k i x), (1)
where " " is the convolution operator, and x is the (unknown) sharp image. The filters k i (i = 1, . . . , n) are typically partial derivative operators, and Ω acts as a regularizer on x, enforcing natural image priors. One often takes Ω(z) = ||z|| 1 (TV-1 model). We propose in this section an end-to-end learnable variant of the method of half-quadratic splitting (or HQS) [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF] to solve Eq. (1). As shown later, a key to the effectiveness of our algorithm is that all linear operations are explicitly represented by convolutions. Let us first introduce notations that will simplify the presentation. Given some linear filters a i and b i (i = 0, . . . , n) with finite support (square matrices), we borrow the Matlab notation for "stacked" linear operators, and denote by A = [a 0 , . . . , a n ] and B = [b 0 ; . . . ; b n ] the (convolution) operators respectively obtained by stacking "horizontally" and "vertically" these filters, whose responses are

A x = [a 0 x, . . . , a n x]; B x = [(b 0 x) , . . . , (b n x) ] ; (2)
We also define A B = n i=0 a i b i and easily verify that (A B) x = A (B x).

A convolutional HQS algorithm

Equation (1) can be rewritten as

min x,z 1 2 ||y -k x|| 2 F + λΩ(z) such that z = F x, (3)
where

F = [k 1 ; . . . ; k n ].
Let us define the energy function

E(x, z, µ) = 1 2 ||y -k x|| 2 F + λΩ(z) + µ 2 ||z -F x|| 2 F . (4)
Given some initial guess x for the sharp image, (e.g. x = y) we can now solve our original problem using the HQS method [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF] with T iterations of the form

z ← argmin z E(x, z, µ); x ← argmin x E(x, z, µ); µ ← µ + δt. (5)
The µ update can vary with iterations but must be positive. We could also use the alternating direction method of multipliers (or ADMM [START_REF] Parikh | Proximal algorithms[END_REF]), for example, but this is left to future work. Note that the update in z has the form

z ← argmin z µ 2 ||z -F x|| 2 F + λΩ(z) = ϕ λ/µ (F x), (6)
where ϕ λ/µ is, by definition, the proximal operator [START_REF] Parikh | Proximal algorithms[END_REF] associated with Ω (a softthresholding function in the case of the 1 norm [START_REF] Elad | Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing[END_REF]) given λ and µ.

The update in x can be written as the solution of a linear least-squares problem:

x ← argmin x 1 2 ||u -L x|| 2 F , (7)
where u = [y;

√ µz] and L = [k;

√ µF ].

Convolutional PCR iterations

Many methods are of course available for solving Eq. [START_REF] Cho | Removing non-uniform motion blur from images[END_REF]. We propose to compute x as the solution of C (u -L x) = 0, where C = [c 0 , . . . , c n ] is composed of n + 1 filters and is used in preconditioned Richardson (or PCR) fixed-point iterations [START_REF] Kelley | Iterative Methods for Linear and Nonlinear Equations[END_REF]. Briefly, in the generic linear case, PCR is an iterative method for solving a square, nonsingular system of linear equations Ax = b. Given some initial estimate x = x 0 of the unknown x, it repeatedly applies the iterations

x ← x -C(Ax -b), (8)
where C is a preconditioning square matrix. When C is an approximate inverse of A, that is, when the spectral radius η of Id -CA is smaller than one, preconditioned Richardson iterations converge to the solution of Ax = b with a linear rate proportional to η [START_REF] Kelley | Iterative Methods for Linear and Nonlinear Equations[END_REF]. When A is an m × n matrix with m ≥ n, and x and b are respectively This method is easily adapted to our context. Since L corresponds to a bank of filters of size w k × w k , it is natural to take C = [c 0 , . . . , c n ] to be another bank of n + 1 linear filters of size w c × w c . Unlike a generic linear preconditioner satisfying CA ≈ Id in matrix form, whose size depends on the square of the image size, C exploits the structure of L and is a linear operator with much fewer parameters, i.e. n + 1 times the size of the c i 's. Thus, C is an approximate inverse filter bank for L, in the sense that

δ ≈ L C = C L = c 0 k + √ µ n i=1 c i k i , (9)
where δ is the Dirac filter. In this setting, C is computed as the solution of

C = argmin C ||δ -L C|| 2 F + ρ n i=0 ||c i || 2 F , (10)
The classical solution using the pseudo inverse of L has cost O (w k + w c -1) 2×3 .

c i = F -1 K * i ρJ + n j=0 | Kj | 2 for i = 0 to n, (11)
using the fast Fourier transform (FFT) with cost O(w 2 c log(w c )) [START_REF] Folland | Fourier Analysis and its Applications[END_REF]. F -1 is the inverse Fourier transform, J is a matrix full of ones, Ki is the Fourier transform of k i (with k 0 = k), K * i is its complex conjugate and the division in the Fourier domain is entrywise. Note that the use of FFT in this context has nothing to do with its use as a deconvolution tool for solving Eq. [START_REF] Cho | Removing non-uniform motion blur from images[END_REF]. Figure 1 shows an example of a blur kernel from [START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF], its approximate inverse when n = 0 and the result of their convolution. Let us define [A] as the linear operator such that [A] B = A B. Indeed, a true inverse filter bank such that equality holds in Eq. ( 9) does not exist in general (e.g., a Gaussian filter cannot be inverted), but all that matters is that the linear operator associated with δ -C L has a spectral radius smaller than one [START_REF] Kelley | Iterative Methods for Linear and Nonlinear Equations[END_REF]. We have the following result.

Lemma 1. The spectral radius of the linear operator Id -[L] [C] , where C is the optimal solution of (10) given by [START_REF] Folland | Fourier Analysis and its Applications[END_REF] is always smaller than 1 when [L] has full rank.

A detailed proof can be found in the supplemental material. We now have our basic nonblind deblurring algorithm, in the form of the Matlab-style CHQS (for convolutional HQS, primary) and CPCR (for convolutional PCR, auxiliary) functions below.

function x = CHQS(y, k, F, µ0) x = y; µ = µ0; for t = 0 : T -1 do u = [y; √ µϕ λ/µ (F x)]; L = [k; √ µF ]; C = argmin C ||δ -C L|| 2 F + ρ n i=0 ||ci|| 2 F ; x = CPCR(L, u, C, x); µ = µ + δt; end for end function function x = CPCR(A, b, C, x0) x = x0; for s = 0 : S -1 do x = x -C (A x -b); end for end function

An end-to-end trainable CHQS algorithm

To improve on this method, we propose to learn the proximal operator ϕ and the preconditioning operator C. The corresponding learnable CHQS (LCHQS) algorithm can now be written as a function with two additional parameters θ and ν as follows.

function x = LCHQS(y, k, F, µ0, θ, ν) x = y; µ = µ0; for t = 0 : T -1 do u = [y; √ µϕ θ λ/µ (F x)]; L = [k; √ µF ]; C = argmin C ||δ -C L|| 2 F + ρ n i=0 ||ci|| 2 F ; x = CPCR(L, u, ψ ν (C), x); µ = µ + δt; end for end function
The function LCHQS has the same structure as CHQS but now uses two parameterized embedding functions ϕ θ τ and ψ ν for the proximal operator and preconditioner. In practice, these functions are CNNs with learnable parameters θ and ν as detailed in Sec. 3. Note that θ actually determines the regularizer through its proximal operator. The function LCHQS is differentiable with respect to both its θ and ν parameters. Given a set of training triplets (x (i) , y (i) , k (i) ) (in i = 1, . . . , N ), the parameters θ and ν can thus be learned end-to-end by minimizing

F (θ, ν) = N i=1 ||x (i) -LCHQS(y (i) , k (i) , F, θ, ν)|| 1 , (12)
with respect to these two parameters by "unrolling" the HQS iterations and using backpropagation, as in [START_REF] Chen | Trainable nonlinear reaction diffusion: A flexible framework for fast and effective image restoration[END_REF][START_REF] Zhang | Learning fully convolutional networks for iterative non-blind deconvolution[END_REF] for example. This can be thought of as the "compilation" of a fully interpretable iterative optimization algorithm into a CNN architecture. Empirically, we have found that the 1 norm gives better results than the 2 norm in Eq. ( 12). 3 Implementation and results

Implementation details

Network architectures. The global architecture of LCHQS shares the same pattern than FCNN [START_REF] Zhang | Learning fully convolutional networks for iterative non-blind deconvolution[END_REF], i.e., n = 2 in Eq. ( 1) with k 1 = [1, -1] and k 2 = k 1 , and the model repeats between 1 and 5 stages alternatively solving the proximal problem [START_REF] Chen | Trainable nonlinear reaction diffusion: A flexible framework for fast and effective image restoration[END_REF] and the linear least-squares problem [START_REF] Cho | Removing non-uniform motion blur from images[END_REF]. The proximal operator ϕ θ is the same as the one introduced in [START_REF] Zhang | Learning fully convolutional networks for iterative non-blind deconvolution[END_REF], and it is composed of 6 convolutional layers with 32 channels and 3 × 3 kernels, followed by ReLU non-linearities, except for the last one. The first layer has 1 input channel and the last layer has 1 output channel. The network ψ ν featured in LCHQS is composed of 6 convolutional layers with 32 channels and 3 × 3 kernels, followed by ReLU non-linearities, except for the last one. The first layer has n + 1 input channels (3 in practice with the setting detailed above) corresponding to the filtered versions of x with the c i 's, and the last layer has 1 output channel. The filters c 1 and c 2 are of size 31 × 31. This size is intentionally made relatively large compared to the sizes of k 1 and k 2 because inverse filters might have infinite support in principle. The size of c 0 is twice the size of the blur kernel k. This choice will be explained in Sec. 3.2. In our implementation, each LCHQS stage has its own θ and ν parameters. The non-learnable CHQS module solves a TV-1 problem; the proximal step implements the soft-thresholding operation ϕ with parameter λ/µ and the leastsquares step implements CPCR. The choice of µ will be detailed below.

Datasets.

The training set for uniform blur is made of 3000 patches of size 180×180 taken from BSD500 dataset and as many random 41 × 41 blur kernels synthesized with the code of [START_REF] Chakrabarti | A neural approach to blind motion deblurring[END_REF]. We compute ahead of time the corresponding inverse filters c i and set the size of c 0 to be 83 × 83 with Eq. [START_REF] Folland | Fourier Analysis and its Applications[END_REF] where ρ is set to 0.05, a value we have chosen with FFT applied on images padded with "edgetaper" as recommended in [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF] and nonpadded images for three image formats; effect of the w c0 /w k ratio on performance; comparison of CG, FFT and CPCR for solving [START_REF] Cho | Removing non-uniform motion blur from images[END_REF].

after cross-validation on a separate test set. We also create a training set for non-uniform motion blur removal made of 3000 180 × 180 images synthesized with the code of [START_REF] Gong | From motion blur to motion flow: A deep learning solution for removing heterogeneous motion blur[END_REF] with a locally linear motion of maximal magnitude of 35 pixels. For both training sets, the validation sets are made of 600 additional samples. In both cases, we add Gaussian noise with standard deviation matching that of the test data. We randomly flip and rotate by 90 • the training samples and take 170 × 170 random crops for data augmentation. Optimization. Following [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF], we train our model in a two-step fashion: First, we supervise the sharp estimate output by each iteration of LCHQS in the manner of [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF] with Eq. ( 12). We use an Adam optimizer with learning rate of 10 -4 and batch size of 1 for 200 epochs. Second, we further train the network by supervising the final output of LCHQS with Eq. ( 12) on the same training dataset with an Adam optimizer and learning rate set to 10 -5 for 100 more epochs without the per-layer supervision. We have obtained better results with this setting than using either of the two steps separately.

Experimental validation of CPCR and CHQS

In this section, we present an experimental sanity check of CPCR for solving [START_REF] Cho | Removing non-uniform motion blur from images[END_REF] and CHQS for solving [START_REF] Aljadaany | Douglas-rachford networks: Learning both the image prior and data fidelity terms for blind image deconvolution[END_REF] in the context of a basic TV-1 problem.

Inverse kernel size. We test different sizes for the w c0 × w c0 approximate inverse filter c 0 associated with a w k ×w k blur kernel k, in the non-penalized case, with λ = 0. We use Eq. ( 11) with ρ set to 0.05. We use 160 images obtained by applying the 8 kernels of [START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF] to 20 images from the Pascal VOC 2012 dataset. As shown in Fig. 3 , the PSNR increases with increasing w c0 /w k ratios, but saturates when the ratio is larger than 2.2. We use a ratio of 2 which is a good compromise between accuracy and speed.

CPCR accuracy. We compare the proposed CPCR method to FFT-based deconvolution (FFT) and conjugate gradient descent (CG), to solve the least-squares problem of Eq. ( 7) in the setting of a TV-1 problem. We follow [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF] and, in order to limit boundary artifacts for FFT, we pad the images to be restored by replicating the pixels on the boundary with a margin of half the size of the blur kernel and then use the "edgetaper" routine. We also run FFT on images padded with the "replicate" strategy consisting in simply replicating the pixels on the boundary. We solve Eq. ( 7) with µ 0 = 0.008, ker-1 ker-2 ker-3 ker-4 ker-5 ker-6 ker-7 ker-8 Aver. Time (s) HQS-FFT (no pad.) 21. [START_REF] Gong | From motion blur to motion flow: A deep learning solution for removing heterogeneous motion blur[END_REF] 1) on 160 synthetic blurry images with 2% white noise. We run all the methods on a GPU. The running times are for a 500 × 375 RGB image. λ = 0.003 and z computed beforehand with Eq. ( 6). The 160 images previously synthesized are degraded with 2% additional white noise. Figure 3 shows the average PSNR scores for the three algorithms optimizing Eq. ( 7). After only 5 iterations, CPCR produces an average PSNR higher than the other methods and converges after 10 iterations.

The "edgetaper" padding is crucial for FFT to compete with CG and CPCR by reducing the amount of border artifacts in the solution. CPCR running time. CPCR relies on convolutions and thus greatly benefits from GPU acceleration. For instance, for small images of size 500 × 375 and a blur kernel of size 55 × 55, 10 iterations of CPCR are in the ballpark of FFT without padding: CPCR runs in 20ms, FFT runs in 3ms and FFT with "edgetaper" padding takes 40ms. For a high-resolution 1280 × 720 image and the same blur kernel, 10 iterations of CPCR run in 22ms, FFT without padding runs in 10ms and "edgetaper" padded FFT in 70ms. Figure 3 compares the running times of CPCR (run for 10 iterations) with padded/nonpadded FFT for three image (resp. kernel) sizes: 500 × 375, 800 × 800 and 1280 × 720 (27 × 27, 55 × 55 and 121 × 121) pixels. Our method is marginally slower than FFT without padding in every configuration (within a margin of 20ms) but becomes much faster than FFT combined to "edgetaper" padding when the size of the kernel increases. FFT with "replicate" padding runs in about the same time as FFT (no pad) and thus is not shown in Fig. 3 . The times have been averaged over 1000 runs.

Running times for computing the inverse kernels with Eq. [START_REF] Folland | Fourier Analysis and its Applications[END_REF]. Computing the inverse kernels c i , with an ratio w c /w k set to 2, takes 1.0 ± 0.2ms for a blur kernel k of size 27 × 27 and 5.4 ± 0.5ms (results averaged in 1000 runs) for a large 121 × 121 kernel. Thus, the time for inverting blur kernels is negligible in the overall pipeline.

CHQS validation. We compare several iterations of HQS using unpadded FFT (HQS-FFT (no pad.)), with "replicate" padding (HQS-FFT (rep. pad)), and the padding strategy proposed in [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF] (HQS-FFT (FDN pad.)), CG (HQS-CG), or CPCR (CHQS) for solving the least-squares problem penalized with the TV-1 regularized in Eq. ( 1) and use the same 160 blurry and noisy images than in previous paragraph as test set. We set the number of HQS iterations T to 10, run CPCR for 5 iterations and CG for at most 100 iterations. We use λ = 0.003 and µ t = 0.008 × 4 t (t = 0, . . . , T -1). Table 1 compares the average PSNR scores obtained with the different HQS algorithms over the test set. As expected, FDN padding greatly improves HQS-FFT results on larger kernels over naive "replicate" padding, i.e. "ker-4" and "ker-8", but overall does not perform FCNN [START_REF] Zhang | Learning fully convolutional networks for iterative non-blind deconvolution[END_REF] EPLL [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] RGCD [START_REF] Gong | Learning deep gradient descent optimization for image deconvolution[END_REF] FDN [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF] Table 2: PSNR scores for Levin [START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF] and Sun [START_REF] Sun | Edge-based blur kernel estimation using patch priors[END_REF] benchmarks, that respectively feature 0.5% and 1% noise. Best results are shown in bold, second-best underlined. The difference may not always be significant between FDN and LCHQS for the Levin dataset. as well as CHQS. For kernels 1, 2, 3 and 5, the four methods yield comparable results (within 0.1 dB of each other). FFT-based methods are significantly worse on the other four, whereas our method gives better results than HQS-CG in general, but is 100 times faster. This large speed-up is explained by the convolutional structure of CPCR whereas CG involves large matrix inversions and multiplications. Figure 2 shows a deblurring example from the test set. HQS-FFT (with FDN padding strategy), even with the refined padding technique of [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF], produces a solution with boundary artifacts. Both HQS-CG and CHQS restore the image with a limited amount of artifacts, but CHQS does it much faster than HQS-CG. This is typical of our experiments in practice.

Discussion. These experiments show that CPCR always gives better results than CG in terms of PSNR, sometimes by a significant margin, and it is about 50 times faster. This suggests that CPCR may, more generally, be preferable to CG for linear leastsquares problems when the linear operator is a convolution. CPCR also dramatically benefits from its convolutional implementation on a GPU with speed similar to FFT and is even faster than FFT with FDN padding for large kernels. These experiments also show that CHQS surpasses, in general, HQS-CG and HQS-FFT for deblurring.

Next, we further improve the accuracy of CHQS using supervised learning, as done in previous works blending within a single model variational methods and learning.

Uniform deblurring

We compare in this section CHQS and its learnable version LCHQS with the non-blind deblurring state of the art, including optimization-based and CNN-based algorithms.

Comparison on standard benchmarks. LCHQS is first trained by using the loss of Eq. ( 12) to supervise the output of each stage of the proposed model and second Fig. 5: Performance of FDN [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF] and LCHQS on the Levin [START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF] (left) and Sun [START_REF] Sun | Edge-based blur kernel estimation using patch priors[END_REF] (right) datasets.

trained by only supervising the output of the final layer, in the manner of [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF]. The model trained in the first regime is named LCHQS G and the one further trained with the second regime is named LCHQS F . The other methods we compare our learnable model to are HQS algorithms solving a TV-1 problem: HQS-FFT (with the padding strategy of [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF]), HQS-CG and CHQS, an HQS algorithm with a prior over patches (EPLL) [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] and the state-of-the-art CNN-based deblurring methods FCNN [START_REF] Zhang | Learning fully convolutional networks for iterative non-blind deconvolution[END_REF] and FDN [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF]. We use the best model provided by the authors of [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF], denoted as FDN 10 T in their paper. Table 2 compares our method with these algorithms on two classical benchmarks. We use 5 HQS iterations and 2 CPCR iterations for CHQS and LCHQS. Except for EPLL that takes about 40 seconds to restore an image of the Levin dataset [START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF], all methods restore a 255 × 255 black and white image in about 0.2 second. The dataset of Sun et al. contains high-resolution images of size around 1000 × 700 pixels. EPLL removes the blur in 20 minutes on a CPU while the other methods, including ours, do it in about 1 second on a GPU. In this case, our learnable method gives comparable results to FDN [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF], outputs globally much sharper results than EPLL [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] and is much faster. As expected, non-trainable CHQS is well behind its learned competitors (Tab. 2).

Number of iterations for LCHQS G and CPCR. We investigate the influence of the number of HQS and CPCR iterations on the performance of LCHQS G on the benchmarks of Levin et al. [START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF] and Sun et al. [START_REF] Sun | Edge-based blur kernel estimation using patch priors[END_REF]. FDN implements 10 HQS iterations parameterized with CNNs but operates in the Fourier domain. Here, we compare LCHQS G to the FDN model trained in a stage-wise manner (denoted as FDN 10 G in [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF]). Figure 5 plots the mean PSNR values for the datasets of Levin et al. and Sun et al. [START_REF] Sun | Edge-based blur kernel estimation using patch priors[END_REF] after each stage. FDN comes in two versions: one trained on a single noise level (green line) and one trained on noise levels within a given interval (blue line). We use up to 5 iterations of our learnable CHQS scheme, but it essentially converges after only 3 steps. When the number of CPCR iterations is set to 1, FDN and our model achieve similar results for the same number of HQS iterations. For 2/3 CPCR iterations, we do better than FDN for the same number of HQS iterations by a margin of +0.4/0.5dB on both benchmarks. For 3 HQS iterations and more, LCHQS saturates but systematically achieves better results than 10 FDN iterations: +0.15dB for [START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF] and +0.26dB for [START_REF] Sun | Edge-based blur kernel estimation using patch priors[END_REF].

Robustness to noise. Table 3 compares our methods for various noise levels on the 160 RGB images introduced previously, dubbed from "PASCAL benchmark". FDN corresponds to the model called FDN 10 T in [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF]. For this experiment, (L)CHQS uses 5 Fig. 6: Example of image deblurring with an additive noise of 3% (better seen on a computer screen). In this example, we obtain better PSNR scores than competitors and better visual results, for example details around the door or the leaves.

HQS iterations and 2 inner CPCR iterations. We add 1%, 3% and 5% Gaussian noise to these images to obtain three different test sets with gradually stronger noise levels. We train each model to deal with a specific noise level (non-blind setting) but also train a single model to handle multiple noise levels (blind setting) on images with 0.5 to 5% of white noise, as done in [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF]. For each level in the non-blind setting, we are marginally above or below FDN results. In terms of average PSNR values, the margins are +0.12dB for 1%, +0.06dB for 3% and -0.05dB for 5% when comparing our models with FDN, but we are above the other competitors by margins between 0.3dB and 2dB. Compared to its noise-dependent version, the network trained in the blind setting yields a loss of 0.2dB for 1% noise, but gains of 0.14 and 0.27dB for 3 and 5% noises, showing its robustness and adaptability to various noises. Figure 6 compares results obtained on a blurry image with 3% noise.

Non-uniform motion blur removal

Typical non-uniform motion blur models assign to each pixel of a blurry image a local uniform kernel [START_REF] Chakrabarti | Analyzing spatially-varying blur[END_REF]. This is equivalent to replacing the uniform convolution in Eq. ( 1) by local convolutions for each overlapping patch in an image, as done by Sun et al. [START_REF] Sun | Learning a convolutional neural network for non-uniform motion blur removal[END_REF] when they adapt the solver of [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF] to the non-uniform case. Note that FDN [START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF] and FCNN [START_REF] Zhang | Learning fully convolutional networks for iterative non-blind deconvolution[END_REF] operate in the Fourier domain and thus cannot be easily adapted to nonuniform deblurring, unlike (L)CHQS operating in the spatial domain. We handle nonuniform blur as follows to avoid computing different inverse filters at each pixel. As in [START_REF] Sun | Learning a convolutional neural network for non-uniform motion blur removal[END_REF], we model a non-uniform motion field with locally linear motions that can well approximate complex global motions such as camera rotations. We discretize the set of the linear motions by considering only those with translations (in pixels) in {1, 3 . . . , 35} and orientations in {0 vance all the 511 35 × 35 local blur kernels and compute their approximate inverses ahead of time. During inference, we simply determine which one best matches the local blur kernel and use its approximate inverse in CPCR. This is a parallelizable operation on a GPU. Table 4 compares our approach (in non-blind setting) to existing methods for locally-linear blur removal on a test set of 100 images from PASCAL dataset nonuniformly blurred with the code of [START_REF] Gong | From motion blur to motion flow: A deep learning solution for removing heterogeneous motion blur[END_REF] and with white noise. For instance for 1% noise, LCHQS G scores +0.99dB higher than CG-based method, and LCHQS F pushes the margin up to +1.13dB while being 200 times faster. Figure 7 shows one non-uniform example from the test set.

Deblurring with approximated blur kernels

In practice one does not have the ground-truth blur kernel but instead an approximate version of it, obtained with methods such as [START_REF] Pan | Deblurring images via dark channel prior[END_REF][START_REF] Whyte | Non-uniform deblurring for shaken images[END_REF]. We show that (L)CHQS works well for approximate and/or large filters, different from the ones used in the training set and without any training or fine-tuning. We show in Figure 8 a deblurred image with an approximate kernel obtained with the code of [START_REF] Pan | Deblurring images via dark channel prior[END_REF] and of support of size 101 × 101 pixels. We obtain with LCHQS F (blind) of Fig. 8: Real-world blurry images deblurred with an 101 × 101 blur kernel estimated with [START_REF] Pan | Deblurring images via dark channel prior[END_REF]. We can restore fine details with approximate, large kernels.

not introduce artifacts as FDN, showing the robustness of CPCR and its embedding in HQS to approximate blur kernels. More results are shown in the supplemental material.

Conclusion

We have presented a new learnable solver for non-blind deblurring. It is based on the HQS algorithm for solving penalized least-squares problems but uses preconditioned iterative fixed-point iterations for the x-update. Without learning, this approach is superior both in terms of speed and accuracy to classical solvers based on the Fourier transform and conjugate gradient descent. When the preconditioner and the proximal operator are learned, we obtain results that are competitive with or better than the state of the art. Our method is easily extended to non-uniform deblurring, and it outperforms the state of the art by a significant margin in this case. We have also demonstrated its robustness to important amounts of white noise. Explicitly accounting for more realistic noise models [START_REF] Foi | Practical poissonian-gaussian noise modeling and fitting for single-image raw-data[END_REF] and other degradations such as downsampling is left for future work.

Acknowledgments. This works was supported in part by the INRIA/NYU collaboration and the Louis Vuitton/ENS chair on artificial intelligence. In addition, this work was funded in part by the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR19-P3IA-0001 (PRAIRIE 3IA Institute). Jian Sun was supported by NSFC under grant numbers 11971373 and U1811461.

Fig. 1 :

 1 Fig. 1: From left ot right: An example of a blur kernel k from the Levin et al. dataset [23]; its approximate inverse kernel c 0 ; the resulting filter resulting from the convolution of k and c 0 (represented as a surface). It gives an approximate Dirac filter δ.

Fig. 2 :

 2 Fig. 2: A blurry image from our test set with 2% white noise and the solutions of Eq. (1) with TV-1 regularization obtained with different HQS-based methods. From the same optimization problem, HQS-FFT displays boundary artifacts. HQS-CG and CHQS produce images with similar visual quality and PSNR values but HQS-CG is much slower.

Fig. 3 :

 3 Fig.3: From left to right: Computation times for CPCR (including computation of C) with FFT applied on images padded with "edgetaper" as recommended in[START_REF] Kruse | Learning to push the limits of efficient FFT-based image deconvolution[END_REF] and nonpadded images for three image formats; effect of the w c0 /w k ratio on performance; comparison of CG, FFT and CPCR for solving[START_REF] Cho | Removing non-uniform motion blur from images[END_REF].

Fig. 4 :

 4 Fig. 4: Comparison of state-of-the-art methods and the proposed LCHQS for one sample of the Levin dataset [23] (better seen on a computer screen). FDN effectively removes the blur but introduces artifacts in flat areas, unlike EPLL and LCHQS.

Fig. 7 :

 7 Fig.7: Non-uniform motion deblurring example with 1 % additive Gaussian noise (better seen on a computer screen). The car and the helmet are sharper with our method than in the images produced by our competitors.

Table 1 :

 1 20.51 22.31 18.21 23.36 20.01 19.93 19.02 20.69 0.07 HQS-FFT (rep. pad.) 26.45 25.39 26.27 22.75 27.64 27.26 24.84 23.54 25.53 0.07 HQS-FFT (FDN pad.) 26.48 25.89 26.27 23.79 27.66 27.23 25.26 25.02 25.96 0.15 HQS-CG 26.39 25.90 26.24 24.88 27.59 27.31 25.39 25.19 26.12 13 CHQS 26.45 25.96 26.26 25.06 27.67 27.51 25.81 25.48 26.27 0.26 Comparison of different methods optimizing the same TV-1 deconvolution model (

  .12 35.11 ± 0.05 35.15 ± 0.04 Sun [37] 32.24 32.46 31.95 32.67 30.36 32.83 ± 0.01 32.93 ± 0.01

					CHQS LCHQS G	LCHQS F
	Levin [23]	33.08	34.82	33.73	35.09 32

  (non-blind) 29.54 ± 0.02 25.59 ± 0.03 23.87 ± 0.06 0.7 LCHQS F (non-blind) 29.53 ± 0.02 25.56 ± 0.03 23.95 ± 0.05 0.7 LCHQS G (blind) 29.22 ± 0.02 25.55 ± 0.03 24.05 ± 0.02 0.7 LCHQS F (blind) 29.35 ± 0.01 25.71 ± 0.02 24.21 ± 0.01 0.7

		1% noise	3% noise	5% noise Time (s)
	HQS-FFT	26.48	23.90	22.15	0.2
	HQS-CG	26.45	23.91	22.27	13
	EPLL [46]	28.83	24.00	22.10	130
	FCNN [44]	29.27	25.07	23.53	0.5
	FDN [22]	29.42	25.53	23.97	0.6
	CHQS	27.08	23.33	22.38	0.3
	LCHQS

• , 6 • , . . . , 174 • }. In this case, we know in ad-G

Table 3 :

 3 Uniform deblurring on 160 test images with 1%, 3% and 5% white noise. Running times are for an 500 × 375 RGB image. The mention "blind" (resp. "nonblind) indicates that a single model handles the three (resp. a specific) noise level(s).

Table 4 :

 4 Table 3 a sharper result than FCNN and do HQS-FFT HQS-CG EPLL [46] CHQS LCHQS G LCHQS F 1% noise 23.49 25.84 25.49 25.11 26.83 ± 0.08 26.98 ± 0.08 3% noise 23.17 24.18 23.78 23.74 24.91 ± 0.05 25.06 ± 0.06 5% noise 22.44 23.10 23.34 22.65 23.97 ± 0.05 24.14 ± 0.05 Non-uniform deblurring on 100 test images with 1%, 3% and 5% white noise. Running times are for an 500 × 375 RGB image.

	Time (s)	13	212	420	0.8	0.9	0.9


Pierre-Franc ¸ois Loos
email: loos@irsamc.ups-tlse.fr

Yann Damour

Anthony Scemama
email: scemama@irsamc.ups-tlse.fr

J P F Leblanc

A E Antipov

F Becca

I W Bulik

-L Chan

C.-M Chung

Y Deng

M Ferrero

T M Henderson

C A Jiménez-Hoyos

M Motta

D M Ceperley

J A Gomez

E Gull

S Guo

T N Lan

J Li

F Ma

K T Williams

Y Yao

L Chen

H Shi

C Niu

U Ray

R J Anderson

The performance of CIPSI on the ground state electronic energy of benzene

we report the performance of the Configuration Interaction using a Perturbative Selection made Iteratively (CIPSI) method on the nonrelativistic frozen-core correlation energy of the benzene molecule in the cc-pVDZ basis. Following our usual protocol, we obtain a correlation energy of -863.4 mE h which agrees with the theoretical estimate of -863 mE h proposed by Eriksen et al. using an extensive array of highly-accurate new electronic structure methods.

Although sometimes decried, one cannot deny the usefulness of benchmark sets and their corresponding reference data for the electronic structure community. These are indeed essential for the validation of existing theoretical models and to bring to light and subsequently understand their strengths and, more importantly, their weaknesses. In that regard, the previous benchmark datasets provided by the Simons Collaboration on the Many-Electron Problem have been extremely valuable. [1][2][3] The same comment applies to the excited-state benchmark set of Thiel and coworkers. [4][5][6][7][8] Following a similar goal, we have recently proposed a large set of highly-accurate vertical transition energies for various types of excited states thanks to the renaissance of selected configuration interaction (SCI) methods [9][10][11] which can now routinely produce near full configuration interaction (FCI) quality excitation energies for small-and medium-sized organic molecules. [12][13][14][15][16] In a recent preprint, 17 Eriksen et al. have proposed a blind test for a particular electronic structure problem inviting several groups around the world to contribute to this endeavour. In addition to coupled cluster theory with singles, doubles, triples, and quadruples (CCSDTQ), 18,19 a large panel of highlyaccurate, emerging electronic structure methods were considered: (i) the many-body expansion FCI (MBE-FCI), [20][21][22][23] (ii) three SCI methods including a second-order perturbative correction [adaptive sampling CI (ASCI), [24][START_REF] Tubman | An efficient deterministic perturbation theory for selected configuration interaction methods[END_REF][START_REF] Tubman | [END_REF] iterative CI (iCI), [27][28][29][30] and semistochastic heat-bath CI (SHCI) [31][32][33] ], (iii) the full coupled-cluster reduction (FCCR) 34,[START_REF] Xu | Towards near-exact solutions of molecular electronic structure: Full coupled-cluster reduction with a second-order perturbative correction[END_REF] which also includes a second-order perturbative correction, (iv) the densitymatrix renornalization group approach (DMRG), [START_REF] White | [END_REF][37][38] and (v) two flavors of FCI quantum Monte Carlo (FCIQMC), 39,40 namely AS-FCIQMC 41 and CAD-FCIQMC. 42,43 We refer the interested reader to Ref. 17 and its supporting information for additional details on each method and the complete list of references. Soon after, Lee et al. reported phaseless auxiliary-field quantum Monte Carlo 44 (ph-AFQMC) correlation energies for the very same problem. 45 The target application is the non-relativistic frozen-core correlation energy of the ground state of the benzene molecule in the cc-pVDZ basis. The geometry of benzene has been optimized at the MP2/6-31G* level 5 and its coordinates can be found in the supporting information of Ref. 17 alongside its nuclear repulsion and Hartree-Fock energies. This corresponds to an active space of 30 electrons and 108 orbitals, i.e., the Hilbert space of benzene is of the order of 10 For the sake of completeness and our very own curiosity, we report in this Note the frozen-core correlation energy obtained with a fourth flavor of SCI known as Configuration Interaction using a Perturbative Selection made Iteratively (CIPSI), 10 which also includes a second-order perturbative (PT2) correction. In short, the CIPSI algorithm belongs to the family of SCI+PT2 methods. The idea behind such methods is to avoid the exponential increase of the size of the CI expansion by retaining the most energetically relevant determinants only, thanks to the use of a second-order energetic criterion to select perturbatively determinants in the FCI space. However, performing SCI calculations rapidly becomes extremely tedious when one increases the system size as one hits the exponential wall inherently linked to these methods.

From a historical point of view, CIPSI is probably one of the oldest SCI algorithm. It was developed in 1973 by Huron, Rancurel, and Malrieu 10 (see also Ref. 46). Recently, the determinant-driven CIPSI algorithm has been efficiently implemented 47,48 in the open-source programming environ-arXiv:2008.11145v3 [physics.chem-ph] 10 Oct 2020 ment quantum package by our group enabling to perform massively parallel computations. [49][50][51] In particular, we were able to compute highly-accurate ground-and excited-state energies for small-and medium-sized molecules (including benzene). [12][13][14][15][16] CIPSI is also frequently used to provide accurate trial wave function for QMC calculations. 47,48,[52][53][START_REF] Caffarel | Using cipsi nodes in diffusion monte carlo[END_REF][START_REF] Scemama | QMC using very large configuration interaction-type expansions[END_REF][START_REF] Scemama | [END_REF][57][58][59][60][61] The particularity of the current implementation is that the selection step and the PT2 correction are computed simultaneously via a hybrid semistochastic algorithm 49,51 (which explains the statistical error associated with the PT2 correction in the following). Moreover, a renormalized version of the PT2 correction (dubbed rPT2 below) has been recently implemented and tested for a more efficient extrapolation to the FCI limit thanks to a partial resummation of the higher-order of perturbation. 51 We refer the interested reader to Ref. 51 where one can find all the details regarding the implementation of the rPT2 correction and the CIPSI algorithm.

Being late to the party, we obviously cannot report blindly our CIPSI results. However, following the philosophy of Eriksen et al. 17 and Lee et al., 45 we will report our results with the most neutral tone, leaving the freedom to the reader to make up his/her mind. We then follow our usual "protocol" [12][13][14][15][16][57][58][59] by performing a preliminary SCI calculation using Hartree-Fock orbitals in order to generate a SCI wave function with at least 10 7 determinants. Natural orbitals are then computed based on this wave function, and a new, larger SCI calculation is performed with this new natural set of orbitals. This has the advantage to produce a smoother and faster convergence of the SCI energy toward the FCI limit. The total SCI energy is defined as the sum of the variational energy E var. (computed via diagonalization of the CI matrix in the reference space) and a second-order perturbative correction E (r)PT2 which takes into account the external determinants, i.e., the determinants which do not belong to the variational space but are linked to the reference space via a nonzero matrix element. The magnitude of E (r)PT2 provides a qualitative idea of the "distance" to the FCI limit. As mentioned above, SCI+PT2 methods rely heavily on extrapolation, especially when one deals with medium-sized systems. We then linearly extrapolate the total SCI energy to E (r)PT2 = 0 (which effectively corresponds to the FCI limit). Note that, unlike excited-state calculations where it is important to enforce that the wave functions are eigenfunctions of the Ŝ 2 spin operator, [START_REF] Applencourt | Spin adaptation with determinant-based selected configuration interaction[END_REF] the present wave functions do not fulfil this property as we aim for the lowest possible energy of a singlet state. We have found that Ŝ 2 is, nonetheless, very close to zero (∼ 5 × 10 -3 a.u.). The corresponding energies are reported in Table II as functions of the number of determinants in the variational space N det .

A second run has been performed with localized orbitals. Starting from the same natural orbitals, a Boys-Foster localization procedure [START_REF] Foster | 64 Indices of molecular orbitals for Boys-Foster localization procedure: core [1-6[END_REF] was performed in several orbital windows: i) core, ii) valence σ, iii) valence π, iv) valence π * , v) valence σ * , vi) the higher-lying σ orbitals, and vii) the higher-lying π orbitals. 64 Like Pipek-Mezey, 65 this choice of orbital windows allows to preserve a strict σ-π separation in planar systems like benzene. As one can see from the energies of Table II, for a given value of N det , the variational energy as well as the PT2-corrected energies are much lower with localized orbitals than with natural orbitals. Indeed, localized orbitals significantly speed up the convergence of SCI calculations by taking benefit of the local character of electron correlation. 17,[66][67][68][69][70] We, therefore, consider these energies more trustworthy, and we will base our best estimate of the correlation energy of benzene on these calculations. The convergence of the CIPSI correlation energy using localized orbitals is illustrated in Fig. 1, where one can see the behavior of the correlation energy, ∆E var. and ∆E var. + E (r)PT2 , as a function of N det (left panel). The right panel of Fig. 1 is more instructive as it shows ∆E var. as a function of E (r)PT2 , and their corresponding four-point linear extrapolation curves that we have used to get our final estimate of the correlation energy. (In other words, the four largest variational wave functions are considered to perform the linear extrapolation.) From this figure, one clearly sees that the rPT2-based correction behaves more linearly than its corresponding PT2 version, and is thus systematically employed in the following.

Our final number are gathered in Table III, where, following the notations of Ref. 17, we report, in addition to the final variational energies ∆E var. , the extrapolation distances, ∆E dist , defined as the difference between the final computed energy, ∆E final , and the extrapolated energy, ∆E extrap. associated with ASCI, iCI, SHCI, DMRS, and CIPSI. The three flavours of SCI fall into an interval ranging from -860.0 mE h (ASCI) to -864.2 mE h (SHCI), while the other non-SCI methods yield correlation energies ranging from -863.7 to -862.8 mE h (see Table I). Our final CIPSI number (obtained with localized orbitals and rPT2 correction via a four-point linear extrapolation) is -863.4(5) mE h , where the error reported in parenthesis represents the fitting error (not the extrapolation error for which it is much harder to provide a theoretically sound estimate). 71 For comparison, the best post blind test SHCI estimate is -863.3 mE h , which agrees almost perfectly with our best CIPSI estimate, while the best post blind test ASCI and iCI correlation energies are -861.3 and -864.15 mE h , respectively (see Table III).

The present calculations have been performed on the AMD partition of GENCI's Irene supercomputer. Each Irene's AMD node is a dual-socket AMD Rome (Epyc) CPU@2.60 GHz with 256GiB of RAM, with a total of 64 physical CPU cores per socket. These nodes are connected via Infiniband HDR100. The first step of the calculation, i.e., performing a CIPSI calculation up to N det ∼ 10 7 with Hartree-Fock orbitals in order to produce natural orbitals, takes roughly 24 hours on a single node, and reaching the same number of determinants with natural orbitals or localized orbitals takes roughly the same amount of time. A second 24-hour run on 10 distributed nodes was performed to push the selection to 8 × 10 7 determinants, and a third distributed run using 40 nodes was used to reach 16 × 10 7 determinants. In total, the present calculation has required 150k core hours, most of it being spent in the last stage of the computation.

  [START_REF] Xu | Towards near-exact solutions of molecular electronic structure: Full coupled-cluster reduction with a second-order perturbative correction[END_REF] Slater determinants. Needless to say that this size of Hilbert space cannot be tackled by exact diagonalization with current architectures. The correlation energies reported in Ref. 17 are gathered in

TABLE I .

 I The frozen-core correlation energy (in mE h ) of benzene in the cc-pVDZ basis set using various methods.

	Method	E c	Ref.
	ASCI	-860.0	17
	iCI	-861.1	17
	CCSDTQ	-862.4	17
	DMRG	-862.8	17
	FCCR	-863.0	17
	MBE-FCI	-863.0	17
	CAD-FCIQMC	-863.4	17
	AS-FCIQMC	-863.7	17
	SHCI	-864.2	17
	ph-AFQMC	-864.3(4)	45
	CIPSI	-863.4	This work
	Table I alongside the best ph-AFQMC estimate from Ref. 45
	based on a CAS(6,6) trial wave function. The outcome of this
	work is nicely summarized in the abstract of Ref. 17: "In our
	assessment, the evaluated high-level methods are all found to
	qualitatively agree on a final correlation energy, with most
	methods yielding an estimate of the FCI value around -863
	mE h . However, we find the root-mean-square deviation of
	the energies from the studied methods to be considerable (1.3
	mE h ), which in light of the acclaimed performance of each of
	the methods for smaller molecular systems clearly displays the
	challenges faced in extending reliable, near-exact correlation
	methods to larger systems."

TABLE II .

 II Convergence of the CIPSI correlation energy of benzene using localized orbitals. Left: ∆E var. , ∆E var. + E PT2 , and ∆E var. + E rPT2 (in mE h ) as functions of the number of determinants in the variational space N det . Right: ∆E var. (in mE h ) as a function of E PT2 or E rPT2 . The four-point linear extrapolation curves (dashed lines) are also reported. The theoretical estimate of -863 mE h from Ref. 17 is marked by a black line for comparison purposes. The statistical error bars associated with E PT2 or E rPT2 (not shown) are of the order of the size of the markers. Variational energy E var. , second-order perturbative correction E PT2 and its renormalized version E rPT2 (in E h ) as a function of the number of determinants N det for the ground-state of the benzene molecule computed in the cc-pVDZ basis set. The statistical error on E (r)PT2 , corresponding to one standard deviation, are reported in parenthesis.The data that support the findings of this study are openly available in Zenodo at http://doi.org/10.5281/zenodo.4075286.

	-400							-600
	-500							-650
								-700
	-600
								-750
	-700
								-800
	-800
								-850
	1000 -900	10 4	10 5	10 6	10 7	10 8	-200	-150	-100	-50	0
	FIG. 1. Natural orbitals				Localized orbitals
	N det		E var.		E var. + E PT2		E var. + E rPT2	E var.		E var. + E PT2		E var. + E rPT2
	1 280	-230.978 056	-231.559 025(212)	-231.463 633(177)	-231.101 676	-231.519 522(149)	-231.472 224(132)
	2 560	-231.043 712	-231.542 344(139)	-231.474 885(120)	-231.161 264	-231.515 577(155)	-231.482 477(140)
	5 120	-231.115 142	-231.534 122(213)	-231.488 815(190)	-231.224 632	-231.516 375(191)	-231.495 022(177)
	10 240	-231.188 813	-231.531 660(516)	-231.502 992(473)	-231.283 295	-231.520 907(271)	-231.507 708(255)
	20 480	-231.260 065	-231.534 172(611)	-231.517 063(573)	-231.330 209	-231.526 433(586)	-231.518 045(561)
	40 960	-231.321 906	-231.538 269(501)	-231.528 301(478)	-231.366 008	-231.532 288(303)	-231.526 639(293)
	81 920	-231.366 895	-231.541 945(813)	-231.535 785(785)	-231.392 888	-231.536 578(614)	-231.532 575(597)
	163 840	-231.392 866	-231.545 499(761)	-231.541 010(739)	-231.414 132	-231.541 400(624)	-231.538 378(609)
	327 680	-231.407 802	-231.548 699(662)	-231.544 980(645)	-231.431 952	-231.545 873(557)	-231.543 532(545)
	655 360	-231.418 752	-231.551 208(661)	-231.548 004(645)	-231.447 007	-231.548 856(498)	-231.547 043(489)
	1 310 720	-231.428 852	-231.552 760(616)	-231.550 006(603)	-231.460 970	-231.552 137(453)	-231.550 723(446)
	2 621 440	-231.439 324	-231.553 845(572)	-231.551 544(560)	-231.473 751	-231.555 261(403)	-231.554 159(397)
	5 242 880	-231.450 156	-231.557 541(534)	-231.555 558(524)	-231.485 829	-231.558 303(362)	-231.557 451(358)
	10 485 760	-231.461 927	-231.559 390(481)	-231.557 796(474)	-231.497 515	-231.562 568(322)	-231.561 901(319)
	20 971 520	-231.474 019	-231.561 315(430)	-231.560 063(424)	-231.508 714	-231.564 707(275)	-231.564 223(273)
	41 943 040	-231.487 978	-231.564 529(382)	-231.563 593(377)	-231.519 122	-231.567 419(240)	-231.567 069(238)
	83 886 080	-231.501 334	-231.566 994(317)	-231.566 325(314)	-231.528 568	-231.570 084(199)	-231.569 832(198)
	167 772 160	-231.514 009	-231.569 939(273)	-231.569 467(271)	-231.536 655	-231.571 981(175)	-231.571 804(174)
	European Union's Horizon 2020 research and innovation pro-
	gramme (Grant agreement No. 863481).

TABLE III .

 III Extrapolation distances, ∆E dist , defined as the difference between the final computed energy, ∆E final , and the extrapolated energy, ∆E extrap. associated with ASCI, iCI, SHCI, DMRG, and CIPSI for the best blind-test and post-blind-test estimates of the correlation energy of benzene in the cc-pVDZ basis. The final variational energies ∆E var. are also reported. See Ref. 17 for more details. All correlation energies are given in mE h .

	Method	∆E var.	∆E final	∆E extrap.	∆E dist
	Best blind-test estimates
	ASCI	-737.1	-835.4	-860.0	-24.6
	iCI	-730.0	-833.7	-861.1	-27.4
	SHCI	-827.2	-852.8	-864.2	-11.4
	DMRG	-859.2	-859.2	-862.8	-3.6
	Best post-blind-test estimates
	ASCI	-772.4	-835.2	-861.3	-26.1
	iCI	-770.7	-842.8	-864.2	-21.3
	SHCI	-835.2	-854.9	-863.6	-8.7
	CIPSI	-814.8	-850.2	-863.4	-13.2

We thank Janus Eriksen and Cyrus Umrigar for useful comments. This work was performed using HPC resources from GENCI-TGCC (2020-gen1738) and from CALMIP (Toulouse) under allocation 2020-18005. PFL and AS have received funding from the European Research Council (ERC) under the