Elias Egho

Nicolas Jay

Chedy Raïssi

Dino Ienco

Pascal Poncelet

Maguelonne Teisseire
email: maguelonne.teisseire@teledetection.frpascal.poncelet@lirmm.fr

Amedeo Napoli
email: amedeo.napoli@loria.frchedy.raissi@inria.fr

A contribution to the discovery of multidimensional patterns in healthcare trajectories

Keywords: complex sequential patterns, multidimensional sequential patterns, data mining, complex data

Sequential pattern mining is aimed at extracting correlations among temporal data. Many different methods were proposed to either enumerate sequences of set valued data (i.e., itemsets) or sequences containing dimensional items. However, in real-world scenarios, data sequences are described as combination of both multidimensional items and itemsets. These heterogeneous descriptions cannot be handled by traditional approaches. In this paper we propose a new approach called MMISP (Mining Multidimensional Itemset Sequential Patterns) to extract patterns from complex sequential database including both multidimensional items and itemsets. The novelties of the proposal lies in: (i) the way in which the data are efficiently compressed; (ii) the ability to reuse and adopt sequential pattern mining algorithms and (iii) the extraction of new kind of patterns. We introduce a case-study on real-world data from a regional healthcare system and we point out the usefulness of the extracted patterns. Additional experiments on synthetic data highlights the efficiency and scalability of the approach MMISP.

Introduction

Real-world databases can be viewed as large and complex sources of information that can be analyzed for discovering new knowledge units or for decision making [START_REF] Yu | Proceedings of the 33rd international conference on Very large data bases[END_REF]. When the temporal dimension is also considered, every bit of information or event is associated with a timeline describing a total order over events. This total ordering introduces complexity in the extraction process. Many efficient approaches were developed to mine patterns depending on time or order, such as PrefixSpan [START_REF] Pei | Prefixspan: Mining sequential patterns by prefixprojected growth[END_REF], SPADE [START_REF] Zaki | Spade: An efficient algorithm for mining frequent sequences[END_REF], ClosSpan [START_REF] Yan | Clospan: Mining closed sequential patterns in large datasets[END_REF],...etc. These approaches focus on a single dimensional sequence dataset. However, there are many situations in which a database can be multidimensional, i.e. several characteristics of data can be ordered over time. Pinto et al. [START_REF] Pinto | Multi-dimensional sequential pattern mining[END_REF], Zhang et al. [START_REF] Zhang | Approxmgmsp: A scalable method of mining approximate multidimensional sequential patterns on distributed system[END_REF] and Yu et al. [START_REF] Yu | Mining sequential patterns from multidimensional sequence data[END_REF] introduced the notion of multidimensionality in a sequence and proposed several algorithms to mine this type of data. In data warehouse environments, a background knowledge is usually available in form of taxonomies, classification or concept hierarchies. Based on that, Plantevit et al. introduced M 3 SP [START_REF] Plantevit | Mining multidimensional and multilevel sequential patterns[END_REF], an algorithm able to incorporate several dimensions and the possible associated posets within the sequential pattern mining process.

The above approaches focus on sequences of homogeneous items. They do not pay attention to real-world complex data described by a vector of heterogeneous elements with different types, i.e. item or itemset. For example, in the healthcare domain, a patient trajectory is defined as a sequence of hospitalizations, where each hospitalization is defined as a vector of three heterogeneous elements: (i) healthcare institution, (ii) diagnosis and (iii) set of medical procedures. The healthcare institution and the diagnosis can be encoded as variables taking values, where values are organized within posets, while medical procedures are not comparable. This example shows that each dimension in data has to be managed in a proper and suitable way.

In this paper, we present an approach for mining multidimensional and heterogeneous patterns from medical patient trajectories, i.e. the trajectory of a patient visiting several hospitals. Our objective is to discover interesting patterns able to characterize patient stays and associated medical procedures. Such patterns can be interpreted by healthcare professionals to better understand patient pathways and improve the organization of care. Such multidimensional and heterogeneous patterns have to be mined by adapted and suitable methods. Accordingly, in this paper, we propose a new method to extract sequential patterns from databases including sequences of heterogeneous vectors. In addition, the approach is able to take into account background knowledge lying in term posets. The approach is original and efficient. An adapted algorithms is proposed which shows a very good behaviour on real-world medical data.

The remainder of this paper is organized as follows, Section 2 describes related work in classical and multidimensional sequential patterns. Section 3 introduces the problem statement and a running example. The algorithm for extracting complex frequent patterns is described in Section 4. Section 5 presents experimental results from both quantitative and qualitative points of views and Section 6 concludes the paper.

Related Work

Agrawal and Srikant [START_REF] Agrawal | Mining sequential patterns[END_REF] introduced the problem of mining sequential patterns over large sequential databases. Formally, given a set of sequences, where each sequence is a list of transactions ordered by time and each transaction is a set of items, the problem amounts to find all frequent subsequences that appear a sufficient number of times with a user-specified minimum support threshold (minsup). Following the work of Agrawal and Srikant many studies have con-tributed to the efficient mining of sequential patterns [START_REF] Mooney | Sequential pattern mining -approaches and algorithms[END_REF]. Most of them are based on the Apriori property, which states that any super pattern of a non-frequent pattern cannot be frequent. The main algorithms are PrefixSpan [START_REF] Pei | Prefixspan: Mining sequential patterns by prefixprojected growth[END_REF], SPADE [START_REF] Zaki | Spade: An efficient algorithm for mining frequent sequences[END_REF], SPAM [START_REF] Ayres | Sequential pattern mining using a bitmap representation[END_REF], PSP [START_REF] Masseglia | The PSP approach for mining sequential patterns[END_REF], DISC [START_REF] Chiu | An efficient algorithm for mining frequent sequences by a new strategy without support counting[END_REF], PAID [START_REF] Yang | Paid: Mining sequential patterns by passed item deduction in large databases[END_REF] and FAST [START_REF] Salvemini | Fast sequence mining based on sparse id-lists[END_REF]. All these algorithms aim at discovering sequential patterns from a set of sequences of itemsets such as customers who frequently buy DVDs of episodes I, II and III of Stars Wars, then buy within 6 months episodes IV, V, VI of the same famous epic space opera.

Many studies about sequential patterns discovery focus on single-dimensional sequences. However, in many situations, the database is multidimensional in the sense that items can be of different nature. For example, a consumer database can hold information such as article price, gender of the customer, location of the store and so on. Pinto et al [START_REF] Pinto | Multi-dimensional sequential pattern mining[END_REF] proposed the first work for mining multidimensional sequential patterns. In this work, a multidimensional sequential database is defined as a schema (ID, D 1 , ..., D m , S), where ID is a unique customer identifier, D 1 , ..., D m are dimensions describing the data and S is the sequence of itemsets. A multidimensional sequence is defined as a vector {d 1 , d 2 , ..., d m }, S 1 , S 2 , ..., S l where d i ∈ D i for (i m) and S 1 , S 2 , ..., S l , are the itemsets of sequence S. For instance, {P aris, M ale}, {mp 1 , mp 2 }, {mp 3 } describes a male patient who underwent procedures mp 1 and mp 2 in Paris and then underwent mp 3 also in Paris.

Here, dimensions remain constant over time, such as the location of the treatment. This means that it is not possible to have a pattern indicating that when the patient underwent procedures mp 1 and mp 2 in Paris then he underwent mp 3 in Nancy. Among other proposals, Yu et al [START_REF] Yu | Mining sequential patterns from multidimensional sequence data[END_REF] proposed two methods Aprior-iMD and PrefixMDSpan for mining multidimensional sequential patterns in the web domain. This study considers pages, sessions and days as dimensions. Actually, these three different dimensions can be projected into a single dimension corresponding to web pages, gathering web pages visited during a same session and ordering sessions w.r.t the day as order.

In real world applications, each dimension can be represented at different levels of granularity, by using a poset. For example, apples in a market basket analysis can be either described as fruits, fresh food or food. The interest lies in the capacity of extracting more or less general/specific multidimensional sequential patterns and overcome problems of excessive granularity and low support. Srikant and Agrawal [START_REF] Srikant | Mining sequential patterns: Generalizations and performance improvements[END_REF] proposed GSP which uses posets for extracting sequential patterns. The basic approach is based on replacing every item with all the ancestors in the poset and then the frequent sequences are generated. This approach is not scalable in a multidimensional context because the size of the database becomes the product of maximum height of the posets and number of dimensions. Plantevit et al [START_REF] Plantevit | Mining multidimensional and multilevel sequential patterns[END_REF] defined a multidimensional sequence as an ordered list of multidimensional items, where a multidimensional item is a tuple (d 1 , ..., d m ) where d i is an item associated with the i th dimension. They proposed M3 SP , an approach taking both aspects into account where each dimension is represented at different levels of granularity, by using a poset. M 3 SP is able to search for sequential patterns with the most appropriate level of granularity. Their approach is based on the extraction of the most specific frequent multidimensional items, which are then used as alphabet to rephrase the original database. Then, M 3 SP uses a standard sequential pattern mining algorithm to extract multidimensional sequential patterns. However, M 3 SP is not adapted to mine sequential databases, where sequences are defined over a combination of sets of items and items lying in a poset. Then it is not possible to have a pattern indicating that when the patient went to uh p for a problem of cancer ca, where he underwent procedures mp 1 and mp 2 , then he went to gh l for the same medical problem ca, where he underwent mp 3 ( i.e, (uh p , ca, {mp 1 , mp 2 }), (gh l , ca, {mp 3 }) ).

Compared to M 3 SP , the main contribution of this article is to generalize the concept of multidimensional sequence by considering multidimensional itemsets instead of multidimensional items. In our approach, an event in a sequence can be seen as a vector of itemsets, whenever M 3 SP represents an event as a tuple of atomic elements. This restriction prevents M 3 SP from extracting condensed and heterogeneous patterns.

Problem Statement

An introductory example

Firstly we propose an example to illustrate the present approach. The French healthcare system called PMSI 3 is a national information system used in France to manage hospital activity with both an economical and a medical points of view [START_REF] Rb Fetter | Case mix definition by diagnosis-related groups[END_REF]. In this system, each patient's stay is a standardized record of administrative and clinical data. Accordingly, each hospitalization can be formalized along three dimensions: (i) healthcare institution, (ii) diagnosis and (iii) medical procedures. The first two dimensions, i.e. healthcare institutions and diagnosis, are considered as variables whose values are lying in a poset. The last dimension is about medical procedures and can be considered as a variable which is setvalued. The basic sets of the healthcare institutions H, the diagnosis DG and the medical procedures MP, are the following:

-H = {T h , uh, gh, uh p , uh n , gh p , gh l }. -DG = {T d , ca, r, r 1 , r 2 , ca 1 , ca 2 , ca 3 }. -MP = {mp 1 , mp 2 , mp 3 , mp 4 }.
Hospital and diagnosis can be described at different levels of granularity through two posets (H, )and (DG, ) which are defined in Figure 1.

In the poset (H, ), uh p and uh n denote university hospitals, uh, with uh p uh and uh n uh; gh p and gh l are general hospitals, gh, with gh p gh and gh l gh.

The hospitalization of a patient is defined as a vector with three components, (h, dg, mp). The component h refers to institution and its value lies in the An elementary sequence e = (e 1 , e 2 , ..., e n ) is defined as a vector of n elements where an element can be either:

an atomic element taken from a partially order set.

a subset of a non ordered set.

Definition 2. (Ordering of elementary sequences)

Given two elementary sequences e = (e Given two sequences S = S 1 , S 2 , ..., S l and T = T 1 , T 2 , ..., T l , S is more specific than T , denoted by S ≤ s T , if there exist a set of indices 1 ≤ i 1 < i 2 < ... < i l ≤ l such that S j ≤ e T ij for all j = 1 . . . l and l l. The most specific sequence is also the longest sequence as l l.

Example 3. S = (uh p , ca A set of m sequences S DB = {S 1 , S 2 , ..., S m } is called hereafter a "sequential database". An example is given in Table 1 where there are four sequences describing four patient trajectories.

Definition 5. (Support of sequence S)

Let S DB = {S 1 , ..., S k } a sequential database. The support of a sequence S, denoted by support s (S) is defined as follows:

support s (S) = |{S i ∈ S DB ; S i ≤ s S}| Definition 6. (Sequential pattern)
Given a positive integer σ as a minsup threshold and a sequential database S DB , the sequence S is called a sequential pattern in S DB iff support s (S) ≥ σ.

Example 4. The sequence S = (uh p , ca, {mp 1 , mp 2 }), (gh l , r, {}) has a support equals to 3 (i.e, support s (S) = 3) in the sequential database S DB (see Table 1). It is a sequential pattern according to minsup threshold equals to 3 (i.e, σ = 3).

Patients

Trajectories

S 1 (uhp, ca1, {mp1, mp2}), (uhp, ca1, {mp1}), (gh l , r1, {mp3}) S 2 (uhn, ca1, {mp4}), (uhp, ca2, {mp1, mp2}), (gh l , r1, {mp2}) S 3 (uhn, ca3, {mp4}), (gh l , r2, {mp3}) S 4
(uhp, ca2, {mp1, mp2}), (ghp, r2, {mp3}), (gh l , r2, {mp2}) Table 1: An example of a database of patient trajectories.

Most specific sequential patterns

In this section, we present the problem of mining the most specific sequential patterns. Given a sequential database, mining all possible frequent patterns results in a huge amount of information that is difficult to manage from the analyst point of view. To overcome this problem, we extract a set of sequential patterns that are not only frequent but also as specific as possible. This second constraint allows to reduce the volume of the final result discarding redundant patterns. An extracted pattern is called a "most specific sequential pattern".

Definition 7. (Most Specific Sequential Pattern (MSSP))

Given a positive integer σ as minsup threshold and a sequential database S DB , a sequence S is a most specific sequential pattern in S DB or MSSP if and only if support s (S) ≥ σ and there does not exist any sequence T such that T ≤ S S with support s (S) = support s (T ) ≥ σ Actually, frequency is monotone (i.e whenever S is frequent, any generalization of S is also frequent). For example, if S = (uh p , c, {mp 1 , mp 2 }) is frequent then T = (uh, c, {mp 1 }) which is more general than S is also frequent. Thus, the most specific sequential patterns are sufficient to retrieve all sequential patterns.

Example 5. Let σ = 3 (i.e. a sequence is frequent if it appears at least three times in S DB ). The sequence S = (uh, ca, {mp 1 }) is frequent but is not the most specific one because the sequence T = (uh p , ca, {mp 1 , mp 2 }) is frequent and verifies T ≤ s S and support s (S) = support s (T ). The sequence T = (uh p , ca, {mp 1 , mp 2 }) is a most specific frequent sequence in S DB : (i) T is frequent (support s (S) ≥ σ) and (ii) there is no other sequence in S DB which is frequent, more specific than T and has the same support of T .

The objective of our approach is to extract a set of frequent sequential patterns that are as specific as possible. In the next section, we present the algorithm MMISP for finding most specific sequential patterns. The basic idea of MMISP consists in transforming the sequential database into an "adapted form" and then to apply a standard algorithm for sequential mining.

The MMISP algorithm

MMISP is based on three steps:

First step (Extraction of frequent elementary sequences)

The algorithm searches for the frequent and specific elementary sequences. It extracts these frequent elementary sequences w.r.t. partial ordering existing between the elements (posets).

Second step (Transformation)

In this step, all frequent elementary sequences extracted in the previous step are mapped to an alternate representation. Then, the sequential database is encoded by using this new representation. Third step (Mining)

In this step, a standard sequential algorithm is applied on the sequential database produced at the preceding step.

Extracting all frequent elementary sequences

Firstly, MMISP considers the elementary sequences in all sequences of S DB . Actually if the elementary sequence is not frequent in a sequential database S DB it does not belong to an extracted sequential pattern. For example, If we have two sequences:

S 1 {s2}{s1}{s3} S 2 {s2}{s4}{s3}
Given a support threshold σ = 2, pattern {s 2 }{s 3 } is frequent. In this example, s 2 and s 3 are frequent items while s 1 and s 4 are not frequent items. In the following, we only consider sequential patterns composed of frequent items.

MMISP extracts all the frequent elementary sequences in a sequential database by taking into account the partial order relation between their elements. The support of an elementary sequence s is defined as follows: Definition 8. (Support of an elementary sequence e)

Let S DB = {S 1 , ..., S m } a sequential database. The support of an elementary sequence e = (e 1 , e 2 , ..., e n ), denoted by support e (e), is defined as follows:

support e (e) = |{S i ∈ S DB ; S i = S i 1 , ..., S i l and ∃j ∈ [1, .., l]; S i j ≤ e e}| Example 6. The support of (gh, r, {mp 3 }) is 3 in S DB as:

-S 1 3 = (gh l , r 1 , {mp 3 }) ≤ e (gh, r, {mp 3 }). -S 3 2 = (gh l , r 2 , {mp 3 }) ≤ e (gh, r, {mp 3 }). -S 4 2 = (gh p , r 2 , {mp 3 }) ≤ e (gh, r, {mp 3 }).
The frequent elementary sequences are ordered over a poset L, denoted by (L, ≤ e ), as follows. Firstly, we generate the most general elementary sequences. In the running example, we consider triples of the form (h, dg, mp) and the most general triple is (T h , T d , {}) where T h and T d denote the most general items in the posets H and D respectively, and the empty set {} denotes the most general item in the set M P . Then, we recursively generate new elementary sequences by starting from the most general one. This generation is done by replacing each element e 1 , e 2 , ..., e n ∈ e with all of its direct specializations.

The set of all direct specializations of an element e i , denoted by desc(e i ), is defined as follows: Definition 9. (Direct specializations of an element e i ) Let e i be an element in the set D i . The direct specializations of e i , denoted by desc(e i ), is defined by:

desc(e i ) = {u ∈ D i ; u ≤ e i and w ∈ D i ; u ≤ w and w ≤ e i } if D i is a poset. {e i ∪ {u}; u ∈ D i \ e i } if D i is a set.
The set of all direct specializations of an elementary sequence e, denoted by desc(e), is defined as follows:

Definition 10. (Direct specializations of an elementary sequence e) Let e = (e 1 , e 2 , ..., e n ) be an elementary sequence. The direct specializations of an elementary sequence e, denoted by desc(s), is defined by: desc(e) = (e 1 , e 2 , ..., e n ) | ∃i ∈ {1, ..., n} (e i ∈ desc(e i )) and (∀j = i)(e j = e j ) .

Example 7. Given the most general elementary sequence s=(T h , T d , {}), the direct specializations of T h are uh and gh, the direct specializations of T d are r and c and the direct specializations of {} are {mp 1 }, {mp 2 }, {mp 3 } and {mp 4 }. Thus, the direct specializations of (T

h , T d , {}) are (uh, T d , {}), (gh, T d , {}), (T h , r, {}), (T h , ca, {}), (T h , T d , {mp 1 }), (T h , T d , {mp 2 }), (T h , T d , {mp 3 }) and (T h , T d , {mp 4 }).
The frequency of an elementary sequence is anti-monotone w.r.t the specificity of elementary sequence, i.e., whenever an elementary sequence e is not frequent, all the specializations of e are also not frequent. For example, (T h , T d , {mp 4 }) is not frequent when minsup threshold is equal to 3, then all the specializations of (T h , T d , {mp 4 }) such as (gh, T d , {mp 4 }) or (gh, r, {mp 1 , mp 4 })... are also not frequent. We use this anti-monotonicity to prune the enumeration space and efficiently build the poset (L, ≤ e ).

Specialization is applied recursively for each new frequent elementary generated sequence. Figure 2 shows an example of generation of a poset of (L, ≤ e ) which is detailed below:

-Considering the most general elementary sequence (T h , T d , {}), our approach generates seven new frequent elementary sequences, which are: (uh, Given a positive integer σ as minsup threshold and a sequential database S DB , an elementary sequence e is a most specific frequent elementary sequence in S DB or MSFES if and only if support e (e) ≥ σ and there does not exist any elementary sequence e such that support e (e ) ≥ σ and e ≤ e e.

T d , {}), (gh, T d , {}), (T h , r, {}), (T h , ca, {}), (T h , T d , {mp 1 }), (T h , T d , {mp 2 }) and (T h , T d , {mp 3 }). -Based on (uh, T d , {}), MMISP generates (uh p , T d , {}), (uh, ca, {}), (uh, T d , {mp 1
Table 2 shows the set of most specific frequent elementary sequences which are extracted from (L, ≤ e ). Algorithm 1 and 2 describe the two steps for extracting all the frequent elementary sequences and the most specific ones. 1 and (ii) the two posets H and DG in Figure 1 and the set M P = {mp 1 , mp 2 , mp 3 , mp 4 } with minsup threshold equals to 3.

id MSFES 1 (uhp, ca, {mp1, mp2}) 2 (gh l , r, {}) 3 (gh, r, {mp3})

Transformation and mining sequences

We now study the temporal relation between the extracted specific frequent elementary sequences as follow. Firstly, we replace each elementary sequence in each sequence of S DB with all its generalizations from M SF ES set. Given a sequence S = S 1 , ..., S n in S DB the replacement consists in substituting each elementary sequence S i in S by several elementary sequences e ∈ M SF ES such that S i e e.

Example 8. The sequence S 1 in S DB , (uh p , ca 1 , {mp 1 , mp 2 }), (uh p , ca 1 , {mp 1 }), (gh l , r 1 , {mp 3 }) is transformed into {(uh p , ca, {mp 1 , mp 2 })}, {(gh l , r, {}), (gh, r, {mp 3 })} : Table 3 shows the transformation of S DB in Table 1 based on the set of all most specific frequent elementary sequences MSFES in Table 2.

-S 1 1 = (uh p , ca 2 , {mp 1 , mp 2 }) is replaced by (uh p , ca, {mp 1 , mp 2 }) from MSFES set in
In a classical sequential pattern mining algorithm, the sequential database to be mined should be a set of pairs (sid, s) where sid is a unique sequence identifier and s is a sequence of itemsets. Thus S DB in Table 3 is transformed as follows:

-Each elementary sequence in the MSFES is assigned a unique id which is used during the mining (see Table 2) . {(uhp, ca, {mp1, mp2})}, {(gh, r, {mp3})}, {(gh l , r, {})} Table 3: Transforming the patient trajectories in Table 1 by using the set of all most specific frequent elementary sequences in Table 2.

-For each sequence S i in the transformed database (see Table 3) and for each elementary sequence T in S i j ; S i j ∈ S i , T is replaced by its id.

Example 9. The sequence S 1 = {(uh p , ca, {mp 1 , mp 2 })}, {(gh l , r, {}), (gh, r, {mp 3 })} in Table 3 is transformed into {1}, {2, 3} such as:

-(uh p , ca, {mp 1 , mp 2 }) in S 1 has id 1 in Table 2.

-(gh l , r, {}) in S 1 has id 2 in Table 2.

-(gh, r, {mp 3 }) in S 1 has id 3 in Table 2.

Table 4 shows the transformation sequential database of Table 3 by using the identifiers of all most specific frequent elementary sequences MSFES in Table 2.

Patients Trajectories S 1 {1}{2, 3} S 2 {1}{2} S 3 {2, 3} S 4 {1}{3}{2}
Table 4: Transformed database in Table 3 Then in MMISP, we use CloSpan [START_REF] Yan | Clospan: Mining closed sequential patterns in large datasets[END_REF] as the sequential pattern mining algorithm. Table 5 displays all sequential patterns in their transformed format and the frequent patient trajectories in which identifiers are replaced with their actual values, with σ = 3.

Experiments

We conduct experiments on both real and synthetic datasets. The MMISP algorithm is implemented in Java and the experiments are carried out on a MacBook Pro with a 2.5GHz Intel Core i5, 4GB of RAM Memory running OS X 10.6.8. Extraction of sequential patterns is based on the public implementation of CloSpan algorithm [START_REF] Yan | Clospan: Mining closed sequential patterns in large datasets[END_REF] supplied by the IlliMine4 toolkit.

Frequent sequential patterns Frequent patient trajectory patterns Support {2} (gh l , r, {}) 4 {3}

(gh, r, {mp3}) 3 {1}{2}

(uhp, ca, {mp1, mp2}), (gh l , r, {}) 3

Table 5: All the most specific sequential patterns extracted from S DB in Table 1 with σ = 3.

Healthcare Trajectory

Mining healthcare trajectories

In order to assess the effectiveness of our approach, we run several experiments on the PMSI system for describing and analyzing patient trajectories. In PMSI, each hospitalization is characterized by the following dimensions: hospital, principal diagnosis and procedures delivered during the stay.

The hospital dimension is associated with a geographical poset of 4 levels: root (France), administrative region, administrative departement and hospital. As illustrated in Figure 4, University Hospital of Nancy (coded as 540002078) is a hospital in Meurthe-et-Moselle and Meurthe-et-Moselle is a department in Lorraine which is a region of France. The number of nodes in this taxonomy is 151 nodes. Principal Diagnosis could be described at 5 levels of the 10 th International classification of Diseases (ICD10): root, chapter, block, 3 characters and 4 characters. As illustrated in Figure 5, chapters such as N eoplasms have specializations: block C30-C39 which is a malignant neoplasms of respiratory, block C50-C50 which is a malignant neoplasms of breast etc. The block C30 -C39 has specializations: malignant neoplasm of larynx (coded as C34), malignant neoplasm of bronchus and lung (coded as C32), etc. C34 (Lung cancer) has specializations: C340 is a cancer of the main bronchus, C341 is a cancer of upper lobe etc. The number of nodes in the disease taxonomy is 1543 nodes.

Procedures were represented by their first CCAM5 code. For example, ZBQK is a chest radiography, GF F A is a pneumonectomy etc.

Our dataset contains 828 patients suffering from lung cancer and living in Lorraine Region, in the East of France. Table 6 shows an example of care trajectories for 3 patients. For example, P atient 1 has two hospitalizations. He was admitted in the University Hospital of Nancy (coded as 540002078), for a Lung cancer (coded as C341), and underwent a chest radiography (coded as ZBQK). Then, he was hospitalized in a private clinic in Metz (coded as 570023630), for a chemotherapy session (coded as Z51) where he had a chest radiography and pneumonectomy (coded as GF F A). As a characterization of the care trajectory database, Figure 6 shows the distribution of the length of healthcare trajectories in our dataset, the median length is of 11 stays. The support value is set to 40 patients (i.e. σ = 5%) for this experiment . MMISP generates 615 different frequent trajectories. Table 7 shows some extracted patterns. Pattern #1 can be interpreted as follows: 5% of patients had a hospitalization in Meurthe et Moselle department for any kind of lung cancer (coded as C34). They underwent three medical procedures: chest radiography (coded as ZBQK) with an electrocardiography (coded as DEQP) and a therapeutic procedure on blood (coded as FELF). Pattern #2 shows that 8% of patients had a hospitalization in Lorraine, because of poisoning.

This kind of information helps healthcare managers and decision makers in planning and organizing healthcare resources at regional level. Besides, sequential patterns can be seen as condensed representations of the care trajectories.

MMISP versus Standard sequential pattern mining method

In this section, we compare MMISP with a standard sequential pattern mining method such as CloSpan. All standard sequential pattern mining algorithms requires the mined dataset to be mined is composed of pairs of the form (id, seq), where id is a sequence identifier and seq is a sequence of itemsets. So, we replace each sequence S i in S DB with an extended-sequence S i . Each elementary sequence of a sequence S i is transformed into a single itemset by replacing each element in the elementary sequence with all its ancestors. For example, with the posets shown in Figure 1, an elementary sequence (uh p , ca 1 , {mp 1 , mp 2 }) would be replaced with {T h , uh, uh p , T d , ca, ca 1 , mp 1 , mp 2 }. Table 8 shows the extended-sequential database S DB of the sequential database S DB in Table 1. Then, we apply CloSpan as a sequential pattern mining algorithm on the extended-sequential database. This way of managing taxonomies has been used in GSP which is proposed by Srikant and Agrawal [START_REF] Srikant | Mining sequential patterns: Generalizations and performance improvements[END_REF].

Patients

Extended sequence S 1 {Th, uh, uhp, Td, ca, ca1, mp1, mp2}, {Th, uh, uhp, Td, ca, ca1, mp1}, {Th, gh, ghl, Td, r, r1, mp3} S 2 {Th, uh, uhn, Td, ca, ca1, mp4}, {Th, uh, uhp, Td, ca, ca2, mp1, mp2}, {Td, gh, ghl, Td, r, r1, mp2} S 3 {Th, uh, uhn, Td, ca, ca3, mp4}, {Th, gh, ghl, Td, r, r2, mp3} S 4 {Th, uh, uhp, Td, ca, ca2, mp1, mp2}, {Th, gh, ghp, Td, r, r2, mp3}, {Th, gh, ghl, Td, r, r2, mp2}

Table 8: An extended sequential database of patient trajectories in Table 1.

Our main goal is to evaluate the quality of the patterns mined with MMISP and its performance compared to the "naive" approach using CloSpan. So, we transform the 828 patients trajectories, then we apply the two approaches using several minimum support threshold ranging from 100 % to 5%. Figure 7 reports the running time and the number of pattens according to different values of support threshold for both CloSpan and MMISP.

Actually, CloSpan cannot finish its calculations for support threshold less than 50 % because the transformation increases number of items in itemsets and generates a large number of similar sequences. By contrast, MMISP runs in acceptable time for support as low as 5 %.

MMISP is able to extract condensed patterns w.r.t. the ones mined by CloSpan. For example, the sequential pattern {T h , uh, uh p , T d , ca, mp 1 , mp 2 } {T h , gh, gh l , T d , r} generated by CloSpan contains redundant information as a hospitalization containing uh p will also contain T h and uh. MMISP is not CloSpan extracts all the sequential patterns while MMISP generates just the more specific ones. For example, if (uh p , ca, {mp 1 , mp 2 }) is a sequential pattern. MMISP does not extract the patterns which are more general such as (uh, T d , {mp 1 }) while CloSpan extracts both the general and specific ones. Figure 7 shows the differences between the number of sequential patterns extracted by CloSpan and MMISP. For example, with a support threshold of 50 %, MMISP extracts 150 sequential patterns while CloSpan extracts 6335683 sequential patterns.

Finally, we may conclude that:

-MMISP is more efficient than CloSpan over extended-sequential database with low support threshold. -The sequential patterns extracted by CloSpan require post processing while this is not the case with MMISP. -MMISP extracts just the most specific sequential patterns while CloSpan extracts both general and specific ones. This means that CloSpan extracts a huge number of sequential patterns. Analyzing all these sequential patterns is not an easy task for healthcare managers and decision makers.

M M ISP versus M 3 SP

Another experiment was carried out for comparing M 3 SP with MMISP. Our main goal is to evaluate the effectiveness of sequential patterns mined by MMISP compared to the ones extracted by M 3 SP . For this purpose, we applied M 3 SP with hospital, diagnosis and medical procedures as analysis dimensions. The support value is set to 40 patients (i.e. σ = 5%). Table 9 reports an example of the extracted patterns with M 3 SP and MMISP.

Firstly, we observe that MMISP is able to extract condensed patterns w.r.t. the ones mined by M 3 SP . For example, 48 sequential patterns, Pattern #1,..., Pattern #48, generated by M 3 SP are summarized by 3 sequential patterns, Pattern #50, Pattern #51 and Pattern #52, extracted by MMISP (see Table 9). This shows that the rigid structure of multidimensional item assumed by M 3 SP limits the expressivity of the results.

Besides that, in M 3 SP , several dimensions can be repeated in the same hospitalization. For example, in M -Some sequential patterns can be extracted by MMISP while they cannot be extracted by M 3 SP . Table 9: Some patterns obtained by M 3 SP and MMISP.

M M ISP +

In our approach, the fundamental step is the first one which is "extraction of frequent elementary sequences", because it provides all elements that will occur in sequences to be mined. In this step, MMISP extracts only the most specific frequent elementary sequences from the poset (L, ≤ e ). As a result, it is impossible to find a sequential pattern that contains an elementary sequence which is comparable with another elementary sequence in the same pattern or in the another one. Formally, the sequential patterns extracted by MMISP have the following property: for any two sequential patterns SP 1 = SP 1 1 , SP For example, in Table 1 the pattern (uh, ca, {}), (gh l , r, {}) is frequent according to minsup threshold equals 3. This pattern does not appear in the results of MMISP because in the first step MMISP extracts (uh p , ca, {mp 1 , mp 2 }) which is more specific than (uh, ca, {}) and still frequent. As a result, MMISP does not extract any pattern which includes elementary sequence more general than (uh p , ca, {mp 1 , mp 2 }).

To solve this problem, we propose an extension to our approach which is called M M ISP + . In M M ISP + , instead of choosing just the most specific elementary sequence in (L, ≤ e ), we choose all the elementary sequences in (L, ≤ e ). Thus, we replace each elementary sequence in each sequence of S DB with all of its generalizations from (L, ≤ e ) and not just with the most specific ones. For example, the elementary sequence (gh p , r 2 , {mp 3 }) in the sequence S 4 in Table 1 is replaced by (gh, r, {mp 3 }), (gh, r, {}), (gh, T d , {mp 3 }), (T h , r, {mp 3 }),(gh, T d , {}), (T h , r, {}), (T h , T d , {mp 3 }) and (T h , T d , {}) from the (L, ≤ e ) set in Figure 3. Then, we use CloSpan as a sequential pattern mining algorithm on the transformed sequential database.

Frequent patient trajectory patterns Support (uh, ca, {}), (gh l , r, {}) 4 (uh, ca, {}), (T h , T d , {}), (gh l , r, {}) 3 (uh, ca, {}), (gh, r, {mp3}) 3 (uhp, ca, {mp1, mp2}), (gh l , r, {}) 3

Table 10: The sequential patterns extracted from S DB with support threshold equals 3.

Applying M M ISP + generates non condensed patterns w.r.t MMISP. For example, the pattern {(T h , T d , {}), (uh, T d , {}), (T h , ca, {}), (uh, ca, {})}{(T h , T d , {}), (gh, T d , {}), (T h , r, {}), (gh l , T d , {}), (gh, r, {}), (gh l , r, {})} generated by M M ISP + contains redundant information. as a patient having the hospitalization (uh, ca, {}) will also have (T h , T d , {}), (uh, T d , {}) and (T h , ca, {}). Redundancy can be removed by post-processing. Table 10 shows all the sequential patterns extracted from S DB in Table 1 by applying M M ISP + . These patterns have been post-processed to remove redundant information.

Figure 8 reports the number of patterns extracted from our dataset (i.e, 828 patient trajectories) according to different values of minsup threshold for applying the two solutions: MMISP and M M ISP + . Actually, M M ISP + cannot finish its calculations for support threshold less than 50 %. This happens because the CloSpan algorithm (i.e, the third step in M M ISP + ) cannot process with support threshold less than 50 %.

Experiments on Synthetic Datasets

In the second experiment, we study the scalability of the MMISP approach. We consider the number of extracted patterns and the running time with respect several parameters:

number of dimensions with an associated subsumption relation.

number of dimensions without any associated subsumption relation.

number of elementary sequences in each sequence (i.e. sequence length).

depth of the poset of elements with an associated subsumption relation.

number of sequences in a sequential database.

In the following, we use the term "set-element" for an element lying in a set and the term "poset-element" for an element lying in a partially ordered set.

The first batch of synthetic data generated contains 1000 sequences defined over 2, 3, 4 and 5 poset-elements and 2 set-elements. Each sequence contains 15 In Figure 11, we study the performance of MMISP by considering several levels of granularity within the posets. We generated 1000 sequences defined over 3 poset-elements and 3 set-elements. Each poset is defined over 3, 4, 5, 6 levels of granularity. Each sequence contains 15 elementary sequences. The number of extracted sequential patterns does not change with each newly added level as MMISP extracts only the most specific sequential patterns.

We study the performance of MMISP and the number of extracted sequential patterns with respect the number of sequences in a sequential database and the length of each sequence. Figure 12 shows the running time and the number of patterns extracted for 1000 sequences with 3 poset-elements and 3 set-elements with varying sequence length. Figure 13 shows the running time and the number of sequential patterns extracted for several number of sequences (1000, 2000, 3000, 4000 and 5000 sequences) also with 3 poset-elements and 3 set-elements. In this experiments, the sequence length is roughly equal to 15.

Figures 9 -13 highlight the fact that MMISP is efficient in terms of runtime for a large panel of sequences when varying different parameters.

Conclusion

In this paper, we propose a new approach to mine a sequential database of heterogeneous sequences. We provide formal definitions and propose a new algorithm MMISP to mine this kind of sequences. The MMISP algorithm relies on external posets to improve the mining process and produces results with appropriate levels of granularity. We conduct experiments on both real-world and synthetic datasets. The method is applied on real-world data where the problem is to mine healthcare patients trajectories and gives potential interesting patterns for healthcare specialists.

For future work, we are planning to use statistical significance tests to evaluate the sequential patterns extracted and choose the most significant ones. On the other hand, proposing a graphical interface to visualize and query the sequential patterns. We are also interested in generalizing our method by considering sequences all elements of which are lying in partially ordered set.

Finally, we are aware that choosing the most specific frequent elementary sequence to mine the sequential patterns prevents us from extracting all of the most specific sequential patterns. Coping with this issue is another interesting extension of the present work that we plan to investigate in the future.

(Fig. 2 :

 2 Fig. 2: The steps of generating the elementary sequences in (L, ≤ e ) with min-sup= 3.

Figure 3

 3 Figure 3 shows the poset (L, ≤ e ) generated with σ = 3. As the objective of MMISP is to extract specific sequential patterns, we retain only the most specific elementary sequences in (L, ≤ e ). The most specific frequent elementary sequences is defined as follows: Definition 11. (Most Specific Frequent Elementary Sequence (MSFES))Given a positive integer σ as minsup threshold and a sequential database S DB , an elementary sequence e is a most specific frequent elementary sequence in S DB or MSFES if and only if support e (e) ≥ σ and there does not exist any elementary sequence e such that support e (e ) ≥ σ and e ≤ e e.

(Fig. 3 :

 3 Fig.3: The poset (L, ≤ e ) is generated by taking into account: (i) the sequential database S DB in Table1and (ii) the two posets H and DG in Figure1and the set M P = {mp 1 , mp 2 , mp 3 , mp 4 } with minsup threshold equals to 3.

Fig. 4 :Fig. 5 :

 45 Fig. 4: A geographical poset of the healthcare institution

Patients Trajectories P atient 1 (

 1 C341, 540002078, {ZBQK}), (Z51, 570023630, {ZBQK, GF F A}), . . . P atient 2 (C770, 100000017, {ZBQK}), (C770, 210780581, {ZZQK, Y Y Y Y }), . . . P atient 3 (H259, 210780110, {Y Y Y Y }), (H259, 210780110, {ZZQK}), . . .

Fig. 6 :

 6 Fig. 6: Distribution length of the patient trajectories

Fig. 7 :

 7 Fig. 7: MMISP versus CloSpan

-

  3 SP , Pattern #48 represents one hospitalization including 9 multidimensional items. Each multidimensional item is associated with the same value of hospital and diagnosis (570000588 and C341) and different values of medical procedures. By contrast, Pattern #50 extracted by MMISP represents the same trajectory as Pattern #48. Pattern #50 has one elementary sequence with three elements: hospital 570000588, diagnosis C341 and a set of medical procedure {ZBQK, DEQP, GF F A, GLLD, GELD, ZZQK, GELE, F CF A, AGLB}. Pattern #50 is much more compact and informative than Pattern #48. Given a minsup threshold, MMISP extracts sequential patterns that are not found by M 3 SP . For instance, Pattern #53 extracted by MMISP is not found by M 3 SP . This is due to the fact that MMISP extracts new frequent elementary sequences not extracted by M 3 SP . For instance e 1 =(Lorraine, Diseases, {GEQE, ACQH, ZCQH}) and e 2 =(Lorraine,Diseases of the respiratory, {ACQH}) are extracted by MMISP. As e 1 and e 2 are frequent and not comparable (i.e. e 1 e e 2 and e 2 e e 1 ), M 3 SP extracts only (Lorraine, Diseases of the respiratory, ACQH) and not (Lorraine, Diseases, ACQH) as (Lorraine, Diseases of the respiratory, ACQH) is more specific than (Lorraine, Diseases, ACQH). From a quantitative point of view, MMISP extracts 803 frequent elementary sequences with 615 sequential patterns while M 3 SP extracts 331 multidimensional items with 470 multidimensional sequential patterns. The execution time of M 3 SP is about 82 seconds while MMISP takes about 98 seconds. Finally, we may conclude that: Several sequential patterns generated by M 3 SP can be summarized by only one sequential pattern mined by MMISP. -Several multidimensional items generated by M 3 SP can be summarized by only one elementary sequence in MMISP. -One elementary sequence in MMISP represents one hospitalization in the trajectory while one multidimensional item in M 3 SP represents only a part of hospitalization in the trajectory.

  C341, GF F A)(570000588, C341, ZZQK)} 2 {(570000588, C341, DEQP )(570000588, C341, GF F A)(570000588, C341, ZZQK)} . . .

Fig. 8 :

 8 Fig. 8: Number of sequential pattern extracted according to different values of support threshold for both MMISP and M M ISP +

  elements & 2 set-elements 3 poset-elements & 2 set-elements 4 poset-elements & 2 set-elements 5 poset-elements & 2 set-elements (a) Runtime sequences over frequency threshold. elements & 2 set-elements 3 poset-elements & 2 set-elements 4 poset-elements & 2 set-elements 5 poset-elements & 2 set-elements (b) Number patterns over frequency threshold.

Fig. 9 :

 9 Fig. 9: Number of sequential pattern extracted and Running time obtained by MMISP with varying in the number of poset-elements.

  elements & 2 set-elements 3 poset-elements & 3 set-elements 3 poset-elements & 4 set-elements 3 poset-elements & 5 set-elements (a) Runtime sequences over frequency threshold. elements & 2 set-elements 3 poset-elements & 3 set-elements 3 poset-elements & 4 set-elements 3 poset-elements & 5 set-elements (b) Number patterns over frequency threshold.

Fig. 10 :

 10 Fig. 10: Number of sequential pattern extracted and Running time for a large panel of sequences when varying over the number of set-elements.

Fig. 11 :

 11 Fig. 11: Number of sequential pattern extracted and Running Time obtained for a large panel of sequences when varying over the levels of granularity between poset-elements.

Fig. 12 :

 12 Fig. 12: Number of sequential pattern extracted and Running time for a large panel of sequences when varying over sequence length.

Fig.

  Fig. Number of sequential pattern extracted and Running time for a large panel of sequences when varying over several number of sequences.

  , ca 1 , {mp 1 , mp 2 }) describes an hospitalisation, where uh p ∈ H, ca 1 ∈ DG and {mp 1 , mp 2 } ⊆ M P .

		T h			T d
	U niversity.H(uh)	General.H(gh)	Respiratory(r)	Cancer(ca)
	uh p	uh n	gh p	gh l	r 1	r 2	ca 1	ca 2	ca 3
	Healthcare Institution (H, h )		Diagnosis (DG, dg )
	Fig. 1: The posets for the healthcare institution set H and the diagnosis set DG
	partially ordered set (H, ). The same thing applies to dg with (DG, ). The
	component mp refers to medical procedures and is a set-valued.
	Then, (uh p 3.2 Basic definitions
	Definition 1. (Elementary sequence)

  1 , e 2 , ..., e n ) and e = (e 1 , e 2 , ..., e n ), e is more specific than e , denoted by e ≤ e e , if for all i = 1...n we have:e i e i where e i , e i are elements of a poset. e i ⊆ e i where e i , e i are sets.Example 1. e = (uh p , ca 1 , {mp 1 , mp 2 }) is an elementary sequence more specific than e = (uh, T d , {mp 1 }) as:uh p uh; uh p , uh ∈ H.-ca 1 T d ; ca 1 , T d ∈ DG. -{mp 1 } ⊆ {mp 1 , mp 2 }; {mp 1 }, {mp 1 , mp 2 } ⊆ M P . = S 1 , S2, ..., S l is a set of elementary sequences S i ordered by the temporal order relation < t , such as S 1 < t S 2 < t S 3 ... < t S l .

	Definition 3. (Sequence)
	A sequence is represented as S Definition 4. (Ordering of sequences)
	Example 2. A patient healthcare trajectory can be considered as a set of hos-
	pitalisations ordered over time. (uh p , ca 1 , {mp 1 , mp 2 }), (gh l , r 1 , {mp 2 }) repre-
	sents a patient trajectory with two hospitalizations. The patient was first admit-
	ted in the university hospital uh p for a lung cancer ca 1 , and underwent proce-
	dures mp 1 and mp 2 . Then he visited the general hospital gh l for a pneumonitis
	r 1 where he underwent procedure mp 2 .

  [START_REF] Agrawal | Mining sequential patterns[END_REF] , {mp 1 , mp 2 }), (gh l , r 1 , {mp 2 }) is a sequence with two elementary sequences S 1 = (uh p , ca 1 , {mp 1 , mp 2 }) and S 2 = (gh l , r 1 , {mp 2 }) where S 1 < t S 2 . The sequence S is more specific than T = (uh, T d , {mp 1 }) , S ≤ s T , as (uh p , ca 1 , {mp 1 , mp 2 }) ≤ e (uh, T d , {mp 1 }).

  }) and (uh, T d , {mp 2 }). -Based on (uh p , T d , {}), MMISP generates (uh p , ca, {}), (uh p , T d , {mp 1 }) and (uh p , T d , {mp 2 }). -Based on (uh p , ca, {}), MMISP generates (uh p , ca, {mp 1 }) and (uh p , ca, {mp 2 }). -Based on (uh p , ca, {mp 1 }), MMISP generates (uh p , ca, {mp 1 , mp 2 }).-By using (uh p , ca, {mp 1 , mp 2 }), no any new frequent elementary sequences can be found, thus the generation stops.

Table 2 :

 2 The most specific frequent elementary sequences extracted from (L, ≤ e ).

	Algorithm 1: Mining All The Most Specific Frequent Elementary Se-
	quence
	input : Sequential Database SDB, the minimum support threshold σ, ground
	sets D1,...,Dn
	output: The set MSFES of all most specific frequent elementary sequences,
	The poset L
	begin
	/* Step1: Generating the most general elementary sequence
	(T1, T2, ..., T k , {}, {}, ...., {})	*/
	for i ← 1 to n do
	e .add(T op(Di));
	M SF ES ← ∅ ;
	L ← ∅ ;
	L ← L ∪ e;
	/* Step2: The recursive generation of the new elementary sequence
	*/
	call get rec msf es(e,σ);
	Algorithm 2: Routine get rec msf es

input : Elementary sequence e, the minimum support threshold σ begin Cand ← {v ∈ desc(e) | supp(v ) σ} ; if Cand = ∅ then M SF ES ← M SF ES ∪ e; else foreach e ∈ Cand do if e ∈ L then L ← L ∪ e; call get rec msf es(e,σ);

Table 2 ,

 2 with (uh p , ca 2 , {mp 1 , mp 2 }) e (uh p , ca, {mp 1 , mp 2 }).

	-S 1 2 = (uh

p , ca 1 , {mp 1 }), none of the elementary sequences in MSFES set l , r 1 , {mp 3 }) is replaced by (gh l , r, {}) and (gh, r, {mp 3 }) from MSFES set.

Table 6 :

 6 Care trajectories of 3 patients

								Trajectory Length
	Number of patients
	1 3 5 7 9 11	14	17	20	23	26	29	32	35	38	41	44	47	50	53	56	59	63	69	76	99
								Trajectory Length

Table 7 :

 7 3-Care trajectories extracted by MMSIP

Programme de Médicalisation des Sytèmes d'Information.

http://illimine.cs.uiuc.edu/

" Classification Commune des Actes Médicaux ": the French classification of medical and surgical procedures.

Acknowledgments

The authors would like to thank the French institution "Institut National du Cancer (INCa)" for its financial support.


Loïc Andolfatto
email: loic.andolfatto@eads.net@lurpa.ens-cachan.fr

François Thiébaut

Claire Lartigue

Marc Douilly

Assisted Decision-Making for Assembly Technique Selection and Geometrical Tolerance Allocation

Keywords: Assembly process planning, assembly technique selection, geometrical tolerance allocation, multi-objective optimization, concurrent engineering

published or not. The documents may come

L'archive ouverte

Introduction 1.Design of an Assembly Process Plan

The assembly of large mechanical structures, such as aeronautical ones, may account for a large share of their total delivery cost. Boothroyd and Dewhurst stressed out the importance for manufacturing companies to assess a product's design by designing assembly process plan as soon as 1 in order to reach the maximum performance [START_REF] Boothroyd | Product design for manufacture and assembly[END_REF]. Computeraided assembly process planning has been the subject of many research works generally aiming at finding the minimum lead time and/or cost. The performance of the assembly process plans is evaluated according to several indicators such as, for example, tooling needs, reorientations of sub-assemblies, technological similarities in consecutive operations and so on [START_REF] Wang | Assembly process planning and its future in collaborative manufacturing: a review[END_REF].

Assembly plan evaluation

Assembly sequence

Assembly techniques

Component geometrical tolerances

Assembly system organisation

Boundary of the study

Fig. 1. Schematic view of the description and evaluation of an assembly process plan with the proposed approach boundary

As summarized in Fig. 1Schematic view of the description and evaluation of an assembly process plan with the proposed approach boundaryfigure.1.1, the description of an assembly process plan can be split into an assembly sequence, the assembly technique selected to make each attachment of the product, the geometrical tolerances allocated to each component and the description of the assembly system 2 . Most of the works presented in the literature only focus on specific aspects of the problem.

Assembly sequence planning, either defined as the part introduction ordering [START_REF] Bourjault | Contribution à une Approche Méthodologique de L'Assemblage Automatisé: Élaboration Automatique des Sequences Opératoires[END_REF][START_REF] De Fazio | Simplified generation of all mechanical assembly sequences[END_REF][START_REF] Dini | Generation of optimized assembly sequences using genetic algorithms[END_REF] or as the joint realization ordering [START_REF] Homem De Mello | A correct and complete algorithm for the generation of mechanical assembly sequences[END_REF][START_REF] Mantripragada | The datum flow chain: A systematic approach to assembly design and modeling[END_REF], has been widely studied in the literature [START_REF] Wang | Assembly process planning and its future in collaborative manufacturing: a review[END_REF]. The joint-based approach is close to the activity-based approach proposed by Cao and Sanderson [START_REF] Cao | Task decomposition and analysis of robotic assembly task plans using petri nets[END_REF], and reflects some of the issues considered for the organisation of the assembly system commonly treated in operational research [START_REF] Becker | A survey on problems and methods in generalized assembly line balancing[END_REF].

But the generation of assembly sequences and the organisation of the assembly system are not the only activities encapsulated in assembly process planning. Assembly technique3 selection [10] is seldom addressed in the literature even though it may have the greatest impact on production cost, according to Abdullah et al. [START_REF] Abdullah | A review of the support tools for the process of assembly method selection and assembly planning[END_REF].

Selecting assembly techniques implies setting the associated geometrical capabilities. The concern of geometrical quality of the assembled product (commonly called tolerancing) is therefore coupled with the assembly process planning. Up to this point, the geometrical tolerances allocated to the components to be assembled may be part of the information contained in a comprehensive assembly process plan [START_REF] Mantripragada | The datum flow chain: A systematic approach to assembly design and modeling[END_REF][START_REF] Camelio | Modeling variation propagation of multi-station assembly systems with compliant parts[END_REF][START_REF] Mounaud | Assembly sequence influence on geometric deviations propagation of compliant parts[END_REF].

Proposed Approach

Massive automation has been the key solution to answer the need for decreasing manufacturing cost and increasing delivery rate in many manufacturing domains. But in the field of aeronautical structure assembly, the complexity of some attachments, the high level of geometrical requirements and the high working area required are several hurdles to an efficient use of automation with respect to manual assembly. Trades-off have to be made.

This paper aims at proposing an original method to select assembly techniques and to allocate component geometrical tolerances in order to minimize the product cost and to maximize the conformity rate associated with the assembly plan. The assembly sequence is assumed to be already defined and the organisation of the assembly system is not considered. Section 2Optimization problemsection.1.2 details how an assembly process plan can be described according to a set of decision variables and how it can be evaluated according to several performance indicators. This formal view of the problem serves the definition and the resolution of a multi objective optimization problem.

A simple use case is presented in Sect. 3Use casesection.1.3 and the results obtained are detailed in Sect. 4Resultssection.1.4. Section 5Conclusionsection.1.5 concludes on the potential benefits of this approach in a concurrent engineering context.

Optimization problem 2.1 Multi-objective Optimization

Considering a given product and a given assembly sequence, the purpose of the presented work is to propose a method to find a set of assembly techniques and a set of geometrical tolerances that minimize the non-conformity rate and the delivery cost of a product at the same time.

A mathematical expression of this problem is given in (1Multi-objective Optimizationequation.1.2.1), where x is a decision variable vector representing an assembly plan, X is the set of feasible assembly plans and f is a function associating a fitness to an assembly plan. The construction of the decision variable vector x ∈ X is later described in Sect. 2.2Parametric Assembly Plansubsection.1.2.2.

inf{f (x) : x ∈ X} (1)
In the present case, the fitness function f must cope with several objectives (nonconformity and cost). Moreover, these objectives are likely to conflict with each other. A common solution to this issue is to combine the objectives in a single-valued fitness function. But it can prove complicated to model the actual objectives into a single performance indicator, which results in making trades-off a priori.

Solving a multi objective optimization problem seems to be more appropriate within the industrial context. From the assembly process planner point of view, an interesting point x * describes an assembly process plan for which one can not find a solution that provides both lower non-conformity rate and lower cost at the same time, as pictured in Fig. 2Non-dominated assembly process plans (circles) among a population (squares) considering the simultaneous minimization of the non-conformity rate and the minimization of the delivery costfigure.1.2.

Assuming f = [NCR(x), C(x)], non-conformity rate function and delivery cost function described in Sect. 2.3Assembly Plan Evaluationsubsection.1.2.3, a solution of the problem given in (1Multi-objective Optimizationequation.

1.2.1) is a point x * satisfying (2Multi-objective Optimizationequation.1.2.2), x * is called a non-dominated point. ∃ x ∈ X, NCR(x) < NCR(x * ) and C(x) < C(x * ) (2)
x *

x * Fig. 2. Non-dominated assembly process plans (circles) among a population (squares) considering the simultaneous minimization of the non-conformity rate and the minimization of the delivery cost

The underlying mathematical problem of this work consists in finding the set of non-dominated points (or a sufficient amount of points if this set is not finite). The assembly process planner finally has to select the plan that provides the best trade-off among the set of non-dominated points a posteriori. The high level, difficult-to-make, decision appears with the rough results in hand instead of during the modelling phase required to build a black-box optimizer. Modelling every single actual objective of the problem alone is likely to be an easier task.

Parametric Assembly Plan

Generalities. Running the optimization problem presented in (1Multi-objective Optimizationequation.1.2.1) requires to convert the assembly process plan description known in technical terms into a mathematical vector x called assembly plan vector in the paper.

A product is a set of components linked together through joints during the assembly process. Each component is subject to geometrical variations bounded by geometrical tolerances. Each joint is made using an assembly technique. The requirement on the assembly plan vector x can be deduced out of those three assertions. It must include elements representing the technique selected for each joint of the product and elements representing the components geometrical deviations.

Assembly techniques. Considering a product with N j joints, the N j first elements of x are dedicated to describe which assembly technique is associated to each of the N j joints. Assembly techniques are stored and indexed in a library with several attributes: index, associated list of assembly operations, list of costing information, geometrical capabilities, etc.

The technique assigned to the joint j is the one the index of which equals the value of the jth element of x, as pictured in Fig. 3Translation of an assembly process plan into a decision variable vector xfigure.1.3.

The assembly process planner has to define the list of techniques suitable to make each joint of the product. Some additional constraints can be set among those N j first elements: two joints can be forced to be made with the same technique for example. This reduces the set X of feasible assembly plans to the technically admissible ones. Technique chosen for the joint j ... ...

Assembly plan vector p i2 p i1

Fig. 3. Translation of an assembly process plan into a decision variable vector x x = ... Geometrical variations and geometrical tolerances. Assuming that the probability distribution of each geometrical variation is known, allocating tolerances amounts to setting these distribution parameters, e.g. the lower bound and upper bound for a uniform probability distribution or the mean value and the standard deviation for a normal distribution. This is illustrated in Fig. 3Translation of an assembly process plan into a decision variable vector xfigure.1.3. If the product's components have N gv geometrical variations, and each of them has k distribution parameters, then the assembly plan vector x also has N vp = Ngv i=1 k i elements to describe the geometrical tolerances allocated.

Intrinsic constraints exist among those N vp elements, such as a lower bound smaller than an upper bound for instance. Extrinsic constraints can also be declared by the user to define other technical limits, such as the minimum size of a tolerance zone. This also reduces the set X of feasible assembly plans.

Assembly Plan Evaluation

Non-Conformity Rate. The conformity of a product can be assessed by verifying that some of its characteristics -called key characteristics (KC) in [START_REF] Mantripragada | The datum flow chain: A systematic approach to assembly design and modeling[END_REF] -are kept within a requirement domain defined during the functional analysis. The Non-Conformity Rate of a product is defined as the probability for a product to have at least one of its KC outside of its requirement domain. In this work, only geometrical KC are considered. Tolerancing studies commonly provide a sensitivity matrix S to link the geometrical deviations δl4 to the KC deviations δKC as in (3Non-Conformity Rateequation.1.2.3) [START_REF] Camelio | Modeling variation propagation of multi-station assembly systems with compliant parts[END_REF][START_REF] Marguet | Geometrical behavior laws for computer aided tolerancing: Anatole a tool for structural assembly tolerance analysis[END_REF][START_REF] Mounaud | Assembly sequence influence on geometric deviations propagation of compliant parts[END_REF].

S • δl = δKC (3)
Samples of δl vectors can be associated to an assembly plan vector using both PDFs of the component geometrical variations and assembly technique capabilities. A Monte Carlo approach is finally used to associate the non-conformity rate N CR(x) to the assembly plan described by x.

Delivery Cost. The delivery cost can be split into two cost sources: the cost associated to the assembly operations and the cost resulting from the allocation of geometrical tolerances.

Each assembly technique is characterized by an interdependent list of assembly operations (an example is given in Table 3Details of the technique Positioning with a robot as tooltable. 1.3). Each operation has a fixed cost c f representing consumables and unitary tool wear for example. It also uses several resources (of hourly rate hr i ) during a certain amount of time t op . The elementary cost of an assembly operation is defined in (4Delivery Costequation.1.2.4):

c op = c f + t op × i hr i (4)
The total activity-based cost C op is the sum of the elementary cost of each operation, the list of which is derived from the list of selected assembly techniques given by x.

The non-recurring cost per product C nr , due to the acquisition cost of all the resources required for the assembly divided by the presumable amount of product to be assembled, is also included to the total delivery cost. As the list of the resources required to assemble the product depends on the assembly technique selection, C nr is a function of x.

The cost associated to a geometrical tolerance allocated is modelled by (5Delivery Costequation.1.2.5) (adapted from [START_REF] Chase | Least cost tolerance allocation for mechanical assemblies with automated process selection[END_REF]). T is the size of the tolerance interval allocated, T lim , a, b, m and k are function parameters identified according to experimental data.

c tol (T ) = a + b • e -m(T -T lim ) •(T -T lim ) -k (5)
The total cost C tol associated to the geometrical tolerances allocated to the product's components is the sum of all the elementary costs evaluated thanks to (5Delivery Costequation.1.2.5) according to the tolerances described by the vector x.

Finally, the sum of C op , C nr and C tol provides the delivery cost C(x) associated to an assembly plan, as depicted in Fig. 4Evaluation scheme and fitness representation of an assembly plan defined by its mathematical representation xfigure.1.4.

3 Use case 3.1 Product to be assembled Assembly sequence: tj1, tj3, tj2, tj4, tj5, J1-J2, tj6, tj7, J3-J4, J5, J6 Fig. 6. Temporary components used during assembly (referred to as Tool k), joints of the product and assembly sequence

KC 1 KC 2 KC 3 KC5 KC6 3 -Splice 4 -Cover KC 4 L 1 L 2 L3 L 2 1 -L_Flank 2 -R_Flank
The method described in the previous section is applied to the assembly of a simple mechanical structure composed of four components (see Fig. 5Details of the product's key characteristics (KC i ), components and their characteristics considered as subject to geometrical variation (L j )figure.1.5). Three additional temporary components are also used during the assembly, as exposed in Fig. 6Temporary components used during assembly (referred to as Tool k), joints of the product and assembly sequencefigure.1.6 in which the assembly sequence is also given (as the order in which the joints are made).

The geometrical variation propagation problem is reduced to a one-dimensional study with six key characteristics and four component's dimensions subject to geometrical variations. The capabilities of the technique selected for the joints tj3 to tj7 are also impacting the conformity rate as these joints are positioning joints involved in the KC values. Deviations on KC 1 and KC 2 must be kept within ±0.6 mm. The value for KC 4 to KC 6 is ±0.3 mm.

The parameters of the cost vs. tolerance law defined by (5Delivery Costequation.1.2.5) are given in Table 1Parameters of the Cost vs. Tolerance law defined by (5Delivery Costequation.1.2.5)table.1.1. The sensitivity matrix S defined in (3Non-Conformity Rateequation.1.2.3) is not detailed in the paper.

Assembly technique library

The assembly technique library is not entirely described in this paper but the Table 2List of the assembly techniques available associated to the joints they can be used fortable.1.2 presents the techniques considered. The Table 3Details of the technique Positioning with a robot as tooltable.1.3 details the information stored in the library for the Positioning with a robot as tool technique. In addition to the feasibility constraints detailed in Table 2List of the assembly techniques available associated to the joints they can be used fortable.1.2, the joints tj3 to tj5 must be done with the same technique. So do the joints tj6 and tj7.

Results

The optimization problem defined in (1Multi-objective Optimizationequation.1.2.1) is solved using the Non-Sorting Genetic Algorithm II (implemented in the Inspyred Python library [START_REF]Inspyred 1.0: Bio-inspired Algorithms in Python[END_REF]) with a population composed of 200 individuals evolving during 20 generations. The sensitivity matrix S is obtained using the tolerance analysis software AnaTole [START_REF] Marguet | Geometrical behavior laws for computer aided tolerancing: Anatole a tool for structural assembly tolerance analysis[END_REF] and the non-conformity rate is evaluated with a Monte Carlo method implemented in OpenTURNS [START_REF]TURNS -Reference Guide -version 1.0[END_REF].

The non-dominated points obtained are displayed in Fig. 7Non-dominated points obtained after optimization, with four separate zones of iso-technique and varying allocated tolerancesfigure.1.7. The assembly process plan corresponding to the square point U(-0.17, 0.16) J4 and J6 Drilling with grid and Manual Fastening

Conclusion

Selecting assembly techniques and allocating geometrical tolerances for a product requires making quality vs. cost trades-off. The method proposed in this paper helps the assembly process planner to find a set of good 5 solutions among which he can select the one that suits the best his interests.

The method is based on an assembly technique library in which the company knowhow can be stored and upon a geometrical variation propagation relation associated to the assembly sequence. A wide range of potential assembly process plans can therefore be investigated and evaluated from an objective point of view. Results can be obtained with various assembly sequence scenarios and even with various product architecture scenarios. Decisions about the product's design and its assembly process plan can be taken from a point of view shared by product designers and manufacturers, enhancing collaborative and concurrent engineering. 5 The non-dominated points, from the mathematical point of view.

Fig. 4 .

 4 Fig. 4. Evaluation scheme and fitness representation of an assembly plan defined by its mathematical representation x.

Fig. 5 .

 5 Fig.5. Details of the product's key characteristics (KCi), components and their characteristics considered as subject to geometrical variation (Lj)

Fig. 7 .

 7 Fig.7. Non-dominated points obtained after optimization, with four separate zones of isotechnique and varying allocated tolerances.

Table 4 .

 4 Assembly process plan corresponding to the square point in Fig.7Non-dominatedpoints obtained after optimization, with four separate zones of iso-technique and varying allocated tolerancesfigure.1.7 Joint Assembly technique selected Dimension Allocated tolerance tj1 and tj2 Back-to-back positioning L1 U(-0.40, 0.30) tj3 to tj7 Positioning with robot as a tool L2 U(-0.40, 0.30) J1 and J2 Traditional bonding L3 U(-0.22, 0.21) J3 and J5 Back-to-back positioning with splitting L4

Table 1 .

 1 Parameters of the Cost vs. Tolerance law defined by (5Delivery Costequation.1.2.5)

		T lim	a	b	m	k
		(mm)	(cost unit) (cost unit / mm) (mm -1 )
	L1 and L2	0.01	0	200	1	1
	L3 and L4	0.01	0	200	1	1

What is called design for manufacturing and assembly (DFMA).

Either spatial or temporal.

Also called assembly process by some authors.

Gathering both geometrical deviations of the components and deviations due to the assembly techniques employed.

Acknowledgment. This research work has been carried out in partnership between EADS Innovation Works and the LURPA-ENS de Cachan in the frame of the GRC-Flexible Assembly of the INNO'CAMPUS program.

in Fig. 7Non-dominated points obtained after optimization, with four separate zones of iso-technique and varying allocated tolerancesfigure.1.7 is detailed in Table 4Assembly process plan corresponding to the square point in Fig. 7Non-dominated points obtained after optimization, with four separate zones of iso-technique and varying allocated tolerancesfigure.1.7table.1.4.

The results illustrate the interest of a multi objective optimization. Considering a single objective for the fitness and the other one as constraint would give an arbitrary boundary between acceptable and non-acceptable solution, leading to an optimal point not necessarily better than other ones. Here, a decision-making team can adapt the final choice with more information in hand.

A deeper analysis of the results displayed in Fig. 7Non-dominated points obtained after optimization, with four separate zones of iso-technique and varying allocated tolerancesfigure.1.7 shows that the points can be classified into four zones for this use case. In each zone, the assembly techniques selected are identical and only tolerances allocated are varying. It is therefore possible to identify the most relevant set of assembly techniques early during the design of the product and to refine the tolerance allocation along its design process.


Audrey Repetti

Emilie Chouzenoux

Jean-Christophe Pesquet

A NONCONVEX REGULARIZED APPROACH FOR PHASE RETRIEVAL

Keywords: Phase retrieval problem, Nonconvex optimization, Nonsmooth optimization, Proximal methods

With the development of new imaging systems delivering large-size data sets, phase retrieval has become recently the focus of much attention. The problem is especially challenging due to its intrinsically nonconvex formulation. In addition, the applicability of many existing solutions may be limited either by their estimation performance or by their computational cost, especially in the case of non-Fourier measurements. In this paper, we propose a novel phase retrieval approach, which is based on a smooth nonconvex approximation of the standard data fidelity term. In addition, the proposed method allows us to employ a wide range of convex separable regularization functions. The optimization process is performed by a block coordinate proximal algorithm which is amenable to solving large-scale problems. An application of this algorithm to an image reconstruction problem shows that it may be very competitive with respect to state-of-the-art methods.

INTRODUCTION

Phase retrieval is one of the oldest while still most challenging problems encountered in imaging. It basically consists of producing an estimate v of an original multivariate signal v from (possibly noisy) magnitude measurements |Hv| where H is a complex-valued linear operator. Such a problem plays a central role in numerous application fields of image processing, in particular in crystallography [START_REF] Harrison | Phase problem in crystallography[END_REF], optical imaging [START_REF] Walther | The question of phase retrieval in optics[END_REF], phase contrast tomography [START_REF] Bauschke | A new generation of iterative transform algorithms for phase contrast tomography[END_REF], and coherent diffraction imaging [START_REF] Shechtman | Efficient coherent diffractive imaging for sparsely varying objects[END_REF].

The most popular methods to estimate v are probably Gerchberg-Saxton (GS) algorithm [START_REF] Gerchberg | A practical algorithm for the determination of phase from image and diffraction plane pictures[END_REF] and its relaxed version Fienup algorithm [START_REF] Fienup | Phase retrieval algorithms : A comparison[END_REF]. Initially introduced for the case when H is a Fourier transform matrix, these algorithms alternate projections onto the range of H and onto the nonconvex set of vectors with magnitude equal to |Hv|. Note that, This work was supported by the MASTODONS project by CNRS (grant 2013MesureHD) since the latter constraint set is nonconvex, these algorithms do not benefit from the convergence guarantees of the famous Projections Onto Convex Set (POCS) algorithm [START_REF] Bauschke | Phase retrieval, error reduction algorithm, and Fienup variants : a view from convex optimization[END_REF]. Convex relaxations of the phase retrieval problem based on Semi-Definite-Programming (SDP) formulations were proposed in [START_REF] Candès | Phaselift : Exact and stable signal recovery from magnitude measurements via convex programming[END_REF][START_REF] Waldspurger | Phase recovery, maxcut and complex semidefinite programming[END_REF], giving rise to the PhaseLift and PhaseCut algorithms, respectively. Both methods have been shown to outperform GS and Fienup, in a certain number of application examples [START_REF] Fogel | Phase retrieval for imaging problems[END_REF].

All the aforementioned methods tend to be sensitive to noise and/or to be less effective in the under-determined case [START_REF] Shechtman | GESPAR : Efficient phase retrieval of sparse signals[END_REF]. Such difficulties are mainly due to the ill-posed character of the phase retrieval problem and they can be alleviated by incorporating some a priori information in the reconstruction process. Several algorithms based on SDP relaxation [START_REF] Jaganathan | Recovery of sparse 1-D signals from the magnitudes of their Fourier transform[END_REF], alternating projection [START_REF] Mukherjee | An iterative algorithm for phase retrieval with sparsity constraints : application to frequency domain optical coherence tomography[END_REF] and greedy pursuit [START_REF] Shechtman | GESPAR : Efficient phase retrieval of sparse signals[END_REF], have been recently introduced, which aim at solving the phase retrieval problem under the prior assumption that the sought signal v has a sparse representation in a (possibly redundant) dictionary. However, up to the best of our knowledge, these approaches may become very computationally intensive as the problem dimension increases. This is particularly critical for SDP and alternating projections methods, when the pseudoinverse of H does not have a closed form, as it usually happens for non-Fourier measurements, or for the greedy method in [START_REF] Shechtman | GESPAR : Efficient phase retrieval of sparse signals[END_REF] when the sparsity degree of the data is not sufficient.

In this paper, we introduce a new phase retrieval strategy relying on the minimization of a penalized criterion. This criterion is composed of (i) a difference of convex functions (DC function) achieving a smooth approximation of the standard nonconvex nonsmooth least-squares data fidelity term, and (ii) a convex non necessarily smooth block separable penalization term. The latter function can be chosen in a flexible manner. For instance, it can be used to incorporate a sparsity assumption on the solution. The resulting nonconvex nonsmooth optimization problem is then solved with the Block Coordinate Variable Metric Forward-Backward (BC-VMFB) method that we have recently introduced in [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF]. This algorithm is ensured to converge to a critical point of the criterion and it takes advantage of the block separable structure of the regularization term. The proposed method thus benefits from a low computational cost which makes it applicable to large-scale problems.

The rest of the paper is organized as follows : Section 2 presents the phase retrieval problem and introduces a new objective function. Section 3 describes the proposed BC-VMFB algorithm and emphasizes its convergence properties. Finally, in Section 4, the good performance of our approach is shown via experiments.

PROBLEM FORMULATION

Observation model

Let v ∈ R M be the original unknown signal. Given H ∈ C S×M an observation matrix with complex-valued elements, the measurements z ∈ [0, +∞) S are assumed to be related to the original signal through the model

z = (z (s) ) 1≤s≤S =|Hv| + w, (1)
where | • | denotes the componentwise modulus operator, and w ∈ [0, +∞) S is a realization of an additive noise.

It is important to note that Model (1) encompasses the case when the target signal is complex since the observed vector can then be expressed under the form

z = [H R + i H I | -H I + i H R ] v R v I + w, (2)
with i 2 = -1 and (•) R (resp. (•) I ) designates the real (resp. imaginary) part of its argument.

Minimization problem

Let W ∈ R M ×N , M ≤ N , be a given frame synthesis operator (e.g. a possibly redundant wavelet representation) [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] such that v = W x. Following a synthesis approach, the estimated signal reads v = W x where the frame coefficient vector x ∈ R N is estimated by minimizing the sum of a data fidelity term F and a regularization function R, i.e.

G = F + R. (3)
In the context of phase retrieval problems, a usual way to define the data fidelity term F is to consider the nonconvex nonsmooth least-squares criterion [START_REF] Fienup | Phase retrieval algorithms : A comparison[END_REF][START_REF] Shechtman | GESPAR : Efficient phase retrieval of sparse signals[END_REF] :

(∀x ∈ R N ) F (x) = S s=1 ϕ (s) ([HW x] (s) ), (4)
where, for every s ∈ {1, . . . , S}, ϕ

(s) = 1 2 (| • | -z (s) ) 2 .
We propose to replace the latter nonsmooth function by a smooth approximation of it, easier to handle, which is defined as the difference of two convex terms, parametrized by a constant δ > 0 :

(∀u ∈ C) ϕ (s) (u) = 1 2 |u| 2 + (z (s) ) 2 -z (s) |u| 2 + δ 2 1/2 .
It can be noticed that, in the limit case when δ = 0, the standard nonsmooth fidelity function is recovered.

Regarding the penalization term, we will focus on the case when R is a block separable function. More precisely, let us define (J j ) 1≤j≤J a partition of the frame coefficient index set {1, . . . , N } into J ≥ 2 subsets with nonzero block-sizes (N j ) 1≤j≤J . In the following, we assume that the regularization function can be expressed as

R(x) = J j=1 R j (x (j) ), (5)
where, for every j ∈ {1, . . . , J}, function R j : R Nj → (-∞, +∞] is proper, lower semicontinuous, convex and continuous on its domain, and x (j) ∈ R Nj denotes the j-th block of components of x.

PROPOSED ALTERNATING OPTIMIZATION METHOD

Optimization tools

We first introduce the notion of proximity operator which will be useful to present our algorithm. Let ψ : R N → (-∞, +∞] be a convex, proper, lower semicontinuous function and let U ∈ R N ×N be a symmetric positive definite (SPD) matrix. The proximity operator of ψ at ξ ∈ R N relative to the metric induced by U ∈ R N ×N is defined as ( [16, Sec. XV.4] and [START_REF] Combettes | Variable metric quasi-Fejér monotonicity[END_REF])

prox U ,ψ (ξ) = argmin ζ∈R N ψ(ζ) + 1 2 ζ -ξ 2 U .
Hereabove, • U denotes the weighted norm defined by

• U = •, U • 1/2
, where •, • is the usual Euclidean scalar product. Note that, if U is the identity matrix, then the usual proximity operator introduced in the pioneering work [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF] is recovered.

Proposed algorithm

Although the proposed data fidelity term ( 4) is a nonconvex function, it presents the advantage of being Lipschitz differentiable. Such a property suggests the use of the Forward-Backward (FB) scheme [START_REF] Tseng | A modified forward-backward splitting method for maximal monotone mappings[END_REF] for minimizing the objective function [START_REF] Bauschke | A new generation of iterative transform algorithms for phase contrast tomography[END_REF]. A major advantage of FB algorithm is that it benefits from proven theoretical convergence guarantees, even in the nonconvex case [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems : proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel methods[END_REF][START_REF] Attouch | On the convergence of the proximal algorithm for nonsmooth functions involving analytic features[END_REF]. However, in practice, this algorithm may suffer from slow convergence so that accelerated versions have been proposed in [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF][START_REF] Chouzenoux | Variable metric forward-backward algorithm for minimizing the sum of a differentiable function and a convex function[END_REF][START_REF] Ochs | iPiano : inertial proximal algorithm for non-convex optimization[END_REF][START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF][START_REF] Tran-Dinh | Composite self-concordant minimization[END_REF][START_REF] Combettes | Variable metric forwardbackward splitting with applications to monotone inclusions in duality[END_REF], still with convergence guarantees in the nonconvex case ( [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF][START_REF] Chouzenoux | Variable metric forward-backward algorithm for minimizing the sum of a differentiable function and a convex function[END_REF][START_REF] Ochs | iPiano : inertial proximal algorithm for non-convex optimization[END_REF][START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]). In this paper, the block separable structure of R is accounted for by adopting the following BC-VMFB method from [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF]. This algorithm combines the advantages of a preconditioning strategy [START_REF] Chouzenoux | Variable metric forward-backward algorithm for minimizing the sum of a differentiable function and a convex function[END_REF] with an alternating minimization scheme :

Algorithm 1 BC-VMFB algorithm.
For every ∈ N, let γ ∈ (0, +∞) and A j (x ) ∈ R Nj ×Nj be a SPD preconditioning matrix. Initialize with x 0 ∈ R N . Iterations :

For = 0, 1, . . .       
Select a block index j ∈ {1, . . . , J}.

x(j ) = x (j ) -γ A j (x ) -1 ∇ j F (x ), x (j ) +1 = prox γ -1 Aj (x ),Rj x(j ) , x ( ) +1 = x ( ) .
In the above algorithm, at each iteration ∈ N, ∇ j F (x ) ∈ R Nj denotes the partial gradient of F with respect to block x (j ) computed at x , and  is the complementary set of j on {1, . . . , J}, i.e.  = {1, . . . , J} \ {j }.

Note that Algorithm 1 is equivalent to the algorithm in [START_REF] Chouzenoux | Variable metric forward-backward algorithm for minimizing the sum of a differentiable function and a convex function[END_REF] for N j ≡ N and reduces to the standard FB algorithm if, additionally, A j (x ) reduces to the identity matrix. The convergence of BC-VMFB requires the design of a proper strategy to adjust the preconditioning matrices (A j (x )) ∈N in an automatic manner. This will be discussed in the next section.

Choice of the preconditioning matrices

At iteration ∈ N, let j be the chosen index in {1, . . . , J} and let x be the -th iterate generated by Algorithm 1. Then, we propose to choose matrix A j (x ) by following the Majorize-Minimize principle [START_REF] Sotthivirat | Image recovery using partitioned-separable paraboloidal surrogate coordinate ascent algorithms[END_REF], i.e. it is such that the quadratic function

Q j (•, x ) = F (x ) + (• -x (j ) ) ∇ j F (x ) + 1 2 • -x (j ) 2 Aj (x ) ,
majorizes the restriction of F to its j -th block on R Nj . We thus have, for every y ∈ R Nj , 1) , . . . , x (j -1) , y, x

F (x (
(j +1) , . . . , x (J) ) ≤ Q j (y, x ).
According to our result in [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF]Prop.4.1], such a majorization condition is fulfilled by the diagonal matrix

A j (x ) = Diag Ω j 1 S , (6)
where 1 S is the unit vector on R S and the elements of Ω j ∈ R S×Nj are given, for every (s, n) ∈ {1, . . . , S} × J j , by

Ω (s,n) = |[HW ] (s,n) R | N n =1 |[HW ] (s,n ) R | + |[HW ] (s,n) I | N n =1 |[HW ] (s,n ) I |. (7)

Convergence result

One appealing feature of Algorithm 1 is that its convergence is guaranteed by the following result deduced from [14, Thm. 3.1] : Theorem 1. Let (x ) ∈N be a sequence generated by Algorithm 1 where (A j (x )) ∈N are defined as in [START_REF] Fienup | Phase retrieval algorithms : A comparison[END_REF]. Assume that (i) G is a coercive function, i.e. lim x →+∞ G(x) = +∞.

(ii) R is a semi-algebraic function. 1(iii) There exists (γ, γ) ∈ (0, +∞) 2 such that, for every ∈ N, γ ≤ γ ≤ 2 -γ. (iv) The blocks are updated according to an essentially cyclic rule, i.e. there exists a constant K ≥ J such that, for every ∈ N, {1, . . . , J} ⊂ {j , . . . , j +K-1 }.

Then, the sequence (x ) ∈N converges to a critical point x of (3). Moreover, G(x ) ∈N monotonically converges to G( x).

EXPERIMENTAL RESULTS

Observation model

We now demonstrate the practical performance of our algorithm in an image reconstruction problem. Starting from a complex-valued original image v ∈ C M (see Fig. 1) with size M = 128 × 128, we generate noisy observations z ∈ R S according to Model (2). The observation matrix H ∈ C S×M results from the composition of two degradation operators, namely a projection matrix modeling S = 23400 Radon projections from 128 parallel acquisition lines and 180 angles regularly distributed on [0, π), and a complex-valued blur operator. For the latter, a 1D convolution kernel with length 3, acting separately on the acquired sinogram for each angle has been used. It can be noticed that such a model is reminiscent of the linearized version of the phase contrast tomography problem addressed in [START_REF] Guigay | Mixed transfer function and transport of intensity approach for phase retrieval in the Fresnel region[END_REF][START_REF] Davidoiu | Nonlinear phase retrieval using projection operator and iterative wavelet thresholding[END_REF] in the 3D case.

Regularization term

We choose the synthesis frame operator W so that each vector x = (x (n) ) 1≤n≤N ∈ R N , N = 8M , is the concatenation of two vectors of R 4M , corresponding to an overcomplete Haar decomposition of v R (resp. v I ) for one resolution level. An example of frame decomposition of the real part of the original image is displayed in Fig. 2.

The regularization function R is defined as the sum, for p ∈ {1, . . . , 4M } of terms of the form : mean value parameter on the sought frame coefficient u (p) . In our experiments, we take ω p = (0.4, 0.6) for the approximation part, and ω p = 0 otherwise. Thus, the resulting criterion G satisfies Assumptions (i)-(iv) of Theorem 1. Note that the proximity operators involved in Algorithm 1 are easy to compute since they have an explicit form [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF].

Simulation results

The magnitude measurement vector is corrupted with an additive real-valued white zero-mean Gaussian noise with variance equal to 0.025 (truncated so as to guarantee the nonnegativity of the observations). In our experiments, parameters ϑ a , ϑ d and δ are adjusted so as to maximize the signalto-noise ratio (SNR) between the original image v and the reconstructed one v. In Algorithm (1), for every j, x (j) corresponds to a vector of given fixed size 8Q which gathers 8 blocks extracted from the approximation and detail subbands of both real and imaginary parts of the image, as depicted in Fig. 2. At each iteration ∈ N, j is then randomly chosen so that each block is updated at least once per J iterations. In practice, the reconstruction time was observed to vary significantly with respect to the block-size parameter Q. In this example, the best choice in terms of convergence speed is reached for Q = 2. Fig. 1 shows the reconstructed image obtained with our algorithm using this optimal blocksize and γ ≡ 1.9. We also present in Fig. 1 the reconstruction result obtained using an enhanced version of the regularized alternating projection approach from [START_REF] Mukherjee | An iterative algorithm for phase retrieval with sparsity constraints : application to frequency domain optical coherence tomography[END_REF], where the non-inversibility of matrix H has been taken into account. Note that due to the large size of the data, it appeared impossible to run SDP approaches from [START_REF] Candès | Phaselift : Exact and stable signal recovery from magnitude measurements via convex programming[END_REF][START_REF] Waldspurger | Phase recovery, maxcut and complex semidefinite programming[END_REF][START_REF] Jaganathan | Recovery of sparse 1-D signals from the magnitudes of their Fourier transform[END_REF], neither the greedy method from [START_REF] Shechtman | GESPAR : Efficient phase retrieval of sparse signals[END_REF]. Moreover, we obtained poor results using the standard GS algorithm.

Fig. 3 illustrates the decay of ( xx / x ) (in log scale) using the proposed algorithm or the non preconditioned variant from [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]. The optimal solution x has been precomputed for each algorithm, using a large number of iterations. One can observe that the variable metric strategy leads to a significant acceleration of the convergence. It must be noticed that the method from [START_REF] Mukherjee | An iterative algorithm for phase retrieval with sparsity constraints : application to frequency domain optical coherence tomography[END_REF] is almost 10 times slower to reach a stabilized SNR value.

Fig. 1 .

 1 Fig. 1. (top) Real and imaginary parts of the original image v and estimated images v using either (middle) Algorithm 1, SNR = 21.27 dB or (bottom) the regularized alternating projection method from [13], SNR = 14.45 dB.

Fig. 2 .

 2 Fig. 2. (left) Example of frame decomposition of the original image. (right) Indices of a block x (j) for Q = 32. Only the coefficients corresponding to the real part of the image are displayed.

  p ∈ E and u (p) = 0, +∞ otherwise, where u (p) ∈ R 2 is the p-th pair of frame coefficients corresponding to the real and imaginary part of the image, and E corresponds to the object background. The sparsity of the detail part of x is promoted by setting κ p = 1 and ϑ p = ϑ d ∈ (0, +∞) for the detail subbands of the decomposition. For the approximation coefficients, we choose κ p = 2 and ϑ p = ϑ a ∈ (0, +∞). Parameter ω p ∈ R 2 corresponds to a

Fig. 3 .

 3 Fig. 3. Convergence profile of Algorithm 1 (solid line) and the non preconditioned variant from [24] (dashed line), Intel(R) Core(TM) i7-3520M @ 2.9GHz using Matlab 7.

Semi-algebraicity is a property satisfied by a wide class of functions, which means that their graph is a finite union of sets defined by a finite number of polynomial inequalities.