Nicolas Cavrel
email: nicolas.cavrel@univ-grenoble-alpes.fr

Damien Pellier
email: damien.pellier@univ-grenoble-alpes.fr

Humbert Fiorino
email: humber.fiorino@univ-grenoble-alpes.fr

An Efficient HTN to STRIPS Encoding for Concurrent Plans

The Hierarchical Task Network (HTN) formalism is used to express a wide variety of planning problems in terms of decompositions of tasks into subtaks. Many techniques have been proposed to solve such hierarchical planning problems. A particular technique is to encode hierarchical planning problems as classical STRIPS planning problems. One advantage of this technique is to benefit directly from the constant improvements made by STRIPS planners. However, there are still few effective and expressive encodings. In this paper, we present a new HTN to STRIPS encoding allowing to generate concurrent plans. We show experimentally that this encoding outperforms previous approaches on hierarchical IPC benchmarks.

Introduction

The Hierarchical Task Network (HTN) formalism [Erol, Hendler, and Nau2003] is used to express a wide variety of planning problems in terms of decompositions of tasks into subtaks. HTN planning is used in many applications as, for instance, in task allocation for robot fleets [Milot et al.2021], video games [Menif, Jacopin, and Cazenave2014] or industrial contexts such as software deployment [Georgievski2020]. One possible explanation is that HTN formalism usually fits better for real-world applications and domain experts' mindset: a HTN planning problem is expressed as a set of tasks to achieve rather than an objective state to reach, and the "processes" achieving these tasks as methods, that is to say task decompositions into "simplier" subtasks. Despite the success of hierarchical planning and the recent revival of this planning technique [Bercher, Alford, and Höller2019], there is comparatively less work on hierarchical planning than in classical STRIPS planning [Fikes and Nilsson1971]. The work of the planning community has been more focused on the development of techniques and heuristics for STRIPS planning, e.g., [Hoffmann2000, [START_REF] Bryce | A Tutorial on Planning Graph Based Reachability Heuristics[END_REF]Kambhampati2007, Hoffmann andNebel2011].

Many techniques are used to solve hierarchical planning problems. Some ad-hoc HTN solvers have been implemented [Bercher, Keen, and Biundo2014, Erol1996, Nau et al.2003]. Another approach consists in encoding HTN problems into SAT problems [Schreiber et al.2019] or into con-straint programming problems [Vidal and Geffner2006].

One particular technique of encoding is to translate hierarchical planning problems as STRIPS problems. Encoding techniques benefit directly from the constant improvements of STRIPS planners. To our best knowledge, two HTN to STRIPS encodings have been published so far [Alford et al.2016,Alford, Kuter, andNau2009] (very recent work have been published concurrently to this paper, which will not be studied here [Behnke et al.2022]). However, they have some limitations on the type of problems they can address and only produce sequential plans. For instance, one of the best current encodings is [Alford et al.2016]. It translates any HTN problem into STRIPS, making it solvable by any STRIPS planner. However, this encoding has three downsides:

1. It depends on a progression bound, which is an integer parameter bounding the size of the task network, meaning that the maximum size of the task network has to be estimated before encoding the HTN problem into STRIPS.

2. The resulting solution plans are sequential, a single action being performed at each time. However, many real-world applications are intrinsically distributed, and need concurrent actions at each time.

3. The encoding generates an high number of actions in the translated STRIPS problem. This makes the actual grounding of the problem difficult. This is particularly true when the makespan increases, as the number of translated actions greatly expands with it.

The contribution of the paper is twofold: (1) we introduce a search procedure called CPFD (Concurrent Partial Forward Decomposition) to generate concurrent plans, and (2) propose an encoding of this search procedure into STRIPS actions.

The rest of this paper is as follows. Section 1 defines the problem statement. Section 2 presents the Concurrent Partial Forward Decomposition procedure (CPFD). Section 3 introduces the concrete STRIPS encoding of CPFD, called Concurrent Task Holders Decomposition encoding (CTHD). In the last section, we compare CTHD with HTN2STRIPS [Alford et al.2016], which is the current state-of-the-art encoding from HTN to STRIPS.

Problem statement STRIPS Planning Problems

A STRIPS planning problem is a tuple P = (L, A, I, G) where L is a finite set of logical propositions, A is a finite set of actions, I ⊆ L is the initial state, and G ⊆ L is the goal.

An A state s is a set of logical propositions. The result of applying an action a to state s is a state s ′ defined by the transition function s ′ = γ(s, a) = (s -del(a)) ∪ add(a) if pre(a) ⊆ s, and undefined otherwise. Let an action layer π be a set of pairwise independent actions, and pre(π) = a∈π pre(a). add(π) and del(π) are defined in the same way. By extension s ′ = γ(s, π) = (s -del(π)) ∪ add(π) if pre(π) ⊆ s, and undefined otherwise. Note that the actions of π can be executed concurrently or in any sequential permutation and still yield exactly the same state s ′ .

A layered plan Π is a sequence of action layers 1).

⟨π 1 , . . . , π n ⟩. Let γ(s, Π) = γ(γ(s, π 1 ), ⟨π 2 , . . . , π n ⟩). π i precedes π j if i < j. Likewise a i ≺ a j if a i ∈ π i , a j ∈ π j , and π i precedes π j . A layered plan Π is a solu- tion to a STRIPS planning problem P = (L, A, I, G) iff G ⊆ γ(s, Π) (see Fig.
In the following, a conditional action will be used, we used the semantic defined by the ADL formalism [Ped-nault1994]. As a regular action, a conditional one is defined by a set of preconditions. Its effects however depend on a set of conditions. For each condition verified, the corresponding effect is applied. This action is used for simplicity reasons but can easily be converted into a set of non condition actions.

HTN Planning Problems

We build on STRIPS planning problem definition to define a HTN planning problem as a tuple P = (L, T , A, M, I, tn) where L is a finite set of logical propositions, T is a finite set of tasks, A is a finite set of actions, M is a finite set of methods, I ⊆ L is the initial state and tn the initial task network. There are two kind of tasks: primitive tasks that can be resolved by a STRIPS action a = (task(a), pre(a), add(a), del(a)) ∈ A , and compound tasks, which can be recursively decomposed into either primitive or compound tasks by a method m ∈ M .

A task network is a tuple tn = (T, ≺, α) such that T is a finite set of tasks symbols, α : T → T maps indexes to tasks in T , and ≺ is a partial order over T representing precedence constraints:

t precedes t ′ if t ≺ t ′ , or equivalently (t, t ′ ) ∈≺. ≺ is transitive. A task α(t), t ∈ T is trailing if ∀t ′ ∈ T , (t ′ , t) /
∈≺ (t has no predecessor in T ). trail(tn) will denote the set of trailing tasks in T . Symmetrically, a task α(t), t ∈ T is a last task if ∀t ′ ∈ T, (t, t ′ ) / ∈≺ (t has no successor in T ).

A method is a tuple m = (task(m), pre(m), tn(m)) where task(m) is the compound task decomposed by the method m, pre(m) is the method's preconditions and tn(m) is a task network. A method m is a resolver of a compound task τ if task(m) = τ . Note that a given compound task can have various methods to resolve it:

task(m) = task(m ′ ) = • • • = τ .
An action a = (task(a), pre(a), add(a), del(a)) can be applied to resolve a primitive task α(t) of the initial task network tn if t is trailing, task(a) = α(t) and pre(a) ⊆ I. The result is a new problem P ′ = (L, T , M, I ′ , tn ′ ) where I ′ = γ(I, a) and tn ′ = (T \ {t}, {(t ′ , t ′′ ) ∈ ≺ | t ′ \ = t}, α\{(t, α(t))}). In a symmetrical manner, a method m = (task(m), pre(m), tn(m)) can be applied to resolve a compound task α(t) of the task network tn if t is trailing, task(m) = α(t) and pre(m) ⊆ I. The result of applying the method m with tn(m) = (T m , ≺ m , α m ) is a new problem P ′ = (L, T , M, I, tn ′ ) where tn ′ = (T ′ , ≺ ′ , α ′ ) and:

T ′ = (T \ {t}) ∪ T m ≺ ′ = {(t ′ , t ′′ ) ∈ ≺ | t ′′ ̸ = t} ∪ ≺ m ∪ {(t ′′ , t ′ ) ∈ T m × T | (t, t ′ ) ∈≺} α ′ = {(t ′ , α(t ′ )), t ′ ∈ T \{t}} ∪ α m
In other words, in ≺ ′ we keep all the precedence constraints of ≺ that does not involve t, add all the precedence constraints in ≺ m , and propagate precedence transitivity between ≺ and ≺ m through t.

Applying either an action a or a method m to resolve a task in a planning problem P is called a progression. If t p ∈ T is a primitive task of P , resolving t p by a is a progression denoted P → a tp P ′ . Similarly, if t c is a compound task of P , decomposing t p using m is a progression denoted P → m tc P ′ .

To conclude, a layered plan Π = ⟨π 1 , . . . , π n ⟩ is a solution for a HTN planning problem P = (L, T , A, M, I, tn) if (1) there exists a sequence of progressions that transforms P into P ′ = (L, T , M, I ′ , (∅, ≺ ′ )) (i.e. all the tasks of P have been resolved), and (2) a i ≺ a j ⇔ task(a i ) ≺ ′ task(a j ) (i.e. action precedence constraints in the layered plan Π are equivalent to the primitive task precedence constraints in P ′ ). In section 2, we show how CPFD builds a layered solution plan by applying progressions on P .

HTN to STRIPS Encoding Problems

The Hierarchical Task Network (HTN) formalism has been shown to be more expressive than STRIPS [Erol, Hendler, and Nau2003]. This means that any STRIPS problem can be formulated as a HTN problem but not the other way around. Therefore, the translation of a HTN problem into a STRIPS problem is not always possible. However, it has been proven by [Alford et al.2016] that this translation is possible if the size of the solution task network can be bounded. Given a HTN problem and a sequential solution plan, the minimum (respectively maximum) bound is the smallest (respectively largest) number of tasks in any task network visited by the sequence of progressions carried out to find this solution plan.

In practice, not all problems have a maximum bound, but all solvable problems have a minimum bound. These bounds are not directly related to the length of a problem solutions, though the minimum progression bound is smaller than the optimal plan length 1 .

Our encoding also assumes the bound existence. In addition, we make two other assumptions on the HTN problem to encode:

1. in the initial HTN problem, T is singleton. Otherwise, it is always possible to add a root dummy-task and a dummymethod to decompose it.

2. every method of the HTN problem has a task network with a single last task. If a method does not have it, a dummy-task with no successor is added to the task network.

3. methods have no preconditions. Otherwise, a trailing dummy action is added to the method task network with no effects and the method's preconditions.

These assumptions are made without loss of generality and will simplify the notations in the following. These assumptions were also made by the current state of the art encoding HTN2STRIPS [Alford et al.2016].

Concurrent Partial-order Forward Decomposition

In this section, we propose a recursive procedure to solve HTN problems called CPFD (Concurrent Partial-order Forward Decomposition) and generate layered plans. CPFD procedure is detailed in Alg. 1. CPDF is an adaptation of PFD (Partial-order Forward Decomposition) [Ghallab, Nau, and Traverso2004] procedure to output layered plans (see Figure 1). A layered plan is solution of a HTN problem if actions resolve all the tasks of the task network, and if the ordering constraints of the actions in the layered plan satisfy the precedence constraints in this task network. CPFD tries to solve recursively the trailing tasks as in the PFD procedure. The difference lies on the resolution of the primitive tasks: while PFD adds actions to a sequential plan, CPFD adds them to layers of independent actions. 1 For more details about the method to compute the bounds of the solution task network see [Alford et al.2016] More precisely, CPFD takes as input four parameters: a HTN problem P = (L, T , A, M, I, (T, ≺)), a layered plan, the index i of the current layer π i and τ the set of primitive tasks resolved by the actions in π i . The initial call of CPFD is CPFD(P, Π = [[]], i = 0, τ = ∅). At each recursive call, CPFD checks if the list of tasks T of the task network is empty, i.e., no more tasks have to be resolved. If this condition is satisfied, the layered plan Π is a solution to P and Π is returned. Otherwise, a task t ∈ T is non deterministically selected among the trailing tasks (tasks without predecessors with respect to precedence constraints), and a resolver is non deterministically chosen. As in the PFD procedure, there are two ways to resolve t depending on whether t is primitive or compound.

Case 1. (Primitive task) The resolvers of t are actions a whose preconditions are satisfied in the current state I and that are independent of all the actions already planned in the current layer π i . Two cases are possible: either t has no resolvers and the current layer π i is empty, meaning no action can solve t in the current state I, and CPFD returns FAILURE. Or t has a resolver but this resolver is not an independent action in the current layer: then CPFD moves to the next layer by applying to the current state all the actions already committed in the current layer. The idea is that t could be resolved by an action in a next state concurrently with other independent actions. Obviously, if t has an independent resolver a, a is added to the current layer π i and t is added to the set of resolved primitive tasks τ . Case 2. (Compound task) CPFD computes all the methods resolving the compound task t, i.e., whose preconditions are satisfied in the current state I. If there is no method, then t cannot be solved, and CPFD returns FAIL-URE. Otherwise, CPFD non deterministically chooses a method m decomposing t, update the task set and the ordering constraints accordingly.

CPFD(P, Π, i, τ ) is then called recursively until the tasks to solve in P are emptied (T = ∅, line 2) or a failure condition is met (line 9 and 24).

Theorem 1 Concurrent HTN is sound and complete.

Proof sketch (Soundness) All produced plans come from a progression of the initial task network, thus there is a sequence of task decomposition that produced the primitive tasks is the solution plan. Furthermore, since a primitive task can be added to a layer if the corresponding node is unconstrained, all tasks preceding the one added have been planned on previous layers. Thus the ordering constraints in ≺ are satisfied in the solution plan. Thus output plans are sounds.

Proof sketch (Completeness) We will show that CPFD is complete based on the demonstration that PFD is complete. Let P = (L, T , A, M, I, tn) be a HTN problem and Π = ⟨π 1 , . . . , π n ⟩ a layered solution plan of P . Let us show that there is a sequence of recursive calls of CPFD outputting Π. First, let us note that given a concurrent Algorithm 1: CPFD(P, Π, i, τ )

1 {P = (L, T , A, M, I, (T, ≺, α)) is the current problem} 2 if T = ∅ then return Π 3 toSolve ← trail(tn) \ τ 4 π i ← get(Π, i) 5 nondeterministically choose t ∈ toSolve 6 if t is primitive then 7 resolvers ← {a ∈ A | task(a) = t,
≺ ′ = {(t ′ , t ′′ ) ∈ ≺ | t ′′ ̸ = t} ∪ ≺ m ∪ 17 {(t", t') ∈ T m × T | (t, t ′ ) ∈≺}
≺← {(t ′ , t ′′ ) ∈ ≺ | t ′′ ̸ = t} ∪ ≺ m ∪ 33 {(t ′′ , t ′ ) ∈ T m × T | (t, t ′ ) ∈≺} 34 end 35 return CPFD(P, Π, i, τ ) layer π = {a 1 , .
. . , a k }, any linearization of that layer (⟨a γ(1) , . . . , a γ(k) ⟩ where γ is a permutation function of {1, . . . , k}) is a sequence of actions which can be applied to the same states as π. This is due to the mutual independence property of the actions within a concurrent layer. From there, any sequential plan produced by linearizing every layer of Π (by taking any permutation function on the layers) is a sound plan that also solves P . Let us consider the linearization Π l defined by the n permutation functions γ 1 , . . . , γ n . Since PFD is a complete algorithm, there is a sequence of recursion of PFD which outputs Π l . We will show that there is an analogous sequence of CPFD recursions outputting Π. At each recursion, PFD and CPFD either solve an unconstrained abstract task, or an unconstrained primitive task.

While they solve abstract tasks the same way, PFD solves a primitive task by a adding an action resolving it to the head of the plan, meanwhile CPFD adds the action resolving the task to the concurrent layer at the head of the plan. If the task can not be resolved, PFD returns a Failure while CPFD tries to add a new concurrent layer to the plan. Thus, for each recursive PFD call resolving an abstract task, the analogous call of CPFD is to solve the same abstract task. Each recursive call of PFD resolving a primitive task is analogous to a CPFD call adding the action to the current concurrent layer. However, CPFD requires extra recursive calls compared to PFD: it needs to select and try to resolve an unsolvable task after each layer. In conclusion, the analogous sequence of recursion of CPFD is the one solving the same abstract and primitive tasks than PFD but with the insertion of recursive calls trying to solve an unsovlable task after resolving the last action of a layer of Π. Thus there is a sequence of CPFD recursions outputting Π and CPFD is complete.

Planning the planning: Translating HTN to STRIPS

We present in this section the concrete STRIPS encoding of the compound algorithm previously presented, called Concurrent Task Holders Decomposition encoding (CTHD).

Taskholder Encoding

The CPFD procedure described previously resolves recursively unconstrained primitive tasks and compound tasks by modifying the initial task network of the problem until the task network contains only an empty set of tasks. To encode this process, we need first to model a task network with STRIPS. To achieve this, we use the concept of taskholder introduced first by [Alford et al.2016]. A taskholder is a STRIPS object that will act as a container for a task. By way of extension, a task network is modeled as a stack of taskholders: static ordering relationship between taskholders defines in which order the taskholders can be allocated to the tasks during the CPFD procedure and fluent relationship define the ordering constraints between the tasks of the task network. The number of taskholders should be fixed before translating the HTN problem. Similarly to HTN2STRIPS, CTHD requires at least as many taskholders as there are tasks in the largest explored task network, thus the number of taskholders in CTHD can be estimated the same way as HTN2STRIPS estimates the number of taskholders.

In addition, we order the taskholders into a stack which defines the order in which the taskholders can be used.

The current layer of the layered plan under construction is represented by a set of propositions, each proposition denoting the fact that an action a ∈ A is planned in the current layer or not.

To fix ideas, let us consider the planning problem defined as a task decomposition graph (see Figure 2), a graphical representation of the initial structures of this problem, i.e., taskholders and current layer, is given Figure 6.

Encoding as Actions

The dynamics of the CPFD procedure (Alg . 1) is defined by three types of STRIPS actions: (1) the first type of ac- tions resolves an unconstrained compound task and update the current task network according to a method decomposition (2) the second type of actions resolves a primitive task and add an action to the current layer; and (3) the last type of actions is the one switching layer, making the algorithm build the next layer of the solution plan. The planning process will choose one action among the three types. Applying one of this action is equivalent to one recursive call of CPFD procedure. The planning process ends when there is no tasks left in the task network, i.e., when all taskholders are empty. These actions are defined as follows:

1. Actions for resolving compound tasks: These actions reflect the decomposition of a compound task into subtasks according to a method. It corresponds to the lines 27 to 30 of CPFD procedure. In order to apply these actions, the following must be verified:

• The taskholder containing the decomposed task is unconstrained. • There is enough taskholders remaining in the stack. It works on our example as represented Figure 3. The substasks of the method M 1 decomposing T 0 are added to the stack of taskholders and the task T 0 is replaced by the last task of the current task network. As M 1 as no last task, the no-op action is used. Finally, the ordering constraints are set over the taskholders.

2. Actions for resolving primitive tasks: These actions resolve an unconstrained primitive task of the task network and add the resolver action a ∈ A into the current layer (lines 21 to 23 of Alg. 1). In CPFD procedure, it translates to adding the action resolving a primitive task contained in an unconstrained taskholder to the current layer, removing the primitive task from the taskholder and marking the taskholder as resolved. These actions can be executed when the following conditions are verified:

• The preconditions of the action are satisfied.

• The taskholder containing the task is unconstrained.

• All actions already planned at the plan step are independent with the task.

An example of application is displayed on Figure 4. The taskholder th2 is unconstrained and contains a primitive task t1 that can be resolved by the action a(t1). The task is resolved by adding a(t1) to the current layer, th2 is emptied and marked as resolved. As in Alg. 1, the constraints implied by the resolved taskholder are not removed yet, they will when going to the next layer. 3. Action for switching of plan layer: This CTHD action corresponds to switching to the next layer. This action empty the current layer, setting up the next one in the solution plan. It also removes the constraint implying the resolved taskholders. This action can be applied in any situation and works as displayed on Figure 5. In this example, only th2 is marked as resolved. After the application of the action, the constraints implying th2 are removed, so the constraint between th2 and th1 is deleted, additionally th2 is unmarked as resolved and set as empty.

Then the current layer is emptied by removing all actions in it.

Concurrent Taskholder Decomposition Encoding

In this section we present the STRIPS formulation of the translated problem. Similarly to HTN2STRIPS, our translation depends on an integer parameter denoted b representing the number of taskholders (i.e. the maximum number of tasks in the task network). We first define the predicates used to encode the problem, then we define the encoding of the three types of actions presented above. We end by presenting the initial state and goal state of the translated problem.

Let P = (L, T , A, M, I, tn) be a HTN problem and let b ∈ N the number of taskholders. CTHD encoding generates a STRIPS problem

CP F D(P, b) = (L ∪ L ′ , A ′ , I ∪ I ′ , G ′ ).
The encoding generates the set propositions L ′ based on the following predicates: • (not constraint ?th1 ?th2 -taskholder) represents the fluent constraints over the taskholders. The predicate is inverted for convenience, so when the proposition (not constraint th1 th2) is false, it means that the task contained in th1 must be planned before the task contained in th2. • (empty ?th -taskholder) represents the fact that a taskholder is empty. The proposition (empty th) is true if the taskholder th does not contain a task. • (in ?t -task ?th -taskholder) is the predicate representing whether or not the task ?t ∈ T is set in the parameter taskholder. The proposition (in t th) is true if t is set in th.

• (not planned ?a -action) is true if a is not planned in the current plan step. • (prec th ?th1 ?th2 -taskholder) is a predicate defining the static relationship between the taskholders and defines the taskholders stack. So the proposition (prec th th1 th2) is true is th1 is above th2 in the stack.

In the following, this order will be fixed and ∀0 ≤ i, j < b the proposition (prec th th i th j ) will be true if and only if i ≤ j. • (resolved ?th -taskholder) is a predicate representing either or not a taskholder have been resolved. So the proposition (resolved th) is true if the taskholder th have been resolved.

The encoding generates the set of actions A ′ = A c ∪ A p ∪ A l where A c is the set of actions for resolving compound tasks, A p the set of actions for resolving primitive tasks and A l the set of actions for switching of plan layer: • Actions for resolving compound tasks:

Let m = (task(m), pre(m), tn(m)) and (task 1 , task 2 , . . . , task k ) the subtasks in tn. We assume that task k is the last task of tn. For all methods m ∈ M there is an action a m ∈ A c with k parameters (?th 1 , ?th 2 , . . . , ?th k ) defined as follows:

-pre(a m ) = (in task(m) ?th 1 ) ∧ b-1 i=0 (not constraint th i ?th 1 ) ∧ k i=2 (prec th ?th i ?th i+1 ) ∧ k i=2 (empty ?th i ) -add(a m ) = k i=2 (in task i ?th i ) ∧ (in task k ?th 1 ) -del(a m ) = k i=2 (empty ?th i ) ∧ taski≺taskj (not constraint ?th i ?th j ) ∧ k i=2 (not constraint ?th i ?th 0 )
Note that the k -1 last taskholder parameters are required to be ordered according to the static relationship defined by the prec th predicate. For instance, if three new taskholders are required to decompose the task in th6, (th6 th2 th4 th7) is a valid combination of parameters, while (th6 th2 th5 th3) is not.

• Actions for resolving primitive tasks:

For all actions a ∈ A there is an action a p ∈ A p with one parameter: a taskholder containing p and denoted ?th.

The action is defined as follows:

pre(a p ) = pre(a) ∧ (in task(a) ?th) ∧ b-1 i=0 (not constraint th i ?th) ∧ t∈nInd(p) (not planned ?a) add(a p ) = add(a) ∧ (empty ?th) ∧ (resolved ?th) del(a p ) = del(a) ∧ (not planned ?a)

• Action for switching layer:

A l is composed of one conditional action a l with no parameter. This action has a non conditional part: emptying the current layer, and a conditional part: unconstraining and making available for a new use the resolved taskholders. It is defined as follows:

pre(a l ) = ∅ add(a l ) = a∈A (not planned ?a) ∧ ∀(?th -taskholder) when (resolved ?th),

b-1 i=0 (not constraint ?th th i ) del(a l ) = ∅ Finally, let us define the encoding for the initial state and goal of the translated problem. The initial state is defined by setting the initial task into the first taskholder. The remaining taskholders are set as empty and ordered in a stack. As there is no constraint over the taskholder yet, all constraints predicate are initialized accordingly.

I ′ = (in task 0 th 0 ) ∧ 1≤i<b (empty th i ) ∧ 1≤i<b-1 i<j<b (prec th th i th j ) ∧ 0≤i,j<b (not constraint th i th j ) ∧ a∈A (not planned a)
The problem is solved when all taskholders are empty, meaning that all task are planned. The goal state is defined accordingly.

G ′ = p-1 i=0 (empty th i )

Differences and Similarities with HTN2STRIPS

If both CTHD and HTN2STRIPS encode HTN solving with taskholders, they aim at finding different kinds of solution plans. While CTHD aims at finding concurrent plans, HTN2STRIPS can only produce sequential ones. In the following, we will experimentally compare the two encodings but one has to keep in mind that CTHD produces plans with higher expressivity than the one produced by HTN2STRIPS.

In addition, while HTN2STRIPS taskholders are unordered, CTHD orders the taskholders into a static stack and imposes an order in which the taskholders parameters are used. The purpose of the stack is to reduce the number of valid operators in the translated problem: let us consider an example with four taskholders and a task T 1 set in the first taskholder th1. T 1 can be decomposed by a method m into four subtasks t 1 , t 2 , t 3 and t 4 . The HTN2STRIPS translation, translates this method into 3! (equivalent) actions (one for each permutation of the three newly used taskholders). Meanwhile CTHD, by constraining the taskholder a fixed order, translates this method into a single action which corresponds to the HTN2STRIPS action where all taskholders are ordered. In the general case, the number of translated actions generated from a method by HTN2STRIPS is factorial over the number of subtasks in the method, while it is quadratic for CTHD.

Experimentation and Results

In this section we will present the experiments and results we used to demonstrate CTHD efficiency. We compare CTHD to the current best HTN to STRIPS translation HTN2STRIPS [Alford et al.2016]. This comparison will be made over IPC HTN benchmarks domains. These domains were chosen because they can express either totally ordered or partially ordered problems. Over each domain, the problems will be divided into totally ordered and partially ordered problems. While the solution plan of a totally ordered problem are necessarily sequential, the solution of a partially ordered problem can be concurrent.

Both HTN2STRIPS and CTHD can use the same progression bound to solve HTN problems. So in order to have a fair comparison, all problems were solved using the same (minimal) progression bound for both encodings. Note however that if there is no benefit for HTN2STRIPS to use a progression bound bigger that the minimal one, CTHD can explore bigger task networks and thus find solution plans with higher concurrency.

The comparison of HTN2STRIPS and CTHD will be made over five metrics:

• Solving time: This include the time spent to instantiate the problem and domain file and the time spent solving the instantiated problem.

• Search time: This corresponds solely to the time spent to solve the instantiated problem.

• Solution makespan: This corresponds to the "length" of the solution plan, meaning the number of actions within the plan for a sequential solution plan or the number of concurrent layers for a concurrent solution plan.

• Number or proposition: This corresponds to the number of proposition generated by instantiating the domain and problem file.

• Number of operators: This corresponds to the number of operators generated by instantiating the domain and problem file.

These five parameters will be evaluated by the IPC scoring metric, which is defined as follows:

IP C(k) = 1 |problems| i∈problems min p∈planners (cost(p)) cost(k)
We ran all experiments on a single core of a Intel Core i7-9850H CPU, using the Fast Downward library [Helmert2011] with the Delfi 1 configuration [Katz et al.2018], With a limit of 8GB of RAM over 600 seconds.

The results of the experiments are displayed on Figure 7. Over the three domains, CTHD consistently obtains a better score on the operator metric. This is due to the difference in scaling between the number of operators in the translated HTN2STRIPS problem and the CTHD one. As stated before, as the number of operators is factorial in HTN2STRIPS, it is quadratic in CTHD.

Proposition wise, both encodings are very similar. The discrepancies between domains can be explained by the number of concurrent actions within the domain: the more concurrency possible between actions, the more proposition are generated by HTN2STRIPS.

For the Makespan metric, on totally ordered problems HTN2STRIPS and CTHD obtained the same maximal score. This is due to the fact that since both encodings use the same minimal progression bound, they both find the shortest When it comes to Search time metric, the results are more mixed. Overall, HTN2STRIPS has a better Search time for the partially ordered problems but on the totally ordered ones, HTN2STRIPS outperforms CTHD only on the Blocksworld domain. These mixed results can be explained by the fact that a planner requires less steps to solve a HTN2STRIPS translated problem than a CTHD translated one: since HTN2STRIPS only produce sequential plans, going to the next layer is implied by resolving a primitive task. On the other hand, going to the next layer is a full solving step in a CTHD translated problem. However, we also saw that HTN2STRIPS generates more operators than CTHD, leading to higher branching factor for the planner solving the translated problem. In the end, there is a trade off between finding a shorter plan with a higher branching factor or finding larger plans with a lower branching factor.

Finally, the Solving time metric represents the total time spent by the planner to go from the parsing of the translated files to the solution plan. It is the sum of the time spent instantiating the files and the Search time. On this metric, CTHD obtains a better score than HTN2STRIPS on all domains except for the partially ordered Satellite one. However on all domains, CTHD improves its score going from Search time to Solving time. Once again, since CTHD requires much less operators, instantiating a CTHD translated domain file is much faster than a HTN2STRIPS one.

Overall, CTHD outperforms HTN2STRIPS both regarding the time spent to solve the translated problems and regarding the makespan of the solution plans.

Conclusion and future work

In this paper, we presented a new HTN procedure, CPFD, solving HTN problem with layered solution plan. We encoded this procedure as a STRIPS problem, thus producing a new HTN to STRIPS translation, CTHD. The translated problem, can be solved by any STRIPS planner. Then we showed experimentally that our translation outperforms the current best one HTN2STRIPS, both in terms of problem representation size, solving efficiency and quality of solution plans. This translation is a new way to solve HTN problems in a concurrent way, offering a new alternative to plan-space HTN algorithms. However, similarly to HTN2STRIPS, our translation still depends on an integer parameter the progression bound, which bounds the size of the explored task networks. In order to improve HTN to STRIPS translation, we feel that this last point is the one to focus on.

  action a is a triplet a = (pre(a), add(a), del(a)) where pre(a) is the action's preconditions, add(a) is its positive effects and del(a) its negative ones, each a set of propositions. Two actions (a, b) are independent iff del(a) ∩ (pre(b) ∪ add(b)) = ∅ and del(b) ∩ (pre(a) ∪ add(a)) = ∅. Note that action independence only depends on action definitions. In the following, for all a ∈ A, nInd(a) will denote the set of actions b ∈ A dependent with a.

Figure 1 :

 1 Figure 1: A layered plan with the successive states resulting from the action layer application. Circles represent states, and rectangles action layers.

T

  pre(a) ⊆ I and (∀b ∈ π i , a independent of b)} ← T \ τ // Update the task network 16

  π i ∪ {a} // Add a to the current layer 23 τ ← τ ∪ {t} // Add t to the resolved tasks set ∈ M | task(m) = t} 28 if resolvers = ∅ then return Failure 29 nondeterministically choose m ∈ resolvers 30 {m = (T m , ≺ m )} 31 T ← (T \ {t}) ∪ T m // Decomposing tn with m 32

Figure 2 :

 2 Figure2: The task decomposition graph of T 0. In this graph, each task node (represented with circles) is linked to the method resolving it (represented with squares). For instance, T 1 can be decomposed into t2 by applying M 2 or into t3 by applying M 3. The ordering constraints are represented with the dotted arrows, so when decomposing T 0 with M 1, three subtasks are generated, T 1, t1 and T 2 where t1 must be planned before T 1.

Figure 3 :

 3 Figure 3: T0 is decomposed into four tasks, two compound ones T1 and T2, and two primitive ones t1 an no-op. Since neither T1, T2 or t1 is a last task, a no-op action is inserted instead of T0. The constraint over the taskholders are represented on the top right graph: each directed edge represents a precedence constraint. So for instance, the task in th2 should planned before the one in th1.

Figure 4 :

 4 Figure 4: The second taskholder was unconstrained, and contains a primitive task. It is added to the plan step by removing the task from the taskholder, adding the action a(t1) resolving t1 to the plan step, and marking the taskholder as resolved (represented by the black dot).

Figure 5 :

 5 Figure 5: The plan step is terminated by emptying it, the constraints implied by the resolved taskholders are removed.

Figure 6 :

 6 Figure 6: The problem is initialized by setting the initial task into the first taskholder.

Figure 7 :

 7 Figure 7: Results for the Satellite, Miconic and Blocksworld domains: TO (Totally ordered) and PO (Partially ordered)


Nk Ko
email: vladimir.solozhenko@univ-paris13.fr

Keywords: boron subarsenide, hardness, bulk modulus

Vickers hardness of boron subarsenide B 12 As 2 has been predicted using three modern theoretical models and experimentally studied by microindentation. The polycrystalline material exhibits hardness of about 31 GPa, and hence B 12 As 2 belongs to a family of (super)hard phases.

Boron subarsenide B 12 As 2 is icosahedral boron-rich compound possessing a fortunate combination of properties (wide energy band gap, high carrier mobility, extraordinary radiation sustainability, etc.) [START_REF] Emin | Unusual properties of icosahedral boron-rich solids[END_REF][START_REF] Slack | Some crystallography, chemistry, physics, and thermodynamics of B 12 O 2 , B 12 P 2 , B 12 As 2 , and related alpha-boron type crystals[END_REF][START_REF] Li | Fundamental studies on the structures and properties of some B 12 -based crystals[END_REF] that makes it a promising material for electronic devices operating under extreme conditions. In the present Letter we report the hardness and other mechanical properties of B 12 As 2 from theoretical and experimental studies.

Vickers hardness (H V ) of boron subarsenide was predicted using three contemporary theoretical models of hardness, i.e. the thermodynamic [START_REF] Mukhanov | The interrelation between hardness and compressibility of substances and their structure and thermodynamic properties[END_REF], the Lyakhov-Oganov [START_REF] Lyakhov | Evolutionary search for superhard materials: Methodology and applications to forms of carbon and TiO 2[END_REF] and the Mazhnik-Oganov [START_REF] Mazhnik | A model of hardness and fracture toughness of solids[END_REF] models. The general thermodynamic model based on crystal structure and thermodynamic properties was used in the formulation developed for the particular case of boron-rich solids [START_REF] Mukhanov | Thermodynamic model of hardness: Particular case of boron-rich solids[END_REF], while within Lyakhov-Oganov approach the strength of covalent bonding, degree of ionicity and topology of the crystal structure were considered. Mazhnik-Oganov model based on elastic constants was also used for evaluation of fracture toughness (K Ic ). The results are summarized in Table . It should be noted that bulk modulus of B 12 As 2 estimated in framework of the thermodynamic model is in good agreement with the experimental value 150(4) GPa [START_REF] Cherednichenko | Equation of state of boron subarsenide B 12 As 2 to 47 GPa[END_REF].

Stoichiometric B 12 As 2 has been synthesized by the method described elsewhere [START_REF] Cherednichenko | Equation of state of boron subarsenide B 12 As 2 to 47 GPa[END_REF]. Polycrystalline bulks for hardness measurements have been produced in a toroid-type apparatus by crystallization from melt at 5.2 GPa. According to X-ray diffraction study (TEXT 3000 Inel) the recovered bulks contain well-crystallized single-phase B 12 As 2 with lattice parameters identical to the literature data [9].

The recovered samples were planar ground with diamond 800 grit and subsequently polished with 9-μm and 1-μm diamond suspensions. Final vibropolishing with 0.04-μm SiO 2 colloidal solution ensured the minimal sample surface damage.

According to scanning electron microscopy study (TESCAN MIRA3 FEG) all samples were homogeneous without residual porosity, and their grain size varied from 0.5 to 3 µm. Vickers hardness measurements have been performed using a Mitutoyo HM-220B Microhardness Testing Machine under loads from 1 to 20 N and 15 seconds dwell time; five indentations were made at each load. Hardness and indentation fracture toughness values were calculated following the standard definitions from the residual imprints upon indentation (for details see [START_REF] Solozhenko | Mechanical properties of boron phosphides[END_REF]).

The measured Vickers hardness of boron subarsenide decreases with the load, and at 6 N reaches the asymptotic value H V = 31(2) GPa that is in perfect agreement with the hardness calculated by the thermodynamic and Mazhnik-Oganov models (see Table ). H V value calculated using Lyakhov-Oganov model is slightly overestimated. Indentation fracture toughness is almost independent of the load, and the experimental value of K Ic = 1.1(3) MPa•m 0.5 is in agreement with the calculated one.

Thus, boron subarsenide B 12 As 2 exhibits experimental Vickers hardness of 31(2) GPa which is higher than that of cubic silicon carbide SiC, a conventional superabrasive. All used theoretical models work properly, however, the thermodynamic model seems to be more reliable in the case of boron-rich solids.

The author thanks Dr. V.A. Mukhanov for fabrication of bulk samples, and Dr. V. Bushlya for sample polishing and SEM study. b Lyakhov-Oganov model [START_REF] Lyakhov | Evolutionary search for superhard materials: Methodology and applications to forms of carbon and TiO 2[END_REF].

H H

V V T T a a L LO O b b M MO O c c e ex xp p. . B B 0 0 G G v v E E   K K I Ic c G GP Pa a - --- M MP Pa a• •m m ½ ½ 3 32
c Mazhnik-Oganov model [START_REF] Mazhnik | A model of hardness and fracture toughness of solids[END_REF].

d Calculated from elastic moduli C ij [START_REF] Lee | Elastic moduli of B 12 and its compounds[END_REF] using Voigt's approach [START_REF] Voigt | Über die Beziehung zwischen den beiden Elasticitätsconstanten isotroper Körper[END_REF].

e Calculated from B 0 and G v values using isotropic approximation.

a

  Thermodynamic model (particular case of boron-rich solids)[START_REF] Mukhanov | The interrelation between hardness and compressibility of substances and their structure and thermodynamic properties[END_REF][START_REF] Mukhanov | Thermodynamic model of hardness: Particular case of boron-rich solids[END_REF].

Table

  Mechanical properties of boron subarsenide: Vickers hardness (H V ), bulk modulus (B 0 ), shear modulus (G v ), Young's modulus (E), Poisson's ratio () and fracture toughness (K Ic ).


A Rakotomamonjy
email: alain.rakotomamonjy@insa-rouen.fr

Litis Ea

Algorithms for Multiple Basis Pursuit Denoising

We address the problem of learning a joint sparse approximation of several signals over a dictionary. We pose the problem as a matrix approximation problem with a row-sparsity inducing penalization on the coefficient matrix. We propose a simple algorithm based on iterative shrinking for solving the problem. At the present time, such a problem is solved either by using a Second-Order Cone programming or by means of a M-Focuss algorithm. While the former algorithm is computationally expensive, the latter is efficient but present some pitfalls like presences of fixed points which are undesiderable when solving a convex problem. By analyzing the optimality conditions of the problem, we derive a simple algorithm.

The algorithm we propose is efficient and is guaranteed to converge to the optimal solution, up to a given tolerance. Furthermore, by means of a reweighted scheme, we are able to improve the sparsity of the solution.

I. INTRODUCTION

In this paper, we consider the problem of simultaneous sparse approximation which can be stated as follows. Suppose that we have measured L signals {s i } L i=1 where each signal is of the form

s i = Φc i + ǫ
where s i ∈ R N , Φ ∈ R N ×M is a matrix of unit-norm elementary functions, c i ∈ R M a weighting vector and ǫ is a noise vector. Φ will be denoted in the sequel as the dictionary matrix. Since we have several signals, the overall measurements can be written as

S = ΦC + E (1)
with

S = [s 1 s 2 • • • s L ] a signal matrix, C = [c 1 c 2 • • • c L ]
and E a noise matrix. Note that in the sequel, we have adopted the following notations. c i,• and c •,j respectively denote the ith row and jth column of matrix C. c i,j is the ith element in the jth column of C.

For the sparse simultaneous approximation problem, the goal is then to recover the matrix C given the signal matrix S and the dictionary Φ under the hypothesis that all signals s i share the same sparsity profile. Such a problem can be formulated as the following optimization problem

min C 1 2 S -ΦC 2 F + λJ p,q (C) (2)
where λ is a user-defined parameter that balances the approximation error and the sparsity-inducing penalty J p,q (C) is

J p,q (C) = i c i,• p q
with typically p ≤ 1 and q ≥ 1. J p,q (C can be interpreted as a relaxed version of a ℓ 0 quasi-norm on the row-sparsity measure of C. Several authors have proposed methods for solving problem [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF]. For instance, Cotter el al. [START_REF] Cotter | Sparse solutions to linear inverse problems with multiple measurement vectors[END_REF] developed an algorithm for solving problem (2) when p ≤ 1 and q = 2, known as M-FOCUSS. Such an algorithm based on factored gradient descent have been proven to converge towards a local or avglobal (when p = 1) minimum of problem (2) if it does not get stuck in a fixed-point. The case p = 1, q = 2, named as M-BP for Multiple Basis Pursuit in the following, is the most natural extension of the so-called Lasso problem [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] or Basis Pursuit Denoising [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], since for L = 1, problem (2) reduced to the Lasso problem. The key point of this case is that it yields to a convex optimization problem and thus it can benefit from all properties resulting from convexity e.g global minimum. Malioutov et al. [START_REF] Malioutov | Sparse signal reconstruction perspective for source localization with sensor arrays[END_REF] have proposed an algorithm based on a second-order cone programming formulation, which at the contrary to M-FOCUSS, always converges to the problem global solution.

In this paper, we develop a simple and efficient algorithm for solving the M-Basis Pursuit problem. We show that by using results from non-smooth optimization theory, we are able to propose an efficient iterative method which only needs some matrix multiplications.

Afterwards, instead of directly deriving a proper algorithm for solving the non-convex optimization problem when p < 1 and q = 2, we introduce an iterative reweighted M-Basis pursuit (IrM-BP) algorithm. We then show that depending on the chosen weights, such an iterative scheme can actually solve problem (2). Our main contribution at this point is then to have translated the non-convex problem (2) into a series of convex problems which are easy to solve with our iterative method for M-BP.

II. ALGORITHMS

This section describes two algorithms for solving the M-BP problem with convex and non-convex penalties J p,q .

A. Iterative shrinking

The M-BP optimization problem is the following

min C W (C) = 1 2 S -ΦC 2 F + λ i c i,• 2 (3)
where the objective function W (C) is a non-smooth but convex function. Since the problem is unconstrained a necessary and sufficient condition for a matrix C ⋆ to be a minimizer

Algorithm 1 Solving M-BP through iterative shrinking C = 0, Loop = 1 while Loop do for i = 1, 2, • • • , M do if c i,• KKT condition is not satisfied then c i,• = 1 -λ Ti + T i
end if end for if all KKT Conditions are satisfied then Loop = 0 end if end while of (3) is that 0 ∈ ∂W (C ⋆ ) where ∂W (C) denotes the subdifferential of our objective value W (C) [START_REF] Bertsekas | Convex Analysis and Optimization[END_REF]. By computing the subdifferential of W (C) with respect to each row c i,• of C, the optimality condition of problem (3) is then

-r i + λg i,• = 0 ∀i where r i = φ t i (S -ΦC) and g i,• is the i-th row of a subdif- ferential matrix G of J 1,2 (C) = i c i,• 2 .
According to this definition of J 1,2 's subdifferential, the optimality condition can be rewritten as

-r i + λ c i,• c i,• 2 = 0 ∀i, c i,• = 0 (4)
r i 2 ≤ λ ∀i, c i,• = 0
A matrix C satisfying these equations can be obtained after the following algebra. Let us expand each r i so that

r i = φ t i (S -ΦC -i ) -φ t i φ i c i,• = T i -c i,• (5)
where C -i is the matrix C with the i-th row being set to 0 and T i = φ t i (S -ΦC -i ). The second equality is obtained by remembering that φ t i φ i =1. Then, equation (4) tells us that if c i,• is non-zero, T i and c i,• have to be collinear. Plugging all these points into equation (4) yields to an optimal solution that can be obtained as :

c i,• = 1 - λ T i + T i ∀i (6)
From this update equation, we can derive a simple algorithm which consists in iteratively applying the update (6) to each row of C. Such an iterative scheme actually performs a blockcoordinate optimization. Although, block-coordinate optimization does not converge in general for non-smooth optimization problem, Tseng [START_REF] Tseng | Convergence of block coordinate descent method for nondifferentiable minimization[END_REF] has shown that for an optimization problem which objective value is the sum of a smooth and convex function and a non-smooth but block-separable convex function, block-coordinate optimization converges towards the global minimum of the problem. Since for M-BP we are considering a quadratic function and a row-separable penalty function, Tseng's results can be directly applied in order to prove convergence of our algorithm.

Our approach, detailed in Algorithm (1), is a simple and efficient algorithm for solving M-BP especially when the dictionary size is large. A similar approach has also been proposed for solving the lasso [START_REF] Friedman | Pathwise coordinate optimization[END_REF], the group lasso [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] and the elastic net [START_REF] Zhou | Regularization and variable selection via the elastic net[END_REF]. Intuitively, we can understand this algorithm as an algorithm which tends to shrink to zero rows of the coefficient matrix that contribute poorly to the approximation. Indeed, T i can be interpreted as the correlation between the residual when row i has been removed and φ i . Hence the smaller the norm of T i is, the less φ i is relevant in the approximation. And according to equation ( 6), the smaller the resulting c i,• is. Insight into this iterative shrinking algorithm can be further obtained by supposing that M ≤ N and that Φ is composed of orthonormal elements of R N , hence Φ t Φ = I. In such situation, we have

T i = φ t i S and T i 2 2 = L k=1 (φ t i s k ) 2
and thus

c i,• =   1 - λ L k (φ t i s k ) 2   + φ t i S
This last equation highlights the relation between the single Basis Pursuit (when L = 1) and the Multiple-Basis Pursuit algorithm presented here. Both algorithms lead to a shrinkage of the coefficient projection. With the inclusion of multiple signals, the shrinking factor becomes more robust to noise since it depends on the correlation of the atom φ i to all signals.

B. Reweighted iterative shrinking algorithm

This subsection introduces an iterative reweighted M-Basis Pursuit (IrM-BP) algorithm which solves problem (2) when p < 1 and q = 2.

Recently, several works have advocated that sparse approximations can be recovered through iterative algorithms based on a reweighted ℓ 1 minimization [START_REF] Zou | The adaptive lasso and its oracle properties[END_REF], [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF], [START_REF] Chartrand | Iteratively reweighted algorithms for compressive sensing[END_REF]. Typically, for a single signal case, the idea consists in iteratively solving the following problem

min c 1 2 s -Φc 2 2 + λ i z i |c i |
where z i are some positive weights, and then to update the positive weights z i according to the solution c ⋆ . Besides, providing empirical evidences that reweighted ℓ 1 minimization yields to sparser solutions than a simple ℓ 1 minimization, the above cited works theoretically support such claims. These results for the single signal approximation case suggest that in the simultaneous sparse approximation problem, reweighted M-Basis Pursuit would lead to sparser solutions than the classical M-Basis Pursuit.

Our iterative reweighted M-Basis Pursuit is defined as follows. We iteratively solve until convergence the optimization problem

min C 1 2 S -ΦC 2 F + λ i z i c i,• 2 (7)
where the positive weight vector z depends on the previous iterate C (n-1) . In our case, we will consider the following weighting scheme

z i = 1 ( c (n-1) i,• 2 + ε) r ∀i (8)
where {c

(n-1) i,•
} is the i-th row of C (n-1) , r a user-defined positive constant and ε a small regularization term that avoids numerical instabilities and prevents from having an infinite regularization term for c i,• as soon as c

(n-1) i,•
vanishes. This is a classical trick that has been used for instance by Candès et al. [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF] . Note that for any positive weight vector z, problem ( 7) is a convex problem that does not present local minima. Furthermore, it can be solved using our iterative shrinking algorithm by simply replacing λ with λ i = λ • z i . Such a scheme is similar to the adaptive lasso algorithm of Zou et al. [START_REF] Zou | The adaptive lasso and its oracle properties[END_REF] but uses several iterations and addresses the simultaneous approximation problem.

The IrM-BP algorithm we proposed above can also be interpreted as an algorithm for solving problem (2) when 0 < p < 1. Indeed, similarly to the reweighted ℓ 1 scheme of Candès et al. [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF] or the one-step reweighted lasso of Zou et al. [START_REF] Zou | One-step sparse estimates in nonconcave penalized likelihood models[END_REF], our algorithm falls in the class of majorizeminimize (MM) algorithms [START_REF] Hunter | A tutorial on MM algorithms[END_REF]. MM algorithms consists in replacing a difficult optimization problem with a easier one, for instance by linearizing the objective function, by solving the resulting optimization problem and by iterating such a procedure.

The connection between MM algorithms and our reweighted scheme can be made through linearization. In effect, in our case, since J p,2 is concave in c i,• for 0 < p < 1, a linear approximation of J p,2 (C) around C (n-1) yields to the following majorizing inequality

J p,2 (C) ≤ J p,2 (C (n-1) )+ i p c (n-1) i,• 1-p 2 ( c i,• -c (n-1) i,•
) then for the minimization step, replacing in problem (2) J p,2 with the above inequality and dropping constant terms lead to our optimization problem [START_REF] Hunter | A tutorial on MM algorithms[END_REF] with appropriately chosen z i and r. Note that for the weights given in equation ( 8), r = 1 corresponds to the linearization of a log penalty i log( c i,• ) whereas setting r = 1p corresponds to a ℓ p penalty (0 < p < 1). According to the convergence properties for MM algorithms towards a local minimum of their objective function [START_REF] Hunter | A tutorial on MM algorithms[END_REF], we can state that our IrM-BP algorithm converges towards a local minimum of problem (2) with p and r being appropriately related.

III. RESULTS

In order to quantify the performance of our algorithms and compare them to other approaches, we have used simulated datasets with different redundancies M N , number k of active elements and number L of signals to approximate. The dictionary is based on M vectors sampled from the unit hypersphere of R N . The true coefficient matrix C has been obtained as follows. The positions of the k non-zero rows in the matrix are randomly drawn. The non-zero coefficients of C are then drawn from a zero-mean unit variance Gaussian distribution. The signal matrix S is obtained as in equation ( 1) with the noise matrix being drawn i.i.d from a zero-mean Gaussian distribution and variance so that the signal-to-noise ratio of each single signal is 10 dB.

We compare performances of M-BP, IrM-BP (with r = 0.5 and r = 1), M-SBL [START_REF] Wipf | An empirical bayesian strategy for solving the simultaneous sparse approximation problem[END_REF] and M-OMP [START_REF] Tropp | Algorithms for simultaneous sparse approximation. part i: Greedy pursuit[END_REF] for different experimental situations. The baseline context is M = 50, N = 25, k = 10 and L = 3. Note that for the M-OMP, we stop the algorithm after exactly k iterations. For this experiment, we did not performed model selection but instead tried several values of λ and σ and chosen the ones that maximize performances.

Figure 1 shows, from left to right, the performance averaged over 50 trials, on sparsity recovery when k increases from 2 to 20, when M goes from 25 to 150 and when L = 2, • • • , 8. We can note that, M-BP performs worse that IrM-BP. This is a result that we could expected in views of the literature [START_REF] Zou | The adaptive lasso and its oracle properties[END_REF], [START_REF] Candès | Enhancing sparsity by reweighted ℓ 1 minimization[END_REF] which compare Lasso and reweighted Lasso, the single signal approximation counterpart of M-BP and IrM-BP.

For all experimental situations, we remark that IrM-BP and M-SBL perform equally well. Again, this similar performances can easily be understood because of the strong relation between reweighted M-BP and M-SBL [START_REF] Rakotomamonjy | Simultaneous sparse approximations : insights and algorithms[END_REF]. When considering M-OMP, although we suppose that k is known, we can see that the M-OMP performance is not as good as those of M-SBL and IrM-BP.

We have also empirically assessed the computational complexity of our algorithms (we used r = 1 for IrM-BP). We varied one of the different parameters (dictionary size M , signal dimensionality N , number of signals L ) while keeping the others fixed. All matrices Φ, C and S are created as described above. Experiments have been run on a Pentium D-3 GHz with 4 GB of RAM using Matlab code. The results, averaged over 50 trials, in Figure 2 show the computational complexity of the different algorithms for different experimental settings. Note that we have also experimented on the M-SBL and M-FOCUSS computational performances owing to the code of Wipf et al. [START_REF] Wipf | An empirical bayesian strategy for solving the simultaneous sparse approximation problem[END_REF]. All algorithms need one hyperparameter to be set, for M-SBL and M-FOCUSS, we were able to choose the optimal one since the hyperparameter is dependent on a known noise level. For our algorithms, the choice of λ is more critical and has been manually set so as to achieve optimal performances. Note that our aim here is not give an exact comparison of computational complexity of the algorithms but just to give an order of magnitude of these complexities. Indeed, careful comparisons are difficult since the different algorithms do not solve the same problem and do not use the same stopping criterion.

We can remark in Figure 2 that with respects to the dictionary size, all algorithms present an empirical exponent between 1.4 and 2.4. Interestingly, we have theoretically evaluate the complexity of M-BP as quadratic whereas we measure a sub-quadratic complexity [START_REF] Rakotomamonjy | Simultaneous sparse approximations : insights and algorithms[END_REF]. We suppose that this happens because at each iteration, only the non-optimal c i,• 's are updated and thus the number of updates drastically reduces along iterations. We can note that among all approaches, M-BP is the less demanding algorithm while IrM-BP is the less efficient one. This is clearly the cost to be paid for trading the resolution of a non-convex problem against several convex ones. Note however, that this complexity can be controlled by reducing the number of iterations while preserving good sparsity recovery. This is the case of many weighted Lasso algorithms which use only two iterations [START_REF] Zou | One-step sparse estimates in nonconcave penalized likelihood models[END_REF], [START_REF] Zou | The adaptive lasso and its oracle properties[END_REF].

The difference between the two top plots in Figure 2 shows that algorithm complexities not only depend on the dictionary size but also on the redundancy of the dictionary. Indeed, on the right plot, signal dimensionality is related to the dictionary size (redundancy is kept fixed) while on the left plot, the signal size is fixed. This results in a non-uniform variation of the complexities which is difficult to understand. It is not clear if it is related to the problem difficulty or is intrinsic to algorithms. Further researches are still needed to clarify this point.

Bottom left plot of Figure 2 depicts the complexity dependency of all algorithms with respects to the number of signal to approximate. The results we obtain is in agreement with theoretical exponents since for M-BP and IrM-BP we have exponents of approximately 1 while the other algorithm complexities do not depend on L. On the bottom right, we have evaluated these exponents with respects to signal dimension.

Here again, we have results in accordance to theoretical expectations : M-BP and IrM-BP have lower complexities than M-SBL and M-FOCUSS. Furthermore, we note that IrM-BP has unexpectedly a very low exponent complexity. We assume that this is due to the fact that as dimension increases, the approximation problem becomes easier and thus needs less M-BP iterations.

IV. CONCLUSIONS

This paper aimed at contributing to simultaneous sparse signal approximation problems on several points. Firstly, we have proposed an algorithm for solving the multiple signal counterpart of Basis Pursuit Denoising named M-BP. The algorithm we introduced is rather efficient and simple and it is based on a soft-threshold operator which only needs matrix multiplications. Then, we have considered the more general non-convex M-FOCUSS problem for which M-BP is a special case. We have shown that M-FOCUSS can also been understood as an ARD approach. Indeed, we have transformed the M-FOCUSS penalty in order to exhibit some weights that automatically influence the importance of each dictionary elements in the approximation. Finally, we have introduced an iterative reweighted M-BP algorithm for solving M-FOCUSS. We also made clear the relationship between M-SBL and such a reweighted algorithm. We also provided some experimental results that show how our algorithms behave and how they compare to other methods dedicated to simultaneous sparse approximation. In terms of performances for sparsity profile recovery, our algorithms does not necessarily perform better than others approaches but they are provided with interesting features such as convexity and convergence guarantees.

Owing to this clear formulation of the problem and its numerically reproducible solution (due to convexity), our perspective on this work is now to theoretically investigate the properties of the M-BP and IrM-BP solutions. We believe that the recent works on the Lasso and related methods can be extended in order to make clear in which situations M-BP and Ir-MBP achieve consistency. Further improvements of algorithm speed can also be interesting so that tackling very large-scale approximation becomes tractable. Fig. 2. Estimation of the empirical exponent of the computational complexity of different algorithms (M-BP, IrM-BP, M-SBL and M-FOCUSS). For the two last algorithms, we have used the code available from Wipf's website. The top plots give the computation time of the algorithms with respects to the dictionary size. On the top left, the signal dimension N has been kept fixed and equal to 128 whereas on the right one, the signal dimension is related to the dictionary size. The bottom plots respectively depict the computational complexity with respects to the number of signals to approximate and the dimensionality of these signals.

Fig. 1 .

 1 Fig. 1. Results comparing performances of different simultaneous sparse algorithms. from left to right, we have varied the number of generating functions k, the dictionary size M and the number of signal to approximate L. The default setting is M = 50, N = 25, k = 10 and L = 3.