Tabea Rebafka
email: tabea.rebafka@upmc.fr.

Céline Lévy-Leduc
email: celine.levy-leduc@telecom-paristech.fr.

Member, IEEE. Maurice Charbit
email: maurice.charbit@telecom-paristech.fr.

OMP-type Algorithm with Structured Sparsity Patterns for Multipath Radar Signals

Keywords: DOA estimation, waveform recovery, multipath propagation, structured sparsity, OMP-type algorithm. I

A transmitted, unknown radar signal is observed at the receiver through more than one path in additive noise. The aim is to recover the waveform of the intercepted signal and to simultaneously estimate the direction of arrival (DOA). We propose an approach exploiting the parsimonious time-frequency representation of the signal by applying a new OMP-type algorithm for structured sparsity patterns. An important issue is the scalability of the proposed algorithm since high-dimensional models shall be used for radar signals. Monte-Carlo simulations for modulated signals illustrate the good performance of the method even for low signal-to-noise ratios and a gain of 20 dB for the DOA estimation compared to some elementary method.

account such as DOAs and relative Times of Arrival (TOA) of reflected paths, relative signal power levels on the different paths and additive noise variance.

In this paper, we shall consider a unique narrowband farfield source propagating through several paths, corrupted with additive Gaussian white noise and observed on a multi-sensor array. This issue is actually very involved since the observed signals on the different sensors are highly correlated and the different replicates are likely to overlap. This prevents us from using methods such as MUSIC, ESPRIT or MVDR which do not work very well with coherent signals [START_REF] Krim | Two decades of array signal processing research: The parametric approach[END_REF]. These approaches, based on subspace methods or on maximum likelihood estimation, have been extended by [START_REF] Zhang | An efficient DOA estimation method in multipath environment[END_REF], [START_REF] Vanderveen | Joint angle and delay estimation (jade) for multipath signals arriving at an antenna array[END_REF], [START_REF] Ma | Detection of coherent signals using weighted subspace smoothing[END_REF] and [START_REF] Shan | On spatial smoothing of estimation of coherent signals[END_REF] to deal with coherent sources. However, to the best of our knowledge, estimating the waveform, which is an infinite-dimensional parameter, in a nonparametric way has never been considered for multipath signals. To solve this problem, our main idea consists in exploiting the sparsity of the time-frequency representation of the waveform, which means that the selection of possibly useful signal components from a huge collection of candidate signals may be driven by a sparsity constraint. Such kind of approach was initially proposed by [START_REF] Fuchs | Multipath time-delay detection and estimation[END_REF] but only for the estimation of the time-delays, which is a finite dimensional parameter, and is extended here for estimating the DOAs and TOAs as well as the waveform.

For dealing with the estimation in sparse linear regression models, the Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] and the greedy orthogonal matching pursuit (OMP) ( [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF], [START_REF] Mallat | Adaptive time-frequency decompositions[END_REF]) have become very popular tools. In our case, an inspection of our model shows that the parameter vector is not sparse in an arbitrary way: its sparsity pattern has a specific structure. In order to include prior information concerning the sparsity structure, different approaches have recently been proposed. On the one hand, there are methods based on composed ℓ 1 /ℓ 2 -penalties (elastic nets [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]; fused Lasso [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF]; group Lasso [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]; composite absolute penalty [START_REF] Zhao | The composite absolute penalties family for grouped and hierarchical variable selection[END_REF];

overlapping groups [START_REF] Jacob | Group lasso with overlap and graph lasso[END_REF], [START_REF] Jenatton | Structured variable selection with sparsity-inducing norms[END_REF]), on the other hand, [START_REF] Lozano | Grouped orthogonal matching pursuit for variable selection and prediction[END_REF] proposed the group-OMP method for nonoverlapping groups, which is an extension of the orthogonal matching pursuit to structured solutions.

In this paper we explain how to extend the group-OMP to overlapping groups and how to deal with the scalability of our algorithm which is a crucial issue in the context of intercepted radar signals since the model dimension is usually very high. A simulation study shows that the proposed method performs well for composed and modulated signals even when the signal-to-noise ratio is low. A gain of more than 20 dB for the DOA estimation is achieved in comparison with some elementary method.

The paper is organized as follows. In Section II the mathematical model for radar signals is introduced, together with a reformulation of the problem in the form of a sparse, partly linear model. In Section III an OMP-type algorithm is developed which is well adapted to our model. Section IV deals with the scalability issue of the algorithm, while Section V presents an appropriate stopping criterion for the March 26, 2011 DRAFT OMP-type algorithm. Finally, the experimental results of Section VI illustrate the performance of the new algorithm and a comparison with some elementary method is provided. Conclusions are made in the final Section VII.

II. MODEL

A. Signal Modelling

Let us denote by s 0 (t) the original radar signal emitted by the source, where t is the time. Due to propagation and reflexion on obstacles, the signal arrives via several pathways on the receiver. With each pathway we associate a direction of arrival, say D u , a time-delay, say T u and a complex constant A u representing the attenuation of the signal on the u-th pathway. The signal is detected by an antenna with C sensors whose array response is denoted by a(•). Then the signal arriving at the sensor array at time t is a C-dimensional vector given by

s(t) = U u=1 A u s 0 (t -T u ) a(D u ) , (1)
where U denotes the number of propagation pathways, which is unknown. As the time lags T u are typically much shorter than the total length of the signal s 0 , the arriving signal s is the superposition of several delayed replicates s 0 with different amplitudes.

In practice, the signal s is observed with some additional noise, say ε, and at a sample frequency, say F s . With ∆ s = 1/F s the observations are hence given by

y m = s(m∆ s ) + ε m , m = 1, . . . , M , (2)
where ε m ∈ C C is a noise vector. Generally one considers circular gaussian independent random vectors with zero mean and covariance σ 2 I C .

B. Linear Regression Model

In practice, all model parameters are unknown, except the array response function a(•), which is determined by the array geometry or by calibration. As the relationship between the parameters is rather involved, we propose to reformulate the estimation problem by linearizing the model. This is obtained, on the one hand, by representing the waveform in an overcomplete basis and, on the other hand, by discretizing the parameter spaces of the other parameters. This leads to a linear model with a highdimensional parameter vector and nonlinear constraints. As this vector is structured and sparse, i.e. most entries are zero, regularization methods with structured sparsity inducing norms can be applied [START_REF] Zhao | The composite absolute penalties family for grouped and hierarchical variable selection[END_REF] .

March 26, 2011 DRAFT

1) Decomposition of the signal :

To estimate the waveform, the radar signal s 0 is represented in some given (overcomplete) basis or dictionary D = {ϕ j , j = 1, . . . , J}, such that

s 0 = J j=1 β j ϕ j , (3)
with appropriate coefficients β j ∈ C. Hence, estimating s 0 becomes recovering the coefficients β j . When using appropriate, large dictionaries, a small number of elements ϕ j is sufficient to reconstruct s 0 . In other words, it is assumed that the dictionary is such that most coefficients β j are zero and hence that the representation of s 0 is sparse.

2) Propagation model: To estimate the parameters A u , D u and T u by discretization, we introduce grids {τ 1 , . . . , τ P } and {θ 1 , . . . , θ Q } of potential values of the delay times T u and the angles of arrival D u respectively, as in [START_REF] Fuchs | Multipath time-delay detection and estimation[END_REF], [START_REF] Malioutov | A sparse signal reconstruction perspective for source localization with sensor arrays[END_REF]. Then s(t) can be rewritten as

s(t) = P p=1 Q q=1 α p,q s 0 (t -τ p ) a(θ q ) , (4)
with

α p,q = U u=1 A u 1{T u = τ p , D u = θ q } ,
where 1 is the indicator function. Clearly, there are only U coefficients α p,q that are nonzero. Hence, when fine grids are used, i.e. P and Q are large, this is another sparse estimation problem. As the number of pathways U is generally unknown, model selection is required, i.e. the estimation of the model order U .

By combining (3) and ( 4), the observation y m reads

y m = J j=1 P p=1 Q q=1 α p,q β j ϕ j (m∆ s -τ p ) a(θ q ) + ε m , (5)
for m = 1, . . . , M . It is convenient to represent the model in matrix notation. Let X be the matrix of size CM × JP Q, where all known quantities are stored. More precisely, the (j, p, q)-th predictor, i.e.

the column x j,p,q of X, is given by

x j,p,q =      ϕ j (∆ s -τ p ) a(θ q ) . . . ϕ j (M ∆ s -τ p ) a(θ q )      ∈ C CM .
In other words, x j,p,q corresponds to the signal received by the antenna (without noise) when the emitted signal is the element ϕ j arriving with delay τ p and angle θ q . Denote the observation vector by Y = (y T 1 , . . . , y T m ) T ∈ C CM and the noise vector by ε = ( ε T 1 , . . . , ε T m ) T ∈ C CM . Finally, let the Kronecker March 26, 2011 DRAFT product β ⊗α of the vectors β = (β 1 , . . . , β J ) T and α = (α 1,1 , . . . , α P,Q ) T be the JP Q-vector w defined as

w =      . . . w j,p,q . . .      =               β 1 α 1,1 β 1 α 1,2
. . . β j α p,q . . .

β J α P,Q               . ( 6
)
Then model ( 5) can be written as

Y = X(β ⊗ α) + ε with α ∈ C P Q , β ∈ C J . (7)
Finally, introduce the following set

W = {w ∈ C JP Q : w = β ⊗ α, α ∈ C P Q , β ∈ C J } ,
which is a subspace of C JP Q that can be parametrized using J + P Q variables. Remark that, given a vector w ∈ W, one can recover the vectors α and β only up to a multiplicative constant, which does not matter in our application. Thus, model ( 5) can be written as the following linear regression model

Y = Xw + ε with w ∈ W, (8)
where the design matrix X is known and the parameter vector w is to be estimated from the observations Y . The particularity of this regression model is that the search space for w is not the entire space C JP Q but the smaller, nonconvex set W.

III. ESTIMATION PROCEDURE

The problem now is to estimate the vectors α and β by fitting model [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], or equivalently, to estimate w in model [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF]. The advantage of ( 8) over [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] is the linearity of the model. Since the dimension of the space W is much higher than that of the observations Y , the minimization problem

min w∈W Y -Xw 2 (9)
is ill-posed. However, we note that w must be a sparse vector, since α and β are assumed sparse. It is well known that sparsity leads to stable estimation procedures based on a regularized fitting criterion with an ℓ 1 -penalty (Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]; basis pursuit [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]; Lars [START_REF] Efron | Least angle regression[END_REF]) or greedy algorithms (as matching pursuit and orthogonal matching pursuit (OMP) [START_REF] Pati | Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[END_REF], [START_REF] Chen | Basis pursuit[END_REF]). A closer inspection of W reveals some structure in the March 26, 2011 DRAFT sparsity patterns, that means, there are constraints on the distribution of the zero entries for all vectors w in W. More precisely, the set of indices S = {(j, p, q) : such that w j,p,q = 0} is called the sparsity pattern of some given vector w. Indeed, all elements w of W verify that α p ′ ,q ′ = 0 implies that the components w j,p ′ ,q ′ = 0 for all j, and likewise, if β j ′ = 0 then w j ′ ,p,q = 0 for all p and q. Thus, w must not be sparse in an arbitrary way, but its sparsity pattern has to respect some structure. Modern estimation procedures force specific structures of the sparse vector. On the one hand, there are methods based on composed ℓ 1 /ℓ 2 -penalties (elastic nets [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]; fused Lasso [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF]; group Lasso [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]; composite absolute penalty [START_REF] Zhao | The composite absolute penalties family for grouped and hierarchical variable selection[END_REF]; overlapping groups [START_REF] Jacob | Group lasso with overlap and graph lasso[END_REF], [START_REF] Jenatton | Structured variable selection with sparsity-inducing norms[END_REF]), on the other hand, the orthogonal matching pursuit can be extended to structured solutions (group-OMP for non-overlapping groups [START_REF] Lozano | Grouped orthogonal matching pursuit for variable selection and prediction[END_REF]).

Our problem can be related to a penalized minimization problem of the form

min w∈C JP Q Y -Xw 2 + Ω(w) , (10)
with a structured ℓ 1 /ℓ 2 -penalty Ω(w) (detailed below) in the sense that the sparsity pattern of a solution of this problem, say w, respects the required structure. However, there is no guarantee that the nonzero components of w satisfy relation [START_REF] Fuchs | Multipath time-delay detection and estimation[END_REF], that is, that there are coefficients α p,q and β j such that wj,p,q = α p,q β j and hence w may not belong to W. The problem lies in the nonzero components of w. To compute the solution of ( 10), different solutions have been proposed in the literature as the active set algorithm in in [START_REF] Jenatton | Active set algorithm for structured sparsity-inducing norms[END_REF] and [START_REF] Jenatton | Structured variable selection with sparsity-inducing norms[END_REF]. In this paper we propose an extension of the group-OMP [START_REF] Lozano | Grouped orthogonal matching pursuit for variable selection and prediction[END_REF] for overlapping groups.

Now, to fit the model given in [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] we propose the following procedure consisting of two steps, that are developed below in more detail.

Step 1. (Model selection) Solve problem [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] via an OMP-type algorithm to obtain a solution w with admissible sparsity structure. Identify the nonzero components wj,p,q and the indices A of the associated coefficients αp,q and βj that must be nonzero.

Step 2. (Estimation in reduced dimension) Reduce the dimension of the regression matrix X by keeping the predictor x j,p,q only if wj,p,q = 0. Then compute the least squares estimator in model ( 7) by the Nelder-Mead simplex method. This problem is now well-posed because of the largely reduced dimension.

A. Model selection step

To specify the penalty Ω(w) in [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF], we describe the structure of the sparsity pattern S of w ∈ W.

To this end, for a given couple (p, q), we denote the set of indices G α p,q = {(j, p, q), j ∈ {1, . . . , J}}, and likewise, for a given j, the set of indices G β j = {(j, p, q), p ∈ {1, . . . , P }, q ∈ {1, . . . , Q}}. It is March 26, 2011 DRAFT clear that the sparsity pattern S of any w in W can be written as the union of all sets G α p,q and G β j with α p,q = 0 and β j = 0.

Given a set of indices G, let w G 2 denote the ℓ 2 -norm of the subvector of w composed of the elements w j,p,q with (j, p, q) ∈ G. Then one can show that the solution, say w, of (10) with the following penalty

Ω(w) = λ 1 p,q w G α p,q 2 + λ 2 j w G β j 2 , (11)
where λ 1 , λ 2 > 0 are regularization parameters, has the required sparsity structure, i.e. the sparsity pattern of w can be written as the union of sets G α p,q and G β j . Note that these groups are possibly overlapping. Indeed, the penalty can be viewed as a mixture of ℓ 1 -and ℓ 2 -norms appropriate for overlapping group structures. More precisely, the penalty is the sum (i.e. the ℓ 1 -norm) of the ℓ 2 -norm of subvectors w G .

Likewise to the standard Lasso, the ℓ 1 -norm entails that for several groups G the terms w G 2 are zero, implying in turn that all entries of w G are exactly zero. By using two regularization parameters λ 1 > 0 and λ 2 > 0 one can use different degrees of sparsity for the vectors α and β.

To solve the problem [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] with penalty Ω(w) given by [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF] we propose an OMP-type algorithm. To describe the algorithm we introduce the collection G = {G g , g ∈ I} of all sets G α p,q and G β j indexed by a common index g with index set I. Note that the cardinality of I is P Q + J. For a set of indices A, denote by X A the matrix made of the corresponding predictors x j,p,q of X such that (j, p, q) ∈ A.

Recall that the principle of OMP consists in adding iteratively the predictor to the current solution w (t) which is the most correlated with the current residual r (t) = Y -Xw (t) . In [START_REF] Lozano | Grouped orthogonal matching pursuit for variable selection and prediction[END_REF] this principle is extended to selecting non-overlapping groups of variables, where G is a set of pairwise disjoint sets G g and one adds all variables of group

G g * if X T G g * r (t) 2 2 = max g∈I X T Gg r (t) 2 2 .
Proceeding in this way guarantees that at every step of the algorithm, the current solution respects the required sparsity structure, which means that the sparsity pattern of the current solution w (t) equals the union of some groups G g ∈ G.

Now this procedure can be extended to our context with overlapping group structures, where the aim remains the same: the sparsity pattern of every intermediate solution w (t) must be the union of some groups G g ∈ G. The difference to the group-OMP of [START_REF] Lozano | Grouped orthogonal matching pursuit for variable selection and prediction[END_REF] is at the level of the selection of variables that may enter the solution. Indeed, due to the overlapping group structure, the sets of variables that may be activated at the next step are not simply the groups G g ∈ G. To be more precise on this issue, we introduce the set Z (t) ⊂ I of indices of zero groups G g associated with the current solution w (t) , such that S (t) = ∪ g∈Z (t) G g is the sparsity pattern of w (t) . Clearly, the indices of the non zero entries of w (t) are given by A (t) = (S (t) ) c , where A c denotes the complementary set of A with respect to G. Now the candidate groups of components that may be activated at the next step are obtained by deleting March 26, 2011 DRAFT one of the indices from the zero groups Z (t) . That means, every potential group of indices has the form

P g0 = (∪ g∈Z (t) \{g0} G g ) c \A (t)
where g 0 ∈ Z (t) . Then, following the philosophy of OMP, one activates the components given by

P g * if X T P g * r (t) 2 2 = max g0∈Z (t) X T Pg 0 r (t) 2 2 .
Once a group P g * is selected, one derives the new set of active indices A (t+1) = A (t) ∪ P g * and updates the current solution by solving the least squares problem of reduced dimension w (t+1)

A (t+1) = arg min w Y -X A (t+1) w 2
2 . We notice that the algorithm can be viewed as a walk through a directed acyclic graph, which also appears in the algorithm proposed in [START_REF] Jenatton | Structured variable selection with sparsity-inducing norms[END_REF]. Formally, our OMP-type algorithm can be described as follows.

Algorithm 1 (OMP-type algorithm for overlapping groups): Initialization. Set t = 0.

Repeat until some stopping criterion is satisfied.

If t = 0 (initial step) Compute (j * , p * , q * ) = arg max (j,p,q)
x T (j,p,q) Y 2 2 . Denote by g * α and g * β the group indices associated with the groups G α p * ,q * and G β j * , respectively Put the set of zero group indices Z (1) = I\{g * α , g * β }. Put the set of active indices A (1) = {(j * , p * , q * )}.

Else

For each g 0 ∈ Z (t) set

P g0 = ∪ g∈Z (t) \{g0} G g c \A (t) .
Compute g * = arg max

g0∈Z (t) (X Pg 0 ) T r (t) 2 2 . Update Z (t+1) = Z (t) \{g * } and A (t+1) = A (t) ∪ P g * .

End

Update the current solution w(t+1)

A (t+1) = arg min w∈C |A (t+1) | Y -X A (t+1) w 2 2 . ( 12
)
Update the residual r (t+1) = Y -X A (t+1) w(t+1) .

Set t = t + 1.

End

Note that in any case the algorithm stops after a finite number of iterations, more exactly after J+P Q-1

steps. Then all groups are activated and the solution w is the ordinary least squares estimator. The important issue of an appropriate stopping criterion for the algorithm is addressed in Section V. Clearly, the aim is to stop the algorithm earlier to provide a relevant, sparse solution.

March 26, 2011 DRAFT

We remark that the proposed algorithm is a natural extension of the orthogonal matching pursuit.

Moreover it is similar to the active set algorithm described in [START_REF] Jenatton | Structured variable selection with sparsity-inducing norms[END_REF], which provides an interesting mathematical analysis of the optimization problem given by [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]. Nevertheless, our OMP-type algorithm provides an approach which is feasible in practice and computationally faster than the active set algorithm.

Indeed, in [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] the solution is updated by the ordinary least squares estimator that is known explicitly, whereas [START_REF] Jenatton | Structured variable selection with sparsity-inducing norms[END_REF] proposes to solve the penalized problem of reduced dimension by a more expensive secondorder cone programming.

B. Solution of the reduced problem

The solution of ( 10) is a sparse vector, say w, whose sparsity pattern serves to identify the indices of the nonzero coefficients α p,q and β j . More precisely, if wj,p,q = 0, this implies that the associated coefficients α p,q and β j are both non zero. Conversely, wj,p,q = 0 implies that at least one of the coefficients is 0.

To obtain estimates of these nonzero coefficients, we reduce the dimension of the regression matrix X by keeping the predictor x j,p,q only if wj,p,q = 0. Denote the resulting matrix by X red . Then compute the least squares estimator by minimizing

L(α, β) = Y -X red ( β ⊗ α) 2 , (13)
where α and β are the associated αand β-vectors of reduced dimension P Q and J, respectively This minimization is a well-posed problem as the dimension of X red is much smaller than that of X.

To compute the solution of (13) one can use the fact that the model is linear in α for fixed β and conversely. That means for fixed β we define the matrix X( β) with P Q columns by

X red ( β) = J j=1 βj (x red j,1,1 , . . . , x red j, P , Q) .
Then the minimum of α → L(α, β) is the ordinary least squares estimator given by

αOLS ( β) = (X red ( β) T X red ( β)) -1 X red ( β) T Y .
Furthermore, as α and β are identifiable only up to some multiplicative constant, one can set β1 = 1, for instance. Finally, it remains to minimize ( β2 , . . . , β J ) → L αOLS ((1, β2 , . . . , β J ) T ), (1, β2 , . . . , β J ) T which can be done by the Nelder-Mead simplex method. This is feasible since the number of unknown variables, i.e. J -1, is quite low.

March 26, 2011 DRAFT

IV. SCALABILITY

To obtain a good representation of the waveform, a large dictionary shall be used. Likewise, to avoid biased estimators of angles and time delays, we may use fine grids. However, large dictionary sizes J and large grid sizes P Q entail a huge regression matrix X having JP Q columns (see Section VI where JP Q ≈ 6.2 10 9 ). This raises computational difficulties, namely concerning the storage of the regression matrix X.

A closer look at the OMP-type algorithm shows that there are no computations involving the entire regression matrix. Essentially, the algorithm considers one by one all possible groups of variables G ∈ G that may enter the solution in the present iteration. For each group of variables G the correlation of the associated predictors with the current residual r (t) = Y -X G w(t) is computed. Hence, instead of the whole matrix X, only the predictors of group G are required.

In short, the storage of X can be avoided by recomputing the required predictors at each iteration.

With such a programming there are almost no limits on the size of the regression matrix, and thus almost arbitrary dictionaries and grids may be used.

To be more precise, let us have a look at the dictionary used in the following simulation study. The dictionary is composed of sinusoids with different frequencies f , different lengths l and different starting points s. The index j is hence a triple index (f, l, s). Now instead of storing the whole vector x j,p,q , we just use the triple index (f, l, s) to reconstruct the dictionary element ϕ j , whenever it is needed. Then further translation and multiplication with the appropriate steering vector yield the required predictor x j,p,q . Finally, we only have to store the steering matrix for all possible angles of arrival, which is a matrix of very moderate size C × P . With this knowledge, the predictor x j,p,q is easily reconstructed at any iteration of the algorithm.

In regression problems it is common to work with a normalized regression matrix, that is all columns have mean zero and standard deviation one. To avoid the computation of the normalization constants for every predictor, they may be computed just once and stored for later usage. Note that the normalization constants do not depend neither on the starting point of the dictionary element nor on the delay on the pathway but only on the frequency f and the length l of the dictionary element and the angle of arrival.

Hence, the number of constants to store remains reasonable.

The most time-consuming step in the algorithm is the scan of all predictors (or more precisely, of all relevant groups of predictors) in every iteration, to identify the one that is the most correlated with the current residual. It is noteworthy that this step can be speeded up by parallelization of the algorithm.

March 26, 2011 DRAFT Indeed, there is no specific order to visit the predictors and the problem can be split in several independent subtasks.

V. STOPPING CRITERION

The OMP-type algorithm requires an appropriate stopping criterion. A common approach to conceive a good stopping criterion is based on cross-validation [START_REF] Zhao | The composite absolute penalties family for grouped and hierarchical variable selection[END_REF]. Note that in our context, we do not have to deal with a single, but with two sparse vectors, namely α and β and their degrees of sparsity may not be equal. Hence, cross-validation would be computationally demanding.

Here we propose a stopping criterion that on the one hand takes into account available prior information on the number of propagation pathways (and hence on the number of α-variables) and on the other hand has a data-driven component to determine the total number of non zero α p,q -and β j -coefficients.

A. Sparsity of α p,q -coefficients

In the context of intercepted radar signals, one generally expects a small number of relevant propagation pathways. Indeed, we are mainly interested in the recovery of the direct path to determine the principal DOA. It is hence convenient to introduce an upper limit for the number of paths U max and to modify the OMP-type algorithm such that at every iteration the number of activated α p,q -coefficients is checked.

If their number has achieved U max , then we stop activating further α p,q -coefficients and concentrate on activating only β j -coefficients.

In simulations we observed that this procedure yields satisfactory results. The only important assumption for a good performance is that the direct path is not too weak compared to the indirect pathways, since the algorithm activates the components by their importance in the signal.

Note that, if U max is smaller than the 'true' number of propagation pathways, then the observed signal Y cannot be completely explained. In other words, the solution obtained by the OMP-type algorithm is not the optimal solution of the minimization problem [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF].

B. Control of the squared error

The degree of sparsity in the β-variables depends on the chosen dictionary. In general, very few prior information is available on the number J of basis elements used to decompose the radar signal.

A simple but useful tool for model selection consists in studying the evolution of the squared error, that is the evolution of the square norm of the current residual L( w(t) ) = Y -X w(t) 2 2 , where w(t)

March 26, 2011 DRAFT denotes the solution at the t-th iteration. Obviously, the squared error L( w(t) ) is monotone decreasing over the iterations.

The OMP-type algorithm always adds the group of variables that are the most correlated with the current residual, that is the components that are likely to diminish the squared error the most. This explains why the algorithm tends to first activating all 'true' variables, and afterwards selecting variables that are not in the true model. Hence the problem of selecting the model order becomes the detection of the point when all 'true' variables are activated and the algorithm starts adding pure noise variables.

In fact, activating a 'true' variable generally improves the squared error a lot, since the new variable March 26, 2011 DRAFT explains a considerable part of the observations. Whereas adding a variable which is not in the true model is less beneficial and the squared error hardly decreases. Consequently, first the squared error L( w(t) )

decreases steeply until the point where all 'true' variables are activated. Then the graph of the squared error becomes very flat. This change of the slope yields a kink in the graph of the squared error L( w(t) ).

Thus, the kink indicates the number of 'true' variables.

An efficient way to determine the kink consists in computing the convex hull of the graph. Then the kink is given by the endpoint of the longest linear piece of the convex hull [START_REF] Lavielle | Using penalized contrasts for the change-point problem[END_REF], [START_REF] Lavielle | Semiparametric estimation of the frequency of unknown periodic functions and its application to laser vibrometry signals[END_REF].

We may hence proceed as follows. The OMP-type algorithm is run with some fixed U max and stopped after a large number of iterations. It is useful to store all intermediate solutions w(t) and the associated square losses. Afterwards we determine the final solution as the intermediate solution w(t) where the squared error presents a kink.

Clearly, the kink is especially pronounced when the signal-to-noise ratio is high. Conversely, when there is too much noise in the data, then it may be difficult to make out a kink. Besides, when the true model contains a relatively 'weak' variable, that is, its contribution to the signal is rather low, then the number of variables may be underestimated as the 'weak' variable can be confounded with noise.

VI. EXPERIMENTAL RESULTS

A. Setting

We illustrate the algorithm by an example. Let us consider a radar signal made of three subsequent sinusoids with different frequencies (100, 140, 110 MHz) and a gap between the two last sinus waves (sinusoids with starting points 0, 150, 300 ns and lengths 150, 100, 120ns), see Figure 1 and (c), where the superposition of the sinusoids results from the multiple propagation pathways. From these figures it is clear that the multiple paths have a strong impact on the signal and hence cannot be neglected in the analysis. The signal is corrupted by additive noise with a signal-to-noise ratio of 5 dB, see Figure 1(d). The signal-to-noise ratio is determined on the direct propagation pathway, defined as the 10 log (energy of the signal on the first path/variance of the noise) as used in [START_REF] Manickam | A least squares algorithm for multipath time-delay estimation[END_REF] or [START_REF] Fuchs | Multipath time-delay detection and estimation[END_REF].

We choose a Fourier dictionary containing sine and cosine functions with frequencies 90, 95,. . . , 150

MHz, signal lengths 80, 85,. . . , 200 ns and starting points 0, 3δ, 6δ,. . . , 255δ ns, where δ = 1/F s ≈ 0.78 ns. The size of the dictionary is thus J = 166 400. For the angles and the delays we use the grids {1 • , 2 • ,. . . , 90 • } and {0, 3δ, 6δ,. . . , 255δ}, respectively Hence, w is a vector of dimension 3.83 • 10 9 .

B. Signal reconstruction by the OMP-type algorithm

In the simulations we set U max = 2, that is, at best we recover two of the three pathways and we stop the algorithm after 9 iterations. Figure 1(e)-(h) illustrate the first four iterations of the OMP-type algorithm. In the first iteration the algorithm selects a single sinusoid. Next, a second pathway is selected leading to two overlapping sinusoids. As now two α-variables are activated and U max = 2, the algorithm can only add further β-variables, that is, further dictionary elements. Thus, in the third iteration another sinusoid is added which appears directly on the two pathways. Finally, another sinusoid is activated.

We note that the graph in (h) is quite similar to the one in (c). Nevertheless, the reconstruction is not perfect, as we allowed the algorithm to find only two pathways, while indeed there are three. In this example the second pathway corresponding to the angle of arrival of 13 • is missing.

C. Model selection device

As mentioned above, we stopped the algorithm after 9 iterations. Then we used the graph of the squared error to determine the number of variables in the model by the point where the curve exhibits a kink.

The kink is determined by computing the convex hull of the graph. Figure 2 illustrates the evolution of the squared loss for a sample with a signal-to-noise ratio of 5 dB. The curve has a slight kink at 5, which is the appropriate number of variables in this case, as the signal is composed of three sinusoids and we set U max = 2. Table I shows the empirical mean and standard deviation of the total number of selected variables, the number of αand β-coefficients for different signal-to-noise ratios over 50 repetitions. Obviously, in this example the model selection tool works very well, as the model selection device successfully distinguishes 'true' variables from noisy variables, even though the algorithm is prevented from selecting all 'true' variables. Hence, the introduction of an upper bound U max for the number of α-coefficients

does not corrupt the model selection.

D. Comparison with an elementary method

A simple, elementary method for estimating the principal DOA in the radar setting consists in first estimating the TOA by thresholding and then using the point of observation at the TOA to fit the angle of arrival. More precisely, to detect the TOA by the elementary method, the absolute value of the observed signal at some time point t is compared to some threshold, here 3σ, where σ 2 is the variance of the noise distribution that is supposed to be known for the elementary method. The first point t where this happens is considered as the TOA. This detection method fails when the signal-to-noise ratio is too low, since then the signal is too weak compared to the noise such that the observed signal never exceeds the threshold. And when no TOA is detected, the elementary method does not provide an estimate of the DOA.

We compare this elementary method to our OMP-type algorithm in the following setting with an FSK-modulated Barker sequence. We consider the 7-bit Barker sequence [1, 1, -1, -1, 1,-

E. Parameter Estimates

It is important to note that in contrast to the elementary method the OMP-type algorithm provides estimates not only of the DOA and the TOA, but of all the other parameters including the waveform.

Indeed, in the setting described above our algorithm identifies several consecutive sinusoids. Table II presents the empirical mean and standard deviation of the estimates obtained by the elementary method and Table III gives the results for the OMP-type method, namely for the estimates of the angles of arrival, the associated delays, the frequencies of the four detected sinusoids, their lengths and their starting points.

Estimates are especially accurate for the angles of arrival and the frequencies, whereas the estimates of the lengths of the sinusoids and their starting points are less accurate. There is a clear tendency to overestimate the length, especially for the shortest sinusoid. Indeed, for high noise levels it is not clear if the OMP-type algorithm detects the shortest (and less energetic) sinusoid. The fourth sinusoid generally has not even the correct frequency, so it is possible that the algorithm adds any non relevant predictor, because the shortest sinusoid is too weak for detection.

We observe that the reconstructed radar signal ŝ has not exactly the form of a modulated signal due to the overlapping sinusoids. It is hence up to the user to make such an interpretation of the signal.

Nevertheless, we emphasize that the algorithm works for any type of modulation (FSK or PSK) without prior knowledge on the modulation and it provides a good idea of the waveform.

Fig. 1 .

 1 Fig. 1. Radar signal and its reconstruction. All signals, except (a), are on the first sensor. (e)-(h) represent the real part.

  (a). The signal is detected by a uniform linear array (ULA) of 5 sensors separated by half a wavelength of the actual narrowband source signals. The sample frequency is F s = 1.28 GHz and the total duration of observation is 0.8 µs. Thus, we have M = 1024 observation points and Y is a vector of length 5120. Furthermore, we consider three pathways (U = 3) with unknown angles of arrivals (9 • , 13 • , 38 • ) and delays (60, 95, 120ns). The signal (without noise) received at the first sensor of the antenna is presented in Figure1(b)

Fig. 2 .

 2 Fig.2. Graph of the squared error of a dataset with a signal-to-noise ratio of 5 dB. The mark indicates the kink, that is the selected number of model variables.

Fig. 3 .

 3 Fig. 3. Comparison of RMSE of the DOA and TOA for the elementary method and the new estimator for different signal-to-noise ratios applied on datasets of size 5 x 1280. The RMSE values are based on 50 repetitions per noise level.

Figure 3

 3 Figure 3 compares the root mean square error (RMSE) of the DOA and TOA estimates obtained by both methods based on 50 simulated datasets for each noise level. Compared to the elementary procedure the OMP-type algorithm achieves a gain of around 20 dB. This gain is mainly due to the fact that the algorithm takes into account all observations and not only a single point. For low signal-to-noise ratios the elementary method completely fails.

TABLE I EMPIRICAL

 I MEAN AND STANDARD DEVIATION OF THE TOTAL NUMBER OF SELECTED VARIABLES OF α-AND β-COEFFICIENTS WITH UMAX = 2.

		9000
		8000
		7000
	squared error	5000 6000
		4000
		3000
		0 2000	1	2	3	4	5	6	7	8	9	10
					number of active variables
				SNR	total nb	nb of α's nb of β's
			5 dB	5 (0)			2 (0)		3 (0)
			-5 dB	5 (0)			2 (0)		3 (0)
			-15 dB 5.1 (0.3)		2 (0)		3.1 (0.3)
	March 26, 2011											DRAFT

  Again, the sample frequency is F s = 1.28 GHz and we use a ULA of 5 sensors separated by half a wavelength of the actual narrowband source signals. The observation time is 1 µs, such that Y is a vector of length 6400. The dictionary contains the sine and cosine functions with frequencies 200, 202,. . . , 240 MHz, lengths 50, 55,. . . , 170 ns and starting points 0, δ, 2δ,. . . , 255δ ns, where δ = 1/F s . The size of the dictionary is thus J = 268 800. We use the grids {1 • , 2 • ,. . . , 90 • } and {0, δ, 2δ,. . . , 255δ} for the angles and the delays, respectively Hence, w is a vector of dimension 6.19 • 10 9 . Data are simulated with different levels of the signal-to-noise ratio. Again we choose U max = 2.

		TABLE II
	EMPIRICAL MEAN AND STANDARD DEVIATION (IN PARENTHESES) OF THE PARAMETER ESTIMATES OBTAINED BY THE
		ELEMENTARY METHOD.
		Elementary Method
		1st angle (DOA) 1st delay (TOA)
	true	6 •	10 ns
	SNR 20 dB 5.99 (0.26)	10.94 (0)
	SNR 0 dB	6.81 (8.43)	22.31 (12.11)
	March 26, 2011		DRAFT

TABLE III EMPIRICAL

 III MEAN AND STANDARD DEVIATION (IN PARENTHESES) OF THE PARAMETER ESTIMATES OBTAINED BY THE OMP-TYPE ALGORITHM.

			OMP-type algorithm
		1st angle	1st delay	2nd angle	2nd delay
	true	6 •	10 ns	42 •	40 ns
	20 dB 6.00 (0)	10.94 (0)	42.0 (0)	40.6 (0)
	0 dB	5.94 (0.24)	10.5 (1.34)	41.9 (0.27)	40.2 (1.26)
		1st frequ.	2nd frequ.	3rd frequ.	4th frequ.
	true	210 MHz	230 MHz	210 MHz	230 MHz
	20 dB 230.0 (0)	210.0 (0)	230.0 (0)	210.0 (0)
	0 dB	230.0 (0)	210.0 (0)	230.1 (0.40) 210.0 (0)
		1st length	2nd length	3rd length	4th length
	true	160 ns	160 ns	80 ns	160 ns
	20 dB 160.0 (0)	160.6 (1.64) 162.4 (12.1) 170.0 (0)
	0 dB	163.1 (4.62) 165.5 (4.20) 162.4 (15.1) 170.0 (0)
		1st start	2nd start	3rd start	4th start
	true	0 ns	160 ns	320 ns	400 ns
	20 dB 0 (0)	158.8 (1.54) 265.3 (12.9) 389.1 (0)
	0 dB	0 (0)	153.1 (5.37) 271.0 (19.7) 389.5 (1.34)
	March 26, 2011				DRAFT

VII. CONCLUSION

The algorithm proposed in this paper is to the best of our knowledge the first method addressing the problem of estimating both the DOA and the unknown waveform of a transmitted radar signal in the presence of multiple propagation pathways. The new method is in the spirit of OMP extended to more sophisticated sparsity structures. The scaling of the OMP-type algorithm has been considered, which is necessary for unbiased estimation. A simulation study illustrates the good performance of the new method and exhibits a considerable improvement with respect to some elementary method.

We would like to emphasize the general character of our method, which is of much interest for the intercepted radar signal setting. First, the method proposed here is not restricted to a specific array geometry. Second and more importantly, the approach can be adapted to a large variety of signals such as sinusoids, chirps, composed signals etc. by using a dictionary where a sparse representation of the signal is possible.

This work was financially supported by the French National Agency: Direction Générale de l'Armement (DGA).


Julie Delon

AND Julien Salomon

Andrei Sobolevski

MINIMUM-WEIGHT PERFECT MATCHING FOR NON-INTRINSIC DISTANCES ON THE LINE

We consider a minimum-weight perfect matching problem on the line and establish a "bottom-up" recursion relation for weights of partial minimum-weight matchings.

Introduction

We start with recalling a few notions from combinatorial optimization on graphs. A matching in an undirected graph is any set of its mutually disjoint edges: no two edges from such set can share a vertex. A matching is called perfect if it involves all vertices of the graph (the number of vertices is then necessarily even).

Depending on the structure of the graph, perfect matchings may be many. Suppose that edges of a graph are endowed with real weights; then it makes sense to look for a perfect matching composed from a set of edges with a minimum sum of weights. In this note we treat a particular case of this minimum-weight perfect matching problem where the graph is complete, all its vertices are located on a line, and edge weights are related to distances between the vertices along the line.

A bipartite version of this problem, in which vertices are divided into two equal classes and edges of the matching must connect vertices of different class, reduces to transport optimization. For a particular class of cost functions (those of "concave type") this problem has been thoroughly treated in the measure-theoretic setting in [START_REF] Mccann | Exact solutions to the transportation problem on the line[END_REF]. Similar problems have also long been considered in the algorithmics literature for the specific case of the distance |x -y|, assuming that the measures are discrete [START_REF] Aggarwal | Efficient minimum cost matching using quadrangle inequality[END_REF][START_REF] Karp | Two special cases of the assignment problem[END_REF][START_REF] Werman | Bipartite graph matching for points on a line or a circle[END_REF]. An algorithm for a general cost function of a concave type has been proposed recently in [START_REF] Delon | Local matching indicators for concave transport costs[END_REF][START_REF]Local matching indicators for transport problems with concave costs[END_REF].

In the present note we consider the minimum-weight perfect matching problem on a complete graph without assuming it to be bipartite. However, for the class of weight functions generated by distances, the non-bipartite problem turns out to be essentially equivalent to a bipartite problem with alternating points due to a "no-crossing" property of the optimal matching. This observation, which we owe to S. Nechaev, allows to employ the theory developed for the bipartite case in [START_REF] Delon | Local matching indicators for concave transport costs[END_REF][START_REF]Local matching indicators for transport problems with concave costs[END_REF] to the non-bipartite matching.

The main contribution of the present note is a specific "bottom-up" recursion relation [START_REF] Edmonds | Maximum matching and a polyhedron with 0, 1-vertices[END_REF] for partial minimum-weight matchings. This recursion follows from a "localization" property of minimum-cost perfect matchings for concave cost functions (Theorem 5) and provides a new perspective on the construction of [START_REF] Delon | Local matching indicators for concave transport costs[END_REF][START_REF]Local matching indicators for transport problems with concave costs[END_REF].

Supported by Agence Nationale de la Recherche under project ANR07-01-0235 OTARIE and the Russian Foundation for Basic Research under project 07-01-92217-CNRSL-a.

The note is organized as follows. Section 2 is a review of the basic construction of metrics on the real line. In Section 3 we recall the problem of minimumweight perfect matching and cite a few useful results about the structure of its solution. Section 4 contains the key technical result, a kind of localization principle for minimum-weight matching. Finally in Section 5 we derive the recursive relation for partial minimum-weight matchings. We also discuss the ensuing algorithm and compare it to modern variants of the Edmonds blossom algorithm for the minimum-weight perfect matching problem.

It is our pleasure to thank the organizers of the conference Optimization and Stochastic Methods for Spatially Distributed Information (St Petersburg, EIMI, May 2010), where an earlier version of this work was presented.

Intrinsic and non-intrinsic distances on the line

Recall the usual axioms for a distance d(•, •) on the real line R: for all x, y, z ∈ R,

(D1) d(x, y) ≥ 0 with d(x, y) = 0 iff x = y; (D2) d(x, y) = d(y, x); (D3) d(x, y) + d(y, z) ≥ d(x, z).
These axioms are satisfied by the distance d(x, y) = |x-y| as well as by any distance of the form (1)

d g (x, y) = g(|x -y|),
where g(•) is a nonnegative concave function defined for all x ≥ 0 such that g(x) = 0 iff x = 0. Here concavity means that g(λx

+ (1 -λ)y) ≥ λg(x) + (1 -λ)g(y)
for all x, y ≥ 0 and all λ ∈ [0, 1]. Note that the distance d g is homogeneous with respect to shifts: (D4) d g (x + t, y + t) = d g (x, y) for all x, y, t ∈ R.

Conversely, any distance satisfying the axioms (D1)-(D4) has the form (1) for a suitable nonnegative concave function g: indeed, take g(x) = d(0, x) for x ≥ 0 and check that concavity follows from (D3).

An example of this construction is the distance |x -y| α with 0 ≤ α ≤ 1, where

|x -y| 0 = 0, x = y, 1, x = y
is the "discrete distance"; here we are mostly interested in the case 0 < α < 1.

An important property of a distance d(•, •) is whether it is intrinsic. To recall the corresponding definition, take two distinct points x, y ∈ R and connect them with parameterized curves taking values x(t) in R, i.e., suppose that 0 ≤ t ≤ 1, x(0) = x, and x(1) = y. By the triangle inequality,

(2) d(x, y) ≤ inf 0≤i<N d(x(t i ), d(x(t i+1 )),
where the infimum is taken over all curves connecting x to y and all meshes 0 = t 2) is an equality for all x, y (see, e.g., [START_REF] Burago | A Course in Metric Geometry[END_REF]).

0 < t 1 < t 2 < • • • < t N = 1 with N ≥ 1. The distance d is called intrinsic if (
The distance |x -y| and its scalar multiples are the only homogeneous intrinsic distances in R. For distances |x -y| α with 0 < α < 1, connecting curves have infinite length unless x = y, and therefore these distances are not intrinsic.

The intuition behind this notion is that whereas the geometry of R equipped with an intrinsic distance is fully determined by the distance, the same line R with a non-intrinsic distance should be viewed as embedded in an auxiliary space of a larger dimension, and geometry on R is induced by that in the embedding space. According to the Assouad embedding theorem (see, e.g., [START_REF] Heinonen | Lectures on Analysis in Metric Spaces[END_REF]), for the distance |x-y| α on R the embedding space has dimension of the order 1/α for small α.

Minimum-weight perfect matchings

Consider an even number of points x 1 < x 2 < • • • < x 2n on the real line R equipped with a distance d and look for a minimum-weight perfect matching in a complete graph K 2n on these points, where the weight of an edge connecting x i to x j is d(x i , x j ). It is convenient to represent R with a horizontal interval and use arcs in the upper halfplane to show the edges of the graph K 2n .

The following lemma, adapted from McCann [10, Lemma 2.1], allows to describe the structure of minimum-weight matchings.

Lemma 1. Suppose the distance has the form d g (1) with a strictly concave function g (i.e., g(λx + (1 -λ)y) > λg(x) + (1 -λ)g(y) for all x, y ≥ 0, x = y, and 0 < λ < 1). Then the inequality

(3) d g (x 1 , y 1 ) + d g (x 2 , y 2 ) ≤ d g (x 1 , y 2 ) + d g (x 2 , y 1 )
implies that the intervals connecting x 1 to y 1 and x 2 to y 2 are either disjoint or one of them is contained in the other.

Proof. Let arcs representing the matching be directed from x's to y's and assume, without loss of generality, that x 1 is leftmost. Then the configuration , where x 1 < y 2 < y 1 < x 2 , cannot satisfy (3) because a strictly concave function g must be strictly growing and therefore d g (x 1 , y 2 ) < d g (x 1 , y 1 ) and d g (x 2 , y 1 ) < d g (x 2 , y 2 ). To rule out configuration , where x 1 < x 2 < y 1 < y 2 , choose 0 < λ < 1 such that

y 1 -x 1 = (1 -λ)(y 2 -x 1 ) + λ(y 1 -x 2 ).
Then, since (y 1 -x 1 ) + (y 2 -x 2 ) = (y 2 -x 1 ) + (y 1 -x 2 ), we obtain

y 2 -x 2 = λ(y 2 -x 1 ) + (1 -λ)(y 1 -x 2 )
and can further use the strict concavity of g to get

d g (x 1 , y 1 ) > (1 -λ)d g (x 1 , y 2 ) + λd g (x 2 , y 1 ), d g (x 2 , y 2 ) > λd g (x 1 , y 2 ) + (1 -λ)d g (x 2 , y 1 ).
The sum of these inequalities contradicts [START_REF] Cook | Computing minimum-weight perfect matchings[END_REF]. All the other configurations where the points x 2 , y 1 , y 2 are located to the right of x 1 , namely , , , , are consistent with (3).

It is possible to show that conversely, for an arbitrary bivariate function d(x, y) the property of Lemma 1 together with homogeneity (D4) imply that d has the form (1) for a suitable concave function g [10, Lemma B4].

We call a matching nested if, for any two arcs (x i , x j ) and (x i ′ , x j ′ ) that are present in the matching, the corresponding intervals in R are either disjoint or one of them is contained in the other.

Theorem 2 ([1, 10]). A minimum-weight matching is nested.

Proof. This statement essentially follows from Lemma 1, because whenever two arcs (x i , x j ) and (x i ′ , x j ′ ) are crossed, the sum of distances corresponding to the "uncrossed" arcs (x i , x j ′ ), (x i ′ , x j ) must be smaller according to [START_REF] Cook | Computing minimum-weight perfect matchings[END_REF].

However uncrossing may introduce new crossings with other arcs, and it remains to be checked that it does decrease the total number of crossings. This is done using an argument adapted from [1, Lemma 1].

Let an arc (x k , x ℓ ) cross either of the arcs (x i , x j ′ ) or (x i ′ , x j ) obtained after uncrossing the arcs (x i , x j ) and (x i ′ , x j ′ ); we would like to prove that before uncrossing the arc (x k , x ℓ ) had the same number of crossings with (x i , x j ) and (x i ′ , x j ′ ). It suffices to consider the following two cases.

Case

x i < x j ′ < x j < x i ′ . Assume that (x k , x ℓ ) crosses (x i , x j ′ ). Then (x k , x ℓ )
, shown as a dotted arc, may be situated in one of the following ways: after uncrossing before uncrossing

x i x j ′ x j x i ′ x i x j ′ x j x i ′
We see that the number of crossings between (x k , x ℓ ) and the other two arcs is the same before and after uncrossing. The case where (

x k , x ℓ ) crosses (x i ′ , x j ) is symmetrical. Case x i < x i ′ < x j < x j ′ .
The following situations of the dotted arc (x k , x ℓ ) are possible: after uncrossing before uncrossing

x i x i ′ x j x j ′ x i x i ′ x j x j ′
Again the number of crossings between (x k , x ℓ ) and the other two arcs is the same before and after uncrossing. It follows that each uncrossing removes exactly one crossing from the matching, and therefore any possible sequence of uncrossings leads in a finite number of steps to a nested matching with a strictly smaller weight. In other words, it suffices to look for the minimum-weight matching only among the nested ones. This result implies that the minimum-weight perfect matching problem is essentially bipartite, and it is indeed the bipartite setting that is considered in [START_REF] Aggarwal | Efficient minimum cost matching using quadrangle inequality[END_REF][START_REF] Delon | Local matching indicators for concave transport costs[END_REF][START_REF]Local matching indicators for transport problems with concave costs[END_REF][START_REF] Karp | Two special cases of the assignment problem[END_REF][START_REF] Mccann | Exact solutions to the transportation problem on the line[END_REF][START_REF] Werman | Bipartite graph matching for points on a line or a circle[END_REF] Corollary 3. In a minimum-weight perfect matching, points with even numbers are matched to points with odd numbers.

Proof. Suppose on the contrary that x i is matched to x j with i and j both even or both odd. The interval connecting x i to x j contains an odd number |i -j| -1 of points, one of which has to be matched to a point outside and thus cause a crossing with (x i , x j ). Since a minimum-weight matching is nested, the result follows.

Observe finally the following simple form of the Bellman optimality principle. Proof. Indeed, if one could rematch these points achieving a smaller total weight, then the full original matching itself would not be minimum-weight: rematching the corresponding subset of arcs (and possibly uncrossing any crossed arcs that might result from the rematch) would give a matching with a strictly smaller weight.

z ′ k ′ y ′ k ′ • • • • • • z ′ 1 y ′ 1 z k y k • • • • • • z 1 y 1 x r x ℓ

Preservation of hidden arcs

Call an arc (x i , x j ) in a nested matching exposed if it is not contained in any other arc, i.e., if there is no arc (x i ′ , x j ′ ) with x i , x j contained between x i ′ and x j ′ . We call all other arcs in a nested matching non-exposed or hidden. Intuitively, exposed arcs are those visible "from above" and hidden arcs are those covered with exposed ones.

Suppose

X = {x i } 1≤i≤2n with x 1 < x 2 < • • • < x 2n and X ′ = {x ′ i ′ } 1≤i ′ ≤2n ′ with x ′ 1 < x ′ 2 < • • • < x ′
2n ′ be two sets such that x 2n < x ′ 1 , i.e., X ′ lies to the right of X. We will refer to minimum-weight perfect matchings on X and X ′ as partial matchings and to the minimum-weight perfect matching on X ∪ X ′ as joint matching. The following result is closely related to properties of "local matching indicators" introduced and studied in [START_REF] Delon | Local matching indicators for concave transport costs[END_REF][START_REF]Local matching indicators for transport problems with concave costs[END_REF].

Theorem 5. Whenever an arc (x i , x j ) [respectively (x ′ i ′ , x ′ j ′ )]
is hidden in the partial matching on X [respectively on X ′ ], it belongs to the joint optimal matching and is hidden there too.

Observe that exposed arcs in partial matchings are generally not preserved in the joint matching: they may disappear altogether or become hidden.

Proof. By contradiction, assume that some of hidden arcs in the partial matching on X do not belong to the joint matching. Then there will be at least one exposed arc (x ℓ , x r ) in the partial matching on X such that some points x i with x ℓ < x i < x r are connected in the joint matching to points outside (x ℓ , x r ).

Indeed, if points inside every exposed arc (x ℓ , x r ) would be matched in the joint matching only among themselves, then by Lemma 4 their matching would be exactly the same as in the partial matching on X, and therefore all hidden arcs between x ℓ and x r would be preserved in the joint matching.

Suppose (x ℓ , x r ) is the leftmost arc of the above kind. Denote all the points in the segment [x ℓ , x r ] that are connected in the joint matching to points on the left of x ℓ by z 1 < z 2 < • • • < z k ; denote the opposite endpoints of the corresponding arcs by y 1 > y 2 > • • • > y k , where the inequalities follow from the fact that the joint matching is nested. Likewise denote those points from [x ℓ , x r ] that are connected in the joint matching to points on the right of x r by z 1). Although k or k ′ may be zero, the number k + k ′ must be positive and even. Indeed, by Corollary 3 the segment [x ℓ , x r ] contains an even number of points and all of them must be matched in a perfect matching; removing from the joint matching all arcs whose ends both lie in [x ℓ , x r ], we are left with an even number of points that are matched outside this segment.

′ 1 > z ′ 2 > • • • > z ′ k ′ and their counterparts in the joint matching by y ′ 1 < y ′ 2 < • • • < y ′ k ′ (fig.
Let us now restrict our attention to the segment [x ℓ , x r ] and consider a matching that consists of the following arcs: those arcs of the joint matching whose both ends belong to [x ℓ , x r ]; the arcs (z 1 , z 2 ), . . . , (z 2κ-1 , z 2κ ), where

1 κ = ⌊k/2⌋; the arcs (z ′ 2 , z ′ 1 ), . . . , (z ′ 2κ ′ , z ′ 2κ ′ -1 ), where κ ′ = ⌊k ′ /2⌋; and (z k , z ′ k ′ ) if both k and k ′ are odd.
Denote by W ′ the weight of this matching. By Lemma 4, it cannot be smaller than the weight W ′ 0 of the restriction of the partial matching on X to [x ℓ , x r ]. For the total weight W of the joint matching on X ∪ X ′ we thus have

(4) W ≥ W -W ′ + W ′ 0 .
The right-hand side of ( 4) is represented in fig. 2 (a) in the case when both k and k ′ are odd. It is a sum of positive terms corresponding to the arcs of the joint matching outside [x ℓ , x r ] (not shown), the arcs of the partial matching on X inside [x ℓ , x r ] (not shown, with exception of (x ℓ , x r ) represented with a solid arc in the upper halfplane), the arcs of the joint matching having one end inside [x ℓ , x r ] and the other end outside this segment (solid arcs in the upper halfplane), and negative terms that come from subtraction of W ′ and correspond to the arcs connecting the z points (solid arcs in the lower halfplane). We now show that by a suitable sequence of uncrossings the right-hand side of (4) can be further reduced to a matching whose weight is strictly less than W .

Step 1. Note that the arcs (z 1 , y 1 ) and (x ℓ , x r ), shown in fig. 2 (a) with thick lines, are crossing. Therefore

d(y 1 , z 1 ) + d(x ℓ , x r ) > d(y 1 , x ℓ ) + d(z 1 , x r ).
Uncrossing these arcs gives the matching represented in fig. 2 (b) and strictly reduces the right-hand side of (4):

W > W -W ′ + W ′ 0 -d(y 1 , z 1 ) -d(x ℓ , x r ) + d(y 1 , x ℓ ) + d(z 1 , x r ).
Now the arcs (y 2 , z 2 ) and (z 1 , x r ) are crossing, so

d(y 2 , z 2 ) + d(z 1 , x r ) -d(z 1 , z 2 ) > d(y 2 , x r )
and therefore

W > W -W ′ +W ′ 0 -d(y 1 , z 1 )-d(y 2 , z 2 )-d(x ℓ , x r )+d(y 1 , x ℓ )+d(z 1 , z 2 )+d(y 2 , x r ).
The right-hand side of this inequality is represented in fig. 2 (c).

1 ⌊ξ⌋ is the largest integer n such that n ≤ ξ.

z ′ k ′ y ′ k ′ • • • • • • z ′ 1 y ′ 1 z k y k • • • • • • z 1 y 1 x r x ℓ (a) z ′ k ′ y ′ k ′ • • • • • • z ′ 1 y ′ 1 z k y k • • • • • • z 1 y 1 x r x ℓ (b) z ′ k ′ y ′ k ′ • • • • • • z ′ 1 y ′ 1 z k y k • • • • • • y 1 x r x ℓ (c) Figure 2.
Step 1 of the proof (see explanation in the text).

Repeating this step κ = ⌊k/2⌋ times gives the inequality

W > W -W ′ + W ′ 0 -d(x ℓ , x r ) - 1≤i≤2κ d(y i , z i ) + 1≤i≤κ d(z 2i-1 , z 2i ) + 1≤i≤κ d(y 2i-1 , y 2i-2 ) + d(y 2κ , x r ),
where in the rightmost sum y 0 is defined to be x ℓ . Note that at this stage all arcs coming to points z 1 , z 2 , . . . from outside [x ℓ , x r ] are eliminated from the matching, except possibly (y k , z k ) if k is odd.

Step 2. It is now clear by symmetry that a similar reduction step can be performed on arcs going from z ′ 1 , z ′ 2 , . . . to the right. Repeating this

κ ′ = ⌊k ′ /2⌋ z ′ k ′ y ′ k ′ • • • • • • x r y ′ 1 z k y k • • • • • • y 1 x ℓ (a) z ′ k ′ y ′ k ′ • • • • • • x r y ′ 1 z k y k • • • • • • y 1 x ℓ (b) z ′ k ′ y ′ k ′ • • • • • • x r y ′ 1 z k y k • • • • • • y 1 x ℓ (c) Figure 3.
Step 3 of the proof. Note that in stage (c) the arc (z k , z ′ k ′ ) gives two contributions with positive and negative signs, which cancel out each other. times gives the inequality

W > W -W ′ + W ′ 0 -d(x ℓ , x r ) - 1≤i≤2κ d(y i , z i ) - 1≤i ′ ≤2κ ′ d(z ′ i ′ , y ′ i ′ ) + 1≤i≤κ d(z 2i-1 , z 2i ) + 1≤i≤κ d(y 2i-1 , y 2i-2 ) + 1≤i ′ ≤κ ′ d(z ′ 2i ′ , z ′ 2i ′ -1 ) + 1≤i ′ ≤κ ′ d(y ′ 2i ′ -2 , y ′ 2i ′ -1 ) + d(y 2κ , y ′ 2κ ′ ),
where y ′ 0 = x r .

Step 3. If k and k ′ are odd, we perform two more uncrossings shown in fig. 3. The final estimate for W has the form

(5) W > W -W ′ + W ′ 0 -d(x ℓ , x r ) - 1≤i≤k d(y i , z i ) - 1≤i ′ ≤k ′ d(z ′ i ′ , y ′ i ′ ) + 1≤i≤κ d(z 2i-1 , z 2i ) + 1≤i ′ ≤κ ′ d(z ′ 2i ′ , z ′ 2i ′ -1 ) + d(z k , z ′ k ′ ) • [k, k ′ are odd] + 1≤i≤κ d(y 2i-1 , y 2i-2 ) + 1≤i ′ ≤κ ′ d(y ′ 2i ′ -2 , y ′ 2i ′ -1 ) + d(y k , y ′ k ′ ) • [k, k ′ are even],
where notation such as [k, k ′ are odd] means 1 if k, k ′ are odd and 0 otherwise.

The right-hand side of ( 5) contains four groups of terms: first,

W - 1≤i≤k d(y i , z i ) - 1≤i ′ ≤k ′ d(z ′ i ′ , y ′ i ′ ),
corresponding to the joint matching without the arcs connecting points inside [x ℓ , x r ] to points outside this segment; second,

W ′ - 1≤i≤κ d(z 2i-1 , z 2i ) - 1≤i ′ ≤κ ′ d(z ′ 2i ′ , z ′ 2i ′ -1 ) -d(z k , z ′ k ′ ) • [k, k ′ are odd],
which comes with a negative sign and corresponds to the arcs of the joint matching with both ends inside [x ℓ , x r ], and cancels them from the total; third,

W ′ 0 -d(x ℓ , x r
), with positive sign, which corresponds to the hidden arcs of the partial matching on X inside the exposed arc (x ℓ , x r ), not including the latter; and finally the terms in the last line of ( 5), corresponding to the arcs matching x ℓ , x r , and points y 1 , . . . , y k , y ′ 1 , . . . , y ′ k ′ , i.e., those points outside [x ℓ , x r ] that were connected in the joint matching to points inside this segment.

Gathering together contributions of these four groups of terms, we observe that all negative terms cancel out and what is left corresponds to a perfect matching with a weight strictly smaller than W , in which all arcs hidden by (x ℓ , x r ) in the partial matching on X are restored. There may still be some crossings caused by terms of the fourth group and not involving the hidden arcs in [x ℓ , x r ]; uncrossing them if necessary gives a nested perfect matching whose weight is strictly less than that of the joint matching. This contradicts the assumption that the latter is the minimum-weight matching on X ∪ X ′ . Therefore all hidden arcs in the partial matching on X (and, by symmetry, those in the partial matching on X ′ ) belong to the joint matching.

Recursion for minimum weights

In this section we show how to apply Theorem 5 to compute the minimumweight perfect matching algorithmically. Let x 1 < x 2 < • • • < x 2n be a set of points on the real line R equipped with a homogeneous distance d of the form (1). For indices i, j of opposite parity and such that i < j, let W i,j be the weight of the minimum-weight perfect matching on the j -i + 1 points

x i < x i+1 < • • • < x j . It is convenient to organize weights W i,j with i < j into a pyramidal table: W 1,2n . . . . . . . . . . . . . . . . . . . W 1,6 W 2,7 W 3,8 W 4,9 . . . W 1,4 W 2,5 W 3,6 W 4,7 W 5,8 . . . W 2n-3,2n W 1,2 W 2,3 W 3,4 W 4,5 W 5,6 W 6,7 . . . W 2n-2,2n-1 W 2n-1,2n
Theorem 6. For all indices i, j of opposite parity with 1 ≤ i < j ≤ 2n, weights W i,j satisfy the following second-order recursion

(6) W i,j = min d(x i , x j ) + W i+1,j-1 , W i,j-2 + W i+2,j -W i+2,j-2
with "initial conditions"

(7) W i,i-1 = 0, W i+2,i-1 = -d(x i , x i+1 ).
Proof. By an abuse of notation, we will refer to the minimum-weight perfect matching on points x r < x r+1 < • • • < x s as the "matching W r,s ." Consider first the matching that consists of the arc (x i , x j ) and all arcs of the matching W i+1,j-1 , and observe that by Lemma 4 its weight d(x i , x j ) + W i+1,j-1 is minimal among all matchings that contain (x i , x j ).

We now examine the meaning of the expression W i,j-2 + W i+2,j -W i+2,j-2 . Denote the point connected in the matching W i,j-2 to x i by x k and the point connected to x i+1 by x ℓ . By Corollary 3, the pairs of indices i, k and i + 1, ℓ both have opposite parity. Assume first that (8)

x i+1 < x ℓ < x k ≤ x j-2 .
Applying Theorem 5 to the sets X = {x i , x i+1 } and X ′ = {x i+2 , . . . , x j-2 } and taking into account parity of k and ℓ, we see that x k and x ℓ (as well as their neighbors x k+1 and x ℓ-1 if they are contained in [x i+2 , x j-2 ]) belong to exposed arcs of the matching W i+2,j-2 . Thus the matching W i,j-2 has the following structure:

, x i x k x i+1 x ℓ W i+2,ℓ-1 • • • • • • W ℓ+1,k-1 • • • W k+1,j-2
where dashed (resp., dotted) arcs correspond to those exposed arcs of the matching W i+2,j-2 that belong (resp., do not belong) to W i,j-2 . Since points x ℓ-1 and x k+1 belong to exposed arcs in the matching W i+2,j-2 , by Lemma 4 we see that the (possibly empty) parts of this matching that correspond to points x i+2 < • • • < x ℓ-1 and x k+1 < • • • < x j-2 coincide with the (possibly empty) matchings W i+2,ℓ-1 and W k+1,j-2 . For the same reason the (possibly empty) part of the matching W i,j-2 supported on

x ℓ+1 < • • • < x k-1 coincides with W ℓ+1,k-1 . Therefore (9) W i,j-2 = d(x i , x k ) + d(x i+1 , x ℓ ) + W i+2,ℓ-1 + W ℓ+1,k-1 + W k+1,j-2 .
Taking into account [START_REF]Paths, trees and flowers[END_REF], observe that in the case k = i + 1 and ℓ = i, which was left out in [START_REF] Heinonen | Lectures on Analysis in Metric Spaces[END_REF], this expression still gives the correct formula W i,j-2 = d(x i , x i+1 ) + W i+2,j-2 . Now assume that in the matching W i+1,j the point x j is connected to x ℓ ′ and the point x j-1 to x k ′ . A similar argument gives [START_REF] Mccann | Exact solutions to the transportation problem on the line[END_REF]

W i+2,j = W i+2,ℓ ′ -1 + W ℓ ′ +1,k ′ -1 + W k ′ +1,j-2 + d(x ℓ ′ , x j ) + d(x k ′ , x j-1 ); in particular, if ℓ ′ = j -1 and k ′ = j, then W i+2,j = W i+2,j-2 + d(x j-1 , x j ).
Suppose that x k < x ℓ ′ . Using again Lemma 4 and taking into account that x k , x k+1 , x ℓ ′ -1 , and x ℓ ′ all belong to exposed arcs in W i+2,j-2 , we can write

(11) W k+1,j-2 = W k+1,ℓ ′ -1 + W ℓ ′ ,j-2 , W i+2,ℓ ′ -1 = W i+2,k + W k+1,ℓ ′ -1 and (12) W i+2,j-2 = W i+2,k + W k+1,ℓ ′ -1 + W ℓ ′ ,j-2 .
Substituting ( 11) into ( 9) and ( 10) and taking into account (12), we obtain

W i,j-2 + W i+2,j -W i+2,j-2 = d(x i , x k ) + d(x i+1 , x ℓ ) + W i+2,ℓ-1 + W ℓ+1,k-1 + W k+1,ℓ ′ -1 + d(x ℓ ′ , x j ) + W ℓ ′ +1,k ′ -1 + d(x k ′ , x j-1 ) + W k ′ +1,j-2 .
The right-hand side of this expression corresponds to a matching that coincides with W i,j-2 on [x i , x k ], with W i+2,j-2 on [x k+1 , x ℓ ′ -1 ], and with W i+1,j on [x ℓ ′ , x j ]. By Lemma 4 this matching cannot be improved on any of these three segments and is therefore optimal among all matchings in which x i and x j belong to different exposed arcs. It follows that under the assumption that x k < x ℓ ′ the expression in the righthand side of (6) gives the minimum weight of all matchings on x i < x i+1 < • • • < x j . the only possible candidates for the optimal matching are those constructed above: one that corresponds to d(x i , x j ) + W i+1,j-1 and one given by the right-hand side of the latter formula.

It remains to consider the case x k ≥ x ℓ ′ . Since x k = x ℓ ′ for parity reasons, it follows that x k > x ℓ ′ ; now a construction similar to the above yields a matching which corresponds to W i,j-2 +W i+2,j -W i+2,j-2 and in which the arcs (x i , x k ) and (x ℓ ′ , x j ) are crossed. Uncrossing them leads to a matching with strictly smaller weight, which contains the arc (x i , x j ) and therefore cannot be better than d(x i , x j ) + W i+1,j-1 . This means that (6) holds in this case too with W i,j = d(x i , x j ) + W i+1,j-1 .

Obviously, recursion (6) can be solved for all 1 ≤ i < j ≤ 2n in O(n 2 ) operations, resulting in computation of weights W ij of all partial optimal matchings. This process is carried out in a "bottom to top" fashion: in the pyramid, weights W i,j with smaller values of j -i are computed first.

To determine the optimal matching on all the points x 1 , x 2 , . . . , x 2n , one should keep track of those pairs (i, j) for which minimum in ( 7) is attained at the first alternative. Indeed, if the minimum is never attained at the first alternative, then it is easy to see that the optimal perfect matching is (x 1 , x 2 ), (x 3 , x 4 ), . . . , (x 2n-1 , x 2n ). Suppose now that the first alternative provides minimum for for some W i0j0 . Then according to Theorem 5 one can retain the matching for the points x i0+1 , . . . , x j0-1 that has been computed by this moment and consider a new, smaller minimumweight perfect matching problem on the points x 1 , x 2 , . . . , x i0 , x j0 , x j0+1 , . . . , x 2n . More precisely, it suffices to remove from the pyramidal table quantities W ij with i, j satisfying at least one of the conditions i 0 < i < j 0 or i 0 < j < j 0 , replace W i0j0 with d(x i0 , x j0 ), stack the cells with either i or j outside (i 0 , j 0 ) into a smaller pyramidal table, and continue solving the recursion. This reduction step is illustrated in fig. 4. The elements to be removed from the table are shown in light gray in the top pane. Assuming that the rows are scanned left to right, at the moment when it is found that the element W i0,j0 = W 3,6 (in brackets) involves the first alternative in [START_REF] Edmonds | Maximum matching and a polyhedron with 0, 1-vertices[END_REF], the elements shown in parentheses have not yet been computed. Those of them that are to be kept in the table (the "black" ones) are then stuck with the already computed elements to form a smaller pyramid shown in the bottom, and the recursion resumes with the element W 3,6 replaced with d(x 3 , x 6 ).

At the reduction step illustrated in fig. 4, the matching is updated with one arc (x 4 , x 5 ), which is guaranteed to belong to the optimal matching by Theorem 5. Generally, every time a reduction step is performed on an element W i0,j0 and results in removal of the points x k1 < x k2 < • • • < x k2m that have been "covered" with (W the arc (i 0 , j 0 ), the matching is updated with arcs (x k1 , x k2 ), . . . , (x k2m-1 , x k2m ). When the process reaches the top element of the pyramid, a similar matching of consecutive pairs is performed on the remaining elements.

A formal description of this procedure is given in [START_REF] Delon | Local matching indicators for concave transport costs[END_REF][START_REF]Local matching indicators for transport problems with concave costs[END_REF] in somewhat different terms (emphasizing "local matching indicators" instead of the quantities W ij ). Since this algorithm requires computing all the weights W ij in the worst case, it takes at most O(n 2 ) operations for the matching on 2n points.

It is interesting to compare performance of our algorithm to the celebrated Edmonds algorithm for minimum-weight perfect matching in a general graph with N vertices and M edges [START_REF] Edmonds | Maximum matching and a polyhedron with 0, 1-vertices[END_REF][START_REF]Paths, trees and flowers[END_REF]. Note that we effectively consider a complete graph on N = 2n vertices, i.e., M = n(2n -1). In this case recent improvements of the Edmonds algorithm run in O(n 3 ) time (see, e.g., the survey [START_REF] Cook | Computing minimum-weight perfect matchings[END_REF]) whereas our specialized algorithm is at most quadratic in n.

Figure 1 .

 1 Figure 1. Notation used in the proof of Theorem 5. Note that in general x ℓ ≤ z 1 and z ′ 1 ≤ x r , but in the figure these pairs of points are shown to be distinct.Lemma 4. Any subset of arcs in a minimum-weight matching is itself a minimumweight perfect matching on the set of endpoints of the arcs that belong to this subset.

Figure 4 .

 4 Figure 4. Reduction of the pyramidal table (i 0 = 3, j 0 = 6, n = 5).

  LTCI CNRS, TELECOM ParisTech, 46 rue Barrault, 75634 Paris, France CEREMADE, Université Paris-Dauphine, Place du Marechal Lattre de Tassigny, 75775 Paris, France Institute for Information Transmission Problems (Kharkevich Institute), 19 B. Karetny per., Moscow 127994, Russia Laboratoire J.-V. Poncelet (UMI 2615 CNRS), 11 B. Vlassievski per., Moscow 119002, Russia E-mail address: sobolevski@iitp.ru

  1,10 ) (W 1,8 ) (W 2,9 ) (W 3,10 ) (W 1,6 ) (W 2,7 ) (W 3,8 ) (W 4,9 ) (W 5,10 ) W 1,4 W 2,5 [W 3,6 ] (W 4,7 ) (W 5,8 ) (W 6,9 ) (W 7,10 ) W 1,2 W 2,3 W 3,4 ) (W 3,10 ) (W 1,6 ) (W 2,7 ) (W 3,8 ) (W 6,9 ) (W 7,10 ) W 1,2 W 2,3 d(x 3 , x 6 )

	W 4,5	W 5,6	W 6,7	W 7,8	W 8,9	W 9,10
		↓
		(W 1,10 )
	(W 1,8 )	(W 2,9 W 6,7	W 7,8	W 8,9	W 9,10