M Guerquin-Kern

M Häberlin

K P Pruessmann

M Unser

A Fast Wavelet-Based Reconstruction Method for Magnetic Resonance Imaging

Keywords: MRI, nonlinear reconstruction, undersampled spiral, non-Cartesian, sparsity, wavelets, total variation, compressed sensing, iterative shrinkage/thresholding algorithm, ISTA, FISTA, FWISTA, thresholded Landweber

In this work, we exploit the fact that wavelets can represent magnetic resonance images well, with relatively few coefficients. We use this property to improve MRI reconstructions from undersampled data with arbitrary k-space trajectories. Reconstruction is posed as an optimization problem that could be solved with the iterative shrinkage/thresholding algorithm (ISTA) which, unfortunately, converges slowly. To make the approach more practical, we propose a variant that combines recent improvements in convex optimization and that can be tuned to a given specific k-space trajectory. We present a mathematical analysis that explains the performance of the algorithms. Using simulated and in vivo data, we show that our nonlinear method is fast, as it accelerates ISTA by almost two orders of magnitude. We also show that it remains competitive with TV regularization in terms of image quality.

I. INTRODUCTION

MRI scanners provide data that are samples of the spatial Fourier transform (a.k.a. k-space) of the object under investigation. The Shannon-Nyquist sampling theory in both spatial and k-space domains suggests that the sampling density should correspond to the field of view (FOV) and that the highest sampled frequency is related to the pixel width of the reconstructed images. However, constraints in the implementation of the k-space trajectory that controls the sampling pattern (e.g., acquisition duration, scheme, smoothness of gradients) may impose locally reduced sampling densities. Insufficient sampling results in reconstructed images with increased noise and artifacts, particularly when applying gridding methods.

The common and generic approach to alleviate the reconstruction problem is to treat the task as an inverse problem [START_REF] Pruessmann | SENSE: Sensitivity encoding for fast MRI[END_REF]. In this framework, ill-posedness due to a reduced sampling density is overcome by introducing proper regularization constraints. They assume and exploit additional knowledge about the object under investigation to robustify the reconstruction.

Earlier techniques used a quadratic regularization term, leading to solutions that exhibit a linear dependence upon the measurements. Unfortunately, in the case of severe undersampling (i.e., locally low sampling density) and depending on the strength of regularization, the reconstructed images still suffer from noise propagation, blurring, ringing, or aliasing This work was supported by the Swiss National Competence Center in Biomedical Imaging (NCCBI). errors. It is well known in signal processing that the blurring of edges can be reduced via the use of nonquadratic regularization. In particular, 1 -wavelet regularization has been found to outperform classical linear algorithms such as Wiener filtering in the deconvolution task [START_REF] Figueiredo | An EM algorithm for waveletbased image restoration[END_REF].

Indicative of this trend as well is the recent advent of Compressed Sensing (CS) techniques in MRI [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF], [START_REF] Gamper | Compressed sensing in dynamic MRI[END_REF]. These let us draw two important conclusions.

• The introduction of randomness in the design of trajectories favors the attenuation of residual aliasing artifacts because they are spread incoherently over the entire image. • Nonlinear reconstructions-more precisely, 1regularization-outperform linear ones because they impose constraints that are better matched to MRI images.

Many recent works in MRI have focused on nonlinear reconstruction via Total Variation (TV) regularization, choosing finite differences as a sparsifying transform [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF], [5], [START_REF] Ying | A statistical approach to SENSE regularization with arbitrary k-space trajectories[END_REF], [START_REF] Liu | Regularized sensitivity encoding (SENSE) reconstruction using Bregman iterations[END_REF]. Nonquadratic wavelet regularization has also received some attention [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF], [START_REF] Lustig | Faster imaging with randomly perturbed, undersampled spirals and |L| 1 reconstruction[END_REF], [START_REF] Liu | Improved spiral sense reconstruction using a multiscale wavelet model[END_REF], [START_REF] Chaâri | Autocalibrated regularized parallel MRI reconstruction in the wavelet domain[END_REF], [START_REF] Guerquin-Kern | Wavelet-regularized reconstruction for rapid MRI[END_REF], but we are not aware of a study that compares the performance of TV against 1 -wavelet regularization.

Various algorithms have been recently proposed for solving general linear inverse problems subject to 1 -regularization. Some of them deal with an approximate reformulation of the 1 regularization term. This approximation facilitates reconstruction sacrificing some accuracy and introducing extra degrees of freedom that make the tuning task laborious. Instead, the iterative shrinkage/thresholding algorithm [START_REF] Figueiredo | An EM algorithm for waveletbased image restoration[END_REF], [START_REF] Bect | A 1 -unified variational framework for image restoration[END_REF], [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] (ISTA) is an elegant and nonparametric method that is mathematically proven to converge. A potential difficulty that needs to be overcome is the slow convergence of the method when the forward model is poorly conditioned (e.g., low sampling density in MRI). This has prompted research in large-scale convex optimization on ways to accelerate ISTA. The efforts so far have followed two main directions:

• generic multistep methods that exploit the result of past iterations to speed up convergence, among them: two-step iterative shrinkage/thresholding [START_REF] Bioucas-Dias | A new twist: Two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF] (TwIST), Nesterov schemes [START_REF] Nestorov | Gradient methods for minimizing composite objective function[END_REF], [START_REF] Weiss | Algorithmes rapides d'optimisation convexe. Application à la restoration d'images et à la détection de changements[END_REF], [START_REF] Becker | NESTA: A fast and accurate first-order method for sparse recovery[END_REF], fast ISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF] (FISTA), and monotonic FISTA [START_REF]Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF] (MFISTA); • methods that optimize wavelet-subband-dependent parameters w.r.t. the reconstruction problem: multilevel thresholded Landweber (MLTL) [START_REF] Vonesch | A fast thresholded Landweber algorithm for wavelet-regularized multidimensional deconvolution[END_REF], [START_REF]A fast multilevel algorithm for wavelet-regularized image restoration[END_REF] and subband adaptive ISTA (SISTA) [START_REF] Bayram | A subband adaptive iterative shrinkage/thresholding algorithm[END_REF].

In this work, we exploit the possibility of combining and tailoring the two generic types of accelerating strategies to come up with a new algorithm that can speedup the convergence of the reconstruction and that can accomodate for every given k-space trajectory. Here, we first consider single-coil reconstructions that do not use sensitivity knowledge. In a second time, we confirm the results with SENSE reconstructions [START_REF] Pruessmann | SENSE: Sensitivity encoding for fast MRI[END_REF].

We propose a practical reconstruction method that turns out to sensibly outperform linear reconstruction methods in terms of reconstruction quality, without incurring the protracted reconstruction times associated with nonlinear methods. This is a crucial step in the practical development of nonlinear algorithms for undersampled MRI, as the problem of fixing the regularization parameter is still open. We also provide a mathematical analysis that justifies our algorithm and facilitates the tuning of the underlying parameters.

This paper is structured as follows: in Section II, we describe the basic data-formation model for MRI and derive the discrete forward model. The representation of the object by wavelet bases is considered in Section II-C2; in particular, we legitimate the use of a wavelet regularization term (which promotes sparsity) to distinguish the solution from other possible candidates. In Section III, we propose a fast algorithm for solving the nonlinear reconstruction problem and present theoretical arguments to explain its superior speed of convergence. Finally, we present in Section IV an experimental protocol to validate and compare our practical method with existing ones. We focus mainly on reconstruction time and signal-to-error ratio (SER) w.r.t. the reference image.

II. MRI AS AN INVERSE PROBLEM

In this section, we present the MR acquisition model and the representation of the signal that is used to specify the reconstruction problem. We motivate the sparsity assumptions in the wavelet domain and the variational approach with 1 regularization that is used to solve the inverse problem of imaging.

A. Model of Data Formation 1) Physics: We consider MRI in two dimensions, in which case a 2D plane is excited. The time-varying magnetic gradient fields that are imposed define a trajectory in the (spatial) Fourier domain that is often referred to as k-space. We denote by k the coordinates in that domain. The excited spins, which behave as radio-frequency emitters, have their precessing frequency and phase modified depending on their positions. The modulated part of the signal received by a homogeneous coil is given by

m(k) = ρ(k) = R 2 ρ(r)e -j k , r dr. (1)
It corresponds to the Fourier transform of the spin density ρ that we refer to as object. The N measurements, concatenated in the vector m = (m 1 , . . . , m N ), correspond to sampled values of this Fourier transform at the frequency locations k n along the k-space trajectory.

2) Model for the Original Data:

x 2 ∈ R + regular quadratic norm x Λ ∈ R + weigthed quadratic norm (x H Λx) x 1 ∈ R + 1 norm ( |x i |)
a) Spatial Discretization of the Object: From here on, we consider that the Fourier domain and, in particular, the sampling points k n , are scaled to make the Nyquist sampling interval unity. This can be done without any loss of generality if the space domain is scaled accordingly. Therefore, we model the object as a linear combination of pixel-domain basis functions ϕ p that are shifted replicates of some generating function ϕ, so that

ρ = p∈Z 2 c[p]ϕ p , with (2)
ϕ p (r) = ϕ (r -p) . (3)
In MRI, the implicit choice for ϕ is often Dirac's delta. Different discretizations have been proposed, for example by Sutton et al. [START_REF] Sutton | Fast, iterative image reconstruction for MRI in the presence of field inhomogeneities[END_REF] with ϕ as a boxcar function or later by Delattre et al. [START_REF] Delattre | Splinebased variational reconstruction of variable density spiral k-space data with automatic parameter adjustment[END_REF] with B-splines. But it has not been worked out in detail how to get back the image for general ϕ that are non-interpolating, which is the case, for instance for B-splines of order greater than 1. The image to be reconstructed (i.e., the sampled version of the object ρ(p)) is obtained by filtering the coefficients c[p] with the discrete filter

P e jω = h∈Z 2 φ (ω + 2πh) , (4)
where φ denotes the Fourier transform of ϕ. b) Wavelet Discretization: In the wavelet formalism, some constraints apply on ϕ. It must be a scaling function that satisfies the properties for a multiresolution [START_REF] Unser | Wavelet theory demystified[END_REF]. In that case, the wavelets can be defined as linear combinations of the ϕ p and the object is equivalently characterized by its coefficients in the orthonormal wavelet basis. We refer to Mallat's book [START_REF] Mallat | A wavelet tour of signal processing[END_REF] for a full review on wavelets. There exists a discrete wavelet transform (DWT) that bijectively maps the coefficients c to the wavelet coefficients w that represent the same object ρ in a continuous wavelet basis. In the rest of the paper, we represent this DWT by the synthesis matrix W . Note that the matrix multiplications c = W w and w = W -1 c have efficient filterbank implementations.

B. Matrix Representation of the Model

Since a finite field of view (FOV) determines a finite number M of coefficients c[p], we handle them as a vector c, keeping the discrete coordinates p as implicit indexing. By simulating the imaging of the object (2), and by evaluating (1) for k = k n , we find that the noise-free measurements are given by

m 0 = Ec, (5)
where E, the MRI system matrix, is decomposed as

E = diag ( φ(k n )) [s 1 . . . s N ] T . (6)
There, s n is a space-domain vector such that s n [p] = e -j ωn , p . A more realistic data-formation model is

m = Ec + b, (7)
or

m = M w + b, (8)
with M = EW and a residual vector b representing the effect on measurements of noise and scanner imprecisions. The inverse problem of MRI is then to recover the M coefficients w (or c) from the N corrupted measurements m. Its degree of difficulty depends on the magnitude of the noise b and the conditioning of the matrix M (or E).

C. Variational Formulation 1) General Framework:

The solution w is defined as the minimizer of a cost function that involves two terms: the data fidelity F(b) and the regularization R(w) that penalizes undesirable solutions. This is summarized as

w = argmin w F(m -M w) + λR(w), (9)
where the regularization parameter λ ≥ 0 balances the two constraints. In MRI, b is usually assumed to be a realization of a white Gaussian process, which justifies the choice

F(b) = b 2
2 as a proper log-likelihood term. The ill-conditioning inherent to undersampled trajectories imposes the use of the suitable regularization term R(w).

Standard Tikhonov regularization corresponds to the quadratic term R(w) = Rw

= (M H M + λR H R) -1 M H m that
is tractable both theoretically and numerically. When the reconstruction problem is sufficiently overdetermined to make noise propagation negligible, regularization is dispensable and the least-squares solution, which corresponds to Tikhonov with λ = 0, is adequate. The approach is statistically optimal if the object can also be considered as a realization of a Gaussian process. Unfortunately, this assumption is hardly justified for typical MR images. The quality of the solution obtained by this means quickly breaks down when undersampling increases.

TV reconstruction is related to the sum of the Euclidean norms of the gradient of the object. In practice, it is defined as R(w) = ∇c 1 , where the operator ∇ returns pixelwise the 2 -norm of finite differences. The use of TV regularization is particularly appropriate for piecewise-constant objects such as the Shepp-Logan (SL) phantom.

2) Sparsity-Promoting Regularization: The main idea in this work is to exploit the fact that the object can be well represented by few non-zero coefficients (sparse representation) in an orthonormal basis of M functions φ p . Formally, we write that

• ∃S ⊂ R 2 , |S|
M (sparse support) and

• ∃a : ρ -p∈S a[p]φ p ρ (small error).
It is well documented that typical MRI images admit sparse representation in bases such as wavelets or block DCT [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF].

The 1 -norm is a good measure of sparsity with interesting mathematical properties (e.g., convexity). Thus, among the candidates that are consistent with the measurements, we favor a solution whose wavelet coefficients have a small 1 -norm. Specifically, the solution is formulated as

w = argmin w C(w), (10)
with

C(w) = m -M w 2 2 + λ w 1 . (11)
This is the general solution for wavelet-regularized inverse problems considered by [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] as well as many other authors.

III. WAVELET REGULARIZATION ALGORITHMS

In this section, we present reconstruction algorithms that handle constraints expressed in the wavelet domain while solving the minimization problem [START_REF] Chaâri | Autocalibrated regularized parallel MRI reconstruction in the wavelet domain[END_REF]. By introducing weighted norms instead of simple Lipschitz constants, we revisit the principle of the standard ISTA algorithm and simplify the derivation and analysis of this class of algorithms. We end up with a novel algorithm that combines different acceleration strategies and we provide a convergence analysis. Finally, we propose an adaptation of the fast algorithm to implement the random-shifting technique that is commonly used to improve results in image restoration.

A. Weighted Norms

Let us first define the weighted norm corresponding to a positive-definite symmetric weight matrix X as w 2 X = w H Xw. The requirement that the eigenvalues of X-denoted λ(X)-must be positive leads to the norm property

0 < λ min (X) w 2 2 ≤ w 2 X ≤ λ max (X) w 2 2 , ∀ w = 0.

B. Principle of ISTA Revisited

An important observation to understand ISTA is to see that the nonlinear shrinkage operation, sometimes called softthresholding, solves a minimization problem [START_REF] Chambolle | Nonlinear wavelet image processing: Variational problems, compression, and noise removal through wavelet shrinkage[END_REF], with

T λ (u) = (|u| -min(λ/2, |u|)) • sgn(u) (12) = argmin w∈C |u -w| 2 + λ|w| .

By separability of norms, this applies component-wise to vectors of C

N : T λ (u) = argmin w u -w 2 2 + λ w 1 .
This means that the 1 -regularized denoising problem (i.e., when M in [START_REF] Guerquin-Kern | Wavelet-regularized reconstruction for rapid MRI[END_REF] is the identity matrix) is precisely solved by a shrinkage operation.

The iterative shrinkage/thresholding algorithm (ISTA) [START_REF] Figueiredo | An EM algorithm for waveletbased image restoration[END_REF], [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF], a.k.a. thresholded Landweber (TL), generates a sequence of estimates w n that converges to the minimizer w of [START_REF] Guerquin-Kern | Wavelet-regularized reconstruction for rapid MRI[END_REF] when it is unique. The idea is to define at each step a new functional C (w, w n ) whose minimizer w n+1 will be the next estimate

w n+1 = argmin w C (w, w n ). (13)
Two constraints must be considered for the definition of C . (i) It is sufficient for the convergence of the algorithm that C (w, w n ) is an upper bound of C(w) and matches it at w = w n ; this guarantees that the sequence {C(w n )} is monotonically decreasing. (ii) The inner minimization [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] should be performed by a simple shrinkage operation to ensure the rapidity and accuracy of the algorithm. In accordance with Constraint (i), C can take the generic quadratically augmented form

C (w, w n ) = C(w) + w -w n 2 Λ-M H M , (14)
with the constraint that (Λ-M H M ) is positive definite, where the weighting matrix Λ plays the role of a tuning parameter. Then, ISTA corresponds to the trivial choice Λ = L 2 I, with the value of L chosen to be greater or equal to the Lipschitz constant of the gradient of M w

2 2 , so that L ≥ 2λ max M H M .
Let us define a = M H m, A = M H M , and

z n = w n + 2 L (a -Aw n ). (15)
Then, using standard linear algebra, we can write

w n+1 = argmin w w -z n 2 2 + 2λ L w 1 (16) = T 2λ L (z n ) . (17)
This shows that Constraint (ii) is automatically satisfied.

Algorithm 1: ISTA Repeat w n+1 ← T 2λ L w n + 2 L (a -Aw n ) ;
Note that both the intermediate variable z n in [START_REF] Nestorov | Gradient methods for minimizing composite objective function[END_REF] and the threshold values will vary depending on L.

Beck and Teboulle [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]Thm. 3.1] showed that this algorithm decreases the cost function in direct proportion to the number of iterations n, so that C(w n ) -C(w ) = O(1/n). Here, we present a slightly extended version of their original result which is valid for any reference point during iterations.

Proposition 1: Let {w n } be the sequence generated by Algorithm 1 with L ≥ 2λ max (A). Then, for any n > n 0 ∈ N,

C(w n ) -C(w ) ≤ L/2 n -n 0 w n0 -w 2 2 . ( 18
)
Proof: [18, Thm. 3.1] gives the result for n 0 = 0. Consider a sequence w n such that w n0 = w 0 . As the iteration does not depend on n, we get w n = w n-n0 . The result follows immediately.

Selecting L as small as possible will clearly favor the speed of convergence. It also raises the importance of a "warm" starting point.

Among the variants of ISTA, FISTA, proposed by Beck and Teboulle [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], ensures state-of-the-art convergence properties while preserving a comparable computational cost. Thanks to a controlled over-relaxation at each step, FISTA quadratically decreases the cost function, with

C(w n ) -C(w ) = O(1/n 2 ). C. Subband Adaptive ISTA (SISTA)
SISTA is an extension of the ISTA that was introduced by Bayram and Selesnick [START_REF] Bayram | A subband adaptive iterative shrinkage/thresholding algorithm[END_REF]. Here, we propose an interpretation of SISTA as a particular case of ( 14) with a weighting matrix that replaces advantageously the step size 2/L in [START_REF] Nestorov | Gradient methods for minimizing composite objective function[END_REF]. The idea is to use a diagonal weighting matrix Λ -1 = diag(τ )-if it is not diagonal, constraint (ii) would not be fulfilled-with coefficients that are constant within a wavelet subband.

1) SISTA: In the same fashion as for ISTA, ( 15) and ( 17) can be adapted to subband-dependent steps and thresholds. Accordingly, SISTA is described in Algorithm 2.

Algorithm 2: SISTA Repeat w n+1 ← T λτ w n + Λ -1 (a -Aw n ) ;
2) Convergence Analysis: By considering the weighted scalar product x H Λy instead of L 2 x H y, we can adapt the convergence proof of ISTA by Beck and Teboulle (see Proposition 1). This result is new, to the best of our knowledge.

Proposition 2: Let {w n } be the sequence generated by Algorithm 2 with (Λ -A) > 0. Then, for any n > n 0 ∈ N,

C(w n ) -C(w ) ≤ 1 n -n 0 w n0 -w 2 Λ . (19)
Proof: We rewrite the cost function [START_REF] Guerquin-Kern | Wavelet-regularized reconstruction for rapid MRI[END_REF] with the change of variable w = Λ 1/2 w. We then apply ISTA to solve that problem with a = Λ -1/2 a, A = Λ -1/2 AΛ -1/2 , and L = 2 (Note that Λ -A is positive-definite iif I -A is positivedefinite). The iteration w n+1 = T λ w n + Λ -1 (a -A w n ) can be rewritten, in terms of the original variable as w n+1 = T λτ (w n + (a -Aw n )). The latter is an iteration of SISTA (see Algorithm 2). According to Proposition 1, we have

C(Λ -1/2 w n ) -C(w ) ≤ w n 0 -Λ 1/2 w 2 n-n0
, which translates directly into the proposed result.

Therefore, by comparing Propositions 1 and 2, an improved convergence is expected. The main point is that, for a "warm" starting point w 0 or after few iterations (n 0 ), the weighted norm in [START_REF]Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF] can yield significantly smaller values than the one weighted by L 2 in (18). 3) Selection of Weights: Bayram and Selesnick [START_REF] Bayram | A subband adaptive iterative shrinkage/thresholding algorithm[END_REF] provide a method to select the values of τ for SISTA. To present this result, let us introduce some notations. We denote by s an index that scans all the S wavelet subbands, coarse scale included, by τ s the corresponding weight constant, and by M s the corresponding block of M . We also define

γ s1,s2 = λ max (M H s2 M s1 M H s1 M s2
). The authors of [START_REF] Bayram | A subband adaptive iterative shrinkage/thresholding algorithm[END_REF] show that, for each subband, the condition

1 τ s > S s =1 γ s,s (20)
is sufficient to impose the positive definiteness of (Λ-M H M ) that is required in [START_REF] Bioucas-Dias | A new twist: Two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF]. In the present context, we propose to compute the values γ s,s by using the power iteration method, once for a given wavelet family and k-space sampling strategy.

D. Best of Two Worlds: Fast Weighted ISTA (FWISTA)

Taking advantage of the ideas developed previously, we derive an algorithm that corresponds to the subband adaptive version of FISTA. In the light of the minimization problem [START_REF] Guerquin-Kern | Wavelet-regularized reconstruction for rapid MRI[END_REF], FWISTA generalizes the FISTA algorithm using a parametric weighted norm. We give its detailed description in Algorithm 3, where the modifications w.r.t. FISTA is the SISTA step in the loop.

Algorithm 3: FWISTA input : A, a, w 0 , and Λ -1 = diag(τ );

Initialization: n = 0, v 0 = w 0 , t 0 = 1; repeat w n+1 ← T λτ v n + Λ -1 (a -Av n ) (SISTA step); t n+1 ← 1 + 1 + 4t 2 n /2; v n+1 ← w n+1 + tn-1 tn+1 (w n+1 -w n ); n ← n + 1; until desired tolerance is reached; return w n+1 ;
In the same fashion as for SISTA, we revisit the convergence results of FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]Thm. 4.4] for FWISTA.

Proposition 3: Let {w n } be the sequence generated by Algorithm 3. Then, for any n ≥ 1,

C(w n ) -C(w ) ≤ 2 n + 1 2 w 0 -w 2 Λ . (21)
Proof: In the spirit of the proof of Proposition 2, we consider the change of variable w = Λ 1/2 w and apply FISTA to solve the new reconstruction problem. The step

w n+1 = T λ (v n + (a -A v n )) is equivalent to w n+1 = T λτ v n + Λ -1 (a -Av n ) .
The convergence results of FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]Thm. 4.4] applies on the sequence {w n }, which

leads to C(Λ -1/2 w n )-C(w ) ≤ 2 n+1 2 w 0 -Λ 1/2 w 2 .
This result shows the clear advantage of FWISTA compared to ISTA (Proposition 1) and SISTA (Proposition 2). Moreover, we note that FWISTA can be simply adapted in order to impose a monotonic decrease of the cost functional value, in the same fashion as MFISTA. The same convergence properties apply [19, Thm. 5.1].

E. Random Shifting

Wavelet bases perform well the compression of signals but can introduce artifacts that can be attributed to their relative lack of shift-invariance. In the case of regularization, this can be avoided by switching to a redundant dictionary. The downside, however, is a significant increase in computational cost. Alternatively, the practical technique referred to as random shifting (RS) [START_REF] Figueiredo | An EM algorithm for waveletbased image restoration[END_REF] can be used. Applying random shifting is much simpler and computationally more efficient than considering redundant transforms and leads to sensibly improved reconstruction.

Here, we propose a variational interpretation that motivates our implementation of FWISTA with RS (see Algorithm 4). We consider the DWT

[W 1 • • • W Ns ]
H , with W i = S i W , where S i represent the different shifting operations required to get a translation-invariant DWT. The desired reconstruction would be defined as the minimizer of

C(c) = m -Ec 2 2 + λ N s [W 1 • • • W Ns ] H c 1 . (22)
In 1D, this formulation includes TV regularization, in other words a single-level undecimated Haar WT without coarsescale thresholding. Rewriting [START_REF] Bayram | A subband adaptive iterative shrinkage/thresholding algorithm[END_REF] in terms of wavelet coefficients, we get

N s C(c) = i C i (W -1 i c), (23)
with

C i (w i ) = m -M S -1 i w i 2 2 + λ w i 1 . (24)
For a current estimate, we select a transform W i and perform a step in the minimization of the cost w.r.t. w i while keeping w j , for j = i fixed. A SISTA step is used, as the minimization subproblem (24) takes the form [START_REF] Guerquin-Kern | Wavelet-regularized reconstruction for rapid MRI[END_REF]. The minimizers of the functionals C i are expected to correspond to images close to each other and to the minimizer of [START_REF] Bayram | A subband adaptive iterative shrinkage/thresholding algorithm[END_REF]. In the first iterations of the algorithm, the minimization steps w.r.t. any w i are functionally equivalent (i.e., the modification is mostly explained by the gradient step). This motivates the use of FWISTA iterations at first and the switch to ISTA steps afterwards.

As the scheme is intrinsically greedy, we do not have a theoretical guarantee of convergence. Yet, in practice, we have observed that the SER stabilizes at a much higher value than it does when using ISTA schemes with no RS (cf. Figure 4).

Algorithm 4: FWISTA with RS

input : A = E H E, a = E H m, c 0 , Λ -1 = diag(τ ), K; Consider the sequence of DWT with RS: W n ; Initialization: n = 0, k = 0, γ = 1, t 0 = 1, C 0 = C(c 0 ), y 0 = c 0 ; repeat w n ← W n -1 y n ; c n+1 ← W n T λτ w n + Λ -1 W -1 n (a -Ay n ) ; C n+1 ← C(c n ) (cost-function evaluation); if C n+1 > C n then k ← k + 1; if k = K then γ ← 0; Λ -1 ← max{τ }I; τ ← max{τ }1; t n+1 ← 1 + 1 + 4t 2 n /2; y n+1 ← c n+1 + γ tn-1 tn+1 (c n+1 -c n ); n ← n + 1;
until stopping condition is met;

IV. EXPERIMENTS

A. Implementation Details

Our implementation uses Matlab 7.9 (Mathworks, Natick). The reconstructions run on a 64-bit 8-core computer, clock rate 2.8GHz, 8GB RAM (DDR2 at 800 MHz), Mac OS X 10.6.5. For all iterative algorithms, a key point is that matrices are not stored in memory. They only represent operations that are performed on vectors (images). In particular, a = M H m is computed once per dataset. Matrixto-vector multiplication with A = M H M , specifically, E H E, have an efficient implementation thanks to the convolution structure of the problem [START_REF] Wajer | Major speedup of reconstruction for sensitivity encoding with arbitrary trajectories[END_REF], [START_REF] Fessler | Toeplitzbased iterative image reconstruction for MRI with correction for magnetic field inhomogeneity[END_REF]. For these Fourier precomputations, we used the NUFFT algorithm [START_REF] Fessler | Nonuniform fast Fourier transforms using min-max interpolation[END_REF] that is made available at http://www.eecs.umich.edu/ fessler/code/. For wavelet transforms, we used the code provided at http://bigwww.epfl.ch/algorithms/mltldeconvolution [START_REF]A fast multilevel algorithm for wavelet-regularized image restoration[END_REF]. This Fourier-domain implementation proved to be faster than Matlab's when considering reconstructed images smaller than 256 × 256 and the Haar wavelet. It must also be noted that the 2D-DFT were performed using the FFTW library which efficiently parallelizes computations.

For Tikhonov regularizations, we implemented the classical conjugate gradient (CG) algorithm, with the identity as the regularization matrix. For TV regularizations, we considered the iteratively reweighted least-squares algorithm (IRLS), which corresponds to the additive form of half-quadratic minimization [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF], [START_REF] Champagnat | A connection between half-quadratic criteria and EM algorithms[END_REF]. We used 15 iterations of CG to solve the linear inner problems, always starting from the current estimate, which is crucial for efficiency. For the weights that permit the quadratic approximation of the TV term, we stabilized the inversion of very small values.

We implemented ISTA, SISTA, and FWISTA as described in Section III, with the additional possibility to use random shifting (see Section IV-C2). For the considered reconstructions using our method, described in Algorithm 4, K = 30 was a reasonable choice. The Haar wavelet transform was used, with 3 decomposition levels when no other values are mentioned. As is usual for wavelet-based reconstructions, the regularization was not applied to the coarse-level coefficients.

Reconstructions were limited to the pixels of the ROI for all algorithms. The regularization parameter λ was systematically adjusted such that the reconstruction mean-squared error (MSE) inside the ROI was minimal. For practical situations where the ground-truth reference is not available, it is possible to adjust λ by considering well-established techniques such as the discrepancy principle, generalized cross validation, or L-curve method [START_REF] Karl | Handbook of image and video processing[END_REF].

B. Spiral MRI Reconstruction

In this section, we focused on the problem of reconstructing images of objects weighted by the receiving channel sensitivity, given undersampled measurements. This problem, which involves single-channel data and hence differs from SENSE, is challenging for classical linear reconstructions as it generates artifacts and propagates noise. We considered spiral trajectories with 50 interleaves, with an interleave sampling density reduced by a factor R = 1.8 compared to Nyquist for the highest frequencies and an oversampling factor 3.5 along the trajectory. Spiral acquisition schemes are attractive because of their versatility and the fact that they can be implemented with smooth gradient switching [START_REF] Glover | Simple analytic spiral k-space algorithm[END_REF], [START_REF] Kim | Simple analytic variable density spiral design[END_REF].

We validate the results with the three sets of data that we present below.

1) MR Scanner Acquisitions: The data were collected on a 3T Achieva system (Philips Medical Systems, Best, The Netherlands). A field camera with 12 probes was used to monitor the actual k-space trajectory [START_REF] Barmet | A transmit/receive system for magnetic field monitoring of in vivo MRI[END_REF]. An array of 8 head coils provided the measurements. We acquired in vivo brain data from a healthy volunteer with parameters T R = 1000 ms and T E = 30 ms. The excitation slice thickness was 3 mm with a flip angle of 30 degrees. The trajectory was designed for a FOV of 25 cm with a pixel size 1.5 mm. It was composed of 100 spiral interleaves. The interleaf distance for the highest sampled frequencies defined a fraction of the Nyquist sampling density (R = 0.9).

The subset used for reconstruction corresponds to half of the 100 interleaves. The corresponding reduction factor, defined as the ratio of the distance of neighboring interleaves with the Nyquist distance, is R = 1.8.

2) Analytical Simulation: We used analytical simulations of the Shepp-Logan (SL) brain phantom with a similar coil sensitivity, following the method described in [START_REF] Guerquin-Kern | Analytical form of Shepp-Logan phantom for parallel MRI[END_REF]. The values of these simulated data were scaled to have the same mean spatial value (i.e., the same central k-space peak) as the brain reference image. A realization of complex Gaussian noise was added to this synthetic k-space data, with a variance corresponding to 40 dB SNR. The 176×176 rasterization of the analytical object provided a reliable reference for comparisons.

3) Simulation of a Textured Object: A second simulation was considered with an object that is more realistic than the SL phantom. We chose a 512 × 512 MR image of a wrist that showed little noise and interesting textures. We simulated acquisitions with the same coil profile and spiral trajectory (176 × 176 reconstruction matrix), in presence of a 40 dB SNR Gaussian complex noise. The height of the central peak was also adjusted to correspond to that of the brain data. The reference image was obtained by sinc-interpolation, by extracting the lowest frequencies in the DFT.

C. Results

In this section, we present the different experiments we conducted. The two main reconstruction performance measures that we considered are • Reconstruction duration, which excludes all aforementioned precomputations and the superfluous monitoring operations.

• Signal to error ratio w.r.t. a reference, defined as SER( ρ,

ρ ref ) = ρ ref 2 / ρ ref -ρ 2 and SERdB( ρ, ρ ref ) = 20 log 10 (SER( ρ, ρ ref )).
Practically, the references are either the ground-truth images or the minimizer of the cost functional. It is known that SER is not a foolproof measure of visual improvement but large SER values are encouraging and generally correlate with good image quality. 1) Convergence Performance of IST-Algorithms: In this first experiment, we compared the convergence properties of the different ISTA-type algorithms, as presented in Section III, with the Haar wavelet transform. The data we considered are those of the MR wrist image. The regularization parameter was adjusted to maximize the reconstruction SER w.r.t. to the ground-truth data. The actual minimizer of the cost functional, which is the common fixed-point of this family of algorithms, was estimated by iterating FWISTA 100 000 times.

The convergence results are shown in Figures 2 and3 for the simulation of the MR wrist image. Similar graphs are obtained using the other sets of data.

For a fixed number of iterations, FISTA schemes (FISTA and FWISTA) require roughly 10% additional time compared to ISTA and SISTA. In spite of this fact, their asymptotic superiority appears clearly in both figures. The slope of the decrease of the cost functional in the log-log plot of Figure 3 reflects the convergence properties in Propositions 1, 2, and 3. When considering the first iterations, which are of greatest practical interest, the algorithms with optimized parameters (SISTA and FWISTA) perform better than ISTA and FISTA (see Figure 3). The times required by each algorithm to reach a 30 dB SER (considered as a threshold value to perceived changes) are 415 s (ISTA), 53 s (SISTA), 12.7 s (FISTA), and 4.4 s (FWISTA). With respect to this criterion, SISTA presents an 8-fold speedup over ISTA, while FWISTA presents a 12-fold speedup over SISTA and nearly a 3-fold speedup over FISTA. It follows that FWISTA is practically close to two orders of magnitude faster than ISTA.

2) Choice of the Wavelet Transform and Use of Random Shifting: The algorithms presented in Section III apply for any orthogonal wavelet basis. For the considered application, we want to study the influence of the basis on performance. In this experiment we considered the Battle-Lemarié spline wavelets [START_REF] Mallat | A wavelet tour of signal processing[END_REF] with increasing degrees, taking into account the necessary postfilter mentioned in [START_REF] Gamper | Compressed sensing in dynamic MRI[END_REF].

We compared the best results for several bases. They were obtained with FWISTA after practical convergence and are reported in Table II. Figure 4 illustrates the time evolution of the SER using ISTA and FWISTA in the case of the SL reconstruction. Similar graphs are obtained with the other experiments.

It is known that the Haar wavelet basis efficiently approximates piecewise-constant objects like the SL phantom, which is consistent with our results. On the other hand, splines of higher degree, which have additional vanishing moments, perform better on the textured images (upper part of Table II).

We present in the lower part of Table II the performances observed when using ISTA with RS. We conclude that, in the case of realistic data, it is crucial to use RS as it improves results by at least 0.7 dB, whatever the wavelet basis is. The remarkable aspect there is that the Haar wavelet transform with RS consistently performs best. Two important things can be seen in Figure 4: FWISTA is particularly efficient during the very first iterations, while SISTA with RS yields the best asymptotic results in terms of SER and stability. Our method combines both advantages.

In Table III, we present the results obtained using different depths of the wavelet decompositions. Our reconstruction method is used together with the Haar wavelet transform and RS. The performances are similar but there seems to be an advantage in using several decomposition levels both in terms of SER and reconstruction speed. The FWISTA scheme seems to recoup the cost of the wavelet transform operations associated to an increase in the depth of decomposition.

3) Practical Performance: We report in Table IV the results obtained for different reconstruction experiments using state-ofthe-art linear reconstruction, TV regularization, and our method. The images obtained when running the different algorithms after approximately 5 s, and after practical convergence as well, are shown in Figure 5. We display in Figure 6 the time evolution of the SER for the different experiments. In each case, we emphasize the time required to reach -0.5 dB of the asymptotic value of SER. Finally, we present in Figure 7 the reconstruction and error maps of the different IST-algorithms at different moments of reconstruction. This was done with the wrist simulated experiment using the Haar wavelet basis and RS.

Firstly, we observe that TV and our method achieve similar SER (Table IV) and image quality (Figure 5). They both clearly outperform linear reconstruction, with a SER improvement from 1.5 to 3 dB, depending on the degree of texture in the original data. Moreover, the pointwise maximal reconstruction error appears to always be smaller with nonlinear reconstructions. Due to the challenging reconstruction task, which significantly undersamples of the k-space, residual artifacts remain in the linear reconstructions and at early stages of the nonlinear ones. Although the k-space trajectory is exactly the same in the three cases, artifacts are less perceived in the in vivo reconstructions, while they stand out for the synthetic experiments. Secondly, it clearly appears that the linear reconstruction, implemented with CG, leads to the fastest convergence, unfortunately with suboptimal quality. For a reconstruction time one order of magnitude longer, our accelerated method provides better reconstructions. This is illustrated in Figure 5 for reconstruction times of the order of five seconds (columns 1, 2, and 4).

Finally, we observe in Figure 7 the superiority of the proposed FWISTA over the other types of algorithms. For the given reconstruction times, it consistently exhibits better image quality as can be seen in both reconstructions and error maps.

D. SENSE MRI Reconstruction

Our reconstruction method is applicable to linear MR imaging modalities. In this section, we report results obtained on a SENSE reconstruction problem. The data were acquired with the same scanner setup as in Section IV-B1. This time, the data from the 8 receiving channels were used for reconstruction, as well as an estimation of the sensitivity maps. An in vivo gradient echo EPI sequence of a brain was performed with T2* contrast. The data was acquired with the following parameters: excitation slice thickness of 4mm, T E = 35ms, T R = 900ms, flip angle of 80 degrees, and trajectory composed of 13 interleaves, supporting a 200 × 200 reconstruction matrix with pixel resolution 1.18mm × 1.18mm. The oversampling ratio along the readout direction was 1.62.

The reference image was obtained using the complete set of data and performing an unregularized CG-SENSE reconstruction. The reconstruction involved 3 of the 13 interleaves, representing a significant undersampling ratio R = 4.33.

The images obtained using regularized linear reconstruction (CG), TV (IRLS), and our method are presented in Figure 8. In Figure 9, the SER evolution with respect to time is shown for the three methods. The times to reach -0.5 dB of the asymptotic SER value are 5.9s (CG), 49.1s (IRLS), and 22.8s (our method).

With this high undersampling, the errors maps show that reconstructions suffer from noise propagation mostly in the center of the image. It appears that TV and our method improve qualitatively and quantitatively image quality over linear reconstruction (cf. Figure 8). As it was observed in Section IV-B with the spiral MRI reconstructions, this in vivo SENSE experiment confirms that our method is competitive with TV. In terms of reconstruction duration, our method proves to converge in a time that is of the same order of magnitude as CG (cf. Figure 9).

V. CONCLUSIONS

We proposed an accelerated algorithm for nonlinear waveletregularized reconstruction that is based on two complementary acceleration strategies: use of adaptive subband thresholds plus multistep update rule. We provided theoretical evidence that this algorithm leads to faster convergence than when using the accelerating techniques independently. In the context of MRI, the proposed strategy can accelerate the reference algorithm up to two orders of magnitude. Moreover, we demonstrated that, by using the Haar wavelet transform with random shifting, we are able to boost the performance of wavelet methods to make them competitive with TV regularization. Using different simulations and in vivo data, we compared the practical performance of our reconstruction method with other linear and nonlinear ones.

The proposed method is proved to be competitive with TV regularization in terms of image quality. It typically converges within five seconds for the single channel problems considered. This brings nonlinear reconstruction forward to an order of magnitude of the time required by the state-of-the-art linear reconstructions, while providing much better quality.

2 2

 2 and leads to the closedform linear solution w

Fig. 1 .

 1 Fig. 1. Reference images from left to right: in vivo brain, SL reference, and wrist.

Fig. 3 .

 3 Fig. 3. Time evolution of the difference in cost function value w.r.t. the minimizer for several ISTA algorithms.

Fig. 4 .

 4 Fig. 4. Time evolution of the SER for several algorithms for the SL simulation using Haar wavelets.

Fig. 6 .

 6 Fig. 6. Evolution of the performance of the algorithms. From top to bottom: SL simulation, wrist simulation, and brain data. Times required to reach -0.5 dB of the asymptotic value are indicated.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Reconstructions (left column) and error maps (right column) for the SENSE EPI experiment using CG (first row), IRLS-TV (second row), and our method (third row). For each reconstruction, the performance in SER w.r.t. the reference (top-left corner), the reconstruction time (top-right corner), the number of iterations (bottom-right corner) and a magnification of the central part (bottom-left) are shown.

TABLE II VALUES

 II OF THE OPTIMAL SER AND CORRESPONDING REGULARIZATION PARAMETERS ARE SHOWN FOR THE DIFFERENT WAVELET BASES.

	Experiment		SL simulation			Wrist simulation
	Wavelet basis	Haar	Spline 2 Spline 4 Spline 6	Haar	Spline 2 Spline 4 Spline 6
	Without RS SER (dB) 12.65	12.16	10.75	9.70	15.93	17.33	17.32	17.07
		λ	1 870	2 510	830	1 460	1 600	946	1 070	1 350
	With RS	SER (dB) 13.38	12.53	11.58	10.38	18.70	18.24	18.05	17.87
		λ	5 650	3 900	7 770	1 370	1 490	850	1 190	1 260

TABLE III RESULTS

 III OF THE PROPOSED WAVELET METHOD FOR DIFFERENT WAVELET DECOMPOSITION DEPTHS. VALUES OF THE REGULARIZATION PARAMETER, THE FINAL SER, THE RELATIVE MAXIMAL SPATIAL DOMAIN ERROR, AND THE TIME TO REACH -0.5dB OF THE FINAL SER.

	Experiment		SL simulation			Wrist simulation			Brain data
	Decomposition depth 1	2	3	4	1	2	3	4	1	2	3	4
	λ opt.	5 330	5 680	5 110	5 200	1 500	1 520	1 570	1 700	9 650	10 800 11 400 11 400
	SER (dB) opt.	13.25 13.34	13.39	13.35	18.54 18.66	18.71	18.70	18.78	19.00	18.99	18.96
	∞ error (%)	50	51	54	53	15.5	16.0	16.0	16.5	17.3	17.5	17.49	17.0
	t -0.5dB (s)	6.56	6.19	6.93	7.59	5.88	5.23	5.33	4.25	8.99	7.03	6.33	7.97

TABLE IV RESULTS

 IV OF THE ALGORITHMS CG (LINEAR), IRLS (TV), AND OUR METHOD (WAVELETS). VALUES OF THE REGULARIZATION PARAMETER, THE FINAL SER, THE RELATIVE MAXIMAL SPATIAL DOMAIN ERROR, AND THE TIME TO REACH -0.5dB OF THE FINAL SER.

				Experiment		SL simulation	Wrist simulation	Brain data
					Method		linear	TV	wavelets linear TV	wavelets	linear TV	wavelets
					λ opt.		0.0247 4 090	6 380	0.436	760	1 620	0.471	6 050	16 800
				SER (dB) opt.	8.46	13.82	13.17	16.14	18.41	18.64	15.81	18.88	18.93
				∞ error (%)	48		49	51	21	16	16	29	12	11
				t -0.5dB (s)	0.286	18.1	5.40	0.209	10.5	4.64	0.205	15.2	6.13
		80
		70								ISTA
	SER (dB)	40 50 60								SISTA FISTA FWISTA
		30
		20
		0	50	100	150	200 time (s)	250	300	350	400

Fig. 2. Time evolution of the SER w.r.t. the minimizer for several ISTA algorithms. Times to reach 35dB are delineated.

ACKNOWLEDGMENTS

We thank J. Fessler and C. Vonesch to make their implementations of NUFFT and multidimensional wavelet transform available. We acknowledge C. Vonesch, İ. Bayram, and D. Van De Ville for fruitful discussions. We also thank P. Thévenaz for useful comments and careful proofreading.


Marcin Jakub Kamiński

Jean-Florent Raymond
email: jean-florent.raymond@mimuw.edu.pl

Théophile Trunck
email: theophile.trunck@ens-lyon.org.

Well-quasi-ordering H-contraction-free graphs *

Introduction

A well-quasi-order is a quasi-order where every decreasing sequence and every collection of incomparable elements (called an antichain) are finite. Well-quasiorders enjoy nice combinatorial properties that can be used in several contexts, from algebra to computational complexity and algorithms. Since its introduction more than sixty years ago, the theory of well-quasi-orders led to major results in Graph Theory and Combinatorics. In particular, Kruskal showed in [START_REF] Kruskal | Well-quasi-ordering, the tree theorem, and Vazsonyi's conjecture[END_REF] that trees are well-quasi-ordered by homeomorphic embedding, and Robertson and Seymour proved that both the minor relation and the immersion relation are well-quasi-orders on the class of finite graphs [START_REF] Robertson | Graph Minors XX. Wagner's conjecture[END_REF][START_REF] Robertson | Graph Minors XXIII. Nash-Williams' immersion conjecture[END_REF]. Most of the usual quasiorders on graphs are not well-quasi-orders in general, though. Given one of these quasi-orders, a natural line of research is to identify the subclasses that are well-quasi-ordered. Our work is motivated by the following results.

Theorem 1 ( [START_REF] Damaschke | Induced subgraphs and well-quasi-ordering[END_REF]). The class of H-induced subgraph-free graphs is well-quasiordered by induced subgraphs iff H is an induced subgraph of P 4 . 1Theorem 2 ( [START_REF] Ding | Subgraphs and well-quasi-ordering[END_REF]). The class of H-subgraph-free graphs is well-quasi-ordered by subgraphs iff H is a subgraph of P n , for some n ∈ N. 1 Theorem 3 ( [START_REF] Chun | Graph Structures and Well-Quasi-Ordering[END_REF]). The class of H-topological minor-free multigraphs is wellquasi-ordered by topological minors iff H is a topological minor of R n , for some n ∈ N. 2Theorem 4 ( [START_REF] Law B Lasiok | Induced minors and well-quasi-ordering[END_REF]). The class of H-induced minor-free graphs is well-quasi-ordered by induced minors iff H is an induced minor of the gem or K4 . 3These results characterize the closed classes defined by one forbidden substructure that are well-quasi-orders. Like the four containment relations on graphs mentioned in the above results, the contraction relation is not a wellquasi-order in general. Let the diamond be the graph obtained from K 4 by deleting an edge. Our main contribution in this direction is the following result.

Theorem 5. The class of connected H-contraction-free graphs is well-quasiordered by contractions iff H is a contraction of the diamond.

The requirement of connectivity in Theorem 5 is necessary in the sense that for every graph H, the class of (not necessarily connected) H-contraction-free graphs contains the infinite antichain {K i , i ∈ N ≥h+1 } (where h = |V (H)|) and therefore is not a well-quasi-order. Theorem 5 can be seen as contraction counterpart of the results mentioned above.

Another line of research when dealing with quasi-orders that are not wellquasi-orders in general is to look at canonical antichains. An antichain is canonical if for every closed subset F of the quasi-order, F is a wqo iff F has a finite intersection with this antichain. Intuitively, a canonical antichain represents all infinite antichains of a quasi-order. As shown by the results below, the question of the presence or absence of a canonical antichain has been studied for several containment relations and graph classes. We give an answer to this question for the containment relation with the following result.

Theorem 10. Under the contraction relation, the class of finite graphs does not have a canonical antichain.

The proof of Theorem 10 relies on the tools introduced in [START_REF] Ding | On canonical antichains[END_REF] that can be used to prove that a quasi-order does not have a canonical antichain.

Organization of the paper. The proof of Theorem 5 contains three parts.

The first one, given in Section 3, is a study of infinite antichains of the contraction relation from which we can deduce that if the class of H-contraction-free graphs is well-quasi-ordered by contractions, then H is a contraction of the diamond. Section 4 contains the second part which is a decomposition theorem for diamond-contraction-free graphs. The last part uses this decomposition to show the well-quasi-ordering result and is presented in Section 5. The proof of Theorem 10 is given in Section 6. Definitions of the terms and notations used are introduced in Section 2.

Preliminaries

We use the notation N ≥k for the set {i ∈ N, i ≥ k}, for every k ∈ N. For every set S, we denote by P(S) the collection of subsets of S.

Graphs. All graphs in this paper are finite, simple, and undirected. We denote by V (G) the vertex set of a graph G and by E(G) its edge set. If X ⊆ V (G), the subgraph of G induced by X, which we write G[X], is the graph with vertex set X and edge set E

(G) ∩ X 2 . Let C be a (not necessarily induced) cycle in a graph G. A pair of vertices {u, v} ⊆ V (C) that are not adjacent in C is a chord of C in G if {u, v} ∈ E(G). Otherwise {u, v} is a non-chord of C in G.
A vertex v of a graph G is a cutvertex if G \ {v} has more connected components than G. A block is a maximal subgraph that has no cutvertex. A clique-cactus graph is a graph whose blocks are cycles and cliques (cf. Figure 1 for an example).

If G is a graph, then G is the graph obtained by replacing all non-edges by edges and vice versa. For every positive integer r we denote by D r the graph

• K 1 ∪ •K r . In particular, D 2 is the diamond. We set D = {D r , r ∈ N} (cf.

Subsets of vertices.

If G is a graph, the degree of a subset S ⊆ V (G) is the number of vertices of V (G) \ S that have a neighbor in S. The subset S is said to be connected if G[S] is connected. We say that S is adjacent to some vertex v (respectively some subset S ⊆ V (G)) if there is an edge from v to a vertex of S (respectively from a vertex of S to a vertex of S ).

Contractions.

: V (H) → P(V (G)) such that: (i) for every v ∈ V (H), ϕ(v) is connected; (ii) {ϕ(v), v ∈ V (H)} is a partition of V (G);
(iii) for every u, v ∈ V (H), the vertices u and v are adjacent in H iff the subsets ϕ(u) and ϕ(v) are adjacent in G.

This definition has several consequences. In particular, the degree of v ∈ V (H) in H is at most the degree of ϕ(v) in G. Also, there is no model of a graph with no dominating vertex in a graph with a dominating vertex.

It is easy to check that H is a contraction of G iff there is a contraction model of H in G. A graph G is said to exclude a graph H as contraction, or to be H-contraction-free, if H is not a contraction of G. We denote the class of connected H-contraction-free graphs by Excl(H).

We say that a graph H is an induced minor of a graph G if it can be obtained from G by deleting vertices and contracting edges.

Sequences and orders. We write s 1 , . . . , s n the sequence containing the elements s 1 , . . . , s n in this order. For every set S, we denote by S the set of all finite sequences over S, including the empty sequence. For any partial order (A, ), we define the relation on A as follows: for every r = r 1 , . . . , r p and s = s 1 , . . . , s q of A , we have r s if there is an increasing function ϕ : {1, . . . , p} → {1, . . . , q} such that for every i ∈ {1, . . . , p} we have r i s ϕ(i) . This generalizes the subsequence relation.

Well-quasi-orders and antichains. Given an order over S, a sequence over S is said to be an antichain of (S, ) if its elements are pairwise incomparable with respect to . A well-quasi-order (wqo for short) is a quasi-order where every decreasing sequence and every antichain is finite.

We will use the two classical results stated below.

Proposition 1 (Folklore). Let (A, ) be a quasi-order and B, C ⊆ A. If both (B, ) and (C, ) are wqo, then so is (B ∪ C, ).

Proposition 2 (Higman's Lemma [START_REF] Higman | Ordering by divisibility in abstract algebras[END_REF]). (A, A ) is a wqo, then so is (A , A ).

If (S, ) is a quasi-order that is not a wqo, a minimal antichain [START_REF] St | On well-quasi-ordering finite trees[END_REF] of (S, ) is an antichain a i i∈N where for every i ∈ N, a i is a minimal element (with respect to ) such that there is an infinite antichain of (S, ) starting with a j j∈{0,...,i} . Observe that every quasi-order that is not a wqo and that has no infinite decreasing sequence has a minimal antichain. For every subset A ⊆ S, we define: Incl(A) = {x ∈ S, ∃y ∈ A, x y and x = y}.

An antichain A is fundamental if (Incl(A), ) is a wqo. A set F ⊆ S is said to be -closed if it satisfies the following property: ∀x ∈ F, ∀y ∈ S, y x ⇒ y ∈ F . An antichain A of a quasi-order (S, ) is canonical if for every -closed subset F ⊆ S, we have

F ∩ A is finite ⇐⇒ (F, ) is a wqo.
Let us end this section by a simple observation.

Observation 1. Every sequence of graphs that is decreasing with respect to contraction is finite. In fact, in such a sequence the number of edges is also decreasing, as every edge contraction decreases the number of edges of a graph by at least one.

A consequence of the above observation is that infinite antichains are the only obstacles for a class of graphs to be well-quasi-ordered by the contraction relation. We deal with them in the next section.

Infinite antichains

A simple but crucial observation in the study of the well-quasi-orderability of classes that are defined by forbidden structures (of any kind) is the following. If none of the graphs of an infinite antichain A contains some graph H, then excluding H does not give a well-quasi-order. Indeed, the class obtained still contains the infinite antichain A. Let us restate this observation in terms of contractions.

Observation 2. Let A be an infinite antichain of the contraction relation. If (Excl(H), ≤ ctr ) is a wqo, then all but finitely many graphs of A contain H as contraction.

For this reason, we deal here with infinite antichains of the contraction relation. The simplest one is certainly the class of complete bipartite graphs with one part of size two:

A K = {K 2,r , r ∈ N ≥2 }.
Lemma 1. For every p, q, p , q ∈ N ≥2 such that p ≤ p and q < q , there is no model of K p,q in K p ,q .

Proof. Let us assume for a contradiction that there is a model ϕ of K p,q in K p ,q . As K p ,q has more vertices than K p,q , there is a vertex v of K p,q such that |ϕ(v)| ≥ 2. Observe that every subset of at least two vertices of K p,q that induced a connected subgraph is dominating. Indeed, such a subset must contain at least a vertex from each part of the bipartition. It follows from the definition of a model that K p ,q has a dominating vertex, a contradiction. Therefore, there is no model of K p,q in K p ,q .

Corollary 1. {K 2,p , p ∈ N ≥2 } is an antichain of ≤ ctr .
Recall that S is the class of stars and that D = {D r , r ∈ N}, where D r is the graph that can be obtained by contracting an edge of K 2,r+1 , for every r ∈ N (cf. Section 2). The following will be useful later.

Observation 3. For every p ∈ N ≥1 , contracting one edge in D p gives either D p-1 , or K 1,p , depending on which edge is contracted.

As we want to identify graphs H such that (Excl(H), ≤ ctr ) is a wqo, we must consider every graph H such that A ∩ Excl(H) is finite, for every antichain A. A first step towards this goal is the following observation.

Lemma 2. Let p ∈ N ≥2 . If H ≤ ctr K 2,p and H = K 2,p , then H ∈ D ∪ S.
Proof. Given that H ≤ ctr K 2,p , there is a sequence of contractions transforming K 2,p into H. If this sequence contains only one contraction, then it is straightforward that H = D p-1 . Therefore in the other cases H is a contraction of D p-1 . We get the result from Observation 3 and the observation that every contraction of a graph of S (i.e. the class of stars) belongs to S. Observation 4. For every positive integers p, q such that p < q, we have

D p ≤ ctr K 2,q .
Indeed, if F is a collection of q -p edges of K 2,q that all incident with the same vertex of degree p, then it is easy to check that contracting F in K 2,q yields D p . An immediate consequence of Observation 4 is that A K ∩ Excl(D p ) is finite for every positive integer p.

From the fact that every graph of D ∪ S is a contraction of D p for some positive integer p, Observation 4 gives. Observation 5. If (Excl(H), ≤ ctr ) is a wqo, then H ≤ ctr D p for some p ∈ N ≥1 However, we will need another antichain in order to find more properties that H must satisfy. Let us consider the antichain of antiholes, which already appeared in [START_REF] Law B Lasiok | Induced minors and well-quasi-ordering[END_REF] in the context of induced minors: 3). This connection with the induced minor relation (where edge contractions and vertex deletions are allowed) is not surprising: as every contraction is an induced minor, every antichain of the induced minor relation is also an antichain of the contraction relation. For completeness, we include the following proof.

A C = {C i , i ∈ N ≥6 } (cf. Fig- ure
Lemma 3 (See also [START_REF] Law B Lasiok | Induced minors and well-quasi-ordering[END_REF]Lemma 1]). A C is an antichain of the contraction relation.

Proof. Towards a contradiction, let us assume that there is a contraction model ϕ of C p in C q for some integers p, q ∈ N ≥3 such that p < q. Recall that the image of ϕ is a partition of

C q . As |C p | < |C q |, there is a vertex v of C p such that |ϕ(v)| ≥ 2.
Observe that for every choice of two vertices of C p there is at most one vertex which is not adjacent to one of them. Therefore, there is at most one set in {ϕ(u), u ∈ V (C p ) \ {v}} that is not adjacent to ϕ(v). This contradicts the fact that ϕ is a model of C p in C q as every vertex of C q is adjacent to all but two vertices. Consequently there is no contraction model of C p in C q , for every integers p, q ∈ N ≥3 , p < q.

Again, we look at graphs H such that Excl(H) ∩ A C is finite. As a wqo must contain none of A K and A C , it is enough to consider graphs such that Excl(H) ∩ A C is finite among those for which Excl(H) ∩ A K is finite.

Lemma 4. If p ≥ 3 then Excl(D p ) ∩ A C is infinite.
Proof. For every p ≥ 3, then graph D p has independence number at least 3. Let q > p. As contracting edges can only decrease the independence number, there is no sequence of contractions transforming C q (which has independent number 2) to D p , for every integer q > p. Therefore C q ∈ Excl(D p ), for every integer q > p.

Corollary 2. If (Excl(H), ≤ ctr ) is a wqo, then H ≤ ctr D 2 .
The next sections are devoted to graphs not containing D 2 as contraction. We will first prove a decomposition theorem for the graphs in this class, which we will use to show that (Excl(D 2 ), ≤ ctr ) is a wqo.

On graphs with no diamond

In this section we show that graphs in Excl(D 2 ) have a simple structure. More precisely, we prove the following lemma. Recall that clique-cactus graphs are the graphs whose blocks are cycles or cliques. Lemma 5. Graphs of Excl(D 2 ) are exactly the connected clique-cactus graphs.

The proof of Lemma 5 will be given after a few lemmas. If C is a cycle of a graph G and {u, v}, {u , v } ⊆ V (C), we say that {u, v} and {u , v } are crossing in C if u, v, u , v are distinct and are appearing in this order on the cycle. Lemma 6. Let G be a graph and let C be a cycle in G. If C has at least one chord and one non-chord in G, then it has one chord and one non-chord that are crossing in C.

Proof. Let {x, x } be a non-chord of C in G and let P and Q be the connected components of C \ {x, x } which obviously are paths. Let us assume that every chord of C in G has both endpoints either in P or in Q (otherwise we are done) and let {y, y } be a chord of C in G, the endpoints of which belong, say, to P . Let z be a vertex of the subpath of P delimited by y and y such that z ∈ {y, y }, and let z be a vertex of Q. If {z, z } is a chord of C in G, then {x, x } and {z, z } are satisfying the required property. Otherwise, {z, z } is a non-chord and now {y, y } and {z, z } are crossing. Proof. Let G be a graph of Excl(D 2 ) and let C be a cycle of G. Towards a contradiction, let us assume that C has at least one chord {u, u } and one nonchord {v, v }. According to Lemma 6 we can assume without loss of generality that they are crossing in C. Let P and Q be the two connected components of C \ {v, v }. Contracting P to a single vertex x and Q to y yields a graph G such that:

• v, x, v , y lie on the cycle in this order;

• {v, v } ∈ E(G ); and
• {x, y} ∈ E(G ) (as {u, u } connects the subgraphs that are respectively contracted to x and y). As G is 2-connected, h i has at least two neighbors on C: let x and y be two of them.

Notice that G [{v, x, v , y}] is isomorphic to D 2 ,
Claim 1. x and y are not adjacent.

Proof. Let us assume that {x, y} ∈ E(G ). Let C be the cycle obtained from C by replacing the edge {x, y} by the path xh i y. This cycle is not induced as x, y are not adjacent in C whereas {x, y} ∈ E(G). It does not induce a clique either since x is not adjacent with the other neighbor of y on C (which is not x as we assume that C has at least 4 vertices). This contradicts Lemma 7 and therefore proves that {x, y} ∈ E(G).

Every pair of distinct vertices of the cycle C defines two subpaths of C meeting only at these vertices. Let u and v be two vertices of C such that h i has at least one neighbor in the interior of each of the subpaths of C defined by u and v, that we will respectively call P and Q. Such vertices exist, as a consequence of Claim 1.

Let us consider the contraction H of G obtained by contracting the interior path of P (respectively Q) to a single vertex w P (respectively w Q ) and then by contracting the edge connecting h i to w P . This edge exists by definition of u and v. Then uw P vw Q is a cycle of H where {w P , w Q } is a chord (because we contracted to w P the vertex h 1 which was adjacent to both w P and w Q ) and {u, v} is a non-chord (as they were non-adjacent vertices of the induced cycle C and that nothing has been contracted to them). According to Lemma 7, the graph H contains D 2 as contraction. As H is a contraction of G, then D 2 ≤ ctr G, a contradiction.

In both cases we reached a contradiction, therefore V (G) \ V (C) is empty: G is a clique or an induced cycle.

We are now ready to prove Lemma 5.

Proof of Lemma 5. The fact that a graph of Excl(D 2 ) is clique-cactus is a straightforward corollary of Lemma 8. It is easy to see that a clique-cactus graph does not contain D 2 as contraction by noticing that D 2 is a contraction of a graph if and only if it is a contraction of one of its 2-connected components. As D 2 is neither a contraction of a cycle, nor of a clique, we get the desired result.

Well-quasi-ordering clique-cactus graphs

We proved in the previous section that graphs of Excl(D 2 ) are exactly the connected clique-cactus graphs. This section contains the last part of the proof of Theorem 5, which is the following lemma. We conclude this section with the proof of Theorem 5.

Lemma 9. Connected clique-cactus graphs are well-quasi-ordered by ≤ ctr .

In this section, we deal with rooted graphs. A rooted graph is a graph which has a distinguished vertex, called root. The contraction relation is extended to the setting of rooted graphs by requiring that a model of a rooted graph H in a rooted graph G maps the root of H to a connected subgraph of G containing the root of G.

Let us denote by C the class of rooted connected clique-cactus graphs. In this class, two isomorphic graphs with a different root are seen as different. It is clear that proving that (C, ≤ ctr ) is a wqo implies Lemma 9. This is what we will do. • stick(G 0 , . . . , G p-1 ) is the graph obtained from U by identifying the vertices root(G 0 ), . . . , root(G p-1 );

Building blocks.

• cycle(G 0 , . . . , G p-1 ) is the graph obtained from U by adding the edges {root(G i ), root(G (i+1) mod p )} for every i ∈ {0, . . . , p -1}; and

• clique(G 0 , . . . , G p-1 ) is the graph obtained from U by adding the edges {root(G i ), root(G j )} for every distinct i, j ∈ {0, . . . , p -1}.

The root of stick(G 0 , . . . , G p-1 ), cycle(G 0 , . . . , G p-1 ) and clique(G 0 , . . . , G p-1 ) is the vertex that is the root of G 0 . These constructors will allow us to encode graphs of C into sequences.

We will now decompose graphs of C along blocks.

For every block B of a graph G, let dec B (G) denote the collection of all the graphs H that can be constructed from some connected component C of G\V (B) by adding a new vertex v adjacent to the vertices of C that are adjacent to a vertex of B in G, and setting root(H) = v.

Observe that as soon as root(G) ∈ V (B), every graph of dec B (G) is a proper contraction of G. Let dec(G) denote the union of the sets dec B (G) for every block B of G containing the root of G. The following observation is a consequence of Lemma 5. Observation 6. For every graph G ∈ C there is a (not necessarily unique) sequence G 0 , . . . , G p-1 ∈ dec(G) (for some p ∈ N) such that G = cycle(stick(G 0 ), . . . , stick(G p-1 )) or G = clique(stick(G 0 ), . . . , stick(G p-1 )). and one graph of {B i } i∈N ≥k . Let us assume that for some i ∈ {0, . . . , π(k) -1} and j ∈ N ≥k we have A i ≤ B j . Then A i ≤ B j ≤ A π(i) , a contradiction with the fact that {A i } i∈N is an antichain. The case B j ≤ A i is not possible by the choice of k. This proves that (B, ≤ ctr ) is a wqo. According to Proposition 2, (B , ≤ ctr ) is also a wqo. Let B = {stick(H), H ∈ B }. Item (iii) of Lemma 10 implies that any antichain in (B , ≤ ctr ) can be translated into an antichain of the same length in (B , ≤ ctr ), hence (B , ≤ ctr ) is a wqo. By the same argument (now using items (i) and (ii) of Lemma 10), we deduce that the quasi-orders ({cycle(H), H ∈ B }, ≤ ctr ) and ({clique(H), H ∈ B }, ≤ ctr ) are well-quasi-orders. Therefore U = {cycle(H), H ∈ B }∪{clique(H), H ∈ B } is well-quasi-ordered by ≤ ctr , as a consequence of Proposition 1. According to Observation 6, we have {A i } i∈N ⊆ U. This contradicts the fact that {A i } i∈N is an infinite antichain. Therefore (C, ≤ ctr ) is a wqo and we are done.

We would like to point out that with a proof similar to Lemma 9, it is in fact possible prove that if a class G of 2-connected graphs is wqo by ≤ ctr , then so is the class of graphs whose blocks belong to G. The interested reader may have a look at [8, Lemma 5] for a result of this flavour, see also [START_REF] Michael | Well quasi orders in subclasses of bounded treewidth graphs and their algorithmic applications[END_REF]Theorem 5] and [START_REF] Law B Lasiok | Induced minors and well-quasi-ordering[END_REF].

The proof of Theorem 5 is now immediate.

Proof of Theorem 5. Let H be a graph such that Excl(H) is a wqo. Then H ≤ ctr D 2 , by Corollary 2. On the other hand, if H ≤ ctr D 2 then Excl(H) ⊆ Excl(D 2 ). Observe that every antichain (respectively decreasing sequence) of (Excl(H), ≤ ctr ) is an antichain (respectively a decreasing sequence) of (Excl(D 2 ), ≤ ctr ). As a consequence of Lemma 9 we get that (Excl(H), ≤ ctr ) is a wqo and we are done.

On canonical antichains

In this section, we will use the following result of Ding in order to prove Theorem 10. Figure 4 illustrates the requirements of the lemma.

Lemma 11 ( [4, Theorem 1.1]). Let (S, ) be a quasi-order, let a i i∈N be a sequence of elements of S and let {W i } i∈N be a sequence of sequences of elements of S. If we have (i) a i i∈N is a fundamental infinite antichain; and

(ii) for every i ∈ N, W i is a fundamental infinite antichain; and

(iii) for every i ∈ N and every H ∈ W i , a i H and there are no other comparable pairs of elements in i∈N W i ∪ {a i } i∈N then (S, ) does not have a canonical antichain. Note that [4, Theorem 1.1] mentions other obstructions to the existence of a canonical antichain, however we will only use that described in Lemma 11. We will now define some sequences of graphs and show that they satisfy the properties of Lemma 11.

For every p, q ∈ N, let W p,q be the graph obtained by adding two non-adjacent dominating vertices to the disjoint union of K p and K 2,q (see Figure 5). These two vertices are called poles, and the two vertices corresponding to the part of K 2,q of size 2 are called semipoles. Observe that the other vertices either have degree two (in which case they are adjacent to the two poles, only), or have degree four (and they are adjacent to both poles and both semipoles). Lemma 12. For every p, p , q, q ∈ N ≥3 , there is a model of W p,q in W p ,q iff (p, q) = (p , q ). Proof. Let us assume that there is a model ϕ of W p,q in W p ,q . Let v be a vertex of W p,q of degree two. By definition of a model, ϕ(v) must be a connected subset of degree 2 of V (W p ,q ). Let us consider the possible choices for this subset.

First case: ϕ(v) is of the form V (W p ,q ) \ {x, y}, for some u, v ∈ V (W p ,q ). Therefore, V (W p ,q ) \ ϕ(v) has two vertices. As W p,q has more than 3 vertices (recall that p, q ≥ 3) which are mapped by ϕ to disjoint subsets of V (W p ,q ), this case is not possible. Second case: ϕ(v) is the subset of vertices inducing the subgraph K 2,q used in the construction of W p ,q . Observe that the poles and semipoles of W p,q (4 vertices in total, as they are distinct from v) all have degree at least 3. Moreover, W p ,q \ ϕ(v) is isomorphic to K 2,p . As every connected subset of degree at least 3 of K 2,p must contain a pole of this graph, there are at most two such subsets that are disjoint. This contradicts the fact that the images by ϕ of the four poles and semipoles of W p,q are disjoint connected subsets of degree at least 3 of W p ,q \ ϕ(v). Therefore, this case is not possible neither. Third case: ϕ(v) = {x} for some vertex x ∈ V (W p ,q ) of degree 2. As we reach this case for every choice of a vertex of degree 2 of W p,q , we deduce p ≤ p . The same argument applied to vertices of degree 4 yields q ≤ q . Let us now consider poles and semipoles.

Let u be a pole. Observe that according to the above remarks, ϕ(u) must be adjacent to vertices of degree two, so it should contain a pole of W p,q . If ϕ(u) contains in addition a vertex of degree 2 or 4 of W p,q , then ϕ(u) is dominating. This is not possible since u is not dominating, therefore ϕ(u) = {v} for some pole v of W p ,q . Let us now assume that u is a semipole of W p,q . As previously, the above remarks imply that ϕ(u) is adjacent to vertices of degree 4 of W p ,q . Hence ϕ(u) contains a semipole of W p ,q (it cannot contain a pole as both belong to the image of poles of W p,q ). Therefore each semipole of W p,q is sent to a subset of V (W p ,q ) containing a semipole. Observe that ϕ(u) cannot contain a vertex of degree two otherwise it would not be connected. Besides, it cannot contain a vertex of degree 4 otherwise it would be adjacent to the image by ϕ of the other semipole of W p,q . Consequently ϕ(u) contains a semipole of W p ,q and no other vertex. We proved that for every u ∈ V (W p,q ), the set ϕ(u) is a singleton. Therefore |V (W p,q )| = |V (W p ,q )|. Given that p ≤ p and q ≤ q (as proved above), this is possible only if p = p and q = q . This concludes the proof. Corollary 3. {W p,q } p,q≥3 is an antichain for ≤ ctr .

For every

i ∈ N ≥3 , let W i = {W i,q } q∈N ≥3 .
Lemma 13. For every p, q, r ∈ N ≥3 , we have K 2,r ≤ ctr W p,q iff r = p + 1.

Proof. Let us consider a contraction model ϕ of K 2,r in W p,q . We call X the vertices of W p,q inducing the K 2,q used in the construction of this graph. Let us consider a vertex u of degree 2 of K 2,r . Exactly as in the proof of Lemma 12, there are three possible choices for ϕ(u). For the same reason as in this proof, the case where ϕ(u) = V (W p,q ) \ {x, y} (for some x, y ∈ V (W p,q )) is not possible. Therefore, either ϕ(u) = X, or ϕ(u) = {x} for some x ∈ V (W p,q ) of degree 2. Since this holds for every vertex of degree 2 of K 2,r , and as W p,q has exactly p vertices of degree 2, we deduce that r ≤ p + 1.

Let v, w denote the poles of K 2,r . Because of the observation above and of the definition of a contraction model, each of ϕ(v) and ϕ(w) must be adjacent to some vertex of degree 2 of W p,q and these sets should not be adjacent. The only possible choice for them is to let ϕ(v) be the singleton containing one pole of W p,q and ϕ(w) be the singleton containing the other pole. Observe that, in the case where p + 1 > r, either one vertex of degree 2 of W p,q or a vertex of X does not belong to the image of ϕ. This contradicts the definition of a model, hence this case is not possible.

The only remaining case is thus r = p + 1. Observe that X induces a connected subgraph. It is not hard to see that contracting X to a single vertex yields K 2,p+1 .

Observation 7. Let p, q ∈ N ≥3 . There is no induced path on four vertices in W p,q , neither in K 2,p .

Then we successively deduce the following consequences.

Corollary 4. For every p, q ∈ N ≥3 , none of the graphs W p,q and K 2,p contains the gem as induced minor.

Corollary 5. No graph of Incl(W i ) and of Incl({K 2,p } p∈N ≥3 ) contains the gem as induced minor, for every i ∈ N ≥3 .

The following observation will allow us to use Lemma 14, which deals with induced minors. Lemma 15. W i is a fundamental antichain, for every i ∈ N ≥3 .

Proof. Let i ∈ N ≥3 . We need to show that (Incl(W i ), ≤ ctr ) is a wqo. Let us call inner edge every edge of W p,q that is not incident with a pole, for every p, q ∈ N ≥3 . Observe that if a graph H ∈ Incl(W i ) has been obtained by contracting at least one edge incident with a pole, then H has a dominating vertex. According to Corollary 7, these graphs are wqo by ≤ ctr , therefore we will here consider graphs of Incl(W i ) that have been obtained by only contracting inner edges. We call I this class.

We first show that I is the union of the two following classes:

• the class I 0 of graphs that can be obtained by adding two non-adjacent dominating vertices to K i + D q for some q ∈ N ≥0 ; and

• the class I 1 of graphs that can be obtained by adding two non-adjacent dominating vertices to K i + S q for some q ∈ N ≥0 .

Again we use the notion of poles to denote the two dominating vertices added to construct graphs of I 0 and I 1 . A semipole is either a dominating vertex of D q (when dealing with graphs of I 0 ), or the dominating vertex of S q (when dealing with graphs of I 1 ).

Contracting an inner edge in W i,q clearly yields a graph of I 0 . Now, observe that any further contraction of an edge connecting a vertex of degree 4 to a semipole gives a graph of I 0 again. If, on the other hand, we contract the edge connecting the two semipoles, then we get a graph of I 1 . On a graph of I 1 , contracting an edge of the star (used in the construction of this graph) still gives a graph of I 1 . Therefore I = I 0 ∪ I 1 .

Let us assume that I is not wqo by ≤ ctr . Therefore it has an infinite antichain. As I = I 0 ∪ I 1 , one of I 0 and I 1 (at least) has an infinite antichain. Let A be such an infinite antichain.

We now look at vertices of graphs of A that are neither poles, nor semipoles, nor have degree 2. These vertices are the vertices of degree 2 of the copy of D q or the vertices of degree one of the copy of S q used in the construction of the graphs of A (depending whether A ⊆ I 0 or A ⊆ I 1 ). We call them inner vertices.

Let A and A be two graphs of A such that A has less inner vertices than A . These graphs exist since the elements of A are distinct. Let q be the number of inner vertices of A and q the one of A .

In both cases A ⊆ I 0 and A ⊆ I 1 we can obtain A from A by contracting q -q inner vertices of A to a semipole. This contradicts the fact that A is an antichain. Therefore (I, ≤ ctr ) is a wqo. This implies that W i is fundamental, as required.

We are now ready to prove Theorem 10.

Proof of Theorem 10. Let A i = K 2,i+1 for every i ∈ N ≥3 .

By the virtue of Corollary 6, {A i } i∈N ≥3 is a fundamental antichain, as well as W i , for every i ∈ N ≥3 , according to Lemma 15. Also, for every i ∈ N ≥3 , we have A i ≤ ctr H for every H ∈ W i (Lemma 13) and there are no other comparable pairs of elements in i∈N ≥3 W i ∪ {A i } i∈N ≥3 (Lemma 12 and Lemma 13).

Hence these sequences of graphs satisfy the requirements of Lemma 11, which implies that there is no canonical antichain for the contraction relation.

Theorem 6 (

 6 [START_REF] Ding | On canonical antichains[END_REF]). Under the subgraph relation, the class of finite graphs has a canonical antichain. Theorem 7 ([4]). Under the induced subgraph relation, the class of finite graphs does not have a canonical antichain. Theorem 8 ([4]). Under the induced subgraph relation, both the class of interval graphs and the class of bipartite permutation graphs have a canonical antichain. Theorem 9 ([8]). Under the multigraph contraction relation, the class of finite (loopless) multigraphs has a canonical antichain.

Figure 1 :

 1 Figure 1: A clique-cactus graph.

Figure 2 )

 2 and S = {K 1,r , r ∈ N}.

Figure 2 :

 2 Figure 2: Graphs of D.

Figure 3 :

 3 Figure 3: Antiholes antichain.

Lemma 7 .

 7 Let G ∈ Excl(D 2 ). Every cycle of G is either an induced cycle, or it induces a clique in G.

  Let us define three graph constructors stick : C → C, cycle : C → C, and clique : C → C. Given a sequence G 0 , . . . , G p-1 ∈ C (for some p ∈ N), if U denote the union of the graphs G 1 , . . . , G p-1 , then we define;

Figure 4 :

 4 Figure 4: The situation described in Lemma 11. Arrows are directed towards larger elements.

Figure 5 :

 5 Figure 5: The graph W 4,3 . Poles are drawn in white and semipoles in gray.

Observation 8 .

 8 Let H and G be two graphs. If both of them have a dominating vertex, then H is a contraction of G iff H is an induced minor of G.Lemma 14 ([1]).Graphs not containing the gem as induced minor are wqo by the induced minor relation.The following corollaries are direct consequences of Lemma 14, Observation 8 and Corollary 5.

Corollary 6 .

 6 Incl({K 2,p } p∈N ≥3 ) is wqo by ≤ ctr .Corollary 7. The graphs of Incl(W i ) with a dominating vertex are wqo by ≤ ctr , for every i ∈ N ≥3 .

  In a graph G, a contraction of the edge {u, v} ∈ E(G) is the operation which adds a new vertex adjacent to the neighbors of u and v and then deletes u and v. We say that a graph H is a contraction of a graph G whenever H can be obtained from G by a sequence of edge contractions, what we denote by H ≤ ctr G.A contraction model of a graph H in a graph G is function ϕ

  By Lemma 7 the cycle C is either an induced cycle, or it induces a clique in G. Let us treat these two cases separately. For contradiction we assume that V (G) \ V (C) is not empty and we call H 1 , . . . , H t the connected components of G \ C, for some t ∈ N ≥1 . Let us consider the graph G where H i , which is connected, has been contracted to a single vertex h First case: C induces a clique in G . Notice that C is then a maximal clique. Let u = h 1 . As C is maximal, there is a vertex v ∈ V (C) such that {u, v} ∈ E(G ). Let x and y be two neighbors of u on C (they exist since G is 2-connected). These vertices define two subpaths of C. Let R be the longest of these paths that contains v. Observe that in this case, R has at least three vertices. The union of {u, x}, {u, y} and R is a cycle of G that we call C . According to Lemma 7, this cycle is either induced or it induces a clique. As {u, v} ∈ E(G ), C cannot induce a clique in G. On the other hand, C is not an induced cycle as every pair of vertices of R are adjacent (and |V (R)| ≥ 3 as mentioned earlier). We reached the contradiction we were looking for. Second case: C is an induced cycle and has at least 4 vertices. Let i ∈ {1, . . . , t}.

however G may also contains other vertices. Let us consider G \ {v, x, v , y}. While G \ {v, v , x, y} contains a connected component adjacent to x or y, we contract it to this vertex (that we keep calling with the same name). Then, while it has a connected component adjacent to v but not v (respectively v but not v), we contract it to v (respectively v ), again keeping the same name for that vertex. Finally, the only remaining connected components (if any) are adjacent to exactly v and v : we contract each of them to a single vertex, adjacent to v and v . Notice that none of these operations create an edge connecting v to v , thus the subset {v, v , x, y} still induces a subgraph isomorphic to D 2 . Let us call G the obtained graph. As observed above, G consists of the subgraph isomorphic to D 2 induced by {v, v , x, y} plus k extra vertices of degree two, z 1 , . . . , z k , each of which is adjacent to v and v . In the case where k = 0, G is isomorphic to D 2 and we reached the contradiction we were looking for. Otherwise, we contract {v , y} (naming the resulting vertex v ), which produces a complete subgraph on vertices v, v , x, and then we contract {v, z i } for every i ∈ {2, . . . , k} (naming the resulting vertex v). The obtained graph is a complete graph on v, v , x, two vertices of which (that are v, v ) are adjacent to an extra vertex, v 1 . This graph is isomorphic to D 2 , as x is not adjacent to v 1 , therefore we reached a contradiction. Therefore C has either no chords or no non-chords in G. It is clear that in the first case C is an induced cycle of G and that in the second case it induces a clique. Lemma 8. Let G ∈ Excl(D 2 ) be a 2-connected graph. Then G is either a cycle, or a clique.

Proof. We assume that |V (G)| > 1, otherwise the result is trivial. Let C be a longest cycle of G. i , for every i ∈ {1, . . . , t}. Observe that G is 2-connected, given that G is 2-connected. Also, G ∈ Excl(D 2 ).

P n is the path on n vertices, for every n ∈ N.

R n is the multigraph obtained by doubling every edge of a path on n edges, for every n ∈ N.

The gem is the graph obtained by adding a dominating vertex to P

and K4 is the graph obtained by adding a vertex of degree 2 to K 4 .

Acknowledgements

The authors thank Jaros law B lasiok for inspiring discussions about the topic of this paper.

* The research was partially supported by the Foundation for Polish Science (Marcin Kamiński), the (Polish) National Science Centre grants SONATA UMO-2012/07/D/ST6/02432 (Marcin Kamiński and Jean-Florent Raymond) PRELUDIUM 2013/11/N/ST6/02706 (Jean-Florent Raymond), and by the Warsaw Center of Mathematics and Computer Science (Jean-Florent Raymond and Théophile Trunck).

From encodings to well-quasi-ordering. The following lemma will allow us to work on sequences in order to show that two graphs are comparable. Proof. Let H = H 1 , . . . , H p and G = G 1 , . . . , G q (for some positive integers p, q), and let H = cycle(H) and G = cycle(G). For the sake of readability we will refer to H i 's (respectively G i 's) either as elements of H (respectively G) or as subgraphs of H (respectively G).

If H ≤ ctr G, then there is, by definition of ≤ ctr , an increasing function ϕ : {1, . . . , p} → {1, . . . , q} such that ∀i ∈ {1, . . . , p} , H i ≤ ctr G ϕ (i). Therefore there is a sequence of edge contractions transforming G ϕ(i) into H i for every i ∈ {1, . . . , p}. Let us perform the following operations on G:

1. for every j ∈ {1, . . . , q} \ {ϕ(i), i ∈ {1, . . . , p}} we contract the subgraph G j to a single vertex v j and we then contract some edge incident with v j ;

2. for every i ∈ {1, . . . , p} we contract the subgraph G i in order to obtain the subgraph H ϕ (i).

Observe that after step 1., we obtain the graph cycle(G -), where G -can be obtained from G be deleting elements of indices in {1, . . . , q}\{ϕ(i), i ∈ {1, . . . , p}}.

Intuitively, we contracted the graphs that do not appear in H and removed their attachment point from the cycle. Then we replace in step 2. every graph of G - by its corresponding contraction of H. Therefore the graph obtained at the end is cycle(H), that is H, as required.

The cases (ii) and (iii) are very similar: H can be obtained from G by following the same operations as above.

Proof of Lemma 9. Let us assume by contradiction that (C, ≤ ctr ) is not a wqo. All decreasing sequences of this quasi-order are finite (as each contraction decreases the number of edges by at least one), therefore (C, ≤ ctr ) contains an infinite antichains. Let us consider a minimal antichain {A i } i∈N of (C, ≤ ctr ). Let B = i∈N dec(A i ), and let us show that (B, ≤ ctr ) is a wqo. For contradiction, let us assume that it is not a wqo and let {B i } i∈N be a minimal antichain of this quasi-order.

By definition of B, for every H ∈ B there is an integer i ∈ N such that H ≤ ctr A i (for instance, an integer i such that H ∈ dec(A i )). Therefore for every i ∈ N there is an integer π(i) such that B i ≤ ctr A π(i) . Let k ∈ N be the integer such that π(k) is minimum. Then the following sequence A = A 0 , . . . , A π(k)-1 , B k , B k+1 , . . . is an infinite antichain of (C, ≤ ctr ). Indeed, as both {A i } i∈N and {B i } i∈N are antichains, every pair of comparable graphs of A involves one graph of {A i } i∈{1,...,π(k)-1}