Mingyuan Jiu

Nelly Pustelnik
email: nelly.pustelnik@ens-lyon.fr.

Stefan Janaqi

Mèriam Chèbre

Lin Qi

Philippe Ricoux

Mériam Chebre

Sparse hierarchical interaction learning with epigraphical projection

come

Introduction

Learning interactions between features is of great interest in data processing. In this work we focus on quadratic interaction effects between the features that extend linear models. Our work focuses both on the regression problem and multiclass SVM, which are still subjects of active research when small dataset are involved (especially in medical [START_REF] Spilka | Sparse Support Vector Machine for Intrapartum Fetal Heart Rate Classification[END_REF] or physics [START_REF] Pascal | Joint Estimation Of Local Variance And Local Regularity For Texture Segmentation. Application To Multiphase Flow Characterization[END_REF] applications).

Feature selection -Our work follows the line of feature selection methods [START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF][START_REF] Blondel | Block coordinate descent algorithms for large-scale sparse multiclass classification[END_REF][START_REF] Chakraborty | Feature selection using a neural framework with controlled redundancy[END_REF][START_REF] Chierchia | A proximal approach for sparse multiclass SVM[END_REF][START_REF] Laporte | Nonconvex regularizations for feature selection in ranking with sparse SVM[END_REF][START_REF] Xu | Semisupervised feature selection based on relevance and redundancy criteria[END_REF], which aim to remove the redundant features and only keep the informative ones. The relevant features are usually correlated and belong to the same group. In this work, additionally to perform feature selection, we learn the interactions between the selected features by making some hierarchy constraints, where the interactions may occur either between the selected features or if only one of both features is active.

Knowledge from the interactions -Features interactions can provide us a new knowledge about the correlation between the subjects. For instance, we can refer to gene interactions aiming to highlight relevant regulatory relationships for genes [START_REF] Pirayre | BRANE Cut: Biologicallyrelated apriori network enhancement with graph cuts for gene regulatory network inference[END_REF] or as mentioned in [START_REF] Bien | A lasso for hierarchical interactions[END_REF] "the co-occurrence of two symptoms may lead a doctor to be confident that a patient has a certain disease, whereas the presence of either symptom without the other would provide only a moderate indication of that disease".

Better discrimination -Dealing with interactions allows us to increase the feature space in order to provide a better discrimination in the learning process [START_REF] Anderson | Quadratic logistic discrimination[END_REF]. The integration of quadratic interactions in the learning problem is not new, for example, discriminant quadratic learning [START_REF] Anderson | Quadratic logistic discrimination[END_REF][START_REF] Ronald | Generalized linear and quadratic discriminant functions using robust estimates[END_REF][START_REF] Witten | Covariance-regularized regression and classification for high dimensional problems[END_REF] learns the covariance matrix in the quadratic term to improve the discrimination ability. We can also refer to efficient multilevel procedures [START_REF] Agarwal | Scalable nonlinear learning with adaptive polynomial expansions[END_REF][START_REF] Hao | Model selection for high-dimensional quadratic regression via regularization[END_REF] starting with fast learning algorithm focus on the linear model and then adding higher-order interaction features, possibly using the learned weights as a guide.

Dimensionality challenge and sparsity -One major challenge when dealing with quadratic interaction learning is that the interaction number quadratically increases with the feature size, for example, 1000 features will have about one million possible interactions, this therefore results in an overfitting problem due to insufficient observation samples in most of the regression problems and for some classification ones. To overcome this weakness, the linear weights and the quadratic interactions can be assumed to be sparse because most features would not contribute to the decision. Sparsity-based regularization is known to be mainly suitable when the feature size is larger than the training samples. Various sparsity regularizations have been proposed and extensively studied in the additive model context, for instance, involving 0 -pseudo-norm [START_REF] Weston | Use of the zero-norm with linear models and kernel methods[END_REF], 1 -norm [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], ∞ -norm [START_REF] Zou | The F∞-norm support vector machine[END_REF], or structural sparsity norm [START_REF] Bach | Structured sparsity through convex optimization[END_REF]. See [START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF] for an exhaustive list of sparse-based regularization in learning and [START_REF] Gui | Feature selection based on structured sparsity : a comprehensive study[END_REF] to refer to structural sparsity. The formalism we consider in this work follows ideas proposed in [START_REF] Bien | A lasso for hierarchical interactions[END_REF][START_REF] Haris | Convex modeling of interactions with strong heredity[END_REF][START_REF] Lim | Learning interactions through hierarchical group-lasso regularization[END_REF][START_REF] She | Group regularized estimation under structural hierarchy[END_REF][START_REF] Zhao | The composite absolute penalties family for groupes and hierarchical variable selection[END_REF]. A detailed comparison with these state-of-the-art methods is provided in next section.

Contributions -Quadratic interaction learning imposes a specific design of the penalization term (this will be formally described in Section 2). Several choices of penalizations have already been proposed in the literature, but each of them relies on a different algorithmic strategy, leading to unfair numerical comparisons. Indeed, do the numerical differences in the learning performance come from the penalization choice or from the algorithmic schemes ? For instance, it is well known that the number of inner iterations (as proposed in [START_REF] Jenatton | Proximal methods for hierarchical sparse coding[END_REF]) is often a parameter tricky to adjust, or that the algorithmic parameters such as the step-size may impact a lot the convergence speed [START_REF] Haris | Convex modeling of interactions with strong heredity[END_REF]. For all these reasons, this work is dedicated to provide a common algorithmic scheme without inner iterations aiming to handle with the most recents state-of-the-art penalizations. Our detailed contributions are listed below:

• Propose a new algorithmic scheme relying of primal-dual algorithm and derive new closed form of epigraphical projections in order to solve the unifying minimization problem (cf. Problem 2.1) proposed in [START_REF] She | Group regularized estimation under structural hierarchy[END_REF] when q = +∞ and z(•) = • r with r = {1, +∞};

• Compare the proposed algorithms to the state-of-the-art methods (FISTA and ADMM), showing a unique and efficient algorithmic scheme for weak and strong formulation and for r = {1, +∞} allowing us to avoid inner iterations.

• Provide the counterpart of the regression minimization Problem 2.1 for classification task.

• Conduct extensive experiments in the regression and classification framework on the simulated data, on two disease applications (i.e. HIV and Parkinson) to analyze the correlations between the factors for the diseases, validating the effectiveness and efficiency of the proposed algorithms.

This work improves over our preliminary contribution [START_REF] Jiu | Multiclass SVM with hierarchical interaction: application to face classification[END_REF] where we restricted our study to classification (cf. section 5 of this work). The algorithm was derived for r = 1, relying only on Proposition 4.3 and for which the proof was not provided. The experiments were conducted on face classification while in this work we focus on synthetic data and diseases applications.

Outline -This paper is organized as follows: Section 2 firstly describes the formal problem in the context of regression and refers to related works. Section 3 derives an epigraphical writing of the objective function. Section 4 presents the primal-dual proximal algorithm iterations and the convergence guarantees. Its specification to our minimization problem with the hierarchical regularization term is specified and new epigraphical projections are provided. Section 5 provides the objective function and the associated algorithm in the context of multiclass SVM with hierarchical interactions. Section 6 evaluates the performance of the proposed strategy compared to FISTA and ADMM formulations proposed in [START_REF] Bien | A lasso for hierarchical interactions[END_REF][START_REF] Jenatton | Proximal methods for hierarchical sparse coding[END_REF] both on synthetic data and real applications. Finally Section 7 gives our conclusion.

Regression model

The regression training set is denoted

R = { y , φ(x ) ∈ R × R N | ∈ {1, .
. . , L}}, where x . It can denote either data without structure or with structure such as a signal, an image or a graph of size M . Then φ(x ) denotes the coefficients associated with the data x , e.g time-frequency coefficients [START_REF] Flandrin | Time-frequency/time-scale analysis[END_REF], scattering coefficients [START_REF] Bruna | Invariant scattering convolution networks[END_REF] or CNN coefficients [START_REF] Lecun | Convolutional networks for images, speech, and time series. The handbook of brain theory and neural networks[END_REF], where generally N M .

Several minimization strategies have been provided in the literature to jointly select discriminating features among the features φ(x ) and identify meaningful interactions between these features. To clarify the state-of-the-art contributions in this context, we propose to recall a general minimization problem derived in [START_REF] She | Group regularized estimation under structural hierarchy[END_REF]:

Problem 2.1 Let R = { y , φ(x ) ∈ R × R N | ∈ {1, .
. . , L}} be the training regression dataset. We aim to estimate the regression weight vector v and the interaction matrix Θ solving:

( v, Θ) ∈ arg min v∈R N Θ∈R N ×N 1 2 L =1 y -φ(x ) v -φ(x ) Θφ(x ) 2 + Ω(v, Θ) (1)
where, for every v = v (i) 1≤i≤N and Θ = Θ (i,j) 1≤i≤N,1≤j≤N ,

Ω(v, Θ) = λ 2 Θ 1 + λ N i=1 v (i) , z Θ (i,•) q + ι C (Θ) (2)
with z : R N → R, 1 ≤ q ≤ +∞, and ι C denotes the indicator function

1 of a closed convex set C ⊂ R N ×N .
Figure 1: Hierarchical feature interaction. The first row is N additive features, and the second row is the interactions between pairs of features.

Choice for the constraint C: Weak/strong hierarchy

The hierarchy penalization Ω creates the relationship between v and Θ, as shown in Fig. 1. Two types of hierarchy constraints have been defined in "hierNet" [START_REF] Bien | A lasso for hierarchical interactions[END_REF]: (i) weak hierarchy where the interactions happen (i.e., Θ (i,j) = 0) only if one of the associated features in v is non-zero (i.e., v (i) = 0 or v (j) = 0) and (ii) strong hierarchy when both associated features are non-zero. These two types of hierarchy are imposed by means of the closed convex set C ⊂ R N ×N . The weak hierarchy is obtained with C = R N ×N and strong hierarchy by imposing a symmetric structure for the matrix of interactions using

C = S = {Θ ∈ R N ×N | Θ = Θ }.

Positioning of Problem 2.1 w.r.t state-of-the-art

Most of the state-of-the-art penalizations may be interpreted as a specific case of Problem 2.1, for instance:

• [24] : z(•) = (•) , q = {2, +∞},
• [START_REF] Bien | A lasso for hierarchical interactions[END_REF][START_REF] Haris | Convex modeling of interactions with strong heredity[END_REF] : z(•) = • 1 and q = +∞,

• [START_REF] She | Group regularized estimation under structural hierarchy[END_REF] : z(•) = (•) and q = {2}.

Note that [START_REF] Jenatton | Proximal methods for hierarchical sparse coding[END_REF] focuses on weak hierarchy using C = R N ×N . [START_REF] Bien | A lasso for hierarchical interactions[END_REF] is encountered under the name "hierNet", while [START_REF] Haris | Convex modeling of interactions with strong heredity[END_REF] is named "FAMILY" and [START_REF] She | Group regularized estimation under structural hierarchy[END_REF] is "Type-A GRESH". The latent overlapping group lasso formulation of [START_REF] Bien | A lasso for hierarchical interactions[END_REF] is known as "glinternet" [START_REF] Lim | Learning interactions through hierarchical group-lasso regularization[END_REF].

From the algorithmic point of view several strategies can be encountered. In [START_REF] Jenatton | Proximal methods for hierarchical sparse coding[END_REF] and [START_REF] She | Group regularized estimation under structural hierarchy[END_REF], the iterations are derived from FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. This procedure appears to be very efficient when C = R N ×N while for C = S it requires inner iterations based on Dykstra's algorithm that significantly slows the convergence [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problem[END_REF]. On the other hand, [START_REF] Bien | A lasso for hierarchical interactions[END_REF] and [START_REF] Haris | Convex modeling of interactions with strong heredity[END_REF] resort to ADMM to deal with C = S.

Using recent developments of non-smooth convex optimization, the objective of this work is to provide a unified framework and an efficient algorithmic strategy based on primal-dual proximal algorithms and epigraphical projection to solve an epigraphical relaxation of (4) when C = R N ×N or C = S, z(•) = • r with r = {1, +∞} and q = +∞. The algorithmic procedure does not resort to inner iterations, thus leading to a faster implementation. 2

Epigraphical formulation

In order to handle with complex optimization problems, the two major solutions encountered in the literature are either to formulate it in the dual or to increase the dimensionality of the problem. Dealing with an epigraphical formulation belongs to this second class and it is mainly adapted to nonlinear constraints. We recall that the epigraph of a function

f : H →] -∞, +∞] is defined as: epif = {(x, ζ) ∈ H × R | f (x) ≤ ζ},
the projection of a point to this epigraph is called epigraphical projection. The utility of such epigraphical formulation can be illustrated when one wants to deal with a constraint of the form x 1 ≤ η, for every x = (x i ) i∈I⊂N and η > 0 is known, whose projection does not have a closed form expression. An alternative is then to replace this constraint with these two constraints: |x i | ≤ ζ i and i∈I ζ i = η where the first one denotes an epigraphical constraint and the second one denotes an hyperplane constraint, both having a closed form expression for the associated projection [START_REF] Chierchia | Epigraphical projection and proximal tools for solving constrained convex optimization problems[END_REF].

The resolution of Problem 2.1 is difficult, since it involves non-smooth convex functions and symmetry constraints. In [START_REF] Jenatton | Proximal methods for hierarchical sparse coding[END_REF], Jenatton et al. solve the problem for r = ∞ in the weak hierarchy formulation by means of a proximal algorithm. However, their algorithm is not well adapted to strong hierarchy. In [START_REF] Bien | A lasso for hierarchical interactions[END_REF], the authors reformulate Problem 2.1 when r = 1 under an epigraphical formulation which represents the problem as a set of hard constraints. This reformulation helps in the interpretation and it highlights a reduction in the shrinkage of certain main effects and an increase in the shrinkage of certain interactions. Its second advantage is to help in the design of the ADMM algorithmic scheme. The following proposition gives the epigraphical formulation of the minimization problem considered in this work.

Proposition 3.1 The minimization problem minimize v,Θ 1 2 L =1 y -φ (x )v -φ (x )Θφ(x ) 2 + Ω(v, Θ) (3)
where, for every v = v (i) 1≤i≤N and Θ = Θ (i,j) 1≤i≤N,1≤j≤N ,

Ω(v, Θ) = λ 2 Θ 1 + λ N i=1 max{v (i) , Θ (i,•) r } + ι C (Θ) (4)
can be written in an epigraphical framework as minimize

(v + ,v -,Θ)∈O×O×C 1 2 L =1 y -φ (x )(v + -v -) -φ (x )Θφ(x ) 2 + λ 2 Θ 1 + λ1 (v + + v -) + N i=1 ι Er v +(i) , v -(i) , Θ (i,•) (5)
where O denotes the positive orthant

, v = v + -v -and E r = {(ω + , ω -, u) ∈ R × R × R N | u r ≤ ω + + ω -} with r = {1, +∞}. Proof 3.1 The proof is given in Appendix A.
The algorithmic strategy designed in next section will focus on this epigraphical formulation (5).

Primal-dual proximal algorithm

Proximal algorithms are derived from two main frameworks that are the forward-backward scheme and the Douglas-Rachford iterations (both being deduced from the Krasnoselskii-Mann scheme) [START_REF] Bauschke | Convex analysis and monotone operator theory in hilbert spaces[END_REF]. FISTA can be presented as an accelerated version of forward-backward iterations [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Chambolle | On the convergence of the iterates of the "Fast Iterative Shrinkage/Thresholding Algorithm[END_REF] while ADMM can be viewed as a Douglas-Rachford procedure in the dual [START_REF] Setzer | Split Bregman algorithm, Douglas-Rachford splitting and frame shrinkage[END_REF]. In the literature dedicated to sparse regression the most encountered algorithm is FISTA when dealing with a sum of two convex functions where one is differentiable with a Lipschitz gradient. When the criterion involves more than two functions, typically an additional constraint such as the constraint S defined previously, is added and most of the works derive an ADMM procedure or compute the proximity operator by means of inner iterations which is often known to leading to an approximate solution even if global convergence can be obtained in specific cases [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problem[END_REF]. Another class of algorithmic procedures allowing to minimize a criterion with more than two functions, possibly including a differentiable function with a Lipschitz gradient, is the class of primal-dual proximal approaches [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF][START_REF] Komodakis | Playing with duality: an overview of recent primaldual approaches for solving large-scale optimization problems[END_REF][START_REF] Vũ | A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF].

Several primal-dual proximal schemes have been derived but one of the most popular is based on forward-backward iterations [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF][START_REF] Vũ | A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF] in order to estimate

w ∈ Argmin w∈H f (w) + g(w) + h(Hw), (6)
where H : H → G denotes a bounded linear operator, H and G being two Hilbert spaces, f : H →] -∞, +∞], g : H →] -∞, +∞] and h : G →] -∞, +∞] denote convex, l.s.c and proper functions. We additionally assume that f is differentiable and ∇f has a Lipschitz constant denoted as β > 0. Iterations are summarized in Algorithm 1 involving the proximity operator defined as

(∀w ∈ H) prox g (w) = arg min u∈H 1 2 u -w 2 2 + g(u)
and h * denotes the Fenchel-Rockafellar conjugate of function h. The proximity operator of the conjugate can be computed according to the Moreau identity that is prox σh * (w) = w -σprox h/σ (w/σ) for σ > 0. The sequence (w [k+1] ) k∈N converges to a minimizer w of ( 6). Moreover, the sequence (u [k+1] ) k∈N converge to a minimizer of the dual formulation of ( 6)

that is u ∈ Argmin u∈G (f + g) * (-H * u) + h * (u).
This algorithmic scheme can be extended for dealing with more than three functions by setting h(Hx)

= K k=1 h k (H k x) involving the computation of prox σh * k .
The interest of this scheme compared to ADMM is twofold: it first makes the possibility to deal with the gradient of the differentiable function and secondly it avoids to invert k H * k H k .

Algorithm 1 Primal-dual splitting algorithm

Parameter settings:

Set τ > 0, σ > 0, such that 1 τ -σ H 2 ≥ β 2 Initialization: (w [0] , u [0] ) ∈ H × G For k = 0, 1, . . . w [k+1] = prox τ g w [k] -τ ∇f (w [k] ) -τ H * u [k] u [k+1] = prox σh * u [k] + σH(2w [k+1] -w [k] )
4.1 Specificity for minimizing [START_REF] Bauschke | Convex analysis and monotone operator theory in hilbert spaces[END_REF] We first provide the iterations of Algorithm 1 specified to the minimization of [START_REF] Bauschke | Convex analysis and monotone operator theory in hilbert spaces[END_REF]. By

setting w = (v + , v -, Θ) ∈ R N × R N × R N ×N
, for weak hierarchy, we can split it as follows:

     f (w) = 1 2 L =1 y -φ (x )(v + -v -) -φ (x )Θφ(x ) 2 + λ1 (v + + v -), g(w) = ι O (v + ) + ι O (v -) + λ 2 Θ 1 , h(w) = N i=1 ι Er v +(i) , v -(i) , Θ (i,•) , (7)
and for strong hierarchy:

         f (w) = 1 2 L =1 y -φ (x )(v + -v -) -φ (x )Θφ(x ) 2 + λ1 (v + + v -), g(w) = ι O (v + ) + ι O (v -) + ι S (Θ), h 1 (w) = λ 2 Θ 1 , h 2 (w) = N i=1 ι Er v +(i) , v -(i) , Θ (i,•) . (8)
The respective iterations are summarized in Algorithm 2 and 4. The projection onto the positive orthant and on S have well known closed form expressions [START_REF] Bauschke | Convex analysis and monotone operator theory in hilbert spaces[END_REF] that are:

P O (•) = max{0, •}, P S (Θ) = Θ+Θ 2 .
The difficulty comes from the computation of P Er whose expressions are provided in the next section.

Algorithm 2 -Weak-PDr -Primal-dual splitting algorithm to solve [START_REF] Bauschke | Convex analysis and monotone operator theory in hilbert spaces[END_REF] when

C = R N ×N Parameter settings: Set β = 2 L =1 |φ(x )| 2 + L =1 |φ(x ) φ(x )| 2 , τ > 0, σ > 0, such that 1 τ -σ ≥ β 2 . Initialization: (v +[0] , v -[0] , Θ [0] , s +[0] , s -[0] , Λ [0] ) ∈ R N × R N × R N ×N × R N × R N × R N ×N . For k = 0, 1, . . .                 b = y -φ (x )(v +[k] -v -[k] ) -φ (x )Θ [k] φ(x ) v +[k+1] = P O v +[k] + τ φ(x )b -τ λ -τ s +[k] v -[k+1] = P O v -[k] -τ φ(x )b -τ λ -τ s -[k] Θ [k+1] = prox τ λ 2 • 1 Θ [k] + τ φ (i) (x )φ (j) (x )b 1≤i≤N,1≤j≤N -τ Λ [k] For i = 1, . . . , N (s +(i)[k+1] , s -(i)[k+1] , Λ (i,•)[k+1] ) = prox σι * Er s +(i)[k] + σ(2v +[k+1] -v +[k] ), s -(i)[k] + σ(2v -(i)[k+1] -v -(i)[k] ), Λ +(i,•)[k] + σ(2Θ (i,•)[k+1] -Θ (i,•)[k] )

New epigraphical projection

We first focus on the derivation of P E∞ .

Algorithm 3 -Strong-PDr -Primal-dual splitting algorithm to solve [START_REF] Bauschke | Convex analysis and monotone operator theory in hilbert spaces[END_REF]

when C = S Parameter settings: Set β = 2 L =1 |φ(x )| 2 + L =1 |φ(x ) φ(x )| 2 , τ > 0, σ > 0, such that 1 τ -2σ ≥ β 2 . Initialization: (v +[0] , v -[0] , Θ [0] , s +[0] , s -[0] , Λ [0] 1 , Λ [0] 2 ) ∈ R N × R N × R N ×N × R N × R N × R N ×N × R N ×N . For k = 0, 1, . . .                    b = y -φ (x )(v +[k] -v -[k] ) -φ (x )Θ [k] φ(x ) v +[k+1] = P O v +[k] + τ φ(x )b -τ λ -τ s +[k] v -[k+1] = P O v -[k] -τ φ(x )b -τ λ -τ s -[k] Θ [k+1] = P S Θ [k] + τ φ (i) (x )φ (j) (x )b 1≤i≤N,1≤j≤N -τ (Λ [k] 1 + Λ [k] 2 ) Λ +[k+1] 1 = prox σλ 2 • * 1 Λ +[k] 1 + σ(2Θ [k+1] -Θ [k] ) For i = 1, . . . , N (s +(i)[k+1] , s -(i)[k+1] , Λ (i,•)[k+1] 2 ) = prox σι * Er s +(i)[k] + σ(2v +[k+1] -v +[k] ), s -(i)[k] + σ(2v -(i)[k+1] -v -(i)[k] ), Λ +(i,•)[k] 2 + σ(2Θ (i,•)[k+1] -Θ (i,•)[k] ) Proposition 4.1 Let u = (u (i) ) 1≤i≤N . The projection of (ω + , ω -, u) ∈ R × R × R N on the epigraphic set E ∞ reads (η + , η -, p) = P E∞ (ω + , ω -, u) with      η -= ω --(N -n+1)(ω + -ω -)+ N i=n ν (i) 1+2(N -n+1) η + = η -+ ω + -ω - p = P [-η + -η -,η + +η -] (u) (9)
where (ν (1) , . . . , ν (N ) ) is an ordered version of (|u (i) |) 1≤i≤N in an ascending order and with ν (0) = -∞, and n ∈ {1, . . . , N } such that ν (n-1) < η

+ + η -≤ ν (n) . Proof 4.1 The proof is given in Appendix B.
The above result is obtained from arguments close to the ones derived in [START_REF] Chierchia | Epigraphical projection and proximal tools for solving constrained convex optimization problems[END_REF] [Proposition 5] for an epigraphical constraint of the form (ω, u) ∈ R×R N | max{τ 1 |u (1) |, . . . , τ N |u (N ) |} ≤ ω where (τ i ) 1≤i≤N denotes positive weights. The next proposition is a preliminary result to derive

P E 1 . Proposition 4.2 Let E + 1 = {(ω + , ω -, u) ∈ R × R × R N | u 1 = ω + + ω -, u ≥ 0}. The projection of (ω + , ω -, u) ∈ R × R × R N on the epigraphic set E + 1 reads (η + , η -, p) = P E + 1 (ω + , ω -, u) with        η -= n i=1 µ (i) -ω + +( n+1)ω - n(1+2/ n) η + = η -+ ω + -ω - p = u -max 0, n i=1 µ (i) -(η + +η -) n (10)
with

n = max{n ∈ {1, . . . , N } |, µ (n) - n i=1 µ (i) -(η + + η -) n > 0} (11)
and where µ = (µ (i) ) 1≤i≤N is an ordered version of u = (u (i) ) 1≤i≤N in a descending order.

Proof 4.2 The proof is given in Appendix C.

Next we get the solution of P E 1 from the ones of projection to E + 1 according to [START_REF] Duchi | Efficient projections onto the l1-ball for learning in high dimensions[END_REF][Lemma 3] and have the following proposition:

Proposition 4.3 The projection of (ω + , ω -, u) ∈ R × R × R N on the epigraphic set E 1 reads (η + , η -, p) = P E 1 (ω + , ω -, u) with (v + , v -, p) = P E + 1 (ω + , ω -, abs(u)) p = sign(u)p ( 12
)
where abs(•) and sign(•) denote the componentwise absolute value and the componentwise sign.

Integrating the projection derived in Proposition 4.

(ω + , ω -, u) = (ω + , ω -, u) -σP Er ω + σ , ω - σ , u σ
ensures the convergence to a minimizer of ( 5) respectively for weak hierarchy and strong hierarchy.

Convergence

The proposed Algorithms 2 and 4 can be viewed as a particular case of the primal-dual algorithm [18, Algorithm 5.1] and thus the convergence guarantees can be derived from [18, Theorem 5.1]).

Theorem 4. [START_REF] Bach | Structured sparsity through convex optimization[END_REF] The sequences (v

+[k+1] , v -[k+1] , Θ [k+1]
) k∈N ) i∈N generated by Algorithm 2 converges to a solution of (5) when C = R N ×N (weak hierarchy).

Theorem 4.5 The sequences (v

+[k+1] , v -[k+1] , Θ [k+1]
) k∈N ) i∈N generated by Algorithm 4 converges to a solution of (5) when C = S (strong hierarchy).

Remark 4.1 In practice, the choice of τ and σ is made as follows: τ = 1.9 β and σ = (1/τ -β/2) for Algorithm 2 (resp. σ = 1/2 * (1/τ -β/2) for Algorithm 4).

Extension to multiclass SVM

Model

Formally, assuming that the training set is denoted

T = {(y , φ(x )) ∈ {1, . . . , K} × R N | ∈ {1, . . . , L}}
with K classes and L training samples. x ∈ R M denotes the data (e.g. an image with M pixels) with the label y . The transform φ : R M → R N maps from the input space onto an arbitrary feature space, for instance, scattering features [START_REF] Bruna | Invariant scattering convolution networks[END_REF].

Our aim being to study quadratic interactions, for the k-th class, we consider a discrimination function taking the form:

f k (x ) = φ(x ) Θ k φ(x ) + v k φ(x ) (13)
where,

v k = v (i) k 1≤i≤N and Θ k = Θ (i,j) k
1≤i≤N,1≤j≤N model respectively the weights of main features and the matrix of interactions for the class k ∈ {1, . . . , K}.

Design of the objective function

A large panel of data term can be encountered in the multiclass SVM literature going from hinge loss to logistic data-term [START_REF] Blondel | Block coordinate descent algorithms for large-scale sparse multiclass classification[END_REF][START_REF] Chierchia | A proximal approach for sparse multiclass SVM[END_REF]. In this work, we focus on the squared hinge loss data-term proposed in [START_REF] Blondel | Block coordinate descent algorithms for large-scale sparse multiclass classification[END_REF] leading to good classification performance and being differentiable with Lipschitz gradient β > 0, whose constant will be specified in Section 5.4. The proposed objective function is written as:

minimize v k ∈R N Θ k ∈R N ×N L =1 k =y max 0, 1 -f y (x ) -f k (x ) 2 + K k=1 Ω(v k , Θ k ) ( 14
)
where the first term denotes the squared hinge loss data-term and the second term penalizes the behavior of the discrimination function. Instead of 1 -norm, making v k and Θ k independent, we can force the hierarchy constraints between v k and Θ k for k-th class to regulate them. In the classification setting the penalization term is written:

Ω(v k , Θ k ) = λ 1 N i=1 max{|v (i) k |, Θ (i,•) k r } + λ 2 Θ k 1 + ι C (Θ k ) ( 15
)
where r = {1, +∞}, Θ

(i,•) k = Θ (i,j) k
1≤j≤N and ι C is defined as in Section 2.

Epigraphical reformulation

The reformulation of problem ( 14) by means of epigraphical constraints3 follows ideas derived in [START_REF] Bien | A lasso for hierarchical interactions[END_REF] (i.e. "hierNet") in the context of regression and for a different penalization Ω.

Proposition 5.1 The minimization problem [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problem[END_REF] with Ω(v k , Θ k ) defined in (15) can be reformulated as

minimize (v + k ,v - k ,Θ k )∈O×O×C L =1 k =y max 0, 1 -f y (x ) -f k (x ) 2 + K k=1 λ 1 1 (v + k + v - k ) + λ 2 N i=1 Θ (i,•) k 1 + N i=1 ι Er v +(i) k , v -(i) k , Θ (i,•) k + ι C (Θ k ) ( 16
)
where O denotes the positive orthant and

E r = {(ω + , ω -, u) ∈ R×R×R N | u r ≤ ω + +ω -} with r = {1, +∞} and f k (x ) = φ(x ) Θ k φ(x ) + (v + k -v - k ) φ(x ).
The variable v k is written as

v k = v + k -v -
k by variable splitting. The proof relies on similar arguments than for Proposition 3.1.

Primal-dual proximal algorithm based on epigraphical projection

The proposed objective function ( 5) is a specific case of ( 6)

when w = (v + k , v - k , Θ k ) ∈ R N × R N × R N ×N , Q = 2 and                f (w) = L =1 k =y max 0, 1 -f y (x ) -f k (x ) 2 + K k=1 λ 1 1 (v + k + v - k ) g(w) = K k=1 ι O (v + k ) + ι O (v - k ) + ι S (Θ k ) h 1 (w) = K k=1 λ 2 Θ k 1 h 2 (w) = K k=1 N i=1 ι Er v +(i) k , v -(i) k , Θ (i,•) k .
where f is differentiable with a Lipschitz constant β = 4(K -1)( φ(x ) 2 + φ(x )φ(x ) 2 ). The primal-dual proximal iterations are displayed in Algorithm 4. It involves four proximity operator computation: the projection onto O and S, the proximal operator of 1 -norm and the epigraphical projection onto E r , whose closed form expression are derived in Section 4.

Algorithm 4 Primal-dual proximal algorithm based on epigraphical projection for classification when C = S

Parameter settings: Set τ > 0, σ > 0, such that 2(1/τ -σ) ≥ β. Initialization: (v

+[0] k , v -[0] k , Θ [0] k , s +[0] k , s -[0] k , Λ [0] k,1 , Λ [0] k,2 ) ∈ R N × R N × R N ×N × R N × R N × R N ×N × R N ×N ∀k = {1, . . . , K}. For t = 0, 1, . . .                          f (v + , v -, Θ) = L =1 k =y max 0, 1 -φ(x ) Θ y φ(x ) + (v + y -v - y ) φ(x )- φ(x ) Θ k φ(x ) -(v + k -v - k ) φ(x ) 2 + K k=1 λ 1 1 (v + k + v - k ); v +[t+1] k = P O v +[t] k -τ ∇ v + k f (v +[t] k , v -[t] k , Θ [t] k ) -τ s +[t] k ∀k = 1, . . . , K v -[t+1] k = P O v -[t] k -τ ∇ v - k f (v +[t] k , v -[t] k , Θ [t] k ) -τ s -[t] k ∀k = 1, . . . , K Θ [t+1] k = P S Θ [t] k -τ ∇ Θ k f (v +[t] k , v -[t] k , Θ [t] k ) -τ (Λ [t] k,1 + Λ [t] k,2 ) ∀k = 1, . . . , K Λ +[t+1] k,1 = prox σλ 2 • * 1 Λ +[t] k,1 + σ(2Θ [t+1] k -Θ [t] k ) ∀k = 1, . . . , K For k = 1, . . . , K, i = 1, . . . , N (s +(i)[t+1] k , s -(i)[t+1] k , Λ (i,•)[t+1] k,2 ) = prox σι * Er s +(i)[t] k + σ(2v +[t+1] k -v +[t] ) k , s -(i)[t] k + σ(2v -(i)[t+1] k -v -(i)[t] k ), Λ +(i,•)[t] k,2 + σ(2Θ (i,•)[t+1] k -Θ (i,•)[t] k ) 6 Experiments
In this section, we firstly provide a comparison between 1 and ∞ proposed approaches with the two closest state-of-the-art procedures, that are "hierNet" [START_REF] Bien | A lasso for hierarchical interactions[END_REF] and Jenatton's framework [START_REF] Jenatton | Proximal methods for hierarchical sparse coding[END_REF], on a simulation dataset into regression framework. Comparisons are of two types, in terms of convergence behavior and in term of performing estimation.

Second, we apply the regression algorithms to HIV disease analysis, which is meaningful to the prevention of the disease. Third, classification experiments are conducted in the context of Parkinson disease classification and face classification. The proposed algorithms and the comparison approaches are implemented in Matlab on a computer with AMD AthlonX4 750 processor.

Simulated data 6.1.1 Dataset

The dataset is created according to [START_REF] Haris | Convex modeling of interactions with strong heredity[END_REF]. It is initially composed of N main features and of N (N -1)/2 interactions (due to symmetry). We denote (v, Θ) ∈ R N × R N ×N the ground truth generated according to strong hierarchy. v denotes a sparse vector where the nonzero values are associated with randomly selected indexes. The value for the non-zero v (i) is randomly selected from the set {-5, -4, . . . , -1, 1, . . . , 5}. Due to strong hierarchical constraint, Θ is only non-zero for those whose main effects are not zero. The values are randomly chosen from the set {-10, -8, . . . , -2, 2, . . . , 8, 10}.

The algorithmic performance is evaluated on two simulated datasets. The first one is composed of N = 30 features, the first 10 main features are non-zeros and the interaction ratio is ρ = 3.45% (Dataset30-345). The second dataset is of size N = 100, the first 30 features are non-zeros and the interaction ratio is ρ = 0.30% (Dataset100-030). The interaction ratio is calculated as the non-zero interaction number divided by the possible interaction number in Θ (i.e. N (N -1)/2).

The datasets are composed of a training, a validation and a testing set with 100 samples each. φ(x ) is randomly generated according to normal distribution N (0, I N ). The observed value for each sample is set by y = φ(x) v + φ(x ) Θφ(x ) + ε , where ε is independent Gaussian noise to make the signal-to-noise ratio approximately equal to 5dB.

Comparison with [24]r = ∞ and weak hierarchy

In order to provide fair comparison with the work in [START_REF] Jenatton | Proximal methods for hierarchical sparse coding[END_REF], we first investigate the performance of proposed primal-dual proximal algorithm (Algorithm 2) when r = ∞ and in the configuration of weak hierarchy (no symmetry constraint is considered, i.e. C = R N ×N ): "Weak-PD-∞ ". The algorithmic procedure designed in [START_REF] Jenatton | Proximal methods for hierarchical sparse coding[END_REF] is based on "FISTA" and it is named "Weak-FISTA-∞ ", where we use the proximity operator of tree-structured variables in "SPAMS" toolbox [START_REF] Jenatton | Proximal methods for hierarchical sparse coding[END_REF], it is worthy mentioning that the core function is written in C++, while the whole implementation of our proposed epigraphical primal-dual algorithm is implemented in MATLAB. Both are compared in two respects: the convergence of the objective function Eq. ( 5) and convergence of the iterates (i.e. w [k] w [+∞] ).

The comparisons are displayed in Fig. 2 (for Dataset30-345) and in Fig. 3 (for Dataset100-030) for the optimal λ (cf. Section 6.1.4). The convergence behaviour are presented both w.r.t. iterations (first row) and time (second row). It is observed that: i) "Weak-PD-∞ " and "Weak-FISTA-∞ " converge to the same value as shown in the first column of both Fig. 2 and3;

ii) From the objective function convergence point-of-view, "Weak-FISTA-∞ " converges faster than 'Weak-PD-∞ " w.r.t. iteration numbers and also time;

iii) From the convergence of the iterates point of view (second column of Fig. 2 and3), for Dataset30-345, "Weak-PD-∞ " converges much faster than "Weak-FISTA-∞ " to the optimal solution, especially from the comparison in terms of time (bottom figures). For Dataset100-030, although it seems that the convergence in iteration for "Weak-PD-∞ " slower, but it costs less time, this may because that the default values for τ and σ are not optimal for Dataset100-030. It shows that the proposed algorithm enables to be closest to the optimal solution with less time; iv) Contrary to"Weak-FISTA-∞ ", our proposed solution can provide a 'Strong-PD-∞ " counterpart of this 'Weak-PD-∞ " without inner iterations.

Further results, especially w.r.t the choice of λ and regression performance on the training/validation/test sets are provided in Section 6.1.4.

Comparison with [7]r = 1 and strong hierarchy

Similar experiments are conducted for r = 1 and strong hierarchy (i.e., C = S). The proposed algorithm (Algorithm 4 with r = 1) is called "Strong-PD-1 " and it is compared with "hierNet" whose iterations are derived from an ADMM scheme, named"Strong-ADMM-1 ", it is noted that the hyperparameter in the inner iteration for ADMM is set to be the default value (i.e. 1) in the following experiments.

The comparisons between these two schemes are displayed in Fig. 4 (for Dataset30-345) and Fig. 5 (for Dataset100-030) for the optimal λ (cf. Section 6.1.4). Both are compared in three respects: the convergence of the objective function Eq. ( 3), the convergence of the iterates (i.e. w [k] w [+∞] ), and the distance to the set S. These quantities are displayed w.r.t. iteration number and time. It is observed that: i) "Strong-PD-1 " and "Strong-ADMM-1 " converge to the same solution, but the convergence of the objective to the optimum solution is very sensitive to the trade-off hyper-parameter for "Strong-ADMM-1 ";

ii) "Strong-PD-1 " is always faster either in iteration number or in time and either in terms of functional or in terms of iterations. The explanation mainly comes from the inner iterations required with "Strong-ADMM-1 " when dealing with C = S;

iii) Third column of Fig. 4 and 5 highlights that the constraint violations (i.e. distance to S) with the proposed method is always smaller than with ADMM;

Further results, especially w.r.t the choice of λ and regression performance on the training/validation/test sets are provided in Section 6.1.4.

Choice of λ

In order to select the optimal λ, we perform 5 simulations and compute the average performance. Hold-out validation method is adopted for the selection of the trade-off parameter λ, where the models with different λ are trained on the training set, and then the best λ is selected based on the performance on the validation set, finally we give the performance of the model on the test set with the selected λ. We also compare the proposed algorithms to other two common methods:

• SVR-2 : it attempts to minimize the criteria with an empirical loss and a 2 -norm on the weights, aiming to learn small weights, which is expressed as:

minimize w∈R N +N ×N λ L =1 y -[φ(x ), φ(x ) φ(x )] w 2 + 1 2 w 2 (17)
We apply LIBSVM library [START_REF] Chang | Libsvm: A library for support vector machines[END_REF] to solve this problem in the dual form with linear kernels on the training samples. In our experiments, λ is validated from 10 to 1000 with a space 10, and the results with the best selected λ are shown in Tab. 1.

• SVR-1 : it integrates the prior knowledge of sparsity to the criteria and aims to learn a set of sparse weights, which can be written as:

minimize w∈R N +N ×N L =1 1 2 y -[φ(x ), φ(x ) φ(x )] w 2 + λ w 1 (18)
We apply forward-backward proximal algorithm [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] to solve it. In our experiments, λ is validated in [10 -6 , 10 -5 , . . . , 1, 2, 4, . . . , 2 5 ], and the results with the best selected λ are shown in Tab. 1.

The average performance in mean square errors (MSE) for different interaction models with different λ on validation set are shown in Fig. 6, and their performance with best selected λ on the three sets are given in Tab. 1. From the performance point of view in the Fig. 6 and Tab. 1, it can be seen that: i) The performance SVR-2 are the worst, which is not surprise because the penalization term aims to learn small values, rather than sparse ones. SVR-1 works better and sparse weights are learned, which meets the sparse property of interactions, but it is still inferior to strong hierarchy with r = 1; ii) "Weak-PD-∞ " and "Weak-FISTA-∞ " almost have the same performance for different λ on both simulation datasets; iii) For Dataset30-345, "Strong-ADMM-1 " have similar performance with "Strong-PD-1 ", but due to sophisticated hyper-parameter selection in ADMM, fluctuations over different λ are observed, especially for large λ values; for Dataset100-030, there are some gaps between different λ, and the performance of "Strong-ADMM-1 " deteriorates seriously when λ becomes large. We conjecture that ADMM would become sensitive to the hyper-parameter selection when the number of features increases.

Discussion regarding the choice of r

From the above algorithmic comparisons, we have observed that the proposed method, either for 1 or ∞ , delivers an accurate solution (as other algorithmic strategies) but faster and with the possibility to include the strong-hierarchy constraint without inner iterations.

In the following analysis, we will thus focus on comparisons between weak/strong 1 or ∞ using the proposed algorithm ("Weak-PD-∞ " , "Weak-PD-1 " , "Strong-PD-∞ " , "Strong-PD-1 "). The average performance with the associated variances are presented in Fig. 6 for different values of the regularization parameter λ. Moreover, the performance obtained by hold-out validation are presented in Tab. 1. It can be observed that:

i) The good behavior of the strong hierarchy constraint clearly appears for both datasets. Indeed, we recall that the data have been created with strong hierarchical structure and we can clearly observe that for both datasets the performances with the strong hierarchical constraint are always better either for r = 1 or r = +∞;

ii) In our set of experiments, the regularization with r = 1 always leads to better performance than with r = +∞;

iii) For "Weak-PD-∞ ", MSE values associated with the Dataset100-030 are larger, probably due to overfitting. This assumption can be validated by the results obtained on the training dataset.

Application to HIV Data

Feature interaction learning is very meaningful to discover the intrinsic correlations between the variables, especially in the real-life applications, for example, genetic diseases and drug analysis. Since the effects usually depends on several genes or factors, the interaction of these associating factors is worth studying. In this section, we apply the proposed algorithms to drug analysis: HIV dataset on the susceptibility of the HIV-1 virus to six nucleoside reverse transcriptase inhibitors (NRTIs). This dataset4 was collected by [START_REF] Rhee | Genotypic predictors of human immunodeficiency virus type 1 drug resistance[END_REF] and it was used to model HIV-1 susceptibility to the drugs because HIV-1 virus can become resistant through genome mutations for different subjects. In the dataset, there are 639 subjects and 240 gene locations. For each observation, the mutation status at each gene location is recorded and also give six drugs (log) susceptibility responses. In this work, we focus on the prediction model for 3TC drug.

The dataset is randomly split into two half sets for training and test, respectively as adopted by [START_REF] Bien | A lasso for hierarchical interactions[END_REF][START_REF] Haris | Convex modeling of interactions with strong heredity[END_REF]. The trade-off parameter λ is selected similarly by hold-out validation method. The performance are measured by the average MSE with the associating variance over 5 random splits on the test set for "Weak-PD-∞ ", "Strong-PD-∞ ", "Weak-PD-1 ", "Strong-PD-1 ", "Weak-FISTA-∞ " and "Strong-ADMM-1 ".

Reduced dataset

Following [START_REF] Haris | Convex modeling of interactions with strong heredity[END_REF], we first work on a reduced dataset over the bins of ten adjacent loci, rather than all 240 genes locations, resulting from the fact that nearby genomes have similar effects to the drug susceptibility, and also leading to less sparse data. So for the first set of experiments we have N = 24 features. The value for each bin is set to 1 if one of genes in that bin undergoes mutation. Fig. 7 (left) shows the performance of six algorithms on the reduced dataset. It can be observed that: i) "Weak-FISTA-∞ has the same performance with "Weak-PD-∞ over different λ; ii) a slight better performance are obtained for "Strong-PD-1 " compared to "Strong-ADMM-1 "; iii) "Strong-PD-∞ " and "Strong-PD-1 " are slightly better than their weaker counterparts, and the best λ for "Strong-PD-∞ " and "Strong-PD-1 " are 12 and 10 respectively. The interactions over one split are visualized in the Fig. 8 (left). In order to better visualize the effect of the penalizations, we also plot the ∞ -norm and 1 -norm over each row of Θ for strong hierarchy, as shown in the middle and right of Fig. 8. Both have the strong effect for 19 th feature, which is consistent with the results observed in [START_REF] Haris | Convex modeling of interactions with strong heredity[END_REF]. We also observe an interesting fact from the visualization of interactions, it can be seen that "Strong-PD-1 " is able to learn a sparser interaction matrix than "Strong-PD-∞ ", however, the maximum strength (i.e. ∞ -norm) and the sum of the strength ( 1 -norm) in Fig. 8 have similar distributions for both cases. The possible reason is that the extra interactions learned by "Strong-PD-∞ " are subtle, it can be also confirmed that both performances are very similar in the Fig. 7 (left).

Full dataset

We further work on the full data without binning operation leading to N = 240 features and each feature gives the indicator of mutation. The data splitting is the same as the one previously described. The average performance of MSE on the test set for six algorithms are shown in Fig. 7 (right). It is clear to demonstrate that: i) on large dataset, there are some performance gaps between FISTA, ADMM and their primal-dual counterparts for weak and strong hierarchy over different λ, especially ADMM deteriorates seriously and the performance becomes worse; ii) strong hierarchy still produces better performance than the weak one. The learned interactions and their feature effects from "Strong-PD-∞ " and "Strong-PD-1 " with best λ = 1 over one split are shown in Fig. 9. It is found that both algorithms can detect that 184 th genes location has the most strong effect. From the maximum strength and sum of strength for each gene in the Fig. 9, we can see the different properties for both algorithms. For "Strong-PD-∞ ", maximum strength for each gene has a more sparse distribution than the sum of strength, which is consistent with its objective that the maximum of absolute values of Θ (i,:) is not larger that v (i) , whereas a more sparse distribution over the sum of strength is obtained for "Strong-PD-1 ", whose objective is to ensure that the sum of absolute values of Θ (i,:) is not larger that v (i) . From Fig. 7(right), it also can be observed that "Strong-PD-1 " is slightly better than "Strong-PD-∞ " because it imposes a stronger penalization for Θ than "Strong-PD-∞ ".

Multiclass learning

In [START_REF] Jiu | Multiclass SVM with hierarchical interaction: application to face classification[END_REF], we performed preliminary validation of the proposed method compared to stateof-the-art in the context of face classification. In this work, we tackle Parkinson disease classification from the speech records of the persons. Thus, x models the speech data (1D) and φ(x ) the time-frequency based features. For all the experiments, we compare the performance in term of accuracy between:

• sparse multiclass SVM without interaction, i.e., f k

(x ) = v k φ(x ) and Ω(v k ) = k v k 1 ,
• sparse multiclass SVM with full interactions:

f k (x ) = φ(x ) Θ k φ(x ) + v k φ(x ) and Ω(v k , Θ k ) = v k 1 + i Θ (i,•) k 1 ,
• proposed approach that is sparse multiclass SVM with strong hierarchy.

The full interaction case differs from the proposed one in the sense that there is no coupling between v k and Θ k in Eq. ( 5). Next we show the results of multiclass SVMs with full interaction and strong hierarchy, taking into account the interactions of features. Fig. 10 shows the performance comparison curves under different combinations of λ 1 and λ 2 for multiclass SVMs with full interactions and strong hierarchy. It is observed that: 1) the accuracy with full interaction on the training set are much higher than multiclass SVMs with 1 -norm and strong hierarchy, but the accuracy on the validation and test set are not reached to as much as on the training set, the possible reason is that by combining interactions of features without further constraints, overfitting problem occurs when the learned weights (i.e. v and Θ) are increased; 2) when strong hierarchy is regularized to v and Θ, the overfitting is alleviated, and the training accuracy are decreased, while the accuracy on the validation and test set are improved compared to the one with full interactions; 3) From Tab. 2 and Fig. 10, the best accuracy on the test set for multiclass SVMs with 1 -norm is 63.98% (λ = 10 -2 ), while the best accuracy for full interaction and strong hierarchy are 65.99% (λ 1 = 10 -4 , λ 2 = 10 -2 ) and 66.57% (λ 1 = 10 -4 , λ 2 = 10 -4 ) respectively, which validates that strong hierarchy can help to improve the classification performance.

Dataset

Conclusion

In this work, we propose a primal-dual proximal algorithm with epigraphical projection for regression and multiclass SVMs with strong hierarchy and apply it to several regression and classification task. Four algorithms are derived ("Weak-PD-∞ " , "Weak-PD-1 " , "Strong-PD-∞ " , "Strong-PD-1 ") allowing to deal with and to properly compare different configurations of the regularization term (4). Compared to other state-of-the-art algorithms, for instance, FISTA or ADMM, the proposed algorithm is computationally more efficient (finding a solution belonging to S), and convergence in terms of iterates is faster. Compared to standard sparse learning strategies, we can integrate the main effects and interaction effects, and then obtain the underlying graphs between the features, as a result, the regularization with strong hierarchy make them more discriminative.

Appendix A: Proof of Proposition 3.1

The proof follows arguments derived in [START_REF] Bien | A lasso for hierarchical interactions[END_REF] for the specific case r = 1. Let v = v + -v -, (5) can be equivalently written as

minimize v,v + ,v -,Θ 1 2 L =1 y -φ (x )v -φ (x )Θφ(x ) 2 + λ 2 Θ 1 + λ1 (2v + -v) subj. to      Θ (i,•) r ≤ 2v +(i) -v (i) (∀i ∈ {1, . . . , N }) v = v + -v - (v + , v -, Θ) ∈ O × O × C. ( 19
)
The constraints involving (v, v + , v -) can be reformulated as (∀i ∈ {1, . . . , N })

     v +(i) ≥ 1 2 ( Θ (i,•) r + v (i) ), v +(i) ≥ v (i) v +(i) ≥ 0.
and that leads to

η -= ω --(N -n+1)(ω + -ω -)+ N i=n ν (i) 1+2(N -n+1) , η + = η -+ ω + -ω -. (25)

Figure 2 :

 2 Figure 2: Comparison between the proposed "Weak-PD-∞ " and "Weak-FISTA-∞ " on Dataset30-345 for λ = 14. (top-left) Objective function in (3) w.r.t. iterations, (bottomleft) Objective function in w.r.t. time, (top-right) Distance w [k] -w [∞] w.r.t. iterations, (bottom-right) Distance w [k] -w [∞] w.r.t. time.

Figure 3 :

 3 Figure 3: Comparison between the proposed "Weak-PD-∞ " and "Weak-FISTA-∞ on Dataset100-030 for λ = 8. (top-left) Objective function (3) w.r.t. iterations, (bottomleft) Objective function (3) w.r.t. time, (top-right) Distance w [k] -w [∞] w.r.t. iterations, (bottom-right) Distance w [k] -w [∞] w.r.t. time.

Figure 4 :

 4 Figure 4: Comparison between the proposed "Strong-PD-1 " and "Strong-ADMM-1 " on Dataset30-345 for λ = 8. (top-left) Objective function (3) w.r.t. iterations, (bottom-left) Objective function (3) w.r.t. time, (top-middle) Distance w [k] -w [∞] w.r.t. iterations, (bottom-middle) Distance w [k] -w [∞] w.r.t. time, (top-right) Distance to the strong hierarchy constraint S w.r.t. iterations, (bottom-right) Distance to the strong hierarchy constraint S w.r.t. time.

Figure 5 :

 5 Figure 5: Comparison between the proposed "Strong-PD-1 " and "Strong-ADMM-1 " on Dataset100-030 for λ = 10. (top-left) Objective function (3) w.r.t. iterations, (bottom-left) Objective function (3) w.r.t. time, (top-middle) Distance w [k] -w [∞] w.r.t iterations, (bottom-middle) Distance w [k] -w [∞] w.r.t. time, (top-right) Distance to the strong hierarchy constraint S w.r.t. iterations, (bottom-right) Distance to the strong hierarchy constraint S w.r.t. time.

Figure 6 :

 6 Figure 6: The performance in MSE on the validation set obtained with the six algorithms ("Weak-FISTA-∞ ", "Weak-PD-∞ ", "Strong-PD-∞ , "Weak-PD-1 ","Strong-ADMM-1 and "Strong-PD-1 ) allowing to compare properly the different configurations of the regularization term (4). (left) Dataset30-345 (right) Dataset100-030).

Figure 7 :

 7 Figure 7: Performance comparisons between "Weak-FISTA-∞ ", "Weak-PD-∞ ", "Strong-PD-∞ , "Weak-PD-1 ","Strong-ADMM-1 and "Strong-PD-1 for different λ values on the test set of HIV dataset when N = 24 (left) and N = 240 (right).

Figure 8 :

 8 Figure 8: The figures in the first row shows the interaction (left column), ∞ -norm over each row of interactions (middle column), 1 -norm over each row of interaction (right column) when N = 24, r = ∞ and λ = 12. The figures in the second row shows the corresponding results when when N = 24, r = 1 and λ = 10.

Figure 9 :

 9 Figure 9: The figures in the first row shows the interaction (left column), ∞ -norm over each row of interactions (middle column), 1 -norm over each row of interactions (right column) when N = 240, r = ∞ and λ = 1. The figures in the second row shows the corresponding results when when N = 240, r = 1 and λ = 1.

Figure 10 :

 10 Figure 10: This figure shows classification accuracy comparison on the training, validation and test for multiclass SVMs with full interactions and strong hierarchy under different combinations of λ 1 and λ 2 on Parkinson Speech dataset.

Table 1 :

 1 The comparison results (MSE) on the train, validation and test set of different algorithms under best selected λ when both datasets are used.

	Data	Approach	λ	TR	VAL	TE
		SVR-2	120 81.698±43.091 1046.957±135.988 1064.482±264.993
		SVR-1	0.01 0.117±0.013	0.520±0.038	0.500±0.036
		Weak-FISTA-∞	14	0.166±0.017	0.642±0.065	0.602±0.057
	Dataset30-347	Weak-PD-∞ Strong-PD-∞	14 16	0.166±0.017 0.190±0.017	0.642±0.065 0.538±0.029	0.602±0.057 0.527±0.048
		weak-PD-1	8	0.198±0.017	0.487±0.032	0.481±0.038
		Strong-ADMM-1	8	0.184±0.023	0.471±0.045	0.471±0.061
		Strong-PD-1	8	0.196±0.017	0.478±0.028	0.473±0.041
		SVR-2	90 96.217±61.627 1843.295±62.330 1870.404±174.944
		SVR-1	0.1	1.496±0.191	2.918±1.596	2.831±1.255
		Weak-FISTA-∞	10	0.032±0.001	2.954±0.649	3.049±0.749
	Dataset100-030	Weak-PD-∞ Strong-PD-∞	2 18	0.001±0.0002 0.101±0.004	2.867±0.808 1.034±0.245	2.999±0.943 1.025±0.285
		weak-PD-1	10	0.109±0.006	1.050±0.217	1.062±0.251
		Strong-ADMM-1	4	0.025±0.001	0.884±0.252	0.814±0.212
		Strong-PD-1	10	0.113±0.006	0.918±0.179	0.929±0.199

  -We first show the Parkinson disease classification by sparse multiclass SVM with strong hierarchy. The experiments are conducted on the Parkinson Speech dataset[START_REF] Sakar | Collection and analysis of a parkinson speech dataset with multiple types of sound recordings[END_REF], which is used to study the underlying relationship between the Parkinson disease and the types of voices of the patients. This dataset consists of training set and test set. For the training set, the data are collected from 20 Parkinson patients (6 female, 14 male) and 20 healthy individuals (10 female, 10 male), for each subject, 26 sound samples including sustained vowels, words, and short sentences are recorded, containing in total 1040 training samples. 26 type of time-frequency based features (such as jitter, shimmer, median pitch and number of pulses etc.) are extracted for each sound samples. For the test set, the data are only collected from another 28 Parkinson patients. They are asked to say particular sustained vowels ('a' and 'o') three times and the same time-frequency based features are extracted, therefore it contains 168 samples. In our experiments, we only use the training set, because the test set does not contain negative samples. We randomly divide the training set to three subsets of equal size for model training, validation and test. The features are preprocessed by Gaussian normalization and 10 6 number of iterations are adopted to guarantee convergence in the learning stage.Performance -We first show the results of multiclass SVM without interaction in Tab. 2. The accuracy on the training, validation and test set with different λ and also the non-zero rate of v are given. From the results, it is seen that a sparse solution is obtained with λ = 10 -2 , as well as it delivers the best test accuracy (63.98%) on the test set. Train acc. 62.72 65.90 68.50 69.36 69.08 69.08 Val. acc. 58.79 60.23 59.94 58.79 58.50 58.50 Test acc. 62.54 63.98 62.82 61.96 61.96 61.96 NZ rate in v 23.08 65.38 92.31 92.31 92.31 100

			lg(λ)
	-1	-2	-3	-4	-5	-6

Table 2 :

 2 The classification accuracy of multiclass SVM with 1 -norm (in %) on Parkinson Speech dataset.

For every x ∈ H, ι C (x) = 0 if x ∈ C and +∞ otherwise.

Matlab codes will be made available at the time of the publication.

The epigraph of a function f :H →] -∞, +∞] is defined as: epif = {(v, ζ) ∈ H × R | f (v) ≤ ζ}and the epigraphical projection onto epif is denoted P epif .

It can be downloaded from https://hivdb.stanford.edu/pages/published analysis/genophenoPNAS2006/

Acknowledgements 7.1 The authors would like to thank Marie-France Benassy, Mickael Bonnard, Laurent Grosset, Georges Oppenheim, Maxime Durot, Leire Oro-Urea from TO-TAL for helpful discussion, Prof. Mark Asch from Université de Picardie Jules Verne for the comments.

This work was supported by a grant from the National Natural Science Foundation of China (No. 61806180)

yielding v +(i) ≥ max{max{0, v (i) }, 1 2 ( Θ (i),• r + v (i) )}. Using |v (i) | = 2 max{0, v (i) } -v (i) , we get the expected result:

Appendix B: Proof of Proposition 4.1

According to the definition of the epigraphical projection, we solve the following minimization problem:

The inner minimization yields a simple projection for u to [-η + -η -, η + + η -], i.e. p = P [-η + -η -,η + +η -] (u), thus Eq. ( 21) becomes:

and thus

where

We now focus on the computation of prox φ (ω + , ω -).

First, we sort (|u (i) |) 1≤i≤N in ascending order to be (ν (1) , . . . , ν (N ) ), such that |ν (1) | ≤ . . . ≤ |ν (N ) |, and also ν (0) = -∞; and second n ∈ [1, N ] can be found such that ν

, so the proximity operator of φ has:

Appendix C: Proof of Proposition 4.2

The Lagrangian associated to the function involved in P E + 1 is:

The KKT conditions are:

From the fourth condition in Eq. ( 27), we know that if p (i) > 0 then ξ (i) = 0 and p (i) = u (i) -α. Thus, if we denote π = (π (i) ) 1≤i≤N is an ordered version of (p (i) ) 1≤n≤N in a decreasing order, we can write

From similar arguments as in [START_REF] Duchi | Efficient projections onto the l1-ball for learning in high dimensions[END_REF][Lemma 2] , n is computed as [START_REF] Chambolle | On the convergence of the iterates of the "Fast Iterative Shrinkage/Thresholding Algorithm[END_REF] and thus

From the second and third equality of the KKT conditions, we have η + = η -+ ω + -ω -. Finally, combining [START_REF] Lim | Learning interactions through hierarchical group-lasso regularization[END_REF] with η -= ω -+ α, yields


Tristan Charrier

Arthur Queffelec

Ocan Sankur

François Schwarzentruber

Franc ¸ois Schwarzentruber

Reachability and Coverage Planning for Connected Agents

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

A number of use cases of planning rose in informationgathering missions from the development of unmanned autonomous vehicles (UAVs). For instance, in search and rescue missions, a fleet of drones can cover a lot of ground in a short amount of time and report any finding to a mission supervisor to narrow the search for the rescue team. Other examples are the terrain analysis for smart farms and for areas in hazardous locations. For this kind of missions, the information gathered is used for decision making at a supervising station. Thus, the robots need to be constantly in communication with the station to report the information gathered during the mission. The use of multiple UAVs to cover an area not only reduces the time required to complete the mission but also enables reaching locations which may not be reachable with a single drone due to connection constraints.

The original multi-agent path finding problem asks for a plan to reach a configuration of agents in a graph [START_REF] Ratner | [END_REF]]. However, an important problem for search and rescue missions or terrain analysis is the coverage of an area. We thus study both the reachability and the coverage problems under a connection constraint over the agents which requires them to be connected to the base either directly or via another agent, who can relay its data. We establish the computational complexity of the connected coverage in its general case and for a practical subclass introduced recently [START_REF] Tateo | [END_REF] in which the UAVs can communicate with others located within one step, called the neighbor-communicable topological graphs. We show that the coverage is PSPACE-complete in the general case, and remains so for neighbor-communicable topological graphs. Thus, restricting to neighbor-communicable graphs does not make the problem feasible, and the relatively high complexity unfortunately remains. Note that this is in line with the PSPACE-completeness of the reachability problem recently reported in [START_REF] Tateo | [END_REF].

Our main result in this paper is the definition of a class of topological graphs which is well adapted and realistic for UAV missions, and for which the coverage and reachability problems admit efficient algorithms. Our subclass, called sight-moveable graphs, is defined assuming that the UAVs cannot communicate through obstacles and are restricted to line-of-sight communication. This class emerged from an ongoing case study for a drone assisted search and rescue project in which the authors take part1 . For this class, we prove that both the reachability and coverage problems are in LOGSPACE. This drastically changes the status of this problem since by LOGSPACE ⊆ NC (this is the class of problems solvable in polylogarithmic time in a parallel machine with a polynomial number of processors), one can build an efficient parallel algorithm [Cook, 1979]. The bounded versions are NP-complete. This means efficient SAT solvers can be used directly to compute bounded executions.

In this work, we consider anonymous agents. Furthermore, we consider the collisions to be handled by the agents themselves, hence are not considered along the results of this paper. We depicted a covering execution of a topological graph by 3 UAVs in Figure 1. In this example, the UAVs need to gather information at each node of the graph while staying connected to the base (red node) during the whole mission.

In Section 2, we present the typical notions used in Multi-Agent Path Planning (MAPP) and their extension for our case and the known results in connected planning. In Sections 3 to 6, we study the complexity of our problems from the general case to the most restrictive one. We describe the related work in Section 7. We conclude in Section 8.

This paper is the follow-up to the extended abstract presented at the 18th International Conference on Autonomous Agents and Multi Agent Systems (AAMAS 2019) [START_REF] Charrier | [END_REF]. An extended version2 is also available.

Preliminaries

In most applications of path planning, the space is discretized in order to generate a graph of movements on which algorithms are executed. For instance, regular grids which decompose the space in square, triangular or hexagonal cells, irregular grids with techniques such as quadtrees [Finkel and Bentley, 1974;Knoll, 2006] or Voronoï diagrams comprehensively discussed in the survey [Aurenhammer, 1991].

Our work is independent of the particular method used to obtain the discretization. We only work under the hypothesis that a feasible plan on the graph generated by the discretization is also feasible in the continuous space.

Topological Graph

Compared to the graphs used in MAPP, we also consider communication edges which specify whether agents at two different locations can communicate. We call graphs with this additional information topological graphs. The formal definition is the following. Definition 1 (Topological graph). A topological graph is a tuple G = V, →,

, with V a finite set of nodes containing a distinguished element B, → ⊆ V × V a set of movement edges and ⊆ V × V a set of undirected communication edges.

The node B symbolizes the supervision base station from which the agents start the mission. A topological graph is undirected if V, → is an undirected graph.

We will now consider three subclasses of interest.

In most situations, if an agent can move to a location in one step, it can also communicate with an agent at that location. This class has been discussed in [START_REF] Tateo | [END_REF]. We call topological graphs satisfying this requirement neighborcommunicable. An example is given in Figure 2b.

Definition 2 (Neighbor-Communicable topological graph). A neighbor-communicable topological graph is a topological graph such that v → v implies v v .

Another class of graphs is that of sight-moveable and is the main one for which we give efficient algorithms. First, this class requires the movement edges to be undirected and reflexive. Second, whenever an agent can communicate with another node, then it can also move to that node while maintaining the communication. This intuitively means that the communication is restricted to line-of-sight and is disallowed through obstacles. The formal definition follows, and an example is depicted in Figure 2c.

Definition 3 (Sight-Moveable topological graph). A sightmoveable topological graph is an undirected neighborcommunicable topological graph in which for all v ∈ V , v → v and whenever v v , there exists a sequence ρ = ρ 1 , . . . , ρ n of nodes such that v = ρ 1 , v = ρ n , v ρ i and ρ i → ρ i+1 for all i ∈ {1, . . . , n}.

Last, we define the complete-communication topological graphs which are simply sight-moveable topological graphs with a complete communication topology. This subclass can model that the communication is not perturbed in the area. An example of such a graph is depicted in Figure 2d, and the formal definition is the following.

Definition 4 (Complete-Communication topological graph).

A complete-communication topological graph is a sightmoveable topological graph such that = V × V .

Observe that complete-communication graphs are reflexive, undirected, connected graphs with = V × V .

Execution

An execution, in MAPP, is a finite sequence of configurations describing the placement of the agents during the mission.

The formal definition of a configuration is the following.

Definition 5 (Configuration). A configuration c of n agents in a topological graph G is an element of V n denoted c = c 1 , . . . , c n in which c i is the location of agent i such that the graph V a , ∩ (V a ×V a ) is connected with V a = {B, c 1 , . . . , c n }. We extend the notation → and denote c → c when c i → c i for all 1 ≤ i ≤ n.
MAPP asks to associate an agent to a specific goal. However, given that we are interested in covering an area with a fleet of agents, the anonymity is useful to get more efficient plans.

Anonymity. In the rest of this paper, agents are anonymous. In other words, a configuration c is equivalent to a configuration c iff c is a permutation of c .

Moreover, an important notion in MAPP is the computation of collision-free plans. In the drone case, in which we are particularly interested, one can place drones at different heights to avoid collisions. Additionally, most drones are, nowadays, equipped with local collision avoidance systems.

Collisions. We do not deal with meet-or head-oncollisions of agents, i.e. we allow two agents to be located in a same node, and to move in opposite directions of an edge within a step.

An execution e of length with n agents in a graph G is a sequence of configuration c 1 , . . . , c such that for c i → c i+1 for all 1 ≤ i < .

A covering execution e = c 1 , . . . c of length with n agents in a graph G is an execution such that c 1 = c = B, . . . , B and for all v ∈ V , there exists i ∈ {1, . . . , } with v appearing in c i .

Decision Problems

We define the MAPP problems, the Reachability problem along with its bounded version, bReachability, for the makespan optimization of the plan. In addition, we define the Coverage problem and the bounded coverage, bCoverage. Definition 6 (Reachability). Given a topological graph G, n ∈ N written in unary and a configuration c of size n, decide if there is an execution c 1 , . . . , c in G such that c 1 = B, . . . , B and c = c. Definition 7 (bReachability). Given a topological graph G, n ∈ N written in unary and a configuration c of size n and ∈ N written in unary, decide if there is an execution c 1 , . . . , c in G s.t. ≤ and c = c. Definition 8 (Coverage). Given a topological graph G and n ∈ N written in unary, decide if there exists a covering execution with n agents. Definition 9 (bCoverage). Given a topological graph G, n, ∈ N written in unary, decide if there exists a covering execution of length such ≤ .

We study the restrictions of the above problems to classes of topological graphs. We denote P C the problem P (one of the four above problems) restricted to a class C of topological graphs (C can either be dir for directed, nc for neighborcommunicable, und for undirected, sm for sight-moveable or cc for complete-communication topological graphs).

Known Results

The complexity of the decision problem associated to the minimization of the makespan with non-anonymous agents and collision, is known to be NP-hard [START_REF] Ratner | [END_REF]. Throughout the study of MAPP, NP-hardness was shown to hold on planar graphs [Yu, 2016] and, later, on 2D grid graphs [Banfi et al., 2017]. Variants of MAPP have been studied such as the package-exchange robot-routing problem [START_REF] Ma | [END_REF] where the robots are anonymous but not the package they exchange, is shown to be NP-hard. A class of grid graphs was shown to be solvable in polynomial time [START_REF] Wang | [END_REF].

The connected version of MAPP was introduced in [START_REF] Hollinger | [END_REF], in which a topological graph discretizes the space and it is proved that the existence of a plan for the reachability of a configuration of agents in a bounded amount of steps is NP-hard: Theorem 10. bReachability restricted to undirected topological graphs is NP-hard [START_REF] Hollinger | [END_REF].

In [START_REF] Tateo | [END_REF], it is shown that deciding the existence of a feasible plan is PSPACE-complete: Theorem 11. Reachability restricted to undirected topological graphs is PSPACE-complete [START_REF] Tateo | [END_REF].

Authors prove this result for graphs with self-loops and a base [START_REF] Tateo | [END_REF] as in our setting (see Discussion following Theorem 1). The only difference with our setting is that the agents start at a specific configuration in [START_REF] Tateo | [END_REF]. Nevertheless, it can be shown that their complexity result holds for our problem by a simple but subtle construction given in the extended version.

In the rest of the paper, we study the upper bounds and the lower bounds complexity of the defined decision problems on the previously defined topological graphs. The following sections present our results, respectively, for the general case, the neighbor-communicable graphs, sight-moveable graphs, and complete-communication graphs.

Directed Topological Graphs

For the bounded versions, we can guess and check a path of bounded length in polynomial-time since the input is encoded in unary: Proposition 12. bCoverage dir and bReachability dir are in NP.

For the unbounded problems, we can design a straightforward NPSPACE algorithm that guesses an execution by keeping in memory the last configuration, and, for Coverage dir , the set of visited regions. We conclude with Savitch's Theorem (NPSPACE=PSPACE) [START_REF] Savitch | [END_REF]: Theorem 13. Coverage dir and Reachability dir are PSPACE-complete.

The lower bound of Reachability dir is given in Theorem 11. We now concentrate on Coverage dir . Lemma 14. Coverage dir is PSPACE-hard.

Proof. The proof is by reduction from Reachability dir in which the base node has a self-loop. As noted in the remark following Theorem 11, this problem remains PSPACE-hard. We map an instance (G, c) of Reachability dir to the instance G of Coverage dir where G is depicted in Fig. 3. Let k denote the number of agents in the instance (G, c). G contains G as a subgraph, plus fresh nodes v 1 , . . . , v k and s 1 , . . . , s k . An agent can move from any node of G to v 1 and back.

Node s 1 can communicate with the base B, and node v k can communicate with all nodes of G . Furthermore, we have the communication edges (s i , s i+1 ) and i ≤ k -1. Now we prove that the k agents can progress to the configuration (c 1 , . . . , c k ) in G if and only if there exists a covering execution in G .

(v i , v i+1 ) for all 1 ≤ B s 1 s 2 s k c 1 c 2 c k v 1 v 2 v k all copy of G
(⇒) If the agents are in the configuration (c 1 , . . . , c k ) then they can progress in one step to configuration (s 1 , . . . , s k ). Then, they have no choice but progress to the configuration (v 1 , . . . , v k ). Once in this configuration, the agent placed on the node v k communicates with the base and with all other agents. This agent stays at v k . Meanwhile the agent placed on the node v 1 will visit all unvisited nodes of G and come back to v 1 while keeping communication to the base through the agent placed on v k . Meanwhile, agents placed on v 2 , . . . , v k-1 come back to B. Finally, when all the nodes have been visited, both agents on v 1 and v k come back to B.

(⇐) If there exists a covering execution of the whole graph G , it means all nodes have been visited. In particular, node s k has been visited and let us consider the first time t s k when s k is visited. Time t s k -1 denotes the time just before t s k .

Fact 15. At time t s k -1, no node v i and no node s i were visited.

Proof. Suppose by contradiction that a node v i was visited by some agent before t s k , then the only possibility such an agent to communicate to the base is that there is also an agent at v k at time t s k . But then, it means that s k was visited strictly before t s k , leading to a contradiction. Thus, no node v i were visited at time t s k (thus at time t s k -1).

As no node v i are visited before t s k , no node s i are visited before t s k -1.

Fact 16. At time t s k -1, the configuration is c 1 , . . . , c k .

Proof. At time t s k , as the agent at s k needs to communicate with the base, the only possibility is that the configuration is s 1 , . . . , s k . Thus, the only possibility is that configuration is c 1 , . . . , c k .

Facts 15 implies that the prefix from time 0 to time t s k -1 of the covering execution is an execution in G. Fact 16 implies that sub-execution reaches c 1 , . . . , c k .

Neighbor-Communicable Topological Graphs

In this subsection, we show that our problems remain hard for neighbor-communicable graphs.

Theorem 17. Coverage nc is PSPACE-complete.

Proof. The upper bound is given by Theorem 13. For the lower bound on Coverage nc , the reduction given in Figure 3 is not adapted for neighbor communicable graphs. Indeed, all nodes may be visited although c 1 , . . . , c k is not reached: v 1 and v k can be reached by two lines of agents connected to the base, making the coverage of the full graph possible. We nevertheless give a similar reduction by adapting the previous reduction.

The details are given in the extended version.

Sight-Moveable Topological Graphs

In this subsection, we show that Reachability sm and Coverage sm are in LOGSPACE while the bounded version bReachability sm is NP-complete.

Upper Bounds

The results of this subsection rely on the problem of checking the connectivity of two nodes s and t in an undirected graph, namely USTCONN. Theorem 18. USTCONN is in LOGSPACE [Reingold, 2008].

Proposition 19. Reachability sm is in LOGSPACE.

Proof. The idea of the proof is to reduce Reachability sm to UCONN, that is the problem of deciding whether an undirected graph is connected. From Theorem 18, we can reduce UCONN to USTCONN by simply looping over all pairs of nodes (s, t) and checking for a path from s to t. Therefore, UCONN is in LOGSPACE. Now we describe the logarithmic space reduction of Reachability sm to UCONN. Let G = V, →, a sight-moveable topological graph and c a configuration. Let V = {c 1 , . . . , c n , B}. The configuration c is reachable iff the restriction of G := (V,

) to the nodes in V is -connected. Indeed, if it is, then c is reachable: each agent follows some →-path from B to c i contained in a -path from B to c i . In other words, (G, c) is a positive Reachability sm -instance iff G is a positive UCONN-instance. The reduction is in logarithmic space: we compute G by enumerating all (u, v)

-edges in G, and we output (u, v) when u, v ∈ V . We recall that we only take into account the working memory for computing G ; the output -G itself -is not taken into account in the used space (see e.g. [Sipser, 1997], Ch. 8, Def. 8.21).

Proposition 20. Coverage sm is in LOGSPACE.

Proof. First we prove that the bounded version of the connectivity in undirected graphs is also in LOGSPACE.

Lemma 21. Bounded-USTCONN, that is the problem, giving an undirected graph G, two nodes s, t, an integer n written in binary, of deciding whether there is a path of length at most n from s to t in G is in LOGSPACE.

Proof. We reduce Bounded-USTCONN to USTCONN in logarithmic space as follows.

From a Bounded-USTCONN instance (G, s, t, n) we construct in logarithmic space a USTCONN instance (G , s , t ): 1. The nodes of G are pairs (v, j) where v is a node of G and j is an integer in {0, n} but smaller than the number of nodes in G ; 2. G contains an edge between (v, j) and (v , j + 1) when there is an edge between v and v in G or when v = v ; 3. s = (s, 0) and t = (t , n).

Now let G = V, →,
be a sight-moveable topological graph and n an integer written in binary. There is a path from any node v to the base B with at most n communication edges iff (G, n) is a positive instance of Coverage. Thus, we test sequentially, for all v, that ((V, ), v, B, n) is a positive instance of Bounded-USTCONN . Hence, we obtain an algorithm in logarithmic space to decide Coverage.

Lower Bounds

We now focus on the NP lower bound of bReachability sm . Proposition 22. bReachability sm is NP-hard for a fixed execution length ≥ 3.

Proof. The proof is by polynomial time reduction from 3-SAT problem (see [Karp, 1972]). Given a 3-SAT instance, set of clauses c 1 , . . . , c m with variables x 1 , . . . , x n , we describe the construction of an instance (G, c) of bReachability sm with k = n + m agents.

The topological graph G = V, →, is constructed as follows. We start by placing the base B from which the agents start their mission.

Please recall that a sight-moveable graph is also a neighbor-communicable graph so all movements edges are also communication edges in the construction below even if not explicitly stated.

For each variable x, we construct a gadget composed of 5 nodes connected to the base depicted in Figure 5a: nodes x, ¬x, staging nodes n x , n ¬x and a goal node g x . We add movement edges from B to n x , from n x to x and from x to g x (resp. from B to n ¬x , from n ¬x to ¬x and from ¬x to g x ). As for the communication, the node x (res. ¬x) communicates with the base.

For each clause c, we construct a gadget composed of 3 nodes depicted in Figure 5b. We create a node c, a staging node n c and a goal node g c . We add movement edges from B to n c , from n c to c and from c to g c . The communication between a clause c and a literal x or ¬x is dictated by the existence of the literal in the clause: c i x j if and only if x j ∈ c i ; and c i ¬x j if and only if ¬x j ∈ c i .

We add movement edges from g xi to g xi+1 , and from g ci to g ci+1 for all 1 ≤ i < n, as well as we from g xn to g c1 . Last, we add a fully connected path containing 3 fresh nodes from g x1 to the base such that g x1 B, in the sense that all nodes of this path have communication edges between them. This translation is polynomial in the number of clauses and variables. The construction is depicted in Figure 4. The snakelike path from g x1 to B is the fully connected path.

From a 3-SAT instance, one can construct the graph G and ask for an execution of length 3 to reach the configuration g x1 , . . . , g xn , g c1 , . . . , g cm .

The rest of the proof is given in the extended version.

From Propositions 12 and 22, we have: Theorem 23. bReachability sm is NP-complete.

Complete-Communication Topological Graphs

The following result relies on the fact that the communication is complete. Proposition 24. bReachability cc is in LOGSPACE.

Proof. From Lemma 21, one can construct an algorithm in LOGSPACE for bReachability cc . Indeed, given a configuration c and ∈ N, the straightforward iteration on the locations c i followed by the verification of a path of at most (given in unary) steps from B to c i yields a sound and complete algorithm for bReachability cc .

Our NP lower bound proof of the bCoverage cc problem is by reduction from the grid Hamiltonian cycle (G-HC) problem which is the Hamiltonian cycle problem restricted to grid graphs and is NP-complete [START_REF] Itai | [END_REF]. Theorem 25. bCoverage cc is NP-complete.

The upper bound follows from Proposition 12. The NPhardness proof is given in the extended version.

Related Work

The coverage planning is an interesting approach to path planning. Indeed, a covering plan can be used for fields such as floor cleaning, lawn mowing, etc. A survey of this field appears in [Choset, 2001]. This multi-agent extension has the ability to reduce the length of the overall mission and also reach parts of the area a single agent would not able to. This problem was studied in [Rekleitis et al., 1997] for two agents. As shown in the survey by Chen et al. [START_REF][END_REF], many coverage problems have been addressed by using analytic techniques. For instance, in [Yanmaz, 2012] and [START_REF] Teacy | [END_REF], they consider UAVs that should cover an area while staying connected to the base, but only empirically study some path planning algorithms without proving formally their soundness and completeness.

We advocate formal methods that give formal guarantees and have already been applied to generate plans for robots and UAVs. Model checking has been applied to robot planning (see [START_REF] Lacerda | [END_REF]) and to UAVs [START_REF] Webster | [END_REF]. Humphrey [Humphrey, 2013] shows how to use LTL (linear-temporal logic) model checking for capturing response and fairness properties in cooperation (for instance, if a task is requested then it is eventually performed).

Bodin et al. [Bodin et al., 2018] treat a similar problem except that the UAVs cover the graph without returning to the base. Without the return-to-the-base constraint, we claim that all our hardness results still hold, except for bCoverage cc . They provide an implementation by describing the problem in Planning Domain Description Language and then run the planner Functional Strips [Francès et al., 2017]. [START_REF] Murano | [END_REF] advocate for a graphtheoretic representations of states, that is, by assigning locations to agents as in Definition 5. In [Aminof et al., 2016 Rubin, 2015], a general formalism is given to specify LTL and monadic second-order logic properties, which are expressive enough to describe the connectivity constraint. They provide an algorithm for parametrized verification in the sense that they a temporal property in a class of graphs. This is relevant for partially-known environments. The algorithm described is non-elementary (i.e. the running time cannot bounded by any tower of exponentials) and therefore not usable in practice. We believe that this is an important problem and our paper identifies an efficient and relevant fragment.

; B x 1 ¬x 1 g x1 x 2 ¬x 2 g x2 x 3 ¬x 3 g x3 c 1 g c1 c 2 g c2

Movement

The multiple traveling salesman problem (mTSP) is a generalization of the traveling salesman problem (TSP) in which multiple salesmen are located at a depot [START_REF] Anbuudayasankar | [END_REF]. mTSP asks for the coverage of all cities so as to minimize the total plan cost by visiting each city exactly once. An overview of TSP and its extensions are presented [START_REF] Matai | [END_REF]. The Coverage problem is related to mTSP, since we use results on Hamiltonian cycle to prove the NP-hardness of bCoverage cc . However, we wish to minimize the length of the execution and not the cost of the execution. Those problems are equivalent on unit graphs, but it is not trivial to use general results on mTSP in order to solve Coverage. Furthermore, to the best of our knowledge, connected versions of mTSP and VRP have not been studied.

Conclusion

Sight-moveable topological graphs we introduced only constrain the communication graph. One can be interested to constrain the movement graph to be planar or a 2D grid given the common usage of grid modelling of the environment. Given the intractability of MAPP on planar graphs [Yu, 2016] and on general 2D grid graphs [Banfi et al., 2017], it is likely that this problem is intractable as well. Furthermore, in [START_REF] Tateo | [END_REF], the decision is proved to stay PSPACE-complete on planar graphs and grids as well. However, one can study this problem on solid grid graphs, given that the Hamiltonian cycle is tractable on such graphs [START_REF] Umans | [END_REF].

Note that our NP lower bounds hold without the anonymity of the agents. Indeed, the bCoverage case is straightforward and for bReachability case, each agent can be associated to a clause or variable, so the reduction would still hold.

We do not know if Coverage remains PSPACE-hard when the →-relation is symmetric (see Figure 6). We think this open issue is important since symmetric →-relations (if UAVs can go from v to v , they can also come back from v to v) are relevant for practical applications. We plan to study the parametrized complexity [START_REF] Downey | [END_REF] of our problems -parameters could be for instance the treewidth of the topological graph or the number of UAVs.

Figure 2 :

 2 Figure 1: Example of a mission execution.

Figure 3 :

 3 Figure 3: Topological graph G constructed from the Reachability dir -instance.

Figure 4 :

 4 Figure 4: Sight-Moveable topological graph computed from the formula (x1 ∨ ¬x2 ∨ x3) ∧ (x2 ∨ x3).

Figure 5 :Figure 6 :

 56 Figure 5: Gadgets in proof of Proposition 22.

EIT Retina project

https://arxiv.org/abs/1903.04300Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence

Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence

Acknowledgements

This work was partially supported by UAV Retina Funded by EIT Digital. Special thanks to Franc ¸ois Bodin for initiating the idea of this work. We thank Eva Soulier for the provided work during her internship.


Tristan Charrier

Arthur Queffelec

Ocan Sankur

François Schwarzentruber

Reachability and Coverage Planning for Connected Agents

Keywords: Multi-Agent Path Planning, Complexity Theory, Connectivity, Coverage, Reachability ACM Reference

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction. Search and rescue missions are difficult to undertake and achieve. In fact, beyond the time constraint of such tasks in order to retrieve the endangered individuals, the place to search can be large as well as hazardous. Thus, unmanned autonomous vehicles (UAV) have caught interest for their ability to be used in unsafe conditions, without (or less) allocation of experts. Furthermore, multiple UAVs can be used to cover a large portion of the place quicker than the experts could in most cases. This approach is used for decision-making and, thus, the UAVs are connected with a base at which an expert can analyze the reported data.

The goal of this work is to study, complexity wise, the problems associated with the usage of connected UAVs. The settings and type of missions we consider are closely linked to the requirements described by the use of UAVs for search and rescue missions. Hence, we examine the complexity of generating reachability and coverage plans of multiple connected agents on topological graphs.

Topological graphs are composed of two types of edges, namely movement and communication edges. The movement edges describe the allowed movements between two nodes. A communication edge enables two agents to communicate from their respective nodes. A plan is a sequence of configurations of agents on the topological graph. For a plan to be valid, it needs to keep the agents connected with the base, which is a specific node of the graph. An agent can be directly connected to the base with a communication edge or it can be connected via other agents who act as relays. Thanks to such relays, agents can cooperatively reach locations behind obstacles which are not directly connected to the base.

Two major types of objectives arise from search and rescue missions. Either an agent is required to reach a specific location in order to sense it, or the agents need to spread over the place to examine several locations. UAVs can, for instance, be used for locating victims of an avalanche by quickly covering as much ground as possible. We explore, in this paper, the complexity of the so-called Reachability problem and the Coveraдe problems along with their bounded versions, bReachability and bCoveraдe.

In Figure 1, we depict a coverage plan execution. The red node denotes the base of the topological graph from which the agents start and need to stay connected with.

The multi-agent path planning (MAPP) with connectivity constraint was introduced in [START_REF] Hollinger | Multirobot Coordination With Periodic Connectivity: Theory and Experiments[END_REF], in which that the existence of a bounded plan for the reachability of a configuration is shown NP-hard. Recently, Tateo et al. showed that the existence of a connected plan for the reachability is PSPACE-complete in [START_REF] Tateo | Multiagent Connected Path Planning: PSPACE-Completeness and How to Deal With It[END_REF]. We extend both of the latter works by considering the connected coverage planning and subclasses of topological graphs on which we can obtain more efficient algorithms for the previously stated problems.

Preliminaries. As mentioned earlier, we study reachability and coverage problems on subclasses of topological graphs. edges and ⊆ V × V is a set of undirected communication edges. The subclasses of topological graphs are named and stated as follows:

A directed topological graph is a tuple G = ⟨V , →, ⟩, V is a set of nodes, →⊆ V × V is a set of directed movement
• A neighbor-communicable topological graph (nc) is a topological graph in which v → v ′ implies v v ′ .
• A sight-moveable topological graph (sm) is a topological graph with → undirected and reflexive and whenever v v ′ there exists a sequence ρ = ⟨ρ 1 , . . . ,

ρ n ⟩ of nodes such that v = ρ 1 , v ′ = ρ n , v ρ i and ρ i → ρ i+1 for all 1 ≤ i < n. • A complete-communication topological graph (cc)
is a topological graph in which → is undirected and reflexive, and

= V × V .
Intuitively, a neighbor-communicable topological graph makes sure that if an agent can move to another node in one step then both nodes are communicating. In sight-moveable ones, the communication is restricted to line-of-sight communication and the agents cannot communicate through obstacles. A complete-communication topological graph can model a situation where a hovering agent connects all agents to the base in the whole area.

On these classes of topological graphs, we study the reachability of a particular configuration of agents in the graph, denoted Reachability. We also study the problem of coverage of the graph, we consider that a plan is covering if and only if every nodes of the graph are visited during the plan by at least one agent, denoted Coverage. Furthermore, we study the bounded versions of the previous problems, respectively denoted bReachability and bCoverage, where a bound on the size of a plan is part of the input, and is encoded in unary. These are the decisions problems associated to the optimization problems for computing minimal plans.

The results stated later consider that the agents are anonymous and, thus, two configurations of agents on the graph are equivalent if one is a reordering of the other. Furthermore, we do not consider head-on-and meet-collisions. Hence, two agents can be located at the same node or they can travel on the same edge in opposite direction during a step of the plan.

Known results. The following theorem is extracted from the original work of Hollinger and Singh in which the bounded connected reachability problem is studied and proved to be NP-hard on undirected topological graph. Thus, we can state this lower bound on directed topological graphs.

THEOREM 1 ([3]

). bReachability dir is NP-hard.

Tateo et al. shown that the existence of a plan for the reachability is PSPACE-complete on directed graphs and discuss that this result extends to directed topological graphs with a base.

THEOREM 2 ([7]

). Reachability dir is PSPACE-complete.

Contributions. We extend previously known results and show that the coverage problems are PSPACE-complete on directed and nc-graphs, and remarkably show that the complexity is reduced to LOGSPACE for sm-and cm-graphs.

Bounded reachability and coverage problems do however remain NP-complete for sm-graphs. All our results are summarized in Fig. 2.

To obtain the LOGSPACE algorithm, we use the result of Reingold who showed that the undirected s-t connectivity problem is in LOGSPACE [START_REF] Reingold | Undirected Connectivity in Log-space[END_REF]. Since LOGSPACE ⊆ NC, it admits polylogarithmic time algorithm on a parallel machine with a polynomial number of processors [START_REF] Cook | Deterministic CFL's Are Accepted Simultaneously in Polynomial Time and Log Squared Space[END_REF]. Although the bounded problems on sight-moveable topological graphs are intractable, it might be possible to use, for instance, SAT solvers, or obtain polynomial-time approximation algorithms.

We lack the complexity result for undirected topological graphs, which is shown by a question mark in the table.

Related work. Coverage planning has been applied to other areas such as lawn mowing and floor cleaning [START_REF] Choset | Coverage for robotics -A survey of recent results[END_REF]. This problem was addressed with analytic techniques [START_REF] Teacy | Maintaining connectivity in UAV swarm sensing[END_REF][START_REF] Yanmaz | Connectivity versus area coverage in unmanned aerial vehicle networks[END_REF]. We advocate for the use of formal methods for provable guarantees on the behaviors; see e.g. [START_REF] Lacerda | Optimal and dynamic planning for Markov decision processes with co-safe LTL specifications[END_REF][START_REF] Webster | Formal methods for the certification of autonomous unmanned aircraft systems[END_REF]. The coverage problem presented is closely related to the traveling salesman problem (TSP) or more precisely, its generalization, multiple TSP. An overview of TSP and its extensions can be found in [START_REF] Matai | Traveling salesman problem: an overview of applications, formulations, and solution approaches[END_REF]. However, to the best of our knowledge, a connected version of mTSP or VRP has not been studied.

Fig. 1 .

 1 Fig. 1. Example of a mission execution.