Voichit ¸a Maxim
email: voichita.maxim@creatis.insa-lyon.fr

Enhancement of Compton camera images reconstructed by inversion of a conical Radon transform

Keywords: image reconstruction, conical Radon transform, Compton transform, analytic algorithms, Compton camera, total variation

We present a new inversion formula for a weighted conical Radon transform modelling Compton camera data. The formula exploits a large proportion of the acquired events and is easy to implement into fast algorithms. We give for it two equivalent formulations relying on known properties of the two-dimensional Radon transform and we test a semi-iterative algorithm for one of them. From a practical point of view, methods robust to measurement noise and to low number of events are required. We show that adding a constraint on the total variation of the final image strongly improves the results. We illustrate our arguments with Monte-Carlo simulated data in both low and realistic noise configurations.

Introduction

The recent interest in the conical Radon transforms was driven by the development of Compton cameras for medical and homeland security imaging applications. The idea of using Compton scattering of γ-rays as a mean to detect their incoming direction and to exploit this information for imaging was first proposed in the early 1970s independently by two groups, Schönfelder et al [START_REF] Schönfelder | A telescope for soft gamma ray astronomy[END_REF] and Todd et al [START_REF] Todd | A proposed Gamma camera[END_REF] in the fields of γ-ray astronomy and nuclear medicine, respectively. This new technology was successfully implemented in the COMPton TELescope COMPTEL, born in 1991 on a satellite Compton Gamma-Ray Observatory for nine years. Other projects of cameras boarded on stratospheric balloon flights were realized since then. The first tests for nuclear medicine application were reported in 1983 by Singh and Doria [START_REF] Singh | An electronically collimated gamma camera for single photon emission computed tomography. Part II: Image reconstruction and preliminary experimental measurements[END_REF]. With the technological progress, medical applications could in a near future benefit from the large acceptance angle and the large energy spectrum allowed by the Compton imaging devices. Nowadays, the application fields of Compton cameras were extended to imaging in homeland security and nuclear decommissioning ( [START_REF] Wulf | Compton imager for detection of special nuclear material[END_REF][START_REF] Lee | A compact Compton camera using scintillators for the investigation of nuclear materials[END_REF][START_REF] Takeda | A portable Si/CdTe Compton camera and its applications to the visualization of radioactive substances[END_REF][START_REF] Wahl | The Polaris-H imaging spectrometer[END_REF]).

The guiding thread of this study is to keep as close as possible to real applications, where the incompleteness of the data set and the acquisition noise reduce the quality of the images. Section 2 is devoted to the introduction of the main concepts employed in the paper. It contains a brief presentation of the Compton camera followed by a short stateof-the-art on conical Radon transforms, an introduction to total variation denoising ( [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]) and to the Simultaneous Iterative Reconstruction Technique (SIRT, [START_REF] Gilbert | Iterative methods for the three-dimensional reconstruction of an object from projections[END_REF]). In section 3 we derive an inversion formula for a class of weighted conical Radon transforms, for cones with vertices in a plane (the scatterer) and arbitrary axis directions and half-opening angles. The filtered backprojection inversion formula we propose is closely related to the one derived in [START_REF] Maxim | Filtered backprojection reconstruction and redundancy in Compton camera imaging[END_REF]. We propose three methods for the numerical computation of the inverse, two of them being analytic and the third semi-iterative. As the noise and the low number of acquired events lead to poor reconstructed images, the three of them are associated with total variation denoising or regularization. Finally we evaluate the methods with Monte-Carlo simulated data reproducing the acquisition of a complex shaped source with a finite-size detector. Two levels of uncertainties on the position and energy measurements are considered in the simulation. Section 4 describes the data, generated with a Monte-Carlo software from nuclear physics. In section 5 we show some numerical examples. We demonstrate that the methods can reconstruct the images from statistically distributed counts, but with a resolution depending on the geometry of acquisition and on the noise level.

Introduction to Compton camera imaging and to total variation denoising

In this section we introduce the principle of Compton camera imaging and we give a concise state-of-the-art on conical Radon transforms and their inverses. We also recall some widely used techniques that will be used in the next sections. They are respectively the Simultaneous Iterative Reconstruction Technique (SIRT, [START_REF] Gilbert | Iterative methods for the three-dimensional reconstruction of an object from projections[END_REF]) that is largely employed for tomographic reconstruction and the total variation (TV) denoising methodology originally introduced by Rudin et al in [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF].

Compton scattering and conical projections

In nuclear medicine imaging Compton cameras might eventually surpass the collimated Anger cameras and allow either reduction of the radio-tracer dose injected to the patient or reduction of the acquisition time ( [START_REF] Fontana | Compton camera study for high efficiency SPECT and benchmark with Anger system[END_REF]). In other medical applications where γ imaging serves to the monitoring of the treatment, the number of emitted photons is too small to allow three dimensional imaging with collimated cameras. Additionally, the collimators hardly cope with a large energy spectrum. For this reasons, applications as treatment monitoring in proton therapy or selective internal radiation therapy are even of more interest but also more challenging (see e.g., [START_REF] Frandes | A tracking Compton-scattering imaging system for hadron therapy monitoring[END_REF][START_REF] Kurosawa | Prompt gamma detection for range verification in proton therapy[END_REF][START_REF] Hilaire | Proton therapy monitoring by Compton imaging: influence of the large energy spectrum of the prompt-γ radiation[END_REF][START_REF] Aldawood | Development of a compton camera for prompt-gamma medical imaging[END_REF][START_REF] Muñoz | Tests of macaco compton telescope with 4.44 mev gamma rays[END_REF]).

Figure 1: Inelastic scattering of a photon on an electron. A part of the energy of the photon is transferred to the electron and the last one is ejected from the atom's shell. In the simplistic model of the electron at rest the two angles, β and ϕ, are uniquelly determined by the energies.

Compton cameras detect γ rays that interact with the detector at least twice. The first interaction should be a Compton scattering, where a part of the energy of the ray is transferred to the recoiling electron (see figure 1). The last one deposits its energy in the detector, making possible to identify its position V 1 and to measure the energy E 1 lost by the photon. The γ particle should further interact at least once at some point V 2 . When the energies E 1 of the recoiled electron and E 0 of the γ ray are known, the scattering angle may be deduced through the formula

cos β = 1 - m e c 2 E 1 (E 0 -E 1 )E 0 , ( 1
)
with m e the mass of the electron at rest and c the speed of the light. In the absence of noise, the point where the gamma ray was emitted lies on the surface of the cone having vertex at V 1 , axis directed by the unitary vector

Ω = --→ V 2 V 1 / --→ V 2 V 1 and half-opening angle β.
Here denotes the usual Euclidean norm of a vector in the three-dimensional space. As a consequence, data acquired with a Compton camera having ideal position and energy resolutions use to be modelled as integrals of the weighted intensity of the source on conical surfaces. These models are gathered under the name of weighted conical Radon transforms, cone transforms or Compton transforms.

The weighted conical Radon transform

More precisely, the weighted conical Radon transform use to be defined as follows. For some given cone, let us denote α ∈ [0, π] and δ ∈ (-π, π] the polar and azimuthal angles of the axis direction Ω, u ∈ R 3 the vector of coordinates of its vertex and β ∈ (0, π) its half-opening angle. The cone may then be defined as the set

C(u, α, δ, β) = {v ∈ R 3 : (v -u).Ω = v -u cos β}, ( 2
)
and the projection associated to it as

C d α,δ,β f (u) = C(u,α,δ,β) f (v) 1 v -u d dv, ( 3
)
Figure 2: Schematic representation of a Compton camera composed of a scatterer (green) and an absorber (blue). The source is represented as the word RADON. Photons emitted from the source are first scattered then absorbed in the detector, and the position where they were emitted can only be determined at this stage as belonging to a conical surface.

where f is the intensity function and d use to be an integer. The weighted conical Radon transform of a function f integrable with respect to the measure vu -d dv over almost all cones with parameters in some manifold

M ∈ R 3 ×[0, π]×(-π, π]×(0, π), is the application (u, α, δ, β) ∈ M → C d α,δ,β f (u). (4)
A first restriction to some M may be imposed by the integrability condition. A second restriction comes from the particularity of the conical Radon transform to be redundant, since the dimension of the data space is larger than the dimension of the image space. Consequently, its inversion requires to select a sub-set from the data and, when possible, to average several inverses in order to reduce the noise.

In the last years some inversion methods accompanied by adapted selection techniques were developed. In [START_REF] Basko | Application of spherical harmonics to image reconstruction for the Compton camera[END_REF], Basko et al have shown that the three-dimensional Radon projections on planes containing a given point can be calculated from the set of projections on cones having the vertex at that point, arbitrary axis direction and given constant half-opening angle. Cone-beam projections (also called divergent-beam projections) can be evaluated from conical projections having the same vertex (see e.g., the works of L. Parra [START_REF] Parra | Reconstruction of cone-beam projections from Compton scattered data[END_REF] and Tomitani and Hirasawa [START_REF] Tomitani | Analytical image reconstruction of cone-beam projections from limited-angle Compton camera data[END_REF]). Detectors placed on a line were considered by Jung and Moon in [START_REF] Jung | Exact inversion of the cone transform arising in an application of a Compton camera consisting of line detectors[END_REF].

The projections on cones with axes perpendicular to some surface provide a complete set. Inversion formulas in such configurations were given by Cree and Bones [START_REF] Cree | Towards direct reconstruction from a gamma camera based on Compton scattering[END_REF] and M. Haltmeier [START_REF] Haltmeier | Exact reconstruction formulas for a Radon transform over cones[END_REF] for the plane, Schiefeneder and Haltmeier [START_REF] Schiefeneder | The Radon transform over cones with vertices on the sphere and orthogonal axes[END_REF] for the sphere, S. Moon [START_REF] Moon | Inversion of the conical radon transform with vertices on a surface of revolution arising in an application of a Compton camera[END_REF] for surfaces of revolution and Moon and Haltmeier [START_REF] Moon | Analytic inversion of a conical Radon transform arising in application of Compton cameras on the cylinder[END_REF] for the cylinder. Inversion methods that restrict the useful set of projections to cones with axes perpendicular on a manifold will generally neglect a large part of the acquired data and thus provide very noisy images. A particular case where they may be of interest is the one of small detectors that move around the patient (see e.g., [START_REF] Kishimoto | First demonstration of multi-color 3-D in vivo imaging using ultra-compact Compton camera[END_REF]). Such detectors can anyway measure projections only for small values of α.

When the values of the projections are known for single axis direction and halfopening angle but for all vertices u ∈ R 3 , the mapping (4) can be seen as a convolution of f with a distribution supported on the cone. Deconvolution formulas were given by Gouia-Zarrad and Ambartsoumian ( [START_REF] Gouia-Zarrad | Exact inversion of the conical Radon transform with a fixed opening angle[END_REF]), F. Terzioglu ( [START_REF] Terzioglu | Some inversion formulas for the cone transform[END_REF]) and V. Palamodov ([29]). The inconvenient of this approach for the application to Compton camera is that measurements are required in all points of the space, including inside the source. This is obviously impossible to realize since the detectors are always placed outside the body. The methods proposed by B. Smith ([30]), Kuchment and Terzioglu ([31,[START_REF] Kuchment | Inversion of weighted divergent beam and cone transforms[END_REF]), Terzioglu et al ( [START_REF] Terzioglu | Compton camera imaging and the cone transform: a brief overview[END_REF]) and one of the methods given by V. Palamodov in [START_REF] Palamodov | Reconstruction from cone integral transforms[END_REF] allow to select the points of some manifold M ⊂ R 3 that do not intersect the source. The manifolds M are characterised by some admissibility (see [START_REF] Kuchment | Three-dimensional image reconstruction from Compton camera data[END_REF][START_REF] Kuchment | Inversion of weighted divergent beam and cone transforms[END_REF]) or completness (see [START_REF] Smith | Reconstruction methods and completeness conditions for two Compton data models[END_REF]) conditions that are very similar to each other, with specificities related to the inversion formula.

Connected to the topic we discuss here is the inversion of the cumulative (or compound) conical Radon transform addressed e.g., in [START_REF] Nguyen | On an integral transform and its inverse in nuclear imaging[END_REF] by Nguyen and Truong. This transform arises in nuclear medicine imaging when one wants to make an image with the photons scattered inside the patient. The differences in the inversion methods are however important and the problem is more complicated in the sense that the projections are not measured individually.

Conical Radon transform for cones with vertices in a plane

A slightly different class of weighted conical transforms was analysed in [START_REF] Maxim | Analytical inversion of the Compton transform using the full set of available projections[END_REF] under the name of Compton transforms. Their specificity is to employ a weight different from (3), weight that accounts for the incident angle of the incoming particle, and to consider the set of cones with vertices in a plane. Two transforms from the class differ either by the polar angle α or by the half-opening angle β and can be averaged together for noise variance reduction as shown in [START_REF] Maxim | Filtered backprojection reconstruction and redundancy in Compton camera imaging[END_REF]. The main result of [START_REF] Maxim | Analytical inversion of the Compton transform using the full set of available projections[END_REF] is the central-slice theorem from section 4.1, proving the invertibility of the transform. Derived afterwards in [START_REF] Maxim | Filtered backprojection reconstruction and redundancy in Compton camera imaging[END_REF], the associated filtered back-projection formula benefits from a relative simplicity and also from the fact that conical projections with axis directions not orthogonal to the detector are allowed.

To ensure uniqueness of the inverse, it is assumed that all the projections on cones with vertices in some given plane are measured. This completeness condition is obviously not verified in practice. The method can still give reasonably good results for Monte-Carlo simulated data, providing addition of some smoothing. In [START_REF] Maxim | Filtered backprojection reconstruction and redundancy in Compton camera imaging[END_REF], this was obtained by frequency damping. In this paper we propose to use total variation either for denoising the solution or as a regularization step associated with an iterative algorithm.

Total variation regularization

Originally proposed by Rudin et al in [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], total variation penalty function is recommended for denoising piecewise constant images, as it promotes smoothness while still conserving sharp edges. Total variation may also improve tomographic reconstruction of images from low-dose acquisitions (see e.g., [START_REF] Sawatzky | Accurate EM-TV algorithm in PET with low SNR[END_REF]), or from limitedangle data (see e.g., [START_REF] Sidky | Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT[END_REF][START_REF] Persson | Total variation norm for three-dimensional iterative reconstruction in limited view angle tomography[END_REF]).

We suppose without loose of generality that the volume is a cube divided in N 3 voxels, rearranged in a vector g, and we note X = R 3N . The volume is transformed through a linear map of matrix R and the measurements p are noisy realizations of Rg. We search for g solution of the regularized problem:

arg min g∈X 1 2 Rg -p 2 2 + µT V (g) , (5)
with µ ≥ 0 the regularization parameter and T V the total variation norm.

For Ω an open subset of R 3 , the total variation of a function g ∈ L 1 (Ω) is defined by:

T V (g) = sup Ω g(x)div (ϕ(x))dx : ϕ ∈ C 1 c (Ω; R 3 ), ϕ(x) ≤ 1 ∀x ∈ Ω . ( 6
) If g has gradient ∇g ∈ L 1 (Ω), then T V (g) = Ω ∇g(x) dx.
The discrete form of the total variation uses to be written as:

T V (g) = 1≤i,j,k≤N (∇g) i,j,k , (7)
and is not differentiable in zero. Let us put Y = (R 3 )

3N . An equivalent expression for the discrete total variation is:

T V (g) = sup { ϕ, ∇g Y : ϕ ∈ Y such that ϕ i,j,k ≤ 1 for all 1 ≤ i, j, k ≤ N } , (8)
where

ϕ, ψ Y = i,j,k ( 3 =1 ϕ i,j,k ψ i,j,k ) is the Euclidean dot product in Y .
With the appropriate definition of the discrete divergence, ϕ, ∇g Y = -div ϕ, g X . In this work we rely on the notations from [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF], and we apply the Chambolle's algorithm proposed therein to solve the denoising problem when R is the identity matrix.

When R is not the identity matrix a splitting approach can be used, with a gradient descent step that decreases the value of the data fidelity term Rgp 2 and a TVdenoising step that smooths the solution. This approach was previously applied by Beck and Teboulle in [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] for deconvolution problems and the gradient descent algorithm was used for the first task. When R stands for the discrete Radon transform, the SIRT algorithm may be applied. SIRT may be considered as a weighted gradient descent algorithm with more stable results.

The SIRT algorithm for the inversion of the Radon transform

In this paper we make use of the iterative inversion of the Radon transform. The Simultaneous Iterative Reconstruction Technique (SIRT, [START_REF] Gilbert | Iterative methods for the three-dimensional reconstruction of an object from projections[END_REF]) is one of the algorithms widely employed for this task. Let R be the system matrix describing the discrete Radon transform in 2D or the X-ray transform in 3D, p the vector of data and g the vector containing all the voxels of the volume we wish to reconstruct. The SIRT algorithm starts from an initial guess g (0) that may be taken as the null vector or the filtered back-projection solution then iterates following the formula:

g (k+1) = g (k) -λ R * Rg (k) -p R1 / (R * 1) , ( 9
)
where λ is a descent step, 1 is a vector of ones and R * is the Hermitian conjugate of R. The divisions of vectors is done component-wise. In all our numerical tests we take λ = 0.2. The SIRT algorithm may give reconstructions more precise than the filtered back-projection, but when the data is very noisy a regularization of the solution is necessary. The best way to introduce this regularization is to add it as a penalty term in the minimization problem, as explained in section 2.4.

Theoretical and numerical inversion of the conical Radon transform

In this section we give an analytic inversion formula for the transform defined in equation ( 3) for d = 1, where the integral of the intensity of the source is weighted by the inverse of the distance from the source point to the vertex of the cone. We propose three equivalent formulations for the inverse and discuss them from a numerical and discrete point of view. One of the methods we propose is semi-iterative, in the sense that a part of the computation, aiming to transform the Compton data in two-dimensional Radon data, is done analytically. Then an iterative algorithm alternating SIRT ( [START_REF] Gilbert | Iterative methods for the three-dimensional reconstruction of an object from projections[END_REF]) iterations and TV denoising is applied in order to calculate a regularized inverse of the Radon transform. The other methods are analytical and based on the commutativity of the Radon transform with differentiation.

Discussion on the Compton camera data modelling

We consider in this work cameras with planar detectors (see figure 2). In the scope of a simpler presentation and without losing generality, hereafter we consider a camera containing a single scatterer layer placed at altitude z s = 0 and an absorber placed at some altitude z a < 0. When the scatterer has several layers, a conical transform should be defined for each layer. The intensity f is then calculated as the mean of their inverses. We suppose that the probabilities for an incident γ ray to be Compton scattered in the first detector and then completely absorbed in the second detector are both equal to one, so that these probabilities do not appear in the model. The differential cross section of a γ photon having initial energy E 0 and scattered by an angle β, is given by the Klein-Nishina formula ( [START_REF] Klein | Über die streuung von strahlung durch freie elektronen nach der neuen relativistischen quantendynamik von dirac[END_REF]), and is denoted by K(β, E 0 ). For a given cone, let u ∈ R 2 × {0} be the position vector of the vertex and Ω = (sin α cos δ, sin α sin δ, cos α) the direction vector of its axis, where α and δ are its polar and azimuthal angles (see figure 3 (a)). When a source of γ particles having intensity distribution represented by a function f : R 3 → R is observed by a Compton camera, the fraction of events recorded with parameters (u, α, δ, β) was considered in [START_REF] Maxim | Analytical inversion of the Compton transform using the full set of available projections[END_REF] proportional to

C α,δ,β f (u) = K(β, E 0 ) cos α C(u,α,δ,β) f (v) cos θdv = K(β, E 0 ) cos α C(u,α,δ,β) f (v) z v -u dv, ( 10
)
where v = (x, y, z) is a vector from R 3 and θ is the polar angle of the vector vu.

The factors cos θ and cos α from (10) account for the arbitrary incidence angle of the γ ray on the scatterer and absorber, respectively. For any given parameters α, δ, β, with α ∈ [0, π], β ∈ (0, π) and δ ∈ (-π, π], the function C α,δ,β f : R 2 × {0} → R was called in [START_REF] Maxim | Analytical inversion of the Compton transform using the full set of available projections[END_REF] a Compton projection of f (see figure 3 (b)). In order to keep the apelation from [START_REF] Maxim | Analytical inversion of the Compton transform using the full set of available projections[END_REF] and to distinguish it from other conical Radon transforms from the literature, the integral transform C operating on functions f : R 3 → R is hereafter referred to as the Compton transform. The intensity of the print is the integral from [START_REF] Maxim | Filtered backprojection reconstruction and redundancy in Compton camera imaging[END_REF].

In figure 4 are represented a deterministic Compton projection and its counterpart from a Monte-Carlo simulation. The source is a small sphere. The deterministic projection from (a) was calculated in each pixel of the scatterer following equation [START_REF] Maxim | Filtered backprojection reconstruction and redundancy in Compton camera imaging[END_REF]. In (b) a Compton camera acquisition with a small spherical source emitting γ photons was simulated with the software MEGAlib (The Medium-Energy Gamma-ray Astronomy library, [START_REF] Zoglauer | MEGAlib-The medium energy gamma-ray astronomy library[END_REF]), software based on the nuclear physics simulation library Geant4 [START_REF] Agostinelli | Geant4-A simulation toolkit[END_REF]. The simulation in figure 4 accounts for Doppler broadening i.e., small errors on the measured energies produced when the photon is scattered on an electron which is not at rest. In this example, the positions of the hits and the sum of energies are exact.

In [START_REF] Maxim | Filtered backprojection reconstruction and redundancy in Compton camera imaging[END_REF] the weight cos θ only accounts for the incidence angle of the γ ray on the scatterer. In a series of papers (see e.g., [START_REF] Wilderman | Fast algorithm for list mode backprojection of Compton scatter camera images[END_REF]), Wilderman et al suggested that the weight should also depend on the distance vu through the solid angle subtended by the element of detector at the source. The model, was recently validated against Monte-Carlo simulations ( [START_REF] Maxim | Probabilistic models and numerical calculation of system matrix and sensitivity in list-mode MLEM 3D reconstruction of Compton camera images[END_REF]). As the source is supposed to lie in the half-space z > 0, the function f is zero for values of v close to u. Thus there is no singularity in the expression under the integral [START_REF] Fontana | Compton camera study for high efficiency SPECT and benchmark with Anger system[END_REF]. As we shall show in section 3.2, the inversion formula given in [START_REF] Maxim | Filtered backprojection reconstruction and redundancy in Compton camera imaging[END_REF] for [START_REF] Maxim | Filtered backprojection reconstruction and redundancy in Compton camera imaging[END_REF] can be adapted easily to d = 1 and the transform

C * α,δ,β f (u) = K(β, E 0 ) C(u,α,δ,β) f (v) cos θ v -u 2 dv, ( 11
) (a) (b)
C 1 α,δ,β f (u) = K(β, E 0 ) C(u,α,δ,β) f (v) 1 v -u dv, ( 12
)
which is closer to the model suggested by the Monte Carlo simulations and in practice seems to give better results.

Inversion of the conical Radon transform

In this section we derive a formula allowing to calculate the source intensity from the conical Radon projections [START_REF] Frandes | A tracking Compton-scattering imaging system for hadron therapy monitoring[END_REF]. To avoid technical difficulties throughout the following, we assume f to be Schwartz function supported in the upper half space.

Let (e 1 , e 2 , e 3 ) be the standard basis of R 3 . The vectors d 1 = (cos δ, sin δ, 0), d 2 = (-sin δ, cos δ, 0) and e 3 also form an orthonormal basis of R 3 , defining a reference frame Otsz obtained by rotation of Oxyz about the Oz axis. Let us consider the cone C(u 0 , α, δ, β) with u 0 = t 0 d 1 + s 0 d 2 . Since the source belongs to the half-space z > 0, for a single-layer scatterer one should always have α ∈ 0, π 2 . No restriction apply on δ, which is taken in the interval (-π, π]. Physically, the Compton angle β can take any value in (0, π). However the parametrization we choose hereafter for the cone constrains to keep only the Compton projections with α + β < π 2 , leading to a selection of the cones that intersect the scatterer in a single point, the vertex. For such admissible values of α and β one has cos 2 α -sin 2 β > 0. We define:

a(α, β) = sin β cos β cos 2 α -sin 2 β , ( 13
) b(α, β) = sin β cos 2 α -sin 2 β , ( 14
) c(α, β) = sin α cos α cos 2 α -sin 2 β . ( 15
)
The intersection of the cone with a horizontal plane at some given altitude z > 0 is an ellipse that may be described by the parametric equation:

t = t 0 + zc(α, β) +za(α, β) cos ϕ s = s 0 +zb(α, β) sin ϕ , ϕ ∈ (-π, π]. ( 16
)
With the cone thought of as a stack of ellipses, its parametric equation may then be readily obtained. The surface integral (12) at a point from the scatterer having coordinates u 0 = t 0 d 1 + s 0 d 2 , may then be expressed as

C 1 α,δ,β f (t 0 d 1 + s 0 d 2 ) = K(β, E 0 )b(α, β) × ∞ 0 2π 0 f ((t 0 + zc(α, β) + za(α, β) cos ϕ)d 1 + (s 0 + zb(α, β) sin ϕ)d 2 + ze 3 )dϕdz. (17)
Hereafter we make use of the two-dimensional Radon transform, that associates to a function g : R 2 → R with bounded support a set indexed on the parameter δ ∈ (-π, π] of line integrals R δ , defined as:

∀ s ∈ R, R δ g(s) = ∞ -∞ g(s cos δ + t sin δ, s sin δ -t cos δ)dt. ( 18
)
We note f z : (x, y) ∈ R 2 → f (x, y, z) ∈ R the restriction of f to some horizontal plane.

For τ ≥ 0 and δ ∈ (-π, π], we define the integral transform given for all s ∈ R by:

P τ,δ f (s) = τ ∞ 0 π -π R δ+ π 2 f z (s + zτ sin ϕ)dϕdz. ( 19
)
Simple calculations starting from ( 17) and ( 18) allow to show that the transforms P and C 1 are related to each other by

P τ,δ = 1 K(β, E 0 ) R δ+ π 2 C 1 α,δ,β , (20)
whenever τ = b(α, β). The transform P produces a data set included in a three dimensional space, indexed by τ > 0, δ ∈ (-π, π] and s ∈ R. The required data set may be further reduced by taking into account the identity

R δ+ π 2 C 1 α,δ+π,β = R δ+ π 2 C 1 α,δ,β , (21)
meaning that from the point of view of the transform P, the projections C 1 α,δ+π,β and C 1 α,δ,β share the same information.

Taking the Fourier transform in [START_REF] Tomitani | Analytical image reconstruction of cone-beam projections from limited-angle Compton camera data[END_REF] leads to:

P τ,δ f (σ) = 2πτ ∞ 0 R δ+ π 2 f z (σ)J 0 (2πzτ σ)dz, ( 22
)
where J 0 is the zero-order Bessel function of the first kind. Then, by the projection-slice theorem applied to the Radon projection and for d 2 = (-sin δ, cos δ, 0) we get:

P τ,δ f (σ) = 2πτ ∞ 0 1 z f z (σd 2 )J 0 (2πzτ σ)zdz. (23)
In [START_REF] Schiefeneder | The Radon transform over cones with vertices on the sphere and orthogonal axes[END_REF], f z = 0 for z close to the origin. Formulas ( 22), ( 23) are valid for all Schwartz functions f supported in the upper half space and both appearing integrals are absolutely convergent. We recognize in [START_REF] Schiefeneder | The Radon transform over cones with vertices on the sphere and orthogonal axes[END_REF] the Hankel transform of f z /z. Finally the inversion of the Hankel transform leads to

f z (σd 2 ) = 2πzσ 2 ∞ 0 P τ,δ f (σ)J 0 (2πzτ σ)dτ, ( 24
)
and to the reconstruction formula:

f (x, y, z) = 2πz π 0 ∞ 0 ∞ -∞ P τ,δ f (σ)J 0 (2πzτ σ)|σ| 3 e 2iπσ(-x sin δ+y cos δ) dσ dτ dδ. ( 25
)
Note that when b(α, β) = τ > 0, either α or β may be calculated from ( 14) as function of the other parameter through the relation:

τ cos α = √ 1 + τ 2 sin β, ( 26
) since both applications β ∈ 0, π 2 -α → b(α, β) ∈ R * + , for some given α ∈ 0, π 2 , and
α ∈ 0, π 2 -β → b(α, β) ∈ R * + , for some given β ∈ 0, π 2
-α , are one-to-one. This means that in [START_REF] Moon | Analytic inversion of a conical Radon transform arising in application of Compton cameras on the cylinder[END_REF], the projections P τ,δ can be calculated as the mean of an infinity of conical projections, one for each polar angle α or one for each Compton angle β.

Numerically, the reconstruction can be done slice by slice (following the direction orthogonal to the camera) as follows. The Radon projections of the Compton projections are filtered in the Fourier domain with a filter depending on the parameter τ , related to the cone pattern generating the projection, and on the altitude z. The filtered projections are then integrated on τ then back-projected on lines in each horizontal slice. The filtering in the Fourier domain relies on the evaluation of the Bessel function J 0 on a grid sufficiently fine to avoid aliasing as much as possible. Consequently, zeropadding of the projections is necessary implying intensive use of the memory when all projections are processed simultaneously.

The evaluation of the reconstruction formula may produce very noisy results as the reconstruction filter |σ| 3 amplifies the noise in the data while emphasizing high frequencies. A practical yet insufficient solution is the damping of the higher frequencies during the reconstruction. Some post-filtering remains necessary and may be obtained by total variation denoising.

Alternative formulations of the inverse Compton transform

In a tempered distribution sense the Bessel function σ → J 0 (2πzτ σ) is the Fourier transform of

k z,τ (s) = 1 2π 2rect s 2zτ (zτ ) 2 -s 2 , ( 27
)
where rect denotes the rectangular function. With ramp the impulse response of the ramp filter σ → |σ|, the filtered projections may be rewritten as:

∞ -∞ P τ,δ f (σ)J 0 (2πzτ σ)|σ| 3 e 2iπσv.d 2 dσ = - 1 (2π) 2 (P τ,δ k z,τ ramp) (v.d 2 ). ( 28
)
With the notation g z,τ,δ = P τ,δ k z,τ we then have

∞ -∞ P τ,δ f (σ)J 0 (2πzτ σ)|σ| 2 e 2iπσ(v.d 2 ) dδ = - 1 (2π) 2 (g z,τ,δ ramp) (v.d 2 ). ( 29
)
The function g z,τ,δ may be calculated either in the temporal or the frequency domain.

The first alternative seems preferable, because the highly oscillatory Bessel function requires a high sampling rate to limit aliasing as much as possible (see figures 5 and 6). This frequency domain sampling may be obtained by zero-padding the projections P τ,δ . Let us now put for each s ∈ R,

p z (s, δ) = z ∞ 0 g z,τ,δ (s)dτ. ( 30
)
From ( 29) we then deduce that

∞ 0 ∞ -∞ P τ,δ f (σ)J 0 (2πzτ σ)|σ| 3 e 2iπσ(v.d 2 ) dδ dτ = - 1 (2π) 2 z (p z (•, δ) ramp) (v.d 2 ),( 31
)
and finally from ( 25) we obtain

f (x, y, z) = - 1 2π R -1 ∂ 2 p z ∂s 2 (x, y), ( 32
)
where R -1 is the two-dimensional inverse Radon transform. We have thus shown that the image of the source may be reconstructed slice by slice as the inverse of the twodimensional Radon transform, from synthetic Radon projections calculated by postprocessing the Compton projections. From a numerical point of view, formula (32) has several advantages:

(i) It gives the possibility to use iterative methods combined with advanced filtering for the inversion of the Radon transform, still conserving a fast method compared to direct use of iterative methods on the raw Compton data; (ii) It has lower memory requirements when the convolution with the filter k is done in the time domain compared to Fourier domain convolution with the highlyoscillating Bessel function;

(iii) It makes use of a more precise calculation of the filtered projections g z,τ,δ which are no more subject to aliasing. An example is shown in figure 6.

In our experiments, only the first point appeared as important. The third point had little impact on the results, but this observation may be biased by the small distance we put between the source and the detector.

Semi-iterative inversion of the Compton transform

The images obtained from (32) may be denoised with a total variation filtering. However, post-filtering increases the projection errors and a better approach is to apply an iterative algorithm to the regularized problem as explained at the end of section 2.4. In this case we choose to solve a problem of type [START_REF] Lee | A compact Compton camera using scintillators for the investigation of nuclear materials[END_REF]. Following the procedure described in section 2.4, SIRT steps are alternated with total variation denoising steps completed with the Chambolle's algorithm from [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. If λ is the descent step of SIRT and µ is the regularization parameter multiplying the TV norm in [START_REF] Lee | A compact Compton camera using scintillators for the investigation of nuclear materials[END_REF], the parameter of the total variation denoising should be taken equal to the product λµ.

Analytic solution written as a second derivative

In [START_REF] Kuchment | Inversion of weighted divergent beam and cone transforms[END_REF] the back-projection is preceded by the differentiation of the projections. The derivative is calculated by finite differences and this operation induces some circular blurring in the reconstructed images.

The commutativity of the Radon transform with derivatives is a classical result (see e.g., [START_REF] Natterer | The mathematics of computerized tomography[END_REF]). For any g : R 2 → R sufficiently smooth,

R δ (∂ 1 g) = (cos δ)(R δ g) and R δ (∂ 2 g) = (sin δ)(R δ g) . and when ∆ = ∂ 11 + ∂ 22 is the Laplace operator, R δ (∆g) = (R δ g) .
We thus may rewrite [START_REF] Kuchment | Inversion of weighted divergent beam and cone transforms[END_REF] as

f (x, y, z) = - 1 2π ∆(R -1 p z )(x, y). ( 33
)
In the continuous case this formula is equivalent to [START_REF] Kuchment | Inversion of weighted divergent beam and cone transforms[END_REF] and thus to [START_REF] Moon | Analytic inversion of a conical Radon transform arising in application of Compton cameras on the cylinder[END_REF]. In the discrete case the equivalence do not hold any more and formula (33) may produce less blurry images. Numerically, the Laplacian can be calculated by finite differences. In this work we used the del2 function from MATLAB.

Description of the data

We generated Compton camera data with the Monte-Carlo simulation software MEGAlib (The Medium-Energy Gamma-ray Astronomy library) described in [START_REF] Zoglauer | MEGAlib-The medium energy gamma-ray astronomy library[END_REF]. This software is composed of several modules, and is based on the Geant4 library [START_REF] Agostinelli | Geant4-A simulation toolkit[END_REF]. The COSIMA module allows to simulate the interactions between the particles and the matter. The output consists in a list of interactions for each individual photon and for the secondary particles produced during these interactions. The ordering of the hits in the sequence is currently available only in simulations, as the time resolution of the detectors is too poor. In a real acquisition, the camera would provide a set of (not ordered) interactions and specific algorithms are required to reconstruct the sequence of hits.

The Monte-Carlo data may then be analysed with one of the two modules SIVAN and REVAN. The first one retrieves the ideal positions of the hits and the exact energies of the secondary particles. Only Doppler broadening may affect the measurements. Additionally the order of the interactions is known in this case and the vertex and axis of the cone are thus the true ones. When the data are processed with REVAN, errors are added both on positions and energies. These errors reverberate on the value of the scattering angle, on the position of the vertex and on the direction of the axis of the cone. The order of the interactions is determined based on kinematics considerations and may be erroneous. The first hit is eventually wrongly identified and leads to a miss-positioning of the vertex of the cone. The same may hold for the second hit, with consequences on the axis direction.

Influence of the errors on the angular parameters

The errors on the angular parameters of the Compton cone have important consequences on the reconstructed images regardless the algorithm. For instance in figure 7 we show iterative Maximum Likelihood Expectation-Maximization (MLEM) reconstructions from deterministic projections of three identical point-like sources. The algorithm was originally adapted to Compton camera imaging by Wilderman et al in [START_REF] Wilderman | Fast algorithm for list mode backprojection of Compton scatter camera images[END_REF] and here we use the implementation from [START_REF] Hilaire | Proton therapy monitoring by Compton imaging: influence of the large energy spectrum of the prompt-γ radiation[END_REF]. In (a) the projections are simulated with α = 2.5°a nd reconstructed as they were acquired at α = 0°. In (b), the projections for β = 10°a nd β = 15°are summed together then reconstructed as they were taken at β = 12.5°. In both cases the intensities of the sources are different after reconstruction and the blur increases as we move further from the detector. Whilst the errors on α produced in this experiment only an elongation of the source, the influence of the errors on β is more critical. Figure 7 shows that the artefacts in the vertical direction are not specific to analytic reconstruction. They are not necessarily related to an erroneous modelling of the data since in this example the model from which the data is generated is perfectly known. Elongation of the source in the direction orthogonal to the camera is usual in Compton camera imaging and is often related to the small solid angle subtended by the detector at a point from the source. As it is produced by a loss of projections, it could be assimilated to a "limited angle" artefact. In this example the camera was designed to be relatively large and close to the sources. We shall thus conclude that data binning strongly reinforces the artefacts in the vertical direction.

Compton camera simulated data

For the purpose of evaluation of the analytic methods we simulated a Compton camera composed of six Silicon scatterer layers and one absorber composed of LaBr 3 crystals. Each layer has dimensions 30 × 30 × 0.2 cm 3 and the absorber is composed of 128 2 crystals arranged in a 73 cm-side square. Each crystal has dimensions 0.5 × 0.5 × 4 cm 3 .

The scatterers are placed with the centres at altitudes z ∈ {-5, -6, . . . , -10} and the absorber at z = -22.

For the noisy data, the spatial resolution in the scatterer is about 1 mm. In the absorber, the hit is always considered at the center of a crystal. The energy errors are randomly drawn from Gaussian distributions. In the scatterer its standard deviation is 1 keV. In the absorber, the standard deviation depends on the measured energy. For instance its value is 11.4 keV at 511 keV. The non-uniform source intensity has the form of the word RADON, with a thickness of 1 cm in the third direction (see figure 8). A total number of 10 9 γparticles was emitted with energies 511 keV. After discarding not Compton and not reconstructible events (see conditions on the angles α and β in section §3.2), about 17 × 10 6 reconstructible events were selected in the ideal case and about 8.1 × 10 6 in the realistic case.

The data were grouped in a five-dimensional array. The dimensions correspond to the layer, the azimuthal angle α, the polar angle δ, the parameter τ and the coordinate s in the Radon projection of the Compton projection P τ,δ . The resulting array had 6 × 45 × 90 × 45 × 170 elements. The higher the sampling density for the angular parameters, the better the reconstructed image should be as the sampling errors produce blur and artefacts that are much more severe than in imaging modalities based on line projections (see figure 7).

Results and discussions

Volumes having 120 × 120 × 33 cubic voxels with side of 2.5 mm were reconstructed from ideal and noisy data. The resolution in the direction perpendicular to the camera is known to be worse than in slices parallel to it. We thus considered two situations where the word was placed either parallel to the camera or orthogonal to it. We also tested a reduction of the number of counts by a factor ten.

We compare three reconstruction methods: analytic stands for equation [START_REF] Kuchment | Inversion of weighted divergent beam and cone transforms[END_REF], Laplacian stands for [START_REF] Terzioglu | Compton camera imaging and the cone transform: a brief overview[END_REF] and semi-iterative stands for the algorithm described in paragraph 3.3.1.

Resolution in slices parallel to the detector

In this test we place the word-shaped source parallel to the detector with the center at altitude z = 0. The central slice is thus at 5 cm from the detector.

The first example is made with ideal data. We recall that ideal data include however some errors on energies induced by Doppler broadening and also statistical errors coming from the random generation of particles in the source. Even in this case, the noise in the reconstructed images is consequent as it may be seen in figure 9. Some denoising is necessary and we suggest that TV denoising would be a good choice in order to preserve the sharp edges of the source. The evolution of the errors for the semi-iterative method is shown in figure 11 for the central slice. The figures of merit are the cost function from equation ( 5) displayed in blue and the root mean square error displayed in red. The iterations are started from the analytic solution.

When a realistic noise is added to the data, in practice by treating the output of the simulation with the module REVAN, the results are largely degraded but the shape of the letters is almost conserved (see figure 12). It should be noted however that the size of the voxel is here of (2.5 mm) 3 which is very low when we know that the spatial Figure 13 shows the evolution of the errors for the semi-iterative method with realistic data.

Figure 13: Evolution of the errors with the number of iterations for the regularized semiiterative method. Realistic data, as produced in the output of the REVAN module.

Low statistics acquisition

The quantity of available data has a strong influence on the noise as it is shown in figure 14. The image is reconstructed from a sub-set of 1.6×10 6 ideal events taken from the previous simulation (thus about one tenth of the sample). Again, imposing low total variation to the reconstructed image helps to strongly improve the result. It may be noted in figure 15 that in this case the methods Laplacian and semi-iterative seem to outweigh the analytic method and to give sharper results.

Resolution in slices perpendicular to the detector

In this test we place the word-shaped source perpendicular to the detector. The bottom of the source is at 3 cm from the up-most scatterer. The central slice refers in this case to y = 0 and this slice is thus perpendicular to the detector, too. A blurring in the vertical direction can be observed. This blurring is usual in Compton camera imaging regardless the reconstruction method and uses to be attributed to projections loss due to the finite size of the detector. At our knowledge, the only theoretical study on the artifact produced by the data incompletness is the one from [START_REF] Palamodov | Reconstruction from cone integral transforms[END_REF], where a "plume"shaped artifact is mentionned but under the convolutional conical Radon transform model. In section 4.1 we showed that the measurement errors and the binning are additional degradation factors. The relative influence the geometry and the errors on the angular parameters may have is an important question not yet elucidated. The first test is done with ideal data and the results are displayed in figure 16. Let us mention that even ideal data are affected by the Doppler broadening that induce errors on the Compton angle. Those errors have a bell-shaped distribution with full width at half maximum of a few degrees, the exact value depending on the material and on the temperature. None of the methods presented here could retrieve the shape of the source when realistic data was used. The result is in all cases blurry and noisy and only the rough location of the word can be established. Again, the results from figure 17 show that the analytic method give the blurriest results. Central slices after TV denoising of the volume. For the semi-iterative method the result is shown after 30 iterations.

Conclusions

The purpose of this work is to push the limit of the attainable resolution in Compton camera imaging. We developed an inverse for a conical Radon transform and we tested three implementations for its numerical calculation. We demonstrated that imposing low total variation to the reconstructed volume strongly improves the result even when the level of the noise in the image is high. The results given by the three methods were very similar for ideal high statistics data. When the quality of the data drops, the Laplacian method, where the derivation is applied to the volume instead of the projections, and the semi-iterative method gave slightly better results in terms of sharpness. The best solution both in time spent in the reconstruction and quality of the image seems to be the Laplacian method. In the less favourable situation where the resolution in the vertical direction is tested with noisy data we did not succeeded to reconstruct an acceptable image of the source. The errors on the axes directions, on the half-opening angle and on the sequence of interactions, cumulated with the approximations induced by the binning, prevent reconstructing a source with small details as the one we chose for our tests. It should be noted however that we did not applied here clever data processing, consisting e.g., to select only events that have deposited in the detector a total energy close to the known initial energy of the photon. Our idea was to explore the limits of the methods, and in some situations the initial energy is unknown.

Figure 3 :

 3 Figure 3: (a) Parameterization of the Compton cone. (b) Construction of the Compton projection C α,δ,β . The Compton projection is the image drawn on the surface of the scatterer by the vertices of the cones with parameters α, δ, β that intersect the source.The intensity of the print is the integral from[START_REF] Maxim | Filtered backprojection reconstruction and redundancy in Compton camera imaging[END_REF].

Figure 4 :

 4 Figure 4: Example of Compton projections of a small spherical source. In (a) the projection is calculated deterministically following (10). In (b) a Monte Carlo simulation of the same detector-source configuration is shown. The model serving for the deterministic projection (a) only approximatelly describes the data in (b).

Figure 5 :

 5 Figure 5: Inverse Fourier transform of the discretized Bessel function (blue line) compared to k z,τ (red line).

Figure 6 :

 6 Figure 6: (a) Example of projection P τ,δ , i.e., Radon projection of the Compton projection, for a source composed of three points. (b) Its filtered versions g z,τ,δ for some arbitrary value z, with the convolution calculated in the frequency domain (blue doted line) and in the time domain (red continuous line).The oscillations due to aliasing that may be observed in the right-hand image are important, but this is not always the case. Aliasing is more visible for large values of zτ . Since they are attenuated by the usual smoothing of projections in the reconstruction process, the influence of the filtering domain may finally have little impact on the images.

Figure 7 :

 7 Figure 7: Image reconstructed with the Maximum Likelihood Expectation-Maximization algorithm (MLEM) from projections with (a) an error of 2.5°on the axis inclination α and (b) an error of 2.5°on the Compton scattering angle β (see text for details). Twenty iterations of the MLEM algorithm were runned.

Figure 8 :

 8 Figure 8: The reference image -one slice from a 1 cm thick source representing the word RADON.

Figure 9 :

 9 Figure 9: Central slices from the word-shaped source placed parallel to the detector. The data are processed with SIVAN, meaning that only mild noise is considered (see text for details). The analytic and Laplacian methods give similar and noisy results.

Figure 10

 10 Figure 10 shows on the left and center the central cuts of the same volumes as in figure 9 but after TV denoising. On the right is shown the TV-regularized SIRT reconstruction. The results for the three methods are very similar although the analytic reconstruction is slightly blurry (see letters R and O).

Figure 10 :

 10 Figure 10: The source is placed parallel to the detector and the collected events are ideal except the Doppler broadening. Central slices from the word-shaped source after TV denoising. For the semi-iterative method the result is shown after 30 iterations.

Figure 11 :

 11 Figure 11: Evolution of the errors with the number of iterations for the regularized semi-iterative method. Ideal data, as produced in the output of the SIVAN module.

Figure 12 :

 12 Figure 12: Realistic data and TV denoising for the word-shaped source placed parallel to the detector. For the semi-iterative method the result is shown after 30 iterations.

Figure 14 :

 14 Figure 14: Low statistics reconstruction from ideal data, with the analytic formula.

Figure 15 :

 15 Figure 15: Low statistics acquisition. Central slices from the word-shaped source placed parallel to the detector after TV denoising. The data are processed with SIVAN. For the semi-iterative method the result is shown after 20 iterations.

  (a) analytic (µ = 0.15) (b) Laplacian (µ = 0.04) (c) semi-iterative (µ = 0.01)

Figure 16 :

 16 Figure 16: Word perpendicular to the camera and ideal data from SIVAN output. Central slices after TV denoising of the volume. For the semi-iterative method the result is shown after 20 iterations.

  (a) analytic (µ = 0.15) (b) Laplacian (µ = 0.04) (c) semi-iterative (µ = 0.01)

Figure 17 :

 17 Figure 17: Word perpendicular to the camera and realistic data from REVAN output. Central slices after TV denoising of the volume. For the semi-iterative method the result is shown after 30 iterations.

Acknowledgements

The author acknowledges the financial support of the French National Research Agency through the ANR project 3DCLEAN (ANR-15-CE09-0009) and the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the ANR. The author would like to thank the anonymous reviewers for careful reading of the manuscript and thoughtful comments.


Andre Manoel

Florent Krzakala

Bertrand Thirion

Gaël Varoquaux

Lenka Zdeborová

Approximate message-passing for convex optimization with non-separable penalties

published or not. The documents may come

Introduction

Proximal algorithms such as FISTA [1] and ADMM [2,3] form these days a large part of the standard convex optimization toolbox, mostly due to the fact that they provide different sorts of theoretical guarantees, such as the convergence to a global minimum of the problem. However, their use is not so interesting in some cases. In particular, for non-separable penalties such as total variation (TV), FISTA requires an inner loop for evaluating the proximal operator, which takes longer and longer to run as the outer loop iterations go by. While variable-splitting algorithms such as ADMM provide a way around this, they are typically very sensitive to the choice of stepsize, and determining an optimal stepsize remains an open problem.

An alternative comes from approximate message-passing (AMP) algorithms [4,[START_REF] Rangan | Generalized approximate message passing for estimation with random linear mixing[END_REF][START_REF] Krzakala | [END_REF], extensions to the belief propagation algorithm [START_REF] Pearl | Probabilistic reasoning in intelligent systems: networks of plausible inference[END_REF] that are capable of dealing with highdimensional inference problems. They have shown remarkable success in applications to compressed sensing, denoising and sparse regression [START_REF] Vila | [END_REF][START_REF] Ziniel | Binary linear classification and feature selection via generalized approximate message passing Information Sciences and Systems[END_REF][START_REF] Guo | [END_REF]11], but despite that remain disfavored compared to classic optimization approaches, most prominently due to convergence issues [START_REF] Rangan | On the convergence of approximate message passing with arbitrary matrices[END_REF][START_REF] Caltagirone | On convergence of approximate message passing[END_REF]. The recently proposed vector approximate message-passing (VAMP) algorithm [START_REF] Rangan | Vector approximate message passing[END_REF] offers a promising, more robust alternative to AMP. For a given class of statistical models, notably generalized linear models (GLMs), different formulations of the VAMP algorithm are capable of approximating both the maximum a posteriori (MAP) and minimum mean-squared error (MMSE) estimates. Moreover, these estimates are conjectured to be asymptotically Bayes-optimal for certain random matrices, and proven to be in specific cases [START_REF] Reeves | The replica-symmetric prediction for compressed sensing with Gaussian matrices[END_REF][START_REF] Barbier | [END_REF][START_REF] Barbier | The Mutual Information in Random Linear Estimation Beyond i[END_REF].

Here, we are interested in the evaluation of the MAP estimate, which is also the solution to a given optimization problem. In particular, we would like to determine whether VAMPlike algorithms can provide a fast alternative to classic optimization approaches in solving a given class of convex optimization problems. VAMP has been typically considered for the minimization of functions consisting of a quadratic loss and a separable penalty, a setting in which coordinate descent algorithms are usually much faster than the alternatives, as discussed in Section 2. We thus look at the more interesting (and challenging) case of penalties that are non-separable-more specifically, separable on disjoint subsets of a linear transformation of the variable being optimized, a notable example being total variation (TV) [START_REF] Beck | [END_REF]19]. For such penalties, proximal algorithms [20] such as FISTA and ADMM are considered state-of-the-art.

While VAMP has already been combined with non-separable denoisers [21,[START_REF] Schniter | Denoising based vector approximate message passing[END_REF], there have been no attempts so far to modify the algorithm to better fit specific types of nonseparable penalties. Moreover, evidence that VAMP has a good performance in real-world applications is still limited. The current work addresses these two issues. In Section 3, we rederive VAMP for a specific class of non-separable penalties, by considering the expectationconsistent (EC) approximation [START_REF] Opper | [END_REF] not on the distribution of the variable of interest, but on that of a linear transform. Moreover, in Section 4, experiments with total variation penalties are performed on benchmark datasets, revealing that the adaptation we propose is competitive with state-of-the-art algorithms.

Motivation and background

Approximate message-passing

In what follows, we consider a typical regression problem and denote by y ∈ χ n the response vector (which can represent e.g. labels, χ = {0, 1}, or continuous value χ = R) and by A ∈ R n×p the feature matrix. Then, given a regularization parameter λ ∈ R and a convex function f , we look at solving x = arg min

x 1 2 y -Ax 2 2 + λ j f (x j ), (1)
The approximate message-passing (AMP) algorithm in [4] aims at solving (1) by means of the following iteration

x t+1 = η λσ t (x t + A T z t ), (2)
z t = y -Ax t + 1 α z t-1 ∇η λσ t (x t + A T z t ) , (3)
where α = n/p, σ t is an empirical estimate on the variance of x, a = 1 na na i=1 a i is the empirical average over the elements of a ∈ R na , and

η λσ t (v) = prox λσ t f (v) = arg min x λσ t f (x) + 1 2 (x -v) 2 , ( 4
)
is a function applied elementwise. In practice, the variance estimate is updated using [START_REF] Donoho | Message passing algorithms for compressed sensing: I. motivation and construction IEEE Information Theory Workshop[END_REF] σ t+1 = 1 + σ t ∇η λσ t /α. AMP's form is close to that of a proximal gradient descent [20], being in fact very similar to the iterative soft thresholding algorithm (ISTA), with, however, an additional term. The extra term in (3), known as Onsager reaction term, comes from a second order correction derived using statistical physics techniques [START_REF] Thouless | [END_REF]26], and has been shown to improve the iteration in synthetic data experiments. An instance of this can be seen in the example of subsection 2.1.1.

An important assumption is made in deriving AMP: that the entries of A are i.i.d., have zero mean, and variance that scales as 1/n. The iteration typically fails to converge whenever this assumption is not met [START_REF] Rangan | On the convergence of approximate message passing with arbitrary matrices[END_REF][START_REF] Caltagirone | On convergence of approximate message passing[END_REF]. A series of works have tried to improve the convergence properties of AMP by means of heuristic modifications [START_REF] Vila | Adaptive damping and mean removal for the generalized approximate message passing algorithm[END_REF][START_REF] Manoel | Swept approximate message passing for sparse estimation International Conference on Machine Learning[END_REF], typically resulting in slower algorithms that provide few theoretical guarantees.

A promising alternative is the vector AMP (VAMP) algorithm, recently introduced by [START_REF] Rangan | Vector approximate message passing[END_REF][START_REF] Ma | [END_REF] as a more robust message-passing scheme. It is shown to converge, in the large n limit, for a wider class of matrices A -namely that of rotationally-invariant random matrices. As AMP, VAMP provides an iteration which upon convergence returns a solution to (1); it uses, however, a different framework than AMP, known as expectation-consistent (EC) approximation [START_REF] Opper | [END_REF], which is briefly described in Appendix A. For problem (1), the VAMP iteration reads

x t = (A T A + ρ t I n ) -1 (A T y + u t ), z t = η λσ t x /(1-σ t x ρ t ) x t -σ t x u t 1 -σ t x ρ t , ( 5
)
u t+1 = u t + (z t /σ t z -x t /σ t x ), ρ t+1 = ρ t + (1/σ t z -1/σ t x ). ( 6
)
While AMP is close to ISTA, the iteration above resembles that of another popular optimization algorithm, the alternating directions method of multipliers (ADMM) [2]. This connection is discussed in subsection 2.2. Comparison between AMP and VAMP for sparse regression on synthetic data, with n = 600, p = 2000 and ρ = 0.1, σ = 10 -5 , λ = 1. Objectives are relative to that of scikit-learn's coordinate descent [30], which runs in a much shorter time (around 5 • 10 -4 s). Left: for Gaussian i.i.d. matrices, VAMP performs worse than AMP, both in terms of convergence rate and preprocessing time. Right: for a matrix which is not Gaussian i.i.d., but instead the product of two Gaussian i.i.d. matrices A = U V T , AMP quickly diverges, as shown in the inset, while VAMP keeps providing a good performance. In this example, both U and V have r = 600 columns.

The σ t

x and σ t z are estimates of the variance of x and of the split variable z, and are obtained from

σ t x = 1 n Tr(A T A + ρ t I n ) -1 , σ t z = σ t x 1 -σ t x ρ t ∇η λσ t x /(1-σ t x ρ t ) x t -σ t x u t 1 -σ t x ρ t . ( 7
)
Iterations ( 5) and ( 7) can be seen as evaluating the means and variances of two distinct distributions, commonly referred to as tilted distributions in the context of EC. By forcing these distributions to concentrate around their modes, the means become proximal operators, whereas the variances come from the curvature of Moreau envelopes of the loss and the penalty. This is detailed in Appendix A. The matrix inversion in ( 5) and ( 7) can be quite costly; however, the iteration can be modified to use the SVD decomposition of A = U ΣV [START_REF] Rangan | Vector approximate message passing[END_REF], in which case

σ t x = 1 n n i=1 1 Σ 2 ii +ρ t .
In particular, if n p, the Woodbury formula can be combined with the SVD to provide an efficient iteration, which is used in our experiments (see Appendix C for more details).

Unless otherwise mentioned, we shall always assume that x 0 , z 0 and u 0 are initialized to all-zero vectors, while ρ 0 = 1 and σ 0 x = σ 0 z = 1/2.

Example: 1 regularization on synthetic data

In order to compare the AMP iteration ( 2)-( 3) to the VAMP one, ( 5)-( 6), we perform a simple experiment on synthetic sparse data. We first sample x from a Bernoulli-Gaussian distribution, P X (x) = j [ρ N (x j ; 0, 1) + (1 -ρ) δ(x j )], then generate a matrix A by doing A ij ∼ N (A ij ; 0, 1/n), to finally obtain measurements as y ∼ N (y; Ax, σ 2 I n ). We try to recover x by using a 1 penalty, f

(x j ) = |x j |. In that case ‡, η λσ t (x) = 1 -λσ t |x| +
x, also known as the soft thresholding operator. In Figure 1 (left), we study how fast each of the schemes performs this task. We look at three different schemes: iterative soft thresholding (ISTA) [31], the original AMP iteration ( 2)-( 3), and the VAMP iteration ( 5)- [START_REF] Pearl | Probabilistic reasoning in intelligent systems: networks of plausible inference[END_REF]. Compared to AMP and ISTA, VAMP has a larger preprocessing time, since it relies on the eigendecomposition of AA T , as detailed in Appendix C. Moreover, in this artificial setting, its convergence rate is smaller than that of AMP. Yet, such Gaussian matrices are precisely those assumed in the derivation of AMP. If a different matrix is used, say A = U V T with both U and V Gaussian, then the AMP iteration quickly diverges, whereas VAMP continues to work (see Figure 1, right). One verifies in practice that VAMP converges for a broad class of matrices, and is actually proven to converge, in the large n limit, for right orthogonally-invariant random matrices [START_REF] Rangan | Vector approximate message passing[END_REF]. This motivates us to perform more detailed experiments with VAMP, as well as to adapt it to our needs.

Also notably, all schemes used in Figure 1 are much slower than coordinate descent, which is currently considered the state-of-the-art for optimization with separable penalties. For this reason, we shall consider the more interesting case of non-separable penalties in Section 3.

Variable splitting methods

We now turn to a class of optimization algorithms known as variable splitting methods. The main very basic idea is to minimize a sum of two functions such as (1) by writing each of them as a function of a different variable, and then enforcing these variables to assume the same value. This is done by means of the following augmented Lagrangian

L(x, z, u) = 1 2 y -Ax 2 2 + λ j f (z j ) -u T (x -z) + ρ 2 x -z 2 2 ρ 2 x-z-u ρ 2 2 -1 2ρ u 2 2 , ( 8
)
One should minimize the expression above with respect to x and z, and maximize its dual L * = min x,z L(x, z, u) with respect to the Lagrange multiplier u, see e.g. [START_REF] Boyd | Convex optimization[END_REF][START_REF] Nocedal | [END_REF]. The alternating direction method of multipliers (ADMM) [2,3] accomplishes this by means of the following iteration

x t = arg min x L(x, z t , u t ) = prox L/ρ (z t + u t /ρ), ( 9
)
z t = arg min z L(x t , z, u t ) = prox R/ρ (x t -u t /ρ), ( 10
)
u t+1 = u t + ρ(z t -x t ). ( 11
)
where we have defined L(x) = 1 2 y -Ax 2 2 and R(z) = λ j f (z j ). Equations ( 9) and (10) perform two sequential gradient descent steps in L, first for x then z, whereas equation ( 11) is effectively a gradient ascent step in L * (u), with stepsize set so as to satisfy the optimality conditions ∇ x,z L = 0.

The strategy employed by ADMM is empirically known as Douglas-Rachford splitting. Different strategies, with the same fixed points but slightly different updates, can be employed. For instance, the Peaceman-Rachford splitting (PRS) consists in the following equations [START_REF] Gabay | Applications of the method of multipliers to variational inequalities[END_REF]35]

x t = prox L/ρ ( u t /ρ), ( 12
)
z t = prox R/ρ (2x t -u t /ρ), ( 13
)
u t+1 = u t + 2ρ(z t -x t ). ( 14
)
While PRS is known to often lead to lower objectives faster than ADMM [START_REF] Gabay | Applications of the method of multipliers to variational inequalities[END_REF]36], it is provably convergent under stronger assumptions, and in practice the objective often does not decrease monotonically.

Writing the proximal operators explicitly in the PRS iteration ( 12) and ( 13) yields

x t = (A T A + ρI n ) -1 (A T y + u t ), z t = η λ/ρ (2x t -u t /ρ). ( 15
)
Interestingly, setting σ t x = σ t z = 1 2ρ in the VAMP iteration ( 5)-( 6) leads to the PRS iteration above; VAMP can thus be seen as a PRS variant where the stepsize is set adaptively. This remark has been first made in [START_REF] Fletcher | Expectation consistent approximate inference: Generalizations and convergence[END_REF].

We thus see that VAMP is more robust than AMP, converging in situations where AMP does not, as shown in Figure 1. It also has a significant advantage over variable-splitting methods, since no stepsize needs to be set. While neither VAMP nor variable-splitting methods are competitive with the state-of-the-art in the particular setting of separable penalties, we now turn to non-sperable penalties, which are more challenging regarding efficiency and convergence control [START_REF] Dohmatob | Benchmarking solvers for tv-1 leastsquares and logistic regression in brain imaging International Workshop on Pattern Recognition in Neuroimaging[END_REF], yet very important in practical applications.

Non-separable penalties

Consider a penalty which is separable not on x but on disjoint subsets of a linear transformation Kx, with K ∈ R r×n x = arg min

x 1 2 y -Ax 2 2 + λ r k=1 f (Kx) k , (16)
for any convex function f . Note that k does not necessarily index a single component of Kx, but more generally a subset of components; moreover, each component belongs to a single subset. An example is the so-called total variation (TV) penalty [START_REF] Beck | [END_REF]19], used to enforce spatial regularity -typically flat regions separated by sharp transitions

R TV (x) = ij (x i,j -x i,j-1 ) 2 + (x i,j -x i-1,j ) 2 , ( 17
)
where we have indexed the vector x by two components i, j such that x k = x i,j for i = k/L and j = k mod L. By defining a gradient operator ∇

(∇x) k = x i,j+1 -x i,j x i+1,j -x i,j , (18)
the 2D TV penalty can be conveniently written as R TV = 2n k=1 (∇x) k 2 ≡ ∇x 2,1 . Different boundary conditions can be used; in our experiments, we assume periodic boundary conditions.

While we define it above for two dimensions, the TV operator is easily generalizated to lattices in three or more dimensions. Notice moreover that we are dealing with the isotropic or rotation-invariant TV operator, instead of the anisotropic one ∇x 1,1 .

Variable splitting algorithms can be easily adapted to that case, by using a split variable z = Kx; the augmented Lagrangian (8) becomes

L(x, z, u) = 1 2 y -Ax 2 2 + λ r k=1 f (z k ) -u T (Kx -z) + ρ 2 Kx -z 2 2 . ( 19
)
and the PRS iteration is now written

x t = (A T A + ρK T K) -1 (A T y + K T u t ), ( 20
)
z t = η λ/ρ (2Kx t -u t /ρ), ( 21
)
u t+1 = u t + 2ρ(z t -Kx t ). (22)
In particular for TV, K = ∇ and the η function is given by the group soft thresholding operator

η λ/ρ (v) = 1 - λ/ρ v 2 + v. ( 23
)

VAMP for non-separable penalties

The following VAMP-like iteration, which we derive in Appendix B, can be used in order to solve problem 16

x t = (A T A + ρ t K T K) -1 (A T y + K T u t ), σ t x = 1 n Tr K(A T A + ρ t K T K) -1 K T , ( 24
)
z t = η λσ t x 1-σ t x ρ t Kx t -σ t x u t 1 -σ t x ρ t , σ t z = σ t x 1 -σ t x ρ t ∇η λσ t x 1-σ t x ρ t Kx t -σ t x u t 1 -σ t x ρ t , ( 25
)
u t+1 = u t + (z t /σ t z -Kx t /σ t x ), ρ t+1 = ρ t + (1/σ t z -1/σ t x ). ( 26
)
Upon convergence, it is verified that z = Kx and σ x = σ z . Once again, if the variances remain fixed throughout the iteration, σ t x = σ t z = 1 2ρ , the PRS iteration ( 20)-( 22) is recovered. Moreover, the Woodbury formula can also be applied, leading to a faster iteration whenever n p. The iteration is particularly fast for the TV matrix K = ∇, for which K T K = ∆ is a circulant matrix and thus diagonal in the Fourier basis, see Appendix C.1.

While we consider convex f functions, we never use this assumption explicitly. In other words, the iteration above could, in principle, also be used when f is non-convex.

Enforcing contractivity

A well-known issue with PRS is that, differently from ADMM, the iteration mapping might not be contractive [START_REF] He | [END_REF]36], which might lead to lack of convergence. The iteration we propose, ( 24)-( 26), suffers from the same issue. A way of enforcing contractivity, done for PRS in [36], is by including a underdetermined relaxation factor 0 < γ < 1 in the update of the Lagrange multiplier u, that is

u t+1 = u t + γ(z t /σ t z -Kx t /σ t x ). ( 27
)
For consistency, we also include the same parameter in the update of the stepsize ρ

ρ t+1 = ρ t + γ(1/σ t z -1/σ t x ). ( 28
)
While this modification improves the convergence properties of the algorithm, it also makes us need to specify an extra parameter γ. This parameter could in principle be set adaptively, as it has been done for other AMP algorithms [START_REF] Vila | Adaptive damping and mean removal for the generalized approximate message passing algorithm[END_REF]. In our experiments, however, we leave this parameter fixed at a small enough value chosen arbitrarily, γ = 0.6.

Numerical experiments

In this section we solve problem ( 16) for a 2D total variation penalty, K = ∇, and different choices of {y, A}. We compare the iteration we propose, ( 24)-( 26), to three other standard algorithms: ADMM with an adaptive stepsize [START_REF] Xu | Adaptive ADMM with spectral penalty parameter selection International Conference on Artificial Intelligence and Statistics[END_REF]; a FISTA [1] adaptation known as FAASTA [START_REF] Varoquaux | [END_REF], as implemented in § the Nilearn package [START_REF] Abraham | [END_REF]; and the Peaceman-Rachford splitting, recovered by leaving variances fixed in our own implementation.

We work in the high-dimensional regime with n p and thus use the numerical tricks described in Appendix C, combining the Woodbury formula to the eigendecompositions of AA T and K T K to speed up the evaluation of the proximal step in VAMP, PRS and ADMM. In all experiments, we set the relaxation factor to γ = 0.6 in VAMP, and to γ = 0.95 in PRS.

All experiments have been performed on a machine with a 12-core Intel i7-8700K processor and 32GB of RAM, using Numpy and a multi-CPU implementation of BLAS (OpenBLAS).

Classification in task fMRI

In our first numerical experiment, we approach a multi-class classification problem, namely a task fMRI one: each row of A ∈ R n×p is given by an fMRI image, while the entries of y ∈ R n tell us what task was being performed by the subject at the time of acquisition. We use the Haxby dataset [43], which is commonly used for benchmarks and contains n = 1452 samples of p = 40 × 64 × 64 = 163840 features each. The classes refer to what image the subject was observing at the time of acquisition. We perform one vs. all classification for three of these classes: "face", "house" and "chair", each associated to 108 of the 1452 samples. While p = 163840, many of the voxels are zero-valued and thus a mask can be applied, leading to a new feature matrix with only 39912 columns. One can pass this smaller matrix to FISTA instead of the full one, as it is capable of dealing with such masks. However, the remaining methods require the full matrix, at least in the way we have implemented them.

Two of the parameters are fixed: the regularization λ = 0.1 and the relaxation γ, which is set to 0.6 for VAMP and to 0.95 for PRS. In this setting, VAMP achieves smaller values of the objective faster than all other methods and for all the three classes considered, as shown in Figure 4. The weight map obtained in "face" classification is displayed in Figure 3.

In Figure 3 (right), we study how the performance of our scheme changes in this experiment as a function of γ. We can see the performance does not change considerably as long as γ is set small enough. While this has been often observed by us, further experiments must be performed to ensure such behavior is indeed recurrent. If that is the case, it should be straightforward to devise an adaptive scheme for γ.

Reconstruction in tomography

A second experiment consists in reconstructing a two-dimensional image from a number of tomographic projections. We use a Shepp-Logan phantom of dimensions 200x200, so that p = 40000, and perform the reconstruction for 10 (n = 2000), 20 (n = 4000) and 50 (n = 10000) equally spaced projections. Moreover, we add a 1% SNR noise to the projections. In other words, the tomographic projections are generated as y = Ax + N (0, σ 2 ), where A is the two-dimensional Radon transform and σ 2 = 0.01 • Ax 2 2 /n. In our experiments, we explicitly generate the Radon matrix, that is, we do not use its operator form.

The regularization parameter is set to λ = 1.0, and the relaxation parameter to γ = 0.6. While for the reconstruction using 20 and 50 projections VAMP is faster than Since there are multiple classes, we perform one vs. all classification for three specific ones: "face", "house" and "chair". All schemes seek to solve the same convex optimization problem ( 16) with a TV penalty, K = ∇, and λ = 0.1. For all three classes considered, VAMP achieves a lower objective faster than all other schemes. VAMP, AADMM PR have a late start due to preprocessing time, during which the eigendecomposition of A T A is evaluated. A TV penalty is used with λ = 1.0. For 10 and 20 projections, VAMP and AADMM are competitive with each other and both faster than FISTA; for 50 projections, VAMP seems to be considerably faster. Meanwhile, PRS performance is highly dependent on the choice of the stepsize ρ. Note that preprocessing time is not discounted for any of the approaches, and that a larger number of projections (i.e. larger n) incurs bigger preprocessing times. the remaining algorithms, for 10 projections it is beaten by adaptive-step ADMM, being, however, still competitive with FISTA. Finally, we show the performance of Peaceman-Rachford for different values of stepsize ρ; while it might eventually perform better than all the other algorithms, its behaviour is highly dependant on the choice of stepsize.

Conclusion

We have proposed a VAMP-like iteration for optimizing functions consisting of linear losses and a given class of non-separable penalties. As the usual VAMP iteration, it can be seen as a Peaceman-Rachford splitting with adaptive stepsizes. Initial experiments reveal that the resulting algorithm is competitive with, and faster than other typically employed optimization schemes, namely ADMM and FISTA. While further experiments should be performed in order to assess the superiority of the algorithm in more general settings, we believe that the results displayed demonstrate its potential for real-world applications. Our analysis also elucidates the connection between message-passing algorithms and other known optimization schemes, expanding on the remark of [START_REF] Fletcher | Expectation consistent approximate inference: Generalizations and convergence[END_REF]. In particular, the Peaceman-Rachford splitting is recovered if variances/stepsizes are kept fixed. Further exploring this connection should lead to improved inference and optimization algorithms, and hopefully also to convergence proofs, which are typically difficult to obtain for messagepassing algorithms.

There are other EC/VAMP adaptations able to deal with constraints of the form z = Kx, but these are typically used on the loss, not the penalty [START_REF] Schniter | Vector approximate message passing for the generalized linear model 50th[END_REF][START_REF] He | Generalized expectation consistent signal recovery for nonlinear measurements[END_REF]. They could be integrated into our framework to study losses other than quadratic. An alternative would be using these adaptations on the penalty instead, by artificially introducing new effective samples that mimic the penalty as an additional loss term; this has been explored for AMP in [START_REF] Borgerding | Generalized approximate message passing for cosparse analysis compressive sensing[END_REF][START_REF] Barbier | Scampi: a robust approximate message-passing framework for compressive imaging[END_REF], and is not as straightforward as the construction we use in this paper.

Yet another perspective is adapting the algorithm for online inference, using the scheme described in [START_REF] Manoel | Streaming bayesian inference: theoretical limits and mini-batch approximate message-passing Allerton Conference on Communication[END_REF]. Interestingly, this is similar to what is already done to adapt ADMM to the online setting [START_REF] Suzuki | Dual averaging and proximal gradient descent for online alternating direction multiplier method[END_REF]. Finally, we leave for future work comparing our approach to other recent promising alternatives, in particular that of [START_REF] Kamilov | [END_REF].

and then insert the eigendecomposition

x t = u t ρ t + A T U U T y -A u t ρ t (d + ρ t ) , (C.2)
where we have used the symbol to denote the Hadamard or elementwise division between two vectors. Note that evaluating the expression above consists only in performing matrixvector products and other vector-vector operations, and can thus be performed in O(np). Finally, iteration [START_REF] Pearl | Probabilistic reasoning in intelligent systems: networks of plausible inference[END_REF] can be replaced by

σ t x = 1 n n i=1 1 d i + ρ t .
(C.3)

Appendix C.1. Further simplifications for TV penalty

For the non-separable penalty ( 16), the same trick can be applied to yield

x t = (K T K) -1 u t ρ t + (K T K) -1 A T U U T (y -A(K T K) -1 u t ρ t (d + ρ t ) , (C.4)
and, by further using the eigendecomposition of

K T K = F T ΛF x t = F T Λ -1    F u t ρ t + ÃT U U T (y -ÃΛ -1 F u t ρ t (d + ρ t )    , (C.5)
for ÃT = F A T . In particular for the TV penalty, K = ∇ is the gradient operator and K T K = ∆ is the Laplacian, which being a circulant matrix has F given by the Fourier transform. We can then replace F by the FFT, so that the expression above can be evaluated almost as efficiently as (C.2). The iteration on σ x is written in this case as

σ t x = 1 ρ t 1 d - 1 nd nd i=1 d i d i + ρ t = 1 nd nd i=1 1 d i + ρ t - d -1 d 1 ρ t , (C.6)
where d is the dimension of the TV operator, i.e. K = ∇ is a nd × n matrix.

Figure 1 .

 1 Figure 1.Comparison between AMP and VAMP for sparse regression on synthetic data, with n = 600, p = 2000 and ρ = 0.1, σ = 10 -5 , λ = 1. Objectives are relative to that of scikit-learn's coordinate descent[30], which runs in a much shorter time (around 5 • 10 -4 s). Left: for Gaussian i.i.d. matrices, VAMP performs worse than AMP, both in terms of convergence rate and preprocessing time. Right: for a matrix which is not Gaussian i.i.d., but instead the product of two Gaussian i.i.d. matrices A = U V T , AMP quickly diverges, as shown in the inset, while VAMP keeps providing a good performance. In this example, both U and V have r = 600 columns.

Figure 2 .

 2 Figure 2. Comparison between VAMP and PRS algorithms for sparse regression on synthetic data, with n = 1200, p = 2000 and ρ = 0.2, σ 2 = 10 -5 , λ = 1. Left: PRS behavior is highly dependent on the value chosen for the stepsize ρ, whereas for VAMP no parameter needs to be set. Center: effective stepsize (top) and variance estimates (bottom) provided by VAMP. Right: difference between PRS and VAMP objectives after 250 iterations.

Figure 3 .

 3 Figure 3. Left: weight map obtained by VAMP in classifying "face" vs. all. Intensities are thresholded at 0.002. Right: objective in classifying "face" vs. all at different numbers of iterations, for different choices of relaxation factor. If the relaxation factor is set small enough, there is no noticeable difference in performance.

Figure 4 .

 4 Figure 4. Classification in task fMRI using the Haxby dataset, n = 1452 and p = 136840.Since there are multiple classes, we perform one vs. all classification for three specific ones: "face", "house" and "chair". All schemes seek to solve the same convex optimization problem (16) with a TV penalty, K = ∇, and λ = 0.1. For all three classes considered, VAMP achieves a lower objective faster than all other schemes. VAMP, AADMM PR have a late start due to preprocessing time, during which the eigendecomposition of A T A is evaluated.

Figure 5 .

 5 Figure 5. Reconstruction of a Sheep-Logan phantom of size 200x200 (p = 40000) using 10 (left), 20 (center) and 50 (right) equally spaced noisy projections (n = 2000, 4000 and 10000, respectively).A TV penalty is used with λ = 1.0. For 10 and 20 projections, VAMP and AADMM are competitive with each other and both faster than FISTA; for 50 projections, VAMP seems to be considerably faster. Meanwhile, PRS performance is highly dependent on the choice of the stepsize ρ. Note that preprocessing time is not discounted for any of the approaches, and that a larger number of projections (i.e. larger n) incurs bigger preprocessing times.

Figure 6 .

 6 Figure 6. Estimates obtained by the different approaches 10 seconds after preprocessing is finished, using 50 noisy projections.

‡ We denote (•) + = max(•, 0).

§ We have performed a small modification to this implementation, so as to allow periodic boundary conditions.

Acknowledgments

This work was supported by LabEx DigiCosme and the DRF Impulsion program. We also acknowledge the support from the ERC under the European Unions FP7 Grant Agreement 307087-SPARCS and the European Union's Horizon 2020 Research and Innovation Program 714608-SMiLe, as well as by the French Agence Nationale de la Recherche under grants ANR-17-CE23-0023-01 PAIL and ANR-17-CE23-0011 Fast-Big.

Appendix A. Expectation-consistent approximation

The expectation-consistent approximation [START_REF] Opper | [END_REF] is defined for a probability distribution of the form P (x) = 1 Z P (x) P r (x). (A.1)

It seeks to approximate (minus) the log-partition functionlog Z = -log dx P (x) P r (x) by means of the following objective, denoted the expectation-consistent free energy

If there are no restrictions on Q, then F is minimized for Q = P , Q r = P r . What we want however is to have a tractable approximation, so the Q are typically replaced by distributions in the exponential family. We use a (unnormalized) Gaussian distribution, written in canonical form as

In our case, P (x) ∝ e -1 2 y-Ax 2 2 and P r (x) = p j=1 e -λf (x j ) . Combined to the explicit form of Q above this gives

There are different ways of looking for local minima of F [START_REF] Opper | [END_REF]; the vector AMP iteration [START_REF] Rangan | Vector approximate message passing[END_REF] is one of them. It consists in the following fixed-point iteration

where we denote by E t ,r the expectation w.r.t. the tilted distributions Qt ,r (x) ∝ P ,r (x)Q t ,r (x), and by Var t r, (x) the variance of these distributions. In particular

and, by defining z(u, ρ) = dxe -λf (x) e -1 2 ρx 2 +ux ,

Note that upon convergence one achieves the stationary conditions ∇F = 0

In other words, the first two moments of Q and Qr must be consistent with those of the approximating posterior Q(x) ∝ Q (x)Q r (x).

Appendix A.1. MAP estimate

In order to go from the equations above to ( 5)-( 7), we perform two additional steps. First, note that one does not need to work with both and r variables; by denoting u ≡ u r , we rewrite the iteration as

and analogously for ρ, thus obtaining [START_REF] Krzakala | [END_REF]. Moreover, since we are interested in the MAP estimate, we use the following trick: we consider instead a distribution P (x) ∝ [P (x) P r (x)] β and take the limit β → ∞, thus making the distribution concentrate around its mode; we then look at 1 β log Z. By rescaling u ,r → βu ,r and ρ ,r → βρ ,r , we ensure that most equations remain the same. The only difference is in the second integral in F, which can be performed using the Laplace method to yield

, we get using (A.7) that

and also that Var t r (x) = ρ -1 ∇η λf /ρ t (u /ρ ) , leading to ( 5) and [START_REF] Pearl | Probabilistic reasoning in intelligent systems: networks of plausible inference[END_REF].

Appendix B. Adapting the EC approximation to TV penalties

In order to adapt iteration ( 5)-( 6) to total variation penalties, we look at the expectationconsistent (EC) approximation [START_REF] Opper | [END_REF] in closer detail. We first define the following distribution

Problem ( 16) is then recovered when considering the MAP estimator for x, given by the mode of P (x). Next we introduce a variable z = Kx, and rewrite the distribution above as

so that the marginal distribution of z becomes

We apply the EC scheme described above for P (z), by setting

thus yielding the following expectation-consistent free energy

We once again parametrize these distributions by canonical-form Gaussians, Q ,r (z; u ,r , ρ ,r ) = e -ρ ,r 2 z T z+u T ,r z , leading to the fixed-point iteration ( 24)- (26).

Appendix C. VAMP iteration for n p

In the n p setting, iteration (5) can be made efficient by combining the Woodbury formula with the eigendecomposition of AA T = U diag(d)U T . We first use the Woodbury formula to replace the iteration on x by