Paul Escande
email: paul.escande@univ-amu.fr

Pierre Weiss
email: pierre.armand.weiss@gmail.com

Fast Wavelet Decomposition of Linear Operators through Product-Convolution Expansions

published or not. The documents may come

Introduction

The efficient computation of linear operators and their inverses is paramount for nearly any scientific computing problem. A large number of numerical approaches have been developed over the years, such as low rank decompositions, product-convolution expansions [START_REF] Escande | Approximation of integral operators using productconvolution expansions[END_REF], hierarchical matrices [START_REF] Hackbusch | Multi-grid methods and applications[END_REF] or wavelet decompositions [START_REF] Meyer | Wavelets, Calderón-Zygmund and Multilinear Operators[END_REF][START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF]. They can yield huge speed-ups for the practical resolution of problems ranging from partial differential equations [START_REF] Candes | Curvelets and fourier integral operators[END_REF][START_REF] Schneider | Wavelet methods in computational fluid dynamics[END_REF] to inverse problems [START_REF] Escande | Sparse wavelet representations of spatially varying blurring operators[END_REF]. We can also expect these ideas to play a growing role in the design of efficient neural networks [START_REF] Fan | Bcr-net: A neural network based on the nonstandard wavelet form[END_REF].

A serious hindrance to their popularization is however the high cost to perform the decomposition in a large scale setting. For a square matrix in R N ×N , the cost of a naive decomposition is typically O(N 3 ), which is incompatible with many current numerical challenges which frequently satisfy N 10 6 . The main objective of this paper is to reduce the computational burden of decomposing an operator in an orthogonal wavelet basis using an alternative decomposition called product-convolution.

Product-convolution expansions Let E denote a finite dimensional vector space of discrete d dimensional functions. The functions in E are defined on Ω = {1, . . . , N 1 } × . . . × {1, . . . , N d } with N = d i=1 N i . Throughout the paper we will work with periodic boundary conditions, i.e. consider that Ω is a discretization of the torus T d . Consider a linear operator G defined for all f ∈ E by G : f → y∈Ω K(•, y)f (y), (1.1) where K ∈ R N ×N is the matrix form (the kernel) of the operator. If all the impulse responses S(•, y) = K(• + y, y) are well approximated by their projections on a lowdimensional subspace U = span (u k , 1 ≤ k ≤ m), it is possible to construct a productconvolution expansion H of G defined as

Hf = m k=1 u k (v k f ), (1.2)
where denotes the convolution operator, u k ∈ E is a convolution filter, denotes the point-wise multiplication and v k ∈ E is a multiplier. The multipliers V = (v 1 , . . . , v m ) can be defined as the projections of the impulse responses on U, i.e. v k (y) = S(•, y), u k for all y ∈ Ω, see [START_REF] Escande | Approximation of integral operators using productconvolution expansions[END_REF]. From a numerical perspective, the expansion (1.2) can be evaluated efficiently using fast Fourier transforms with a complexity O(mN log(N )).

Product-convolution expansions have appeared at least 3 decades ago and have been given different names in various fields, see e.g. [START_REF] Busby | Product-convolution operators and mixed-norm spaces[END_REF][START_REF] Nagy | Restoring images degraded by spatially variant blur[END_REF][START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF][START_REF] Gilad | A fast algorithm for convolution integrals with space and time variant kernels[END_REF][START_REF] Denis | Fast approximations of shift-variant blur[END_REF][START_REF] Alger | Scalable matrix-free adaptive product-convolution approximation for locally translation-invariant operators[END_REF]. A remarkable aspect of these expansions is their ability to interpolate linear operators from scattered impulse responses: the bases U and V can be evaluated efficiently for operators with slowly varying impulse responses. Assuming that a few impulse responses S(•, y i ) are known at scattered locations y i ∈ Ω, it is possible to evaluate the basis U using a principal component analysis and to interpolate the projection coefficients over the domain Ω to obtain the multipliers V. This process not only provides a minimax estimate of the operator G [START_REF] Bigot | Estimation of linear operators from scattered impulse responses[END_REF], but also a compact representation compatible with fast numerical algorithms. We refer the reader to [START_REF] Escande | Approximation of integral operators using productconvolution expansions[END_REF] and Section 3.3 of this paper for more insight on the construction and approximation properties of this decomposition.

Wavelet decompositions On the other hand, multi-scale representations of integral operators have led to practical breakthroughs in the theoretical and numerical analysis of partial differential equations and inverse problems. Despite their considerable impact for the compact representation of functions, especially in the field of signal processing, the role of wavelets for coding operators still seems marginal, at least at an industrial level.

A possible explanation is the high cost related to the computation of the wavelet representation of an operator. Let (ψ λ ) λ∈Λ denote an orthogonal wavelet basis and Ψ * : E → R N denote the associated forward wavelet transform. The wavelet representation of G in the wavelet basis is the matrix Θ

= (θ[λ, µ]) λ,µ∈Λ of coefficients θ[λ, µ] = Gψ λ , ψ µ .
For an arbitrary operator G, computing these coefficients has a complexity of order O(N 3 ) operations since the operator G needs to be applied to each of the N wavelets. This cost is prohibitive for large N .

The contribution We address this issue by providing a fast decomposition algorithm for product-convolution operators. Its complexity roughly scales as O(mN log 2 2 N η -α ) for any arbitrary precision η > 0, where α > 0 is a quantity that depends on the smoothness of the operator and on the number of vanishing moments of the wavelet basis.

One may wonder what is the interest of using a wavelet decomposition if a productconvolution is already available. The reason comes from the higher numerical efficiency of wavelet decompositions when incorporated in iterative methods. For instance, we showed in [START_REF] Escande | Sparse wavelet representations of spatially varying blurring operators[END_REF][START_REF] Escande | Accelerating 1 -2 deblurring using wavelet expansions of operators[END_REF] that the simultaneous sparsity of operators and images in the same orthogonal wavelet basis could be leveraged to accelerate the resolution of inverse problems by one or two orders of magnitude. The gain comes from two complementary facts:

• the method can operate in the wavelet domain only and avoid to continuously swap between different domains such as the spatial domain, the Fourier domain and the wavelet domain.

• In addition, efficient diagonal preconditioners can be designed since the matrix Θ is concentrated mostly along its diagonal. This idea is often referred to as a multi-level preconditioner [START_REF] Cohen | Numerical Analysis of Wavelet Methods[END_REF].

We will illustrate the power of these ideas on a practical image deblurring problem.

Preliminaries

Let E denote the linear space of d-dimensional discretized functions of N samples defined on Ω = {1, . . . , N 1 }×. . .×{1, . . . , N d } with d i=1 N i = N . This space E will be often identified with R N through a bijective mapping ω : Ω → {1, . . . , N }. To simplify the discussion, we assume that the numbers N 1 = . . . = N d = 2 J and we use circular boundary conditions. Assuming that the number of pixels is large enough, methods and analyses proposed in this work can be extended to any boundary conditions by using boundary wavelets.

We let f (x) denote the value of a function f at x, v[i] denote the i-th coefficient of a vector v and A[i, j] denote the (i, j)-th element of a matrix A.

Let A be an N × N matrix, we let supp A denote its support i.e. the set of indexes associated to a non-zero coefficient: supp A = (i, j) ∈ {1 . . . N } 2 | A[i, j] = 0 . The approximation rates stated in this paper will be expressed with respect to the spectral norm

• 2→2 defined by A 2→2 = sup f ∈R N , f 2 =1
Af 2 .

One dimensional orthogonal wavelet bases

We first recall the construction of wavelets on the continuous interval [0, 1]. Let φ and ψ denote the scaling and mother wavelets. Translated and dilated versions of the wavelets are defined, for j ≥ 0, as follows

φ j,n = 2 j/2 φ 2 j • -n , ψ j,n = 2 j/2 ψ 2 j • -n ,
with n ∈ {0, . . . , 2 j -1}. A mother wavelet ψ is said to have M vanishing moments when

∀0 ≤ m < M, [0,1] t m ψ(t)dt = 0.
We now detail the construction of a 1D discrete orthogonal wavelet transform. A function f : [0, 1] → R can be discretized leading to a vector v ∈ R 2 J samples. The Fast Wavelet Transform is derived from the observation that wavelets φ and ψ are associated to a filter bank (h, g) [START_REF] Mallat | A Wavelet Tour of Signal Processing -The Sparse Way[END_REF]Theorem 7.7,p. 348]. From these filters (h, g) the discrete wavelets φ j and ψ j can be defined recursively using convolutions and sub-sampling. Setting

φ J,n = 1 {n} , we get for 0 ≤ j < J φ j,l = n∈Z h[n -2l]φ j+1,n , and ψ j,l = n∈Z g[n -2l]φ j+1,n .
We assume that the wavelets are compactly supported i.e. supp(ψ) = [-δ + 1, δ]. Note that δ is related to the number of vanishing moments of the mother wavelet. Let M be this number, we have δ ≥ M , with equality for Daubechies wavelets, see e.g. [START_REF] Mallat | A Wavelet Tour of Signal Processing -The Sparse Way[END_REF]Theorem 7.9,p. 294]. Together with [START_REF] Mallat | A Wavelet Tour of Signal Processing -The Sparse Way[END_REF]Theorem 7.5,p. 286] and g[i] = (-1) 1-i h[1 -i], we deduce that supp h = supp φ = [0, 2δ -1] and supp g = [-2(δ -1), 1]. Therefore for 0 ≤ j ≤ J, supp φ j,0 = 0, (2 J-j -1)(2δ -1) and supp ψ j,0 = -(2 J-j -1)2(δ -1), (2 J-j -1) .

Orthogonal wavelet bases on E

In dimension d, we use isotropic separable wavelet bases, see, e.g., [START_REF] Mallat | A Wavelet Tour of Signal Processing -The Sparse Way[END_REF]Theorem 7.26,p. 348]. Let l = (l 1 , . . . , l d ). Define ρ 0 j,n = φ j,n and ρ 1 j,n = ψ j,n . Let e = (e 1 , . . . , e d ) ∈ {0, 1} d and T j = {0, . . . , 2 j -1} d . For the ease of reading, we will use the shorthand notation λ = (j, e, l) and |λ| = j. We also let J = 1 d log 2 N and

Λ = (j, e, l) | 0 ≤ j ≤ J -1, l ∈ T j , e ∈ {0, 1} d . (2.1)
The wavelet ψ λ is defined by ψ λ (x 1 , . . . , x d ) = ψ e j,l (x 1 , . . . , x d ) = ρ e 1 j,l 1 (x 1 ) . . . ρ e d j,l d (x d ). With these definitions, every signal f ∈ E can be written as

u = u, ψ 0 0,0 ψ 0 0,0 + e∈{0,1} d \{0} J-1 j=0 l∈T j u, ψ e j,l ψ e j,l = λ∈Λ u, ψ λ ψ λ .
Finally, we let Ψ * : E → R N denote the discrete forward wavelet transform and Ψ its inverse. We refer to [START_REF] Mallat | A Wavelet Tour of Signal Processing -The Sparse Way[END_REF][START_REF] Daubechies | Ten Lectures on Wavelets[END_REF][START_REF] Cohen | Wavelets on the interval and fast wavelet transforms[END_REF] for more details on the construction of wavelet bases.

Main results

Assumptions

The main result will be stated under mild regularity and decay assumptions on the kernels u k defined below. Definition 3.1 (Smoothness class A α [7, p. 281]). A convolution kernel u k is said to be αasymptotically smooth if there exists constants C ≥ 0, α > 0 and a compactly supported wavelet basis with α vanishing moments such that the following inequality holds:

| u k ψ λ , ψ µ | ≤ C2 -(d/2+α)||λ|-|µ|| ϑ(λ, µ) -(d+α) , (3.1)
where ϑ(λ, µ)

:= 1 + 2 min(|λ|,|µ|) dist(supp ψ λ , supp ψ µ )
measures the distance between the wavelets supports.

Remark 3.1 (Asymptotic smoothness in a continuous setting ). In a continuous setting, it can be shown that sufficient conditions for the kernels u k to be α-asymptotically smooth is u k ∈ W α,∞ (R d ) with the following decay property

|∂ p u k (x)| ≤ C k (1 + x 2 ) -(|p|+d) , (3.2)
for all multi-indexes p with |p| ≤ α. Hence, the α smoothness property mostly boils down to regularity properties of the kernel. The larger α, the more regular the kernels.

Our main results will be stated under the following assumptions.

Assumption 3.1 (Regularity of the kernel).

• the operator H has the form (1.2).

• there exists α > 0 such that u k ∈ A α for all 1 ≤ k ≤ m.

• the mother wavelet is compactly supported on a hypercube of sidelength 2δ.

Overall the set of assumptions in Assumption 3.1 describes a fairly large variety of operators comprising blurring operators [START_REF] Escande | Sparse wavelet representations of spatially varying blurring operators[END_REF], singular integral operators and pseudo-differential operators [START_REF] Cohen | Numerical Analysis of Wavelet Methods[END_REF].

The algorithm and its guarantees

The main objective of this paper is to efficiently compute a sparse approximation Θ of the wavelet representation Θ of H defined by Θ = Ψ * HΨ.

(3.

3)

The proposed algorithm heavily relies on the peculiar structure of product-convolution expansions. In what follows, we let U k denote the convolution operator with

u k , V k = diag(v k ) denote the k-th multiplier, A k = Ψ * U k Ψ denote the wavelet representation of U k and B k = Ψ * V k Ψ denote
O(N log 2 N ) Threshold A k to get an k approximation A k O(N log N ) Construct B k (using a sparse cascade algorithm) O(N log N ) Compute C k = A k B k O(N log N -d/t k ) Accumulate Θ ← Θ + C k O(N log N -d/t k ) end for return Θ an η-approximation of Θ
Our main theoretical result reads as follows.

Theorem 3.1. Under Assumption 3.1, Algorithm 1 produces a matrix Θ satisfying Θ - Θ 2→2 ≤ η in O δ d+1 m max(1,d/α) N log 2 2 N η -d/α operations.
Furthermore, the number of coefficients in Θ can also be bounded above by O δ d N log 2 2 N η -d/α . Remark 3.2. In many applications, the kernels of integral operators are smooth. In particular, for α ≥ d the term max(1, d/α) = 1. The complexity is thus linear in m. Remark 3.1 shows that this is satisfied when the kernels u k are of class W d,∞ .

Remark 3.3. The bound on the number of coefficients in Θ provided by Theorem 3.1 allows to bound the complexity of matrix-vector products with Θ. If we assume that the operator G belongs to the class A α , a direct application of [START_REF] Beylkin | Fast wavelet transforms and numerical algorithms I[END_REF][START_REF] Escande | Sparse wavelet representations of spatially varying blurring operators[END_REF] (see also Theorem 5.1), shows that one can construct an η approximation G of G, containing no more than O(δ d N η -d/α ) coefficients. The result of Theorem 3.1 is therefore optimal up to the log 2 N factor. We believe that it could be discarded by assuming further regularity of the multipliers v k . We could then use an additional approximation in place of B k , in the exact same fashion as for the A k 's.

Getting a product-convolution expansion

While the last two items in Assumption 3.1 are now well established, the combination with the first one may seem problematic. We review some numerical approaches to design product-convolution expansions below.

Naive interpolation of impulse responses

The simplest and probably the most widespread product-convolution expansions in image processing are based on the assumption of slow variations of the impulse responses of the operator in the domain Ω. Under this assumption, it is possible to set

u k = S(•, y k )
i.e. (u k ) m k=1 are the impulse responses of the operator at some locations (y k ) m k=1 . By setting the locations (y k ) m k=1 as a coarse Euclidean grid [START_REF] Nagy | Restoring images degraded by spatially variant blur[END_REF], it is possible to keep a small value for m. The functions (v k ) m k=1 are then used to interpolate the impulse responses i.e. they define a partition of unity with the constraints v k (y k ) = 1 for all 1 ≤ k ≤ m. The choice of these v k are discussed in many works, without particular guarantees on the approximation obtained. We refer the interested reader to [START_REF] Denis | Fast approximations of shift-variant blur[END_REF] for a nice overview and [START_REF] Escande | Approximation of integral operators using productconvolution expansions[END_REF] for approximation rates.

In the context of PDEs, product-convolution expansions are encountered in Schur complement methods for solving PDEs, Dirichlet-to-Neumann maps and in PDE-constrained optimization as Hessians. In [START_REF] Alger | Scalable matrix-free adaptive product-convolution approximation for locally translation-invariant operators[END_REF], the functions (u k , v k ) m k=1 are chosen as above, but instead of being defined on a Euclidean grid, locations (y k ) m k=1 are adaptively sampled to best approximate the operator with a fixed number m.

Interpolating expansions of impulse responses

Under mild regularity assumptions of the impulse responses, product-convolution expansions can be designed from scattered impulse responses in a more efficient fashion than what was just described [START_REF] Flicker | Anisoplanatic deconvolution of adaptive optics images[END_REF][START_REF] Gentile | Interpolating point spread function anisotropy[END_REF][START_REF] Bigot | Estimation of linear operators from scattered impulse responses[END_REF].

The main idea is to use a principal component analysis of the observed impulse responses S(•, y k ) to obtain a low-dimensional orthogonal basis (u k ). The coefficient maps (v k ) can then be computed by interpolating the known projection coefficients at the positions (y k ). We showed in [START_REF] Bigot | Estimation of linear operators from scattered impulse responses[END_REF] that this approach was optimal from an approximation theoretic point of view.

Singular value decomposition

When the operator G in (1.1) is fully known, the optimal way (in Frobenius norm) for fixed m to get a product-convolution expansion of an operator (1.1) is to construct a lowrank approximation of the SVIR S (and not the kernel K). For many practical problems, the SVIR S has a much lower rank than the kernel K. A typical example is the identity operator which has a kernel of rank N and an SVIR of rank 1. We refer to [START_REF] Escande | Approximation of integral operators using productconvolution expansions[END_REF] for more examples and mathematical insights.

The low-rank decomposition of space varying impulse response can be achieved by performing a singular value decomposition (SVD) on S. For most problems, this approach is however intractable due to the size of matrices S.

Alternatively, we may just assume that the action of an operator with a matrix form S can be evaluated on any vector. The couples (u k , v k ) m k=1 can now be obtained using a randomized SVD of the matrix S [START_REF] Halko | Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions[END_REF]. The complexity of such algorithms scales as O(m 2 N ), where m is the number of matrix-product evaluations.

4 Numerical experiments -Spatially varying deblurring

Description of the operators

To assess the performance of the algorithm, we propose to perform numerical experiments on a large scale image deblurring problem. We synthesize two different blurring operators:

• The first one, in Figure 1a, is made of isotropic Gaussian impulse responses with a variance σ(y 1 , y 2 ) increasing along the vertical direction only (i.e. σ(y 1 , y 2 ) = 3y 2 for (y 1 , y 2 ) ∈ [0, 1] 2 ). The impulse responses are then truncated out, so that 99% of the Gaussian's energy is kept. This allows to encode the operator as a sparse matrix to accelerate the explicit computation of Θ. Notice that this is however unnecessary to apply our algorithms. As an indication, the largest support has a size 19 × 19 pixels.

• The second operator, in Figure 1b models realistic degradations appearing in optical systems [START_REF] Simpkins | Parameterized modeling of spatially varying optical blur[END_REF]. The impulse responses are rotated and skewed Gaussian, with parameters depending on their location with respect to the center of the field. As an indication, the support of the impulse responses are of size 25 × 25.

We compute the product-convolution expansions of each of these operators using a randomized SVD of the matrix S as explained in Section 3.3.3. The order m of the expansion has been minimized so as to obtain a fixed deblurring performance, see paragraph 4.3.1. The first operator -Figure 1a -is decomposed with m = 5 coefficients, while the second one -Figure 1b -with m = 25 coefficients.

Comparing decomposition timings

In this section, we compare Algorithm 1 to standard decomposition methods. The direct algorithm to compute Θ from a product-convolution expansion is given in Algorithm 2. It consists in using a matrix-vector product and a wavelet transform for each column. Its complexity is therefore in O(N 3 ) operations (this complexity is reduced here to O(N 2 ) since we use the fact that the impulse responses are compactly supported). In this comparison, we used parallel implementations running on 12 cores on a workstation with 256GB of RAM in double precision. Figure 2 compares the computing times using 3 Algorithms:

• Algorithm 2 with a direct computation with a sparse matrix (Alg. 2 -exact).

• Algorithm 2 with a product-convolution expansion (Alg. 2 -PC).

• The proposed Algorithm 1 (Alg. 1).

We use square images containing N = 2 n × 2 n pixels with n ∈ {7, . . . , 12}. The maximal image size is therefore 4096 × 4096. We run the algorithms with a Symmlet wavelet basis of order 6. This choice is explained in Remark 4.1. The precision η = 5.10 -4 has been chosen so as to obtain good performance for a deblurring experiment. This choice is further detailed in Section 4.3.1.

In this setting, Algorithm 1 is much faster than any instance of Algorithm 2. For n = 12, corresponding to images of 16 millions pixels, the speed-up from Alg. 2 -productconvolution to Alg. 1 is 4680 for the operator on Figure 1a and 2330 for the operator on Figure 1b.

We also evaluate the decomposition times for various precisions with N = 256 × 256 being fixed, see Figure 3. Assuming that the clock time is proportional to the number of operations, we deduce that the number of operations is bounded by a quantity proportional to η -d/α . We fit the curves with a line to obtain an approximation of the smoothness α of each operator. Blur Fig. 2b η -2/(2.38) In [START_REF] Escande | Sparse wavelet representations of spatially varying blurring operators[END_REF] we performed extensive numerical experiments of the approximation properties of a wavelet basis w.r.t. the number of vanishing moments. As expected, taking as many moments as possible was preferable for the approximation rate. However, increasing it too much led to insignificant approximation quality while deteriorating the numerical complexity significantly. This is mostly due to the constants of Theorem 5.1 which increase with the number of vanishing moments.

We performed additional numerical simulations -not reported here -to assess the behavior of the complexity of Algorithm 1 w.r.t. the sidelength δ of the wavelet basis. This sidelength is related to the number of vanishing moments M since δ ≥ M [23, Theorem 5.7, p. 286]. Numerically, we observed that the algorithm actually performs better than the rate δ d+1 in Theorem 3.1.

A deblurring experiment

The setting We assume that an observed image f 0 ∈ E reads

f 0 = Hf + b,
where f ∈ E is the clean image to recover, b ∼ N (0, σ 2 Id N ) is a white Gaussian noise of standard deviation σ and H ∈ R N ×N is the blurring operator. The image f used in this experiment is shown in Figure 4a and Figure 4 contains the degraded version f 0 with the two operators considered.

The optimization problem In this experiment we will seek to recover the image f by solving the following variational problem arg min

f ∈E E(f ) = 1 2 Hf -f 0 2 2 + Ψ * f 1,w , (4.1)
where z 1,w is a weighted 1 -norm with weights w ∈ R N and defined by The solution of this variational problem does not correspond to the state-of-the-art in terms of image quality since the prior is simple, but it is known to perform well in short computation times. Figure 5 shows the solutions of (4.1) for the two blurs, with Ψ being a Symmlet wavelet basis of order 6, since it offers a good compromise between computing times and visual quality of the results [START_REF] Escande | Accelerating 1 -2 deblurring using wavelet expansions of operators[END_REF]. We use this basis in all the experiments. The proposed algorithms Many algorithms were designed to solve non-smooth convex problems of the form (4.1). The two predominant algorithms are the (accelerated) proximal gradient descent [START_REF] Nesterov | Gradient methods for minimizing composite functions[END_REF] also known as FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] and the alternating descent method of multipliers (ADMM) [START_REF] Glowinski | Lectures on Numerical Methods for Non-linear Variational Problems[END_REF]. As shown in [START_REF] Escande | Accelerating 1 -2 deblurring using wavelet expansions of operators[END_REF], a very efficient version of FISTA can be derived by approximating (4.1) in the wavelet domain i.e. arg min

z 1,w = λ∈Λ w[λ] |z[λ]| .
z∈R N 1 2 Θ L z -z 0 2 2 + z 1,w , (4.2)
where z 0 = Ψ * f 0 , and Θ L is an L-sparse approximation of Θ = Ψ * HΨ. This idea allows to avoid constantly swapping between the wavelet and spatial domains, where most of the time is spent in traditional solvers. In addition, the matrix Θ * L Θ L -which appears in the gradient of the quadratic term -is mostly concentrated around its diagonal that decreases across sub-bands. An efficient diagonal preconditioner can be designed exploiting this property and further reduce the number of iterations by a factor roughly equal to 2 without compromising image quality.

The ADMM on its side requires the inversion of symmetric positive definite linear systems with matrix Θ * L Θ L + τ I N . This step could also be accelerated using the same ideas, but we observed that it was less efficient in practice and will not further report about this.

Table 1 compares the performance of various algorithms detailed in more details in Appendix A.

Algorithm Description Complexity

MVP FFT FWT Diag

FISTA-Spa

FISTA with the exact operator

2 - 2 - FISTA-PC
FISTA using a PC expansion

- 2m 2 2m FISTA-W FISTA with sparse Θ L 2 - - - FISTA-WP
FISTA with Θ L and Jacobi preconditioner

2 - - 1 ADMM-PC See Appendix A.2 - 3m 3 6m+2
Table 1: List of algorithms compared in the numerical experiments with their complexity per iteration in terms of matrix vector products (MVP), fast Fourier transforms (FFT), wavelet transforms (WT) and product with diagonal matrices (Diag).

Choice of m and η

The approximation properties of product-convolution expansions are investigated in [START_REF] Escande | Approximation of integral operators using productconvolution expansions[END_REF] in terms of Frobenius distance on the operators. We are dealing here with an inverse problem thus the approximation needs are different. Because of the regularization term in (4.1), fine approximations of the operator are unnecessary, allowing to reduce the computational burden, see Figure 6. We choose the order of expansion m in a such a way that the pSNR of the deblurred images does not differ from more than 0.01dB from the pSNR of the solution of (4.1). From Figure 6, we obtain m = 5 for the operator on Figure 1a and m = 25 for the operator on Figure 1b. Note that m is larger for the second blur since the impulse response variations are more complex.

In a similar fashion, we select the precision η in such a way that the pSNR of the solution of (4.2) does not differ from more that 0.04dB from the pSNR of the solution of (4.1). We do not report the plots, but we found that η = 5.10 -4 was sufficient to ensure this requirement.

Computing times

In this paragraph, we precisely compare the different methods designed to solve (4.1). The methodology consists in:

• finding the number L for which the solution of (4.2) has a decrease of pSNR of less than 0.2dB w.r.t. the pSNR of the solution of (4.1), • for each algorithm, finding a number of iteration N it leading to a precision.

E(f (N it) ) -E(f (0) ) ≤ 10 -3 E(f (0) ). ( 4.3)
The cost functions are reported in Figure 7, while the timings are reported in Tables 2 and3. Deblurring were performed on Matlab with automatic multi-threading disabled (using the -singleCompThreah option), to avoid uncontrolled behaviors while timing. To further demonstrate the efficiency of FISTA-WP, we implement it in CUDA and run it on a GeForce GTX Titan X in double precision. 1b for the algorithms in Table 1 and in Figure4a (N = 1024 × 1024). 1 for an image of size 1024 × 1024. The cost functions are displayed with respect to the number of iterations on the left column and the time on the right column. The first row corresponds to the blur on Figure 1a and the second one to the blur on Figure 1b.

This example confirms that solving problem (4.1) by an approximation in the wavelet domain (4.2) leads to dramatic reduction of computation times (×16 and ×62 speedup factors) with little loss in the quality (less than 0.2dB). Furthermore, the use of the structure of Θ L allows to derive efficient diagonal preconditioners which lead to an extra ×2 speed-up.

Proof of the main result

The proof exploits the peculiar structure of product-convolution expansions:

• the convolution operator U k = u k • can be decomposed in the wavelet basis and yields a matrix

A k = Ψ * U k Ψ in O(N log 2 2 N ) operations.
As long as U k belongs to the smoothness class described in Definition 3.1, the matrix A k can be thresholded to obtain a sparse approximation A k of A k .

• the multipliers V k = v k • are diagonal. Therefore the sparsity pattern of B k = Ψ * V k Ψ is included in the set of wavelet with supports intersecting the diagonal and contain at most O(N log 2 N ) non-zero coefficients, with known locations. They can be computed in O(N log 2 N ) operations using a cascade algorithm.

The general approach is therefore to compute each couple (A k , B k ), compute the products A k B k and accumulate the sum to get an approximation Θ = m k=1 A k B k of Θ. The last two steps are also critical since in the general case, the number of coefficients can explode in the product A k B k or when summing them all. However, relying again on the peculiar structure of the matrices A k and B k , one can prove that the complexity of the algorithm stays quasi-linear. In what follows we precisely describe each of these steps.

Structure of A k

The matrix A k = ΨU k Ψ * is the wavelet representations of a convolution operator U k . As such, it inherits a peculiar structure. As observed in [START_REF] Escande | Accelerating 1 -2 deblurring using wavelet expansions of operators[END_REF] it has circulant sub-bands. To make it precise, we recall that a wavelet representation Θ can be decomposed into its wavelet sub-bands: Θ = Θ e,e j,j e,e ∈{0,1} d j,j ≤J-1

, with Θ e,e j,j [l, l ] = Θ[λ, µ].

for l ∈ T j , l ∈ T j and λ = (j, e, l), µ = (j , e , l ). For instance on 1D signals with J = 2, the sub-bands of Θ are shown in Figure 8. The sub-bands Θ e,e j,k themselves have a specific structure captured by the following definition.

Definition 5.1 (Rectangular circulant matrices). Let P ∈ R 2 jd ×2 kd denote a rectangular matrix and τ a be the translation operator by vector a ∈ Ω. The matrix P is called circulant if and only if

• When k ≥ j : there exists p ∈ R 2 kd such that P [l, :] = τ 2 k-j l p for all l ∈ T j .

• When k < j : there exists p ∈ R 2 jd such that P [:, l] = τ 2 j-k l p for all l ∈ T k . Furthermore, the matrices A k are assumed to belong to the smoothness class A α governing the decay of its coefficients. This properties allows to obtain efficient approximations of A k as described by Theorem 5.1.

Theorem 5.1 ([7, Theorem 4.6.2]). Let A ∈ A α . For all L ≥ 0, one can construct a matrix A L with at most L non-zero entries of A per row and column such that

A -A L 2→2 L -α/d . (5.1)
As a summary, Lemma 5.1 together with Theorem 5.1 ensure that the computation of the matrices A k and their approximations at any precision are tractable in large dimensions.

Remark 5.1. Note that the approximation step is crucial to obtain a quasi-linear complexity of the overall algorithm. If no approximation were made, the complexity of the product A k B k could become quasi-quadratic.

Remark 5.2. The proof of Theorem 5.1 gives an explicit way of building A L from A. We recall it here for sake of completeness. The matrix A L is obtained by zeroing all entries satisfying either (i) ||λ| -|µ|| > J(L) or (ii) ϑ(λ, µ) > 2 J(L) where J(L) is linked to L through the relation L J(L)2 dJ(L) . A straightforward algorithm building A L consists in looping along all rows λ and rows µ and discard elements satisfying either (i) or (ii). This algorithm requires O(N 2 ) operations in the general case. Due to the particular structure of matrix A -only log 2 (N ) coefficient are non-zero per row (Lemma 5.1) -this algorithm has O(N log 2 N ) complexity. In practice, many authors propose to simply threshold A. This comes from the observation that keeping coefficients that are larger than 2 -(d/2+α)J do not satisfy (i) and (ii) (from Definition 3.1). Again, since A is encoded with O(N log 2 N ) coefficients, thresholding A costs O(N log 2 N ) instead of O(N 2 ) in the general case.

Structure of B k

We now turn on discussing the structures of the matrices B k = ΨV k Ψ * . Since the operators V k are diagonal, the coefficients B k [λ, µ] = v k ψ λ , ψ µ necessarily vanish as long as supp ψ λ and supp ψ µ do not intersect. Define

I = {(λ, µ) ∈ Λ × Λ | supp ψ λ ∩ supp ψ µ = ∅} ,
we thus have supp B k ⊆ I for all 1 ≤ k ≤ m. This property allows to derive upper bounds for the number of coefficients in a column and a row of B k (Lemma 5.2) as well as for the cardinality of supp B k (Lemma 5.3). Lemma 5.2. Assume that the mother wavelet is compactly supported on a hypercube of sidelength 2δ. Let λ ∈ Λ and j = |λ|, there are at most (2δ) d (2 d -1)j + 2 d(J-j) indexes µ ∈ Λ such that supp ψ λ ∩ supp ψ µ = ∅.

Proof. Let λ = (j, e, l) ∈ Λ. From the definition of ψ λ we get that supp ψ λ ⊆ -(2 J-j -1)(δ -1) + 2 J-j l, (2 J-j -1)δ + 2 J-j l .

Consider now another µ = (j , e , l ) ∈ Λ. Necessary conditions for the intersections of the supports of ψ λ and ψ µ are -(2 J-j -1)(δ -1) + 2 J-j l ≤ (2 J-j -1)δ + 2 J-j l -(2 J-j -1)(δ -1) + 2 J-j l ≤ (2 J-j -1)δ + 2 J-j l , where the ≤ has to be understood as a component-wise inequality for d > 1. Without loss of generality, we assume that l = 0.

We let N j,j be the set of l ∈ Λ satisfying the above conditions i.e.

   l ≥ -2 j -J 2 J-j + 2 J-j -2 δ + 2 j -J (2 J-j -1) l ≤ 2 j -J 2 J-j + 2 J-j -2 δ -2 j -J (2 J-j -1)
The volume of N j,j is therefore

#N j,j = 2 j -J 2 J-j + 2 J-j -2 (2δ -1) + 1 d ≤ 2 j 2 -j + 2 -j (2δ -1) + 1 d ≤ 2 j 2 -j + 2 -j 2δ + 1 -2 j -j -1 d ≤ 2 j 2 -j + 2 -j 2δ d Therefore, # {µ ∈ Λ | supp ψ µ ∩ supp ψ λ = ∅} ≤ J-1 j =0 e∈{0,1} d \{0} l ∈T j 1 N j,j (l ) = (2 d -1) J-1 j =0 2 j 2 -j + 2 -j 2δ d ≤ (2 d -1)(2δ) d   j-1 j =0 1 + J-1 j =j 2 d(j -j)   = (2 d -1)(2δ) d j + 2 -jd 2 dJ -2 dj 2 d -1 ≤ (2δ) d (2 d -1)j + 2 d(J-j)
which gives the upper-bound on the number of intersections. Proof.

#I = # {(λ, µ) ∈ Λ × Λ | supp ψ µ ∩ supp ψ λ = ∅} ≤ J-1 j=0 (2 d -1)2 jd 2 d δ d (2 d -1)j + 2 d(J-j) ≤ (2 d -1)2 d δ d   (2 d -1) J-1 j=0 j2 jd + J-1 j=0 2 dJ   ≤ (2 d -1)2 d δ d (2 d -1)J 2 dJ -1 2 d -1 + J2 dJ ≤ 2(2 d -1)2 d δ d J2 dJ
The numerical computation of B k = Ψ * V k Ψ is conceptually simple. It consists in two steps:

Step 1 Compute the wavelet transform of each column of

V k i.e. build a matrix V k such that V k [λ, i] = V k [•, i],
ψ λ for all 1 ≤ i ≤ N and λ ∈ Λ. Each column of V k is a discrete Dirac mass, i.e. contains only one element. Therefore, its wavelet transform can be computed with a minimal number of operations using a cascade algorithm taking advantage of the signals sparsity. A cascade algorithm consists in applying discrete convolutions with filters h and g and subsampling the outputs in a recursive fashion. It is derived from the multiresolution structure of a wavelet transform. We refer to [START_REF] Mallat | A Wavelet Tour of Signal Processing -The Sparse Way[END_REF]Section 7.3] for more details.

Step 2 Compute the wavelet transform of each row of

V k i.e. B k [λ, µ] = V k [λ,
•], ψ µ for all λ, µ ∈ Λ. Similarly, the λ-th row of V k is sparse: it contains at most (2 J-j 2δ) d contiguous coefficients. The same sparse cascade algorithm can be used to compute its wavelet transform.

We illustrate the sparsity structure of the matrices V k and B k in Figure 9 for the case of 1D signals. Observe that each row of V k contains only a few contiguous nonzero indexes.

To control the complexity of each step, we will use the following result.

Lemma 5.4. Let f ∈ E be supported on the hypercube of Ω defined by the lower vertex q 0 ∈ Ω and the side length κ. We call q 1 ∈ Ω the upper vertex of the hypercube i.e. q 1 = q 0 + κ. Let c = Ψ * f be its wavelet transform and c j,e its (j, e)-sub-band i.e. c j,e [l] = c[λ] = f, ψ λ for all l ∈ T j and λ = (j, e, l). We have

supp c j,e = 2 j-J q 0 -(1 -2 -J ), 2 j-J q 1 + (1 -2 -J )2(δ -1) # supp c j,e ≤ 2 d 2 d(j-J) κ d + (2δ) d
Furthermore, the number of operations to compute c is bounded above by Proof. Let λ = (j, e, l) ∈ Λ. From the definition of ψ λ we get that supp ψ λ ⊆ -(2 J-j 1)(δ -1) + 2 J-j l, (2 J-j -1)δ + 2 J-j l .

2dδ2 d κ d + (2 d -1)(2δ) d J .
Necessary conditions for the intersections of the supports of ψ λ and f are

q 0 ≤ (2 J-j -1)δ + 2 J-j l q 1 ≥ -(2 J-j -1)(δ -1) + 2 J-j l.
where the inequalities have to be understood component-wise. These conditions are equivalent to 2 j-J q 0 -(2 J-j -1) ≤ l ≤ 2 j-J q 1 + (1 -2 J )2(δ -1) .

The volume of such a set is # supp c j,e = 2 j-J q 1 -2 j-J q 0 + (1 -

2 -J )2(δ -1) + (1 -2 -J ) + 1 d = 2 j-J (q 1 -q 0 + 1) + (1 -2 -J )2δ d ≤ 2 d 2 d(j-J) κ d + (2δ) d
Since filters g and h have the same length, the details and average coefficients of the wavelet transforms have the same size in each sub-band -albeit not the same support. Therefore we only analyze the detail coefficients. The computation of one c j,e is obtained from d convolutions of average coefficients at scale j + 1 with wavelet filters of size 2δ. Hence the number of operations to compute one c j,e is d2δ# supp c j,e and the total number of operations needed to compute c is

J-1 j=0 (2 d -1)2dδ# supp c j,e ≤ 2dδ(2 d -1)2 d J-1 j=0 2 d(j-J) κ d + (2δ) d ≤ 2dδ(2 d -1)2 d κ d 1 2 d -1 + (2δ) d J (5.2)
which ends the proof.

Gathering everything

We will now gather all the ingredients to prove Theorem 3.1. The first part of the proof certifies the accuracy of the approximation.

Θ -Θ 2→2 ≤ m k=1 A k B k -A k B k 2→2 ≤ m k=1 B k 2→2 A k -A k 2→2 ≤ m k=1 V k 2→2 k = m k=1 v k ∞ k = η.
We used the fact that B k 2→2 = V k 2→2 = v k ∞ since Ψ is an orthogonal wavelet basis and V k is a diagonal matrix. The last equality is derived from the definition of

k = η m v k ∞ .
The second part deals with the overall complexity.

Proof of Theorem 3.1. The overall complexity of computing Θ is bounded above by a quantity proportional to

c(d)    m k=1 N log 2 2 N Lemma 5.1 +N log 2 N + δ d+1 N log 2 2 N Lemma 5.6 +2 δ d N log 2 N -d/α k Lemma 5.7    ≤ c(d) 3δ d+1 mN log 2 2 N + 2N log 2 N m k=1 η m v k ∞ -d/α ≤ c(d) δ d+1 mN log 2 2 N + N log 2 N η -d/α m d/α m k=1 v k d/α ∞ The constant c(d) depends only on d. The term m k=1 v k d/α ∞ depends on the considered operator H. Finally, Θ is computed in O δ d+1 m max(1,d/α) N log 2 2 N η -d/α
. The bound on the number of coefficients in Θ is obtained using Lemma 5.9. Since there exists a support S such that all supp C k ⊆ S we get that supp Θ ⊆ S.

Remark 5.3. The complexity analysis does not take into account that v k may have some smoothness properties. In this case, the associated B k might belong to a class A α and could therefore be thresholded to further speed up the algorithm.

The Forward-Backward algorithm applied to the right hand-side of (A.1) starts with an initial guess z (0) = y (1) ∈ R N and iterates for i ≥ 1 as

     z (i) = Prox • 1,w y (i) -τ Ψ * H * (HΨy (i) -f 0 ) y (i+1) = z (i) + i -1 i + 2 z (i) -z (i-1)
with a step-size τ ≤ H -2 2→2 . This algorithm can be applied verbatim with the exact spatial operator (coined as FISTA-Spa). In the same fashion, we design another algorithm, named FISTA-PC, using the product-convolution expansion.

One can also use FISTA to solve (4.2), the approximation of (4.1) in the wavelet domain. It starts by setting z (0) = y (1) ∈ R N and iterates for i ≥ 1 as

     z (i) = Prox Σ • 1,w y (i) -τ ΣΘ * L (Θ L y (i) -z 0 ) y (i+1) = z (i) + i -1 i + 2 z (i) -z (i-1)
where Σ is a diagonal preconditioner equal to the identity or to the Jacobi preconditioner diag(Θ * L Θ L ). Note that in this alternative, one iteration costs only 2 matrix-vector products with the sparse matrix Θ L . We will refer to this algorithm as FISTA-W when Σ = I N and FISTA-WP when Σ = diag(Θ * L Θ L ) (see [START_REF] Escande | Accelerating 1 -2 deblurring using wavelet expansions of operators[END_REF] for more details).

A.2 ADMM for product-convolution

The resolution of (4.1) using an ADMM or a Douglas-Rachford is expensive due to the resolution of non trivial linear systems at each iteration. We can adapt ideas proposed in [START_REF] O'connor | Total variation image deblurring with spacevarying kernel[END_REF] and exploit the product-convolution structure to accelerate the resolution of the linear systems and improve its efficiency. We describe this idea below. A product-convolution operator H can be decomposed as: The augmented Lagrangian associated to this problem reads

Hf = m k=1 u k v k f = UVf, with V = V 1 . . . V m T ∈ R
L(f, g, z, ξ, ζ) = 1 2 Ug-f 0 2 2 + z 1,w + ξ, g -Vf + ζ, z -Ψ * f + β 1 2 g-Vf 2 2,γ + β 2 2 z-Ψ * f 2 2 .
The norm • 2,γ is just a weighted l 2 -norm for the augmented term of g = Uf . Weighting the influence of each U k allows accelerating the convergence of the ADMM, i.e. it acts as a preconditioner. This norm is defined by The ADMM amounts to iterating:

• g (i+1) = arg min g∈R mN L f (i) , g, z (i) , ξ (i) , ζ (i) .

The solution of this sub-problem is given by

g (i+1) = (U * U + β 1 D γ ) -1 U * f 0 -ξ (i) + β 1 D γ Vf (i) .
The matrix U * U contains m × m blocks of size N × N with the (k, l)-th block populated with U * k U l . The direct inversion of this matrix will be inefficient even in the Fourier domain. However, using the Woodbury matrix identity [START_REF] Hager | Updating the inverse of a matrix[END_REF], the solution can be expressed as

g (i+1) = D -1 γ -D -1 γ U * (Id N + UD -1 γ U * ) -1 UD -1 γ U * f 0 -ξ (i) + β 1 D γ Vx (i) . Matrix (Id N + UD -1 γ U * ) = Id N + m k=1 γ k U k U * k ∈ R N ×N
is diagonal in the Fourier domain and can be efficiently inverted. In this step, the computation of 3m FFTs and 4m + 1 products with diagonal matrices are involved.

• z (i+1) = arg min

z∈R N L f (i) , g (i+1) , z, ξ (i) , ζ (i) = Prox 1 β 2 • 1,w Ψ * f (i) -1 β 2 ζ (i) .
This step costs one wavelet tranform plus a thresholding operation.

• f (i+1) = arg min i+1) , z (i+1) , ξ (i) , ζ (i) .

f ∈R N L f, g (
The solution of this sub-problem reads

f (i+1) = (β 1 V * D γ V + β 2 Id N ) -1 V * (ξ (i) + β 1 D γ g (i+1)
) + Ψ(β 2 z (i+1) + ζ (i) ) .

Matrix

β 1 V * D γ V + β 2 Id N = β 1 m i=1 γ i V i V i + β 2
Id N is N × N diagonal matrix and can be efficiently inverted.

In this step, one wavelet transform and m + 1 products with diagonal matrices are involved.

• Update Lagrange multipliers

ξ (i+1) = ξ (i) + β 1 D γ (g (i+1) -Vf (i+1) ) ζ (i+1) = ζ (i) + β 2 (z (i+1) -Ψ * f (i+1) ).
This step requires m products with diagonal matrices and one wavelet transform.

To summarize, one iteration of the ADMM costs 3m FFTs, 3 wavelet transforms and 6m + 2 products with diagonal matrices. This is substantially more than one interation of FISTA with the product-convolution expansion i.e. 2m FFTs, 2 wavelet transforms and m products with diagonal matrices.

Algorithm 1

 1 the wavelet representation of the multiplier. The proposed methodology is described in Algorithm (1). Decomposition of product-convolution in an orthogonal wavelet basis Require: A precision η > 0, the Set k = η m v k ∞ for all 1 ≤ k ≤ m Set Θ an empty sparse matrix for all k = 1 → m do Compute A k (by exploiting the convolution structure)

Algorithm 2

 2 Naive decomposition of matrices in an orthogonal wavelet basis for all λ ∈ Λ do Get ψ λ and compute w λ = Hψ λ Compute Ψ * w λ to obtain ( Hψ λ , ψ µ ) µ∈Λ . Set it as the λ-th column of Θ end for (a) Operator 1 (applied to a 1024 × 1024 image). (b) Operator 2 (applied to a 1024 × 1024 image)

Figure 1 :

 1 Figure 1: An illustration of the spatially varying blurs used in the numerical experiments. Blurring operators are applied to a Dirac comb to obtain the impulse responses at various locations.

Figure 2 :

 2 Figure 2: Running times for various number of pixels N = 2 n × 2 n with n ∈ {7, . . . , 12}. Left: for the blurring operator on Figure 1a and Right: for the one on Figure 1b.

Remark 4 . 1 .

 41 The choice of the wavelet basis is of crucial importance as it has to simultaneously provide a good representation for the image and for the operator. A key parameter is the number of vanishing moments of the wavelet basis. As suggested by Definition 3.1, the number of vanishing moments should be at least equal to the smoothness α of the convolution kernels u k in order to get the best possible approximation.

Figure 3 :

 3 Figure 3: Decomposition times for various precisions and N = 256 × 256.

( a )

 a Original image f -1024 × 1024 pixels (b) Blurred with blur Figure 1a -SNR = 25.30 dB (c) Blurred with blur Figure 1b -SNR = 22.72 dB

Figure 4 :

 4 Figure 4: Images f 0 , degraded versions of f for the two operators considered. Noise level σ = 5.10 -3

( a )

 a Restored image -SNR = 27.66 dB (b) Restored image -SNR = 29.71 dB

Figure 5 :

 5 Figure 5: Solutions of (4.1) for the two operators considered. The exact operator is used with w[λ] = 2 • 10 -2 • |λ|.

Figure 6 :

 6 Figure 6: pSNR of the solution of (4.1) obtained with the product-convolution expansion of order m (blue) and the exact operator (red). Left: for the operator on Figure 1a. Right: for the operator on Figure 1b.

Figure 7 :

 7 Figure 7: Performance of the deblurring methods in Table1for an image of size 1024 × 1024. The cost functions are displayed with respect to the number of iterations on the left column and the time on the right column. The first row corresponds to the blur on Figure1aand the second one to the blur on Figure1b.

Figure 8 :

 8 Figure 8: Sub-band structure of convolution matrices in the wavelet domain. Here on 1D signals with J = 2.

Lemma 5 . 1 (

 51 [START_REF] Escande | Accelerating 1 -2 deblurring using wavelet expansions of operators[END_REF] Theorem 3]). The sub-bands of matrix A k are circulant. Therefore, only O (2 d -1)N log 2 N coefficients are needed to encode A k . They can be computed with no more than O (2 d -1)N log2 2 N operations.

Lemma 5 . 3 .

 53 We have #I ≤ c(d)δ d J2 dJ , where the constant c(d) = 2(2 d -1)2 d . Therefore since J = 1 d log 2 (N ), the total number of intersections is at most c(d) d δ d N log 2 (N ).

( a )Figure 9 :

 a9 Figure 9: Wavelet representations of a diagonal matrix. Here we illustrate the case of a 1D wavelet transform of size 1024 × 1024 with J = 3 decomposition levels.

Lemma 5 . 9 .

 59 There exists S ⊂ Λ 2 s.t. supp A k B k ⊂ S for all 1 ≤ k ≤ m. Moreover, #S = O c(d)δ d N log 2 N (min k k ) -d/α .Proof. Let k 0 be the index such that k 0 = min k k . Since A k 0 belongs to the smoothness class A α , there exists a F ⊂ Λ 2 s.t. supp A k 0 ⊆ F. Furthermore all A k belong to the same A α we have that supp A k ⊆ supp A k 0 .On the other hand, since all V k are diagonal matrices supp B k ⊆ I for all k. Define I (resp. F ) a N × N matrix with ones on I (resp. F) and zeros outside. Let S = supp F I. Then obviously supp A k B k ⊆ S. Moreover, from Lemma 5.7, the number of non-zero coefficients in F I is bounded above by c(d)δ d N log 2 N

Lemma 5 . 10 .

 510 The matrix Θ returned by Algorithm 1 satisfies Θ -Θ 2→2 ≤ η Proof. Using triangular inequality and since • 2→2 is multiplicative we get

2 Ug -f 0 2 2

 22 mN ×N and U = U 1 . . . U m ∈ R N ×mN . This decomposition allows to split the problem (A.1) as min f ∈R N ,g∈R mN ,z∈R N z=Ψ * f,g=Vf 1 + z 1,w

• 2 2 ,

 2 γ = D γ •, • R mN , with D γ = mN ×mN a block diagonal matrix where γ k > 0 and m k=1 γ k = 1. The choice γ k = U k 2→2 led to good practical behaviors.

Table 2 :

 2 Timings and iterations number to reach criterion (4.3) for the blur in Figure1afor the algorithms in Table1and in Figure4a (N = 1024 × 1024).

		FISTA-Spa FISTA-PC ADMM-PC FISTA-W FISTA-WP FISTA-WP (GPU)
	# iterations	84	80	32	84	38	38
	Time (s)	58.3	44.4	37.6	3.6	1.9	0.087
	Speed-up	-	1.3	1.5	16.2	30.7	670
		FISTA-Spa FISTA-PC ADMM-PC FISTA-W FISTA-WP FISTA-WP (GPU)
	# iterations	58	54	34	54	28	28
	Time (s)	616	119.1	172.6	9.9	5.1	0.164
	Speed-up	-	5.2	3.6	62	120	4089

Table 3 :

 3 Timings and iterations number to reach criterion 4.3 for the blur in Figure

Acknowledgments

P. Weiss is supported by the ANR JCJC Optimization on Measures Spaces ANR-17-CE23-0013-01 and the ANR-3IA Artificial and Natural Intelligence Toulouse Institute. Both authors thank the GDR ISIS its support.

We are now ready to bound the complexity of step 1. Proof. The i-th column of V is supported on {ω -1 (i)}. Using Lemma 5.4 together with J = 1 d log 2 N give the result. The cost of step 2 is controlled by the following Lemma.

Lemma 5.6. The computation of B k from V k requires at most

Proof. To bound the number of operations needed to compute B from V , we need to know the support of each row of V . Let λ = (j, e, l). From Lemma 5.4, the necessary condition for V [λ, i] to be non zero is

Equivalently the condition becomes

. Using again Lemma 5.4 we conclude that the number of operations required to compute B from V is bounded above by

where we used J = 1 d log 2 N to conclude the proof. Overall, the computation of one matrix B k can be achieved in O δ d+1 N log 2 N operations. As seen in this two-step algorithm, the crucial ingredient of its efficiency is the sparse cascade algorithm. Its implementation is very similar to the standard algorithm but loops on the non-zero elements of the input signal in E. Furthermore, the indexes of each non-zero coefficient at the next scale have to be tracked and reused at the next filtering stage.

Products A k B k

Two questions are important while investigating the product A k B k . First, how many operations are needed for its computation and how many non-zero coefficients are in the product matrix. These questions are answered in the following Lemma. The product A k B k could be implemented by building the matrix A k from the set of circulant vectors and then rely on a standard implementation of sparse matrix-matrix products. However, the building step could be memory-intensive and long time-wise. This step is avoided by taking advantage of the sub-band circulant structure of A k : locations of non-zero coefficients can be predicted from the elements of the circulant vectors.

Before proving Lemma 5.7, we present a useful Lemma.

Lemma 5.8 ( [START_REF] Gustavson | Two fast algorithms for sparse matrices: Multiplication and permuted transposition[END_REF][START_REF] Yuster | Fast sparse matrix multiplication[END_REF]). Let A and B two N × N matrices. Let α i be the number of nonzero entries in the i-th column of A and β i be the number of non-zeros entries in the i-th row of B. Then, the number of operations required to perform AB is only N i=1 α i β i . Furthermore, the number of non-zero coefficients in AB is bounded above by the same quantity.

We are now ready to prove Lemma 5.7.

Proof of Lemma 5.7. From Theorem 5.1 we know that the number of non-zero coefficients in the λ-th column of A is bounded by

Similarly, Lemma 5.2 shows that the number of coefficients in the λ-th row of B is bounded by

Therefore the total number of operations is bounded above by:

Furthermore, L is of the order

Accumulation of A k B k

The last step of the algorithm consists in summing all the products C k = A k B k into Θ. The following results show that the number of coefficients in Θ does not explode while adding the C k and that the returned Θ k reaches the prescribed accuracy.

A Deblurring algorithms

We briefly describe the deblurring algorithms used in the numerical experiments.

A.1 FISTA

Let Ψ : R N → R N denote an orthogonal wavelet transform and H : R N → R N denote a linear operator. We then have arg min

Hff 0


MS Augustin Belkadi

MS Céline Dietrich

MSc François Machavoine

Jefferson Victor

PhD Maria Leite-De-Moraes
email: maria.leite-de-moraes@parisdecartes.fr

PhD Jefferson Russo Victor

γδ T cells amplify Blomia tropicalis -induced allergic airway disease

Keywords: Blomia tropicalis, gd T cells, atopy, asthma, airway hyperreactivity, lung, IL-4, patients

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

To the Editor:

Allergic asthma is a chronic and heterogeneous immunological disease characterized by airway inflammation, hyper-IgE production, mucus hypersecretion and airway hyperreactivity (AHR). The most studied house dust mites (HDM) respiratory allergens belong to Dermatophagoides pteronyssuns and Dermatophagoides farinae mite species. However, the storage mite Blomia tropicalis represents an important allergenic source and sensitization to its allergens are commonly associated with typical clinical manifestations of allergy, including asthma (1)(2)(3). Although these mites may coexist in tropical and subtropical regions, B. tropicalis is the major allergen in some countries (4,5). A limited number of reports have focused on the allergic processes involved in the immune responses triggered by B. tropicalis (3,6,7). Here, we used allergic asthma murine models and ex vivo analysis of B. tropicalis allergic subjects to better explore the role of immune cells in allergic responses obtained in the context of B. tropicalis exposure.

Mice received B. tropicalis allergens exclusively administered intranasally (without adjuvant) for both sensitization and challenge phases (Appendix S1). This model reproduced the major allergic asthma hallmarks, namely enhanced BALF cell counts accompanied by airway eosinophilia and airway hyperresponsiveness (AHR) (Figure 1A andB). We addressed the involvement of IL-4, as a major Th2 cytokine by comparing lung inflammation in IL-4-deficient (IL-4 -/-) and wild-type (WT) mice. The lack of IL-4 resulted in the absence of the typical airway eosinophilia and AHR observed in WT mice following B. tropicalis sensitization and challenge (Figure S1A andB). These results show that IL-4 is required for B. tropicalis allergen-induced asthma in our model. Distinct immune cells, such as iNKT cells and gd T cells, may be capable to secrete IL-4 and favour asthmatic inflammatory Th2-cytokine-dependent responses [START_REF] Lisbonne | Cutting edge: invariant V alpha 14 NKT cells are required for allergen-induced airway inflammation and hyperreactivity in an experimental asthma model[END_REF][START_REF] Hahn | Different potentials of gamma delta T cell subsets in regulating airway responsiveness: V gamma 1+ cells, but not V gamma 4+ cells, promote airway hyperreactivity, Th2 cytokines, and airway inflammation[END_REF]. As previously described, iNKT cells can be divided into three major subsets, iNKT1 (IFN-g and IL-4 producers), iNKT2 (IL-4 producers) and iNKT17 (IL-17 producers) [START_REF] Michel | Identification of an IL-17-producing NK1.1(neg) iNKT cell population involved in airway neutrophilia[END_REF]. These subsets were present in the lung of B. tropicalis treated mice (Figure S2A). No major differences were observed concerning their ability to secrete IL-4, IL-13 or IFNg but the frequency of IL-17-producing iNKT cells was enhanced when compared to controls (Figure S2A). However, this enhanced frequency was not critical for airway inflammation since iNKT cell-deficient mice presented similar airway eosinophilia than WT mice (Figure S2B).

Concerning gd T cells, including the Vg1.1 + (hereafter mentioned as Vg1 + ) gd T cell subset, we found that their number were higher in the lung of B. tropicalis-treated than in NaCl-treated mice (Figure S3). It is noteworthy that airway eosinophilia and AHR were significantly reduced in B. tropicalis-sensitized and -challenged TCRd -/-mice when compared to WT controls (Figure 1C, 1D). Further, IgE levels in the serum, Mucin-5 mRNA expression and mucus deposition in the lung were also impaired in TCRd -/-compared to WT mice (Figure S4). Previous reports have shown that the influence of gd T cells on allergic responses may be complex [START_REF] Hahn | Different potentials of gamma delta T cell subsets in regulating airway responsiveness: V gamma 1+ cells, but not V gamma 4+ cells, promote airway hyperreactivity, Th2 cytokines, and airway inflammation[END_REF]. Here, we focused our analysis on Vg1 + and Vg1 -gd T cell subsets (Figure S5A). Among gd T cells, IL-4 and IL-13 were mainly produced by the Vg1 + gd T cell subset in the lungs of the asthmatic mice (Figure S5B andC). Of note, the frequency of IL-4 + and IL-13 + among gated Vg1 + gd T cells from B. tropicalis-treated mice was higher than in controls (Figure S5B). These cells also secreted higher levels of IFNg and IL-17 in asthmatic mice (Figure S5B). In contrast, the percentage of IL-4 + , IL-13 + , IL-17 + and IFNg + among gated Vg1 -subset was similar in mice treated with B. tropicalis and controls (Figure S5B andC). These findings indicating that cytokine-producing Vg1 + gd T cells may be critical for this inflammatory allergic response.

The enhanced frequency of IL-4-producing Vg1 + gd T cells in the lung of B. tropicalistreated mice (Figure S5B) led us to the question whether airway eosinophilia and AHR could be restored by transferring fully competent Vg1 + gd T cells to B. tropicalis-sensitized andchallenged TCRd -/-mice. This was the case since the adoptive transfer of 25 000 Vg1 + gd T cells from WT mice one hour before the first challenge restored airway inflammation, particularly for airway eosinophilia, and AHR in TCRd -/-mice (Figure 1E andF). In contrast, adoptive transfer of Vg1 + gd T cells from IL-4 -/-mice had no effect (Figure 1E andF). The failure to develop an effective AHR response in TCRd -/-mice and the inability of Vg1 + gd T cells from IL-4 -/-mice to restore AHR can be ascribed to the lack of IL-4 production by Vg1 + gd T cells, which compromises their contribution to asthmatic symptoms. Thus, IL-4producing Vg1 + gd T cells are required for airway inflammation and AHR in our B. tropicalis asthma model.

We further addressed the relevance of gd T cells to atopic immune process observed in B. tropicalis allergic subjects. We obtained blood samples from subjects clinically classified as atopic or non-atopic (no allergic responses to any tested allergen) to B. tropicalis allergens (Table S1). The frequency of peripheral blood gd T cells was significantly lower in B. tropicalis atopic than non-atopic individuals (Figure 2A andS6). However, these cells were more activated in atopic subjects, as they spontaneously produced higher levels of IFNg and IL-4 than those of non-atopic donors (Figure 2B, 2C and S6). The frequency of IL-17producing gd T cells was similar in atopic and non-atopic subjects (Figure 2D andS6).

Overall, these observations suggest that gd T cells may be involved in immune allergic processes in B. tropicalis atopic patients. Further studies are required to determine the major blood gd T cell subsets activated in these patients.

This study has limitations since the precise nature of the stimuli (B. tropicalis antigens or cytokines) required to induce IL-4 production by gd T cells and the analysis of these lymphocytes in BALF from subjects susceptible to B. tropicalis remain to be determined. On-line repository

METHODS

Mice and reagents

Eight to ten-week-old specific pathogen-free C57BL/6J, IL-4 -/-and TCRd -/-mice were bred in our facility. All animal experiments were carried out according to the guidelines for care and use of animals approved by the French Institutional Committee (APAFIS#4105-

201511171831592).

Airway allergen sensitization and challenge model.

To induce allergic airway inflammation, we adapted a previously described protocol (1). In brief, mice were sensitized and challenged by intranasal administrations of B. tropicalis extracts (Greer laboratories, USA) at day 0 (D0) with 100µg/mouse then at D7, D9, D11, D14

and D16 with 50µg/mouse.

Twenty-four hours after the last challenge, mice were anesthetized with a mixture of ketamine (150 mg/kg) and xylasine (400 µg/kg) and their tracheas were cannulated (tracheostomy with ligation). FlexiVent apparatus (SCIREQ) was used to access airwayspecific resistance (Rn, tidal volume of 10 ml/kg at a respiratory rate of 150 breath/min in response to increasing doses of aerosolized acetyl-b-methylcholine chloride (methacholine;

Sigma-Aldrich). Assessments were performed at least three times and the maximum R value obtained after each dose of methacholine was used for the measure.

Airway inflammation was assessed in cytospin preparations of cells in bronchoalveolar lavage fluid (BALF, 3 x 0.5 mL washes with PBS) that were stained with May-Grünwald/Giemsa (Merck). For some experiments, BALF cells were also analysed by flow cytometry.

Adoptive transfer

In some experiments, 25 x 10 3 electronically sorted Vg1 + or Vg1 -gd T cells were adoptively transferred intravenously (i.v.) to B. tropicalis immunized mice 1h before the first challenge (D7). Mice were then further challenged at at D7, D9, D11, D14 and D15 with 50µg/mouse, as described above.

Leucocytes from lung tissues

Lung tissues were cut into pieces using a GentleMACS Dissociator (Miltenyi Biotec) and treated with collagenase type 4 (Thermo Fischer Scientific) plus DNAse I (Roche). The lymphocyte-enriched fraction was collected at the 35-70% interface of Percoll gradients (GE Healthcare). Cells were immediately stained or stimulated for 4h with 10 -8 M PMA and 1 μg ml -1 ionomycin, in the presence of 10 μg ml -1 brefeldin A (all from Sigma-Aldrich).

Immunofluorescence.

Murine cells were incubated with CD1d-PBS57-APC tetramers, anti-TCRb-APCeFluor780, anti-Vg1.1-PE and anti-TCRgd-BV605 (Biolegend). Fixable viability dye was used to exclude dead cells (ThermoFischer Scientific). For intracellular staining, cells were further fixed with 4% PFA, washed, and permeabilized with 0.5% saponin (Sigma-Aldrich), and then incubated with anti-IL-17-PerCPCy5.5, anti-IL-4-PE-Cy7, anti-IFN-g-V450, anti-anti-IL-13-PE-eF610

or isotype control (eBiosciences). The cells were washed and fluorescence was detected using a LSRFortessa (Becton Dickinson) and further analysed using the FlowJo 10.4.1 software (Tree Star).

gd T cell sorting.

Splenocytes from WT or IL-4 -/-mice, isolated as previously described (2, 3), were labelled with anti-Vg1-PE and anti-TCRgd-BV605, as described above, and then sorted using a FACSAria cell sorter (Becton Dickinson).

mRNA expression.

RNAs were extracted using the RNeasy Plus Minikit (Qiagen) including a DNase treatment.

Then RNA was reverse transcribed using the High Capacity RNA-to-cDNA Kit (ThermoFisher Scientific), according to the manufacturer's instructions. Primers and probes for real-time PCR were provided by ThermoFisher Scientific under references: Mucin 5b (Muc5b): Mm00466391_m1 and HPRT: Mm01545399_m1. All reactions were performed in triplicate with TaqMan® Fast Advanced Master Mix according to the supplier's instructions for a Step One Plus apparatus (ThermoFisher Scientific). All data were normalized to the internal standard, namely HPRT expression in each sample, and expressed as relative expression using the DDCt method versus the reference sample.

Histology of the airways.

Lungs were fixed with 10% formalin via the trachea, removed and stored in 10% formalin.

Lung tissues were embedded into paraffin and 3µm sections were stained with periodic acid Schiff (PAS).

Patient samples

Blood samples were collected from subjects who were previously clinically classified as atopic or non-atopic individuals and voluntarily subjected to a skin prick test (SPT) to confirm the atopic state. These individuals were classified into two groups: as Bt-atopic individuals (clinically allergic and reactive to Blomia tropicalis extract, with co-reaction to at least one allergen), and non-atopic individuals (without any clinical allergy symptoms and not reactive to any tested allergen). Additional information about these individuals is shown in table 1.

Each sample of peripheral blood mononuclear cells (PBMCs) was provided from a different donor. The ethics committees at the HCor and the School of Medicine at the University of São Paulo approved this study (CAAE: 15507613.4.0000.0060).

Skin prick test (SPT) and blood sample collection

The SPT were performed in accordance with European standards (4) with an adapted panel of allergens that included the profile of Brazilian allergens (i.e., Blomia tropicalis, Canis familiaris, Periplaneta americana, Aspergillus fumigatus, Penicillium notatum, Alternaria alternata, Cladosporium herbarum, Dermatophagoides pteronyssinus, Dermatophagoides farinae, and Felis domesticus -IPI ASAC), as previously described (5).

Briefly, one drop of each allergen extract or controls was applied to the volar forearm and a superficial skin puncture was made using a hypodermic needle (Alko, Brazil). After 15 minutes, the results were considered positive when wheals measure reach a diameter of 3 mm greater than that of the negative control. We excluded patients who used antihistamines, glucocorticosteroids or certain other systemic drugs that can influence the SPT results within 15 days before the test or with severe eczema or dermographism.

Human flow cytometry

Flow cytometry using human PBMCs to evaluate phenotype and spontaneous intracellular cytokine production were performed as previously described (6,7). Briefly, PBMC separations were performed using Ficoll-Paque Plus (GE Healthcare) after centrifugation and suspensions of PBMCs were phenotypic evaluated ex vivo or cultured to intracellular cytokine production evaluation. PBMCs were stained with mouse anti-human gdTCR-FITC, CD161-PECy5 or isotype control antibodies (BD Pharmingen).

To evaluate intracellular cytokine production, PBMCs were cultured in RPMI 1640 medium containing 10% FBS and 1 µg/mL of Brefeldin A (Sigma-Aldrich) for 24 hours without additional stimulus. Cultured PBMCs were then washed, fixed with formaldehyde and stained with mouse anti-human gdTCR-FITC, IL-17A-Alexa700, IL-4-PECy7, IFN-g-HorizonV450 or isotype control antibodies (BD Pharmingen) in 100 µL of PBS containing 0.05% saponin. Thirty minutes later, PBMCs were washed, fixed, acquired using a

LSRFortessa cytometer (BD Biosciences, USA), and analysis was performed using FlowJo

Figure 1 :

 1 Figure 1: gd T cells were required for asthma symptoms induced by B. tropicalis. (A, C, E) The number of total cells (Total), monocyte (Mono), eosinophil (Eosino), neutrophil (Neutro) and lymphocyte (Lympho) were determined in BALF of wild-type (WT) or TCRd -/-asthmatic (B. tropicalis) or non-asthmatic (NaCl) mice (n=7 to 20). (B, D, F) Lung resistance (Rtr) was measured 24h after the last B. tropicalis challenge or controls (NaCl) (n=7 to 15). (E, F) BALF cell counts (E) and lung resistance (Rtr) were measured 24h after the last B. tropicalis challenge of WT, TCRd -/-, TCRd -/-mice adoptively transferred with sorted Vg1 + gd T cells from WT or IL-4 -/-mice (n=5 to 10). * p<0.05; **p<0.01; ****p<0.0001.Figure 2. gd T cells were activated in B. tropicalis atopic patients. (A) Percentages of gd T cells in PBMC from non-atopic versus B. tropicalis atopic subjects. (B, C, D) Frequency of spontaneous IFNg + (B), IL-4 + (C) and IL-17A + (D) cells among gated gd T cells from non-atopic versus B. tropicalis atopic subjects.

Figure 2 .

 2 Figure 1: gd T cells were required for asthma symptoms induced by B. tropicalis. (A, C, E) The number of total cells (Total), monocyte (Mono), eosinophil (Eosino), neutrophil (Neutro) and lymphocyte (Lympho) were determined in BALF of wild-type (WT) or TCRd -/-asthmatic (B. tropicalis) or non-asthmatic (NaCl) mice (n=7 to 20). (B, D, F) Lung resistance (Rtr) was measured 24h after the last B. tropicalis challenge or controls (NaCl) (n=7 to 15). (E, F) BALF cell counts (E) and lung resistance (Rtr) were measured 24h after the last B. tropicalis challenge of WT, TCRd -/-, TCRd -/-mice adoptively transferred with sorted Vg1 + gd T cells from WT or IL-4 -/-mice (n=5 to 10). * p<0.05; **p<0.01; ****p<0.0001.Figure 2. gd T cells were activated in B. tropicalis atopic patients. (A) Percentages of gd T cells in PBMC from non-atopic versus B. tropicalis atopic subjects. (B, C, D) Frequency of spontaneous IFNg + (B), IL-4 + (C) and IL-17A + (D) cells among gated gd T cells from non-atopic versus B. tropicalis atopic subjects.

software 10 . 1 (LEGENDSFigure S1 :

 101S1 Figure S1: IL-4-producing cells were required for asthma symptoms induced by B. tropicalis. (A) The number of total cells (Total), monocyte (Mono), eosinophil (Eosino), neutrophil (Neutro) and lymphocyte (Lympho) were determined in BALF of wild-type (WT) or IL-4 -/- asthmatic (B. tropicalis) mice. (B) Lung resistance (Rtr) was measured 24h after the last B. tropicalis challenge or controls (NaCl) (n=7 to 15).

Figure S2 :

 S2 Figure S2: iNKT cells were not required for asthma symptoms induced by B. tropicalis. (A) A representative Facs profile showing the gating strategy used to identify iNKT (CD1d tetramer + TCRb + ) cells and the percentage of IFNg, IL-4, IL-13 and IL-17A positive cells, and their respective isotype controls, among gated iNKT cells in the lungs of WT asthmatic (B. tropicalis) or non-asthmatic (NaCl) mice. (B) The number of total cells (Total), monocyte (Mono), eosinophil (Eosino), neutrophil (Neutro) and lymphocyte (Lympho) were determined in BALF of wild-type (WT) or iNKT (Ja18 -/-)-deficient asthmatic (B. tropicalis) mice. (n=5 to 7).

Figure S3 :

 S3 Figure S3: gd T cells and the Vg1 + gd T cell subset were enhanced in the lung of B. tropicalis asthmatic mice. (A) A representative Facs profile showing the gating strategy used to identify TCRgd + and Vg1 + gd T cells. (B) The number of total gd T cells and of the Vg1 + gd T cell subset were determined in the lungs of wild-type asthmatic (B. tropicalis) or non-asthmatic (NaCl) mice (n=5 to 10).

Figure S4 :

 S4 Figure S4: Impaired seric IgE and mucus production in the lung of TCRd -/-B. tropicalis asthmatic mice. (A) Total IgE levels in the serum of WT and TCRd -/-mice following B. tropicalis or NaCl administration (B) Mucin-5 mRNA expression assessed by quantitative RT-PCR in the lung of WT and TCRd -/-mice following B. tropicalis or NaCl administration.

Figure S5 .

 S5 Figure S5. Distinct cytokine production by lung Vg1 + and Vg1 -gd T cell subsets. (A) Flow cytometric gating strategy used to identify TCRd + Vg1 + and TCRd + Vg1 -gd T cell subsets and the frequency of these cells in the lung of B. tropicalis-versus NaCl-treated mice. (B, C) A representative Facs profile showing the percentage of IFNg, IL-4, IL-13 and IL-17A positive cells among gated TCRd + Vg1 + (B) and TCRd + Vg1 -(C) gd T cell subsets and their frequency in the lungs of WT asthmatic (B. tropicalis) or non-asthmatic (NaCl) mice. * p< 0.05. (D) Facs profile represents the respective IL-13 versus IL-4 and IL-17A versus IFN-g isotype controls.

Figure S6 .

 S6 Figure S6. gd T cells were activated in B. tropicalis atopic patients. (A) Gate strategy to obtain viable lymphocytes. (B) A representative Facs profile showing the gating strategy used to identify TCRgd + cells and the percentage of IFNg, IL-4, IL-13 and IL-17A positive cells among gated TCRgd + cells.

Table S1 :

 S1 Characteristics of the adults included in the study

Non-atopic Bt-Atopic p

  Augustin Belkadi, MS a, *, Céline Dietrich, MS a, *, François Machavoine, MSc a , Jefferson Russo Victor, PhD b,c,d and Maria Leite-de-Moraes, PhD a

	Number	15	13
	Age, years (mean ± SE)	30.8±2.1	29.6±2.8 0.41
	Sex (male/female)	7/8	7/6	0.83
	SPT results (n/%)
	Blomia tropicalis	0/0	13/100
	Dermatophagoides pteronyssinus	0/0	9/69
	Dermatophagoides farinae	0/0	8/61
	Aspergillus fumigatus	0/0	5/38
	Canis familiaris	0/0	5/38
	Alternaria alternata	0/0	4/30
	Penicillium notatum	0/0	4/30
	Felis domesticus	0/0	2/15
	Cladosporium herbarum	0/0	1/8
	Periplaneta americana	0/0	0/0
	Others
	Clinical manifestation (past 6 mo)	0/0	13/100
	Regular use of antihistamines	0/0	13/100
	SPT -Skin prick test; mo -months

We are indebted to the technical assistance of E. Panafieu and Rachel Rignaud and to the help of Clélia Comte with the histology figures. Our thanks go also to the Cytology and Histology platforms of the INEM. This work was supported by grants from CNRS (Centre National de la Recherche Scientifique), INSERM (Institut National de la Santé et de la Descartes and Legs Poix (Chancellerie


Andriniaina Harry Solo
email: andriniaina.harry@laplace-univ.tlse.fr

Malyk Benmouffok

Pierre Freton

Jean Jacques Gonzalez

STOICHIOMETRY AIR -CH4 MIXTURE: COMPOSITION, THERMODYNAMIC PROPERTIES AND TRANSPORT COEFFICIENTS

Keywords: air-CH 4, plasma composition, constant pressure, constant mass density, thermodynamic properties, transports coefficients

HAL is

Introduction

The CH 4 -air gas mixture is generally used in the ignition studies of thermal engines [START_REF] Maly | Initiation and propagation of flame fronts in lean methane/air mixtures by the three modes of the ignition spark[END_REF][START_REF] Kravchik | Numerical modeling of spark ignition and flame initiation in a quiescent methane-air mixture[END_REF][START_REF] Zaepffel | Experimental study of an electrical discharge used in reactive media ignition[END_REF]. The used mixture proportions change according to the literature works, and are sometimes subject of parametric studies [START_REF] Maly | Initiation and propagation of flame fronts in lean methane/air mixtures by the three modes of the ignition spark[END_REF][START_REF] Zaepffel | Experimental study of an electrical discharge used in reactive media ignition[END_REF]. However, it is difficult to find the data for the desired conditions (temperature range, pressure, or mixture proportions). So mixing laws are sometimes used, whereas they are inappropriate for reactive gases. Other authors consider the proportion of methane in the mixture to be small (5-10%) and use air data [START_REF] Kravchik | Numerical modeling of spark ignition and flame initiation in a quiescent methane-air mixture[END_REF][START_REF] Zaepffel | Etude expérimentale et numérique d'une décharge électrique appliquée à l'allumage d'un milieu réactif[END_REF]. These assumptions may lead to erroneous results, or at least far from real behavior. Our goal in this paper consists to determine a database for a mixture of 90.5% air and 9.5% CH 4 in order to be able to implement in the future a plasma modelling after the spark phase of an ignition process. We present our calculation methods and the equation system under the hypothesis of local thermodynamic equilibrium (LTE). The results of the air-methane mixture will be compared with those of air, and with results from the literature.

Theory

Partition function

The determination of the data need generally the internal partition functions Q int of species and their derivatives. They are calculated from equations (1 or 2). The input parameters used come from the NIST [START_REF] Kramida | Nist atomic spectra database (ver. 5.6.1)[END_REF] and JANAF [START_REF] Chase | Nist-janaf thermochemical tables (ver. 1.0)[END_REF] sites. These equations are solved for a temperature step of 1K and the derivatives passed in a spline to avoid oscillations in thermodynamic properties as well as transport coefficients.

Atomic Species

Q int (T ) = imax i g i .e - E i k B .T (1)
with g i the degeneration of the electronic level i given by g i = 2.J + 1, where J is the quantum number of angular moment, k B the Boltzmann constant, T the temperature, E i the energy of electronic level, and i max the limit of the summation such as:

E imax ≤ E ion -E (2)
E ion is the ionization energy and E the lowering of ionization potential [START_REF] Harry Solo | The virial effect -application for SF6andCH4 thermal plasmas[END_REF].

Molecular Species vol.

no. / air -CH 4 mixture: Transport and thermodynamic properties

Laws Equations

Mass action

N i=1 n vi i = N i=1 Q vi tot,i (4)
Conservation of atomic nucleus

j N i n i .C i,k = k N i n i .C i,j (5)
Conservation of electrical neutrality

N i Z i .n i = 0 (6) Dalton P = N i n i .k B .T (7)
Conservation of mass density ρ

= N i n i .m i ( 8
)
Table 1. Basis equations used to determine the chemical composition.

Q int (T ) = P r k B .T h 2 2.π.m.k B .T 3 2 exp 1 N a .k B [H(0) -H(T r )] JAN AF T - G -H(T r ) T JAN AF (3)
with P r and T r respectively the pressure and reference temperature in the JANAF table [START_REF] Chase | Nist-janaf thermochemical tables (ver. 1.0)[END_REF], h the Planck constant, m the mass of the species, Na the Avogadro number, H(0) -H(T r ) the enthalpy evaluated at zero kelvin, and [G -H(T r )]/T the Gibbs energy dependent of the temperature.

Chemical composition

The calculation of chemical composition was developed using the mass action law including the chemical base concept [START_REF] Godin | A robust and efficient method for the computation of equilibrium composition in gaseous mixtures[END_REF]. This method allows us to generalize the equations for different types of pure gases and gas mixtures. If calculation is at constant pressure, Dalton's law is used with the corrections of Debye-Hückel δP Debye and virial δP virial [START_REF] Harry Solo | The virial effect -application for SF6andCH4 thermal plasmas[END_REF][START_REF] Harry Solo | Compositions chimiques et propriétés thermodynamiques à l'etl d'un mélange air -CH4[END_REF]. For a calculation at constant density, the mass density equation is used. The basic equations for the determination of chemical composition are summarized in Table 1.

In Table 1 , N represents the total number of chemical species,n i the density of the species i, v i the stoichiometry coefficients of the reaction,Q vol (tot,i) the total partition function by volume [START_REF] Harry Solo | The virial effect -application for SF6andCH4 thermal plasmas[END_REF], the atomic proportion of the species j or k, C i,k et C i,j the elements of the composition matrix [START_REF] Harry Solo | The virial effect -application for SF6andCH4 thermal plasmas[END_REF], and Z i the charge of the species.

Thermodynamic properties

By knowing the densities of each species, as well as their internal partition functions, the thermodynamic properties can be calculated (such as mass density ρ, enthalpy H or specific heat at constant pressure C p ). These properties are generally presented for a given pressure [START_REF] Harry Solo | Compositions chimiques et propriétés thermodynamiques à l'etl d'un mélange air -CH4[END_REF]. Here we will look at internal energy U and specific heat at constant volume C v . These data are used particularly in variable pressure and constant volume systems.

Their equations are illustrated in Table 2 with N i is the particle density of the i species given by

N i = ni ρ
, and

E ref i
the reference energy.

Thermodynamic properties Equations

Internal energy

U = 3 2 .k B .T. N i N i + k B .T 2 . N i N i . δQ int i δT + N i N i .E ref i (9) Specific heat at constant volume C v = δU δT V ( 10
)

Transport coefficients

The conventional formulations of the literature based on approximation of Chapman-Enskog method are used for the determination of transport coefficients [START_REF] André | Transport coefficients of plasmas consisting of insulator vapours: Application to pe, pom, pmma pa66 and pc[END_REF][START_REF] Cressault | Thermodynamic properties and transport coefficients in Ar -H2 -Cu plasmas[END_REF]. The equations to be solved require the knowledge of additional input data that are the collision integrals. In our code, they are calculated in first approximation from the hard sphere potentials. However, following available data, the appropriate formalism according to the type of collision of the species (Table 3) is considered.

Case of Charged-Charged interaction: The Coulomb potential limited by the Debye's length is used. It depends of the species charge, its temperature and the Debye's length. The collision integrals are so calculated from Mason [START_REF] Mason | Transport properties of gases obeying a modified buckingham (exp-six) potential[END_REF] and Devoto [START_REF] Devoto | Transport coefficients of partially ionized argon[END_REF][START_REF] Devoto | Simplified expressions for the transport properties of ionized monatomic gases[END_REF][START_REF] Devoto | Transport coefficients of ionized argon[END_REF] tables.

Case of Neutral-Neutral interaction: Lennard-Jones potential (atoms or molecules). The reduced temperature T reduced is used. This formalism depends on the maximal energy attraction. Tabulated values for the Lennard-Jones potential (12-6) are used.

Case of Neutral-Charged: Electron/Neutral and Charged/Neutral (Input: charge, polarity and T reduced ).

The parameters used to calculate the collision integrals come from Capitelli et al [START_REF] Capitelli | Collision integrals of high temperature air species[END_REF] for air, and from Bar et al [START_REF] Baronnet | Transport coefficient of hydrogen-methane thermal plasma[END_REF] as well as Sanon [START_REF] Sanon | Contribution au calcul des propriétés thermodynamiques et des coefficients de transport de plasmas thermiques de mélangesAr[END_REF][START_REF] Sanon | Transport coefficients of Ar -C -H -O -N systems thermal plasma at atmospheric pressure[END_REF] for methane.

Results and Discussion

We present the results obtained for a mixture of 90.5% air and 9.5% CH 4 for a temperature range from 300K to 30kK. 51 species were considered in our calculations:

N , N + , N +2 , N 2 , N + 2 , N O, N O + , N 2 O, N H, CN , CN O, O, O + , O +2 , O -, O 2 , O + 2 , O - 2 , OH, OH + , CO 2 , HO 2 , N O 2 , C , C + , C +2 , C 2 , CO, CO + , CH, CHN , CHO, C 2 H, C 2 H 2 , H, H + , H -, H 2 , H 2 O, H 2 O 2 , HN O, N H 2 , N H 3 , CH 2 , CH 3 , CH 4 , CH 2 O, Ar, Ar + , Ar +2 , e -.
A comparison is made between the thermodynamic properties and the transport coefficients of the air-methane mixture and the pure air. For validation, we compare our results with those available in the literature.

Chemical composition

The chemical compositions for a pressure of 1 bar and 10 bars and at a density of 0.1481 kg.m -3 and 1.111kg.m -3 are shown respectively in Figure 1 and Figure 2. At a given pressure ( Figure 1), we observe that neutral species dominate the medium at low temperature (T<10kk). The densities of molecular species gradually decrease with temperature. Then these molecules dissociate to give neutral atoms. The various processes of electron formation, recombination and ionization of species then occur. At high temperatures the ionized species become predominant and neutrals negligible. Increasing pressure changes species behavior. The processes of species transformation are then shifted to higher temperatures. Chemical compositions (Figure 2) with constant density have characteristics similar to those at a given pressure. The species N 2 , H 2 O and CO 2 are the most predominant at low temperature. Increasing the temperature, they are replaced by neutral atoms, then electrons as well as ionized species. Species densities increase with the value of the mass density. For a given value of mass density, densities of dominant species vary slightly with the temperature. Variations occur only when a certain energy is reached (dissociation or ionization temperature) allowing the species to dissociate or ionize according to the type.

Thermodynamic properties

In Figures 3 and4 the internal energy and the constant volume specific heat of the air-methane mixture and air are compared. Figure 3 shows an increasing evolution of internal energy with temperature. Small waves are observed that correspond to the various chemical reactions occurring. The value of internal energy is inversely proportional to the mass density. For the gases compared here, above around 4kK, the mixture of 90.5% air and 9.5% CH 4 allows a better heat exchange between the system and the external medium compared to pure air because its internal energy is more predominant. With regard to Figure 4 representing specific heat at constant volume, three main peaks, more or less pronounced are present. They correspond to the dissociation of the molecule O 2 , then that of N 2 , and the ionization of N that we find on the case of air. The peaks that correspond to the different processes move towards higher temperatures with the increase in mass density as we found on chemical compositions (Figure 2). For 9.5% of methane in air, the results show significant differences compared to pure air.

Transport coefficients

Figures 5, 6 and 7 respectively show the thermal conductivity, the electrical conductivity and the viscosity of 90.5% air and 9.5% CH 4 mixture and the pure air at a pressure of 1 and 10 bars. In order to validate our developments, we compare these data with the available literature results.

Figure 5 describes the evolution of thermal conductivity. They result from the contribution of three terms including internal thermal conductivity, translation and reaction. Below 20kK, the medium is dominated mainly by thermal conductivity of reaction, and by small contributions of the internal component as well as the translation of heavy particles. This explains the presence of peaks series similar to the specific heat curves reflecting the chemical reactions that occur. Through the formation of electrons and ionization processes at high-temperature, neutral species begin to disappear. From there, the electron translation component represents an essential part of the total thermal conductivity. We also notice in this case the dependence of the thermal conductivity to the pressure. The peaks also shift to the high temperatures and their amplitudes decreases when the pressure increases. However, when the electronic translation component becomes dominant, the phenomenon reverses and the thermal conductivity increases with the pressure. In comparison with air, the air -CH 4 mixture allows overall a better ability to drive heat. For electrical conductivity, Figure 6 compares our work with two literature results. There is a strong increase for temperatures up to 10kK, followed by a quasi-linear variation and finished with a constant profile where the electron density begins to stabilize as can be seen in the case at 1 bar above 25kK (Figure 6 left). The electrical conductivity depends heavily on the electronic density, the differences observed between the two gases studied are in principle at the level of the strong gradient (T<10kk). In this area of temperature, electrons vary as a result of ionization phenomena such as molecule NO. In the case of the air-CH 4 mixture, the formation of this molecule results from several species, and as the proportion of air is reduced in the mixture, its quantity in the medium also decreases. As a result, the ionization of NO is delayed compared to that of air, hence a relatively low electronic density.

In Figure 7 the viscosities of 90.5% air and 9.5% CH 4 mixture, as well as the air at a pressure of 1 and 10 bars are confronted with the literature. We find the typical bell shape, with a maximum who reflects the transition from neutral or weakly ionized plasma to a fully ionized medium. As the molar fraction of the species is involved in the viscosity equation, it also depends on the pressure. We then see an increase values and the movement of the maximum toward higher temperatures when the pressure increases. For the gap observed between the two gases presented here, they can be explained by the difference in the values of molar masses and the number of collisions between particles, especially neutrals. For the transport coefficients compared with the literature presented in this work, we find generally a good agreement. The differences may be due to the lack of consideration of some species or the choice of collision integrals.

Conclusions

The objective of this study was to determine the LTE data for the mixture of 90.5% air and 9.5% CH 4 (chemical composition, thermodynamic properties and transport coefficients). For this purpose, tools have been developed allowing access to these quantities for a large temperature range. We first showed the ability of our developments to calculate the chemical composition for a given pressure or density value.The results illustrate fairly similar behavior on the different species transformation processes between the case at P=1 or 10 bars and at ρ=0.1481 or 1.111 kg.m -3 . However, for a given mass density, species evolve differently with the particularity of keeping a linear variation up to their respective dissociation or ionization temperatures. We then presented the thermodynamic properties and transport coefficients of the air -CH 4 mixture compared with those of air. For a proportion of 9.5% of CH 4 in air, the results show an important difference in the properties of these two gases. If we want a better prediction of the real plasma behavior in the models, these differences cannot be ignored and the data of 90.5% of air and 9.5% of CH 4 mixture must be used. [START_REF] Capitelli | Transport properties of high temperature air in local thermodynamic equilibrium[END_REF][START_REF] Boulos | Thermal plasmas: Fundamentals and applications[END_REF]

Figure 1 .

 1 Figure 1. Chemical composition: 90.5% air and 9.5% CH4 -1 bar and 10 bars -LTE.

Figure 2 .

 2 Figure 2. Chemical composition: 90.5% air and 9.5% CH4 for 0.1481 and 1.111 kg.m -3 at LTE.

HAL

Figure 3 .

 3 Figure 3. Internal energy: 90.5% air and 9.5% CH4 for 0.1481 and 1.111 kg.m -3 at LTE.

Figure 4 .

 4 Figure 4. LTE Specific heat capacity at constant volume: 90.5% air and 9.5% CH4 for 0.1481 and 1.111 kg.m -3 at LTE.

Figure 5 .Figure 6 .

 56 Figure 5. Thermal conductivity: 90.5% air and 9.5% CH4 and Air for 1 and 10 bars at LTE. Comparison with literature[START_REF] Sourd | High temperature transport coefficients in e -C -H -N -Omixtures[END_REF][START_REF] Capitelli | Transport properties of high temperature air in local thermodynamic equilibrium[END_REF][START_REF] Angola | Thermodynamic and transport properties in equilibrium air plasmas in a wide pressure and temperature range[END_REF]

Table 2 .

 2 Equations to calculate thermodynamic properties.

	Species(i)/Species(j) Electron	Neutral	Charged
	Electron	Coulomb	Polarity	Coulomb
	Neutral	Reciprocity Lennard Jones Polarity
	Charged	Reciprocity Reciprocity	Coulomb

Table 3 .

 3 Used formalism according to the species collision type.