Matthias Roth

Jean-Jacques Lesage

Lothar Litza

The concept of residuals for fault localization in discrete event systems

Keywords:

In this paper an approach for fault localization in closed-loop Discrete Event Systems is proposed. The presented diagnosis method allows fault localization using a fault-free system model to describe the expected system behavior. Via a systematic comparison of the observed and the expected behavior, a fault can be detected and a set of fault candidates is determined. Inspired by residuals known from diagnosis in continuous systems, different set operations are introduced to generate the fault candidate set. After fault detection and a first fault localization. a procedure is given to render the fault localization more precisely by an analysis of the further observed system behavior. Special emphasis is given to the use of identified models for the fault-free system behavior. The approach is explained using a laboratory manufacturing facility.

t. Introduction

Diagnosis in Discrete Event Systems (DES) has been a vital research area in the last 15 years. The class of technical systems that can be modeled as DES includes production, transportation and communication systems [START_REF] Cassandras | Introduction to discrete event systems[END_REF]. One of the main purposes of diagnosis methods for technical systems is to increase their availability by helping the system operator to find fault sources as quickly as possible. An increased availability leads to an improved system dependability and has positive effects on economic key issues such as productivity.

An important class of DES diagnosis approaches is model based. The idea behind is to compare the modeled and the observed system behavior in order to detect and to localize faults. Model-based approaches can be divided in two groups. The first group considers models that contain fault-free behavior as well as system behavior for given faults. The second possibility uses models of fault-free system behavior only.

A prominent example for DES-approaches with automaton models including the faulty system behavior is given in Sampath, [START_REF] Sam Path | Failure diagnosis using discrete-event models[END_REF]. One of the main features ofthis "diagnoser" approach is the possibility to give guarantees concerning the diagnosability of faults that are considered in the underlying model. If certain conditions hold, faults considered in the model can be precisely localized. Various improvements of this approach have been proposed in the last E-mail addresses:mroth@eit.uni-kl.de (M. Roth), lesage@lurpa.ens-cachanJr U.-J. Lesage), Iitz@eit.uni-kl.de (I.. Litz), years mainly dealing with adapting the method for distributed diagnosis e.g. [START_REF] Sayed-Mouchaweh | Decentra lized diagnosis based on Boolean discrete event models: Application on manufacturing systems[END_REF] or incorporating timing information e.g. Hashtrudi Zad, [START_REF] Hashtrudi Zad | Fault diagnosis in discrete-event systems: Incorporating timing information[END_REF].

Apart from using automata as system models for DES, Petri nets are a widely spread modeling formalism. In [START_REF] Fanti | Fault diagnosis and identification of discrete event systems using Petri nets[END_REF] an overview of existing diagnosis approaches using Petri nets is given. In the most recent works it is proposed to perform diagnosis based on integer linear programming, see e.g. [START_REF] Dotoli | On-line fault detection in discrete event systems by Petri nets and integer linear program ming[END_REF]. The idea is to check if an observed behavior is consistent with some modeled faulty beha vior. An inherent disadvantage of approaches relying on fault models is that only faults explicitly considered in the system model can be detected and localized.

Diagnosis methods only working with a model of the nominal fault-free system behavior avoid this disadvantage. Faults that lead to an inconsistency between the nominal and the observed behavior can be detected when comparing modeled behavior and measure ments without having to be considered explicitly (Cordier et aI., 2004;[START_REF] Reiter | A theory of diagnosis from first principles[END_REF]. Another advantage of this class of methods is that model building is straightforward since no special knowledge of faulty system behavior is necessary. A drawback of this second class of approaches is that fault localization is more ambitious than in the case of methods relying on fault models, since the models have been built using less knowledge. A second drawback is that diagnosability of given faults can usually not be guaranteed.

In this paper, fault localization is studied based on a fault-free system model. The approach has been developed for a widely used class of closed-loop DES consisting of a plant and a controller. Although the proposed method is not restricted to identified models, special emphasis will be given to the implications of models obtained by identification. The paper is structured as follows: In Section 2 an overview of the proposed diagnosis approach is given. Section 3 sketches the identification method that has been used to obtain the fault-free system model. It also defines the semantics of the used model class. Section 4 represents the kernel of the paper: It is explained how to detect faults and to define appropriate residuals for fault localization using the model defined in Section 3. In order to show the relevance of the proposed method for real systems, a case study is treated in Section 5.

Overview of the proposed diagnosis approach

According to [START_REF] Isermann | Trends in the application of model based fault detection and diagnosis of technical processes[END_REF] the term fault diag nosis is defined as follows.

Diagnosis: Determination of the time of detection, kind, size and location of a fault. Diagnosis follows detection, includes fault isolation and identification.

The definition of [START_REF] Isermann | Trends in the application of model based fault detection and diagnosis of technical processes[END_REF] is proposed with the background of time continuous systems. Since this paper deals with a different class of systems and different models, not each aspect of this definition applies to the presented approach. Systems that are considered in this paper are DES consisting of a closed-loop of plant and controller (Fig. 1). The controller has a certain number of binary outputs connected with actuators in the plant. Sensors of the plant are connected to binary controller inputs. Possible faults in this scenario are broken actuators and sensors as well as damaged plant hardware. In the presented approach the only information used to perform the diagnosis task are the signals exchanged between controller and plant. The aim is to detect if a sensor or an actuator fails or if some plant hardware gets damaged which leads to an abnormal behavior. Recognizing such an abnormal behavior will be referred to as fault detection. In order to allow systematic repair actions, it is also helpful to give information of where a detected fault is located. This process will be referred to as fault localization. The aim of fault localization in the context of this paper is to determine controller I/Os (inputs and outputs) which are related to a faulty component. Based on the binary signals exchanged between controller and plant it is not possible to gather any information concerning the "size" of a fault. Hence, in this paper the diagnosis task is considered to be completed after fault localization.

A model-based diagnosis method that is well established for continuous systems is to use residuals. [START_REF] Isermann | Trends in the application of model based fault detection and diagnosis of technical processes[END_REF] define residuals as follows:

A residual is a fault indicator, based on a deviation between measurements and model-based computation.

In the upper part of Fig. 2 this principle is depicted for continuous systems. The deviation between a measured signal evolution (dashed line) and the signal evolution computed using a model (continuous line) is quantified. Analyzing this quantification it is possible to detect and localize faults in the observed system. Fig. 2 also shows measurements of a DES and an automaton that models the con sidered system. Typically, observations of a DES are event sequences as shown in the lower part of the figure. The expected event trace can be derived from the automaton modeling the system. The observed sequence a, b can be reproduced by the state trajectory 0, 4, 5. When event c is observed, the sequence is no longer reproducible with the given automaton and a fault can be detected. The approach presented in this paper focuses on the analysis of the difference between measured and expected DES behavior. Two generic fault symptoms will be considered in Section 4: Faults leading to observed but unexpected events and faults that lead to missed events. The proposed fault detection and localization principle for closed loop DES is depicted in Fig. 3. The closed-loop DES is observed by collecting the signals exchanged between controller and plant. This measured system output is given to an evaluator which contains an automaton modeling the fault free system behavior. The evaluator tries to reproduce the observed system output using the embedded model. This evaluation creates the model output that is compared to the measured output. If measured and modeled output differ, a fault is detected and fault localization starts in order to determine where the fault is located.

Model identification

Data collection and arrangement

As explained in the former section. model-based diagnosis is driven by the comparison of modeled and observed behavior. The closer the model to the real behavior is. the more accurate results model-based diagnosis can yield. The usual way for manual building of systems consisting of controller and plant (see Fig. 1) is to build models for the plant GPlanc and for the controller algorithm GConrcollec (Sampath et al., 1996). The system model G is obtained by parallel composition of the plant and controller models: G = GPlancllGConrcoller' Usually. GConrcoller only models the idealized behavior of the con troller algorithm. Hardware specific aspects of the controller are usually not considered. An example for such an aspect are industrial controllers where several controller outputs can change their value within one controller cycle. An example for effects usually not considered during manual model building of the plant model GPlanc are chattering sensors. As a consequence. manually built models often represent the specified system behavior rather than the real one. If the real behavior is compared with an idealized model. the results can either be false alerts or non-detectable faults.

The main advantage of working with an identified model is that it represents the real behavior of the closed loop system with respect to the identification data base. As identification data base. it is possible to use behavior observed during normal system operation. This behavior includes consequences implied by the specific controller hardware as well as non-ideal (but not faulty) plant behaviors like chattering sensors.

For model identification and for diagnosis itself it is necessary to collect the signals exchanged between controller and plant when the system is in operation. The identification ofthe model is based on sequences of observed I/O vectors. For the identification it is necessary to observe a certain number of system evolutions and to collect the according I/O vectors. The sequence of I/O vectors that is exhibited during a given system evolution is built by I/O vectors in the order of their appearance: Definition 2 (//0 vector sequence). If during the h-th system evolution lh I/O vectors Uh have been observed. the sequence is denoted as a(h) = (Uh(1).Uh(2) •...•Uh(lh))'

Hence lh denotes the length of the h-th I/O vector sequence. We assume that for two successive I/O vectors u(t) #-u(t-1) holds. In order to translate the observed I/O vector sequences into a model. a formal definition of the observed system behavior is necessary. We define the set of observed words with length q observed during p different system evolutions: Definition 3 (Observed word set and language). The observed words of length q are denoted as

W6bs = _. •U j (j+ q-1))) .
i~1 (/;:Q\Ui(j),Ui(j+ 1)•.

With the observed word set we can define the observed language of length n of the system as

n L3bs = U Whbs' j=1
The observed language of length n consists of the I/O vector sequences up to length n that have been observed to get the necessary data base for identification.

Model definition

Since the observed system language has to be translated into an appropriate DES model, we have to define an automaton that reflects the characteristics ofthe closed-loop DES. The closed-loop DES from Fig. 1 is an autonomous event generator that exhibits a system output. A well programmed controller can be considered as deterministic, whereas the physical plant must generally be considered as non-deterministic. Hence the coupled system of plant and controller must be considered as non-deterministic. An automaton that reflects these properties is the Non-Deterministic Autonomous Automaton with Output that is defined by Definition 4 (NDAAO). NDAAO = (X,QJ).•xo) with X finite set of states. Q output alphabet. f :X---> 2 x non-deterministic transition function. A. : X---> Q output function and Xo the initial state. This automaton will be presented as a digraph in the usual manner [START_REF] Cassandras | Introduction to discrete event systems[END_REF]. The words and language generated by this automaton are given by Definition 5 (Words and Language of the NDAAO). The set of words of length n generated from a state x(i) is defined as W;(.) I = (we QI : w = A.(x(i))} and W;(.) 1 = {we Qn : w = (A.(x(i».A.(x(i + 1)).

... •A.(x(i+n-1))): x(j+ 1) ef(x(j))Vi:s;j:s; i+n-2}.

The language generated by the NDMO is given by n LrdenC = U UW~. i=lx€X Ifthe output alphabet Q consists of I/O vectors. it is possible to reproduce the observed language by performing state trajectories in the automaton.

Identification and model properties

For identification of DES various methods exist (see e.g. [START_REF] Fanti | Fault diagnosis and identification of discrete event systems using Petri nets[END_REF] for an overview of identification of Petri nets). An identification algorithm that has especially been designed to deliver appropriate models for fault detection purposes is given in [START_REF] Klein | Identijication of discrete event systems for fault detection purposes[END_REF]. It works on the basis of words of the parametric length k and yields a non-deterministic autonomous automaton with output (NDMO) as introduced in Definition 4. Basically. states representing observed words of length k are connected in the order the words have been observed. The algorithm delivers an automaton that is k+ 1-complete [START_REF] Moor | Supervisory control of hybrid systems via I-complete approximations[END_REF]. This means that the automaton generates exactly the set of observed words of length k+ 1: L~tsl =Lfd~~t. This characteristic is advantageous for fault detection: If the observed language has been collected during fault-free system behavior only, it can be made sure that each fault leading to an abnormal word of length k+ 1 can be detected since the model (that is intended to be fault free) cannot erroneously reproduce it. Faulty I/O vector sequences (words) that are reproduced by the model represent undetectable faults.

A second important model property is shown in Roth, Jean Jacques, and Litz (201Oa): If it is possible to state L~t/ =L~-:;i (with L~-:;i denoting the fault-free language of the considered system) for an arbitrary k, L?d~n~ ~ L~~k holds: If the system language of length k+ 1 (L~-:;;) has been completely observed. the identified model simulates the complete original system language of length larger than k (L~~k). This reduces the number of false alerts when using the identified model for diagnosis purposes: If L~tsl =L~-:;i holds for a given k, the model is able to reproduce any fault-free word which can be exhibited by the system. Hence, none of the fault-free words will lead to a false alert. Since both model properties relate the identified and the observed system language, they help validating the model against the observed behavior.

Since for the second property it is necessary to state ~t, 1 =~.

the most important question for the identification approach is to determine a k for which this assumption holds. Fig. 4 gives a typical evolution of the observed word set for a closed loop system over time. The data has been captured for the case study system that will be described in detail in Section 5. The word sets have been collected during fault-free system evolutions. It can be seen that 1Lb b , 1 con verges to a stable level after about 8 system evolutions. This allows concluding that probably each word of length one (single I/Ovectors) that is exhibited during normal system behavior has been observed.

It is reasonable to state that this also holds IL~,I and IL6bsi (I/O vector sequences of length two and three). If word sets containing longer sequences are considered, it can be seen that the according curves do not perfectly converge. For the data depicted in Fig. 4 it can be concluded that ~tsl =~ holds for k=2 due to the convergence of the according curves. A model identified with k=2 will thus not lead to false alerts during online diagnosis. If k is chosen to a larger value, it is likely that the system exhibits new fault-free words of length k+ 1 after the identification data has been collected since the according languages do not converge to a stable level. Since these new words cannot be reproduced by the model, false alerts would be the consequence.

Note that problems similar to the choice of k must be solved if the model is obtained from manual model building. The engineer creating the model has always to decide upon the level of abstraction and necessary model accuracy. It is up to his model ing experience to ensure that the model that is supposed to represent the fault-free system behavior contains any possible fault-free behavior and does not erroneously contain faulty sequences. In general, neither working with identified nor work ing with manually created models guarantees having a perfect system representation.

Fault detection and localization

Definition of I/O behavior

In Section 3.2 it has been explained that the NOMO can be used to reproduce I/O vector sequences of the considered system if its output alphabet consists of I/O vectors. In case of a fault it is not sufficient to know that an observed I/O vector sequence cannot be represented by the model (fault detection). It is necessary to determine I/Os which are possibly responsible of the malfunctioning (fault localization). Hence, it should be possible to determine the behavior of single I/Os based on the observed and modeled I/O vectors. If the system produces a new I/O vector, it is the result of at least one I/O changing its value. If a binary I/O changes its value a rising or falling edge can be observed (see Fig. In order to determine edges appearing if two arbitrary I/O vectors are considered, an edge function is defined: Definition 7 (Edge function). Let10m, 1O,{k) be the i-th controller I/O in the j-and k-th I/O vector.

lOi_1 if IOj (j) = 0 and IOj(k) = 1, Edge(lOj(J),lO j(k» = 10j_0 if IOj(j) = 1 and IOj(k) = 0, { IO;_B if IOj(j) =IOj(k)
delivers the resulting edge of the i-th I/O from the comparison of its values from the I/O vectors uU) and u(k).

Instead of Edge(IO,U), IOj(k» it is also written Edge(uU)[i], u(k)[i])

. Since more than one I/O can change their value when new I/O vectors are produced, the evolution set that summarizes the edges 'between' two I/O vectors is defined: Definition 8 (Evolution set). ES(u(j)

.u(k)) = {Edge(u(j)[I1.u(k)[l1) E EI Edge(u(j)[i],u(k)[i]) #

The NDAAO as part of the evaluator

As explained in Section 2, the first step in diagnosis is fault detection. In the presented approach fault detection is performed using the NDMO as fault-free system model. Fig. 3 shows that the model is used in the evaluator to reproduce the observed system output.

A more detailed view of this procedure is given in Fig. 7. The part with dashed lines will be introduced in the next section. The current I/O vector u(t) that is exhibited by the closed-loop DES is given to the evaluator. Following Algorithm 1, the evaluator tries to determine the current state of the NDMO. The main idea of this algorithm is to determine a current state estimation follow ing the observed evolution set. ).(x) determines the output of an NDMO state according to Definition 4 and ES is the evolution set according to Definition 8. With X t -1 and X t the state estimations after the occurrence of the t-th and t -l-th I/O vector are denoted. Before the algorithm is started, both sets are empty.

Algorithm t. Evaluator algorithm.

Require New observed I/O vector u(t) and former state

estimation X C -1 1: if IX t -1 1> 0: then 2: Xc <-Ix E XI3Xpre E X t-1 : X Ef(x pre ) 1\ ES(),(xpre),).(x» = ES().(xpre),u(t»} 3: else 4: X t <-(x E XI).(x) = u(t)} 5: end if 6: X t-1 <-Xt 7: return X t
The algorithm checks if the former state estimation X C -1 contains at least one state. If the observation is being initialized or after a fault has been detected, this set is empty. In case of an empty estimation IX C -1 1= 0, the evaluator determines each NDMO state with the observed I/O vector as output as possible current state (see line 4 of Algorithm 1) and adds it to the current state estimation Xt. If the former state estimation X t-1 was not empty (line 2), the algorithm checks the observed evolution set of I/O edges in order to determine the current state estimation. Each NDMO state that can be reached by reproducing the observed I/O edges starting in one of the states from the former state estimate Evaluator

[~J=U(l)

X C -1 is added to the current state estimation Xc. At the end of the algorithm, the former state estimation X t -1 for the next run ofthe algorithm is prepared by copying Xc to it.

The complexity of the evaluator algorithm with respect to the number of I/Os in the system is linear since the number of [/Os defines the number of comparisons carried out in line 2 or line 4. Practical experience with the algorithm showed that it is able to process a newly observed I/O vector even with a large model within a few milliseconds which is sufficient for most industrial applications.

The resulting state estimation X t of the evaluator as well as the observed I/O vector u(t) are given to the analyzer where the fault detection policy as well as the fault localization operations are implemented. The fault detection policy is given by

-_ {faUlt if IXtl #-1, (1) FD(Xtl-OK if IXtl = 1.
A fault is detected if the evaluator cannot determine an unam biguous state estimation. Fig. 8 shows an example of the state estimation and fault detection process. In the lower part of the figure, the NDMO is depicted. In the part "observation" an observed I/O vector sequence is shown. Between NDMO states and between I/O vectors, the resulting edges are given. It is supposed that the algorithm starts without a defined initial NDMO state Xo like it is the case if the start of the diagnosis procedure is not synchronized with the start of the closed-loop DES. After the observation of the first I/O vector, the evaluator algorithm goes to line 4 since IX c -1 1= 0. The evaluator determines the NDMO states with

).(x) = u(t) and adds them to Xc. State Xl and state X3 have the first observed I/O vector as output. The fault detection policy returns "fault" since the state estimation is not unambiguous. This changes as soon as the next I/O vector is observed which leads to 10 1 _ °. I/O vectors by determining the state trajectory X3, X4. When the last I/O vector is observed, the resulting evolution set {lOt_a} cannot be produced by leaving state X4. The only evolution set that can be generated by leaving X4 is {lOt_O, 102_0} which differs from the observed one. Hence the state estimation from line 2 in Algorithm 1 results in an empty set which leads to fault detection. During the initial phase of the observation a fault is detected until the state estimation becomes unambiguous. This can be avoided if an initial NOMO state Xo is given and the start of the diagnosis process and the closed-loop OES are synchronized.

Definition of the residuals

After a fault has been detected, the next step is to determine which sensor, actuator or hardware part of the plant is possibly affected. Since sensors and actuators are directly connected to controller l/Os, it is helpful to give a certain number of I/Os that could be related to the fault. In the following two subsections four residuals will be introduced that formalize generic fault symp toms. The idea of the residuals follows the definition of [START_REF] Isermann | Trends in the application of model based fault detection and diagnosis of technical processes[END_REF] where they are defined as fault indicators based on a deviation between measurements and model-based computa tion (see Section 2). In order to calculate the residuals the state estimation before a fault was detected is supposed to be unam biguous IXt-tl = 1. Hence, only one NOMO state was formerly considered as possible current state. If this condition does not hold, the residuals are not calculated.

Unexpected behavior

The first class of residuals has the aim to localize faults that led to an observed behavior that was unexpected in the given context. The current context is defined by the last estimated actual state x in the automaton. The first residual is

Res1(x,u(t)) = ES(A,(x),u(t))\ U ES(A,(x),A,(x')).

(2)

'Ix' ef(x)
With ES(A,(x),u(t)) the rising and falling edges are determined that are observed when comparing the I/O vector of the last estimated current state and the I/O vector that led to fault detection. This set represents what actually happened when the fault was detected and refers to the notion measurements in the residual definition. A,(x')) represents the union of the sets of rising and falling edges when the last estimated current state and each of its direct successor states are considered. It represents the expected behavior. This refers to the model-based computation in the residual definition. The set difference of the observed (ES(A,(x), u(t))) and the expected (Uvx'ef(xlES(A,(x),A,(x'))) behavior is built in the residual equation and thus represents the unexpected beha vior. In Res1, the expected behavior is given by the union of each possible following behavior of the last estimated current state. A stricter formulation of the expected behavior is used in the second residual:

Uvx' ef(x)ES(A,(X),

Res2(X,u(t)) = ES(A,(x),u(t))\ n ES(A,(x),A,(x')).

(3)

Vx' ef(x)

Instead of a union over the expected behavior of the possible following states, an intersection is used. The intersection delivers the I/O edges that must be observed no matter which following state in the model is taken. It is obvious that Res1 c;;Res2 since Uvx'eflX)ES(A,(X),A,(x')) ;2 flvx' ef(x)ES(A,(x),A,(x')).

Example. Fig. 9 shows an example for unexpected behavior that led to fault detection (instead ofthe I/O vectors only the I/O edges between two states are given): from the estimated actual NOMO state Xt it is not possible to take a transition that has exactly the observed falling edges 10 3 _ °and 10 4 _ °.Hence a fault is detected. The result of Res1 (with x= Xt ) is

Res1(x,u(t)) = {I03_0,104_0}\({I0t_O,102_1,104_0} U {I0t_O,102_1}) = {I03_0}.
This result means that 10 3 _°was an unexpected event which implies that the system operator should check the sensor or actuator that is connected with 10 3 • However, it is possible that the fault cannot be found at this component. If this is the case, Res2 should be calculated in order to use a stricter formulation of the expected behavior. In Res1 each possible following behavior is subtracted from the observation. Using Res2 only the behavior that must occur no matter which regular following behavior is considered:

Res2(x,u(t)) = {I03_0,104_0}\({I0t_O,102_1,104_0} n (lO t _O,10 2 _1}) = {I°3_0,10 4_O}
This result implies that the occurrence of a change in value of 10 4 is not always expected and thus another possible fault localization.

Missed behavior

In contrast to an observed but unexpected behavior it is also possible that a fault can be localized by determining a missed event. Set operations that help to localize an expected but unobserved behavior are given by the third and the fourth residual:

Res3(x,u(t)) = n ES(A,(x),A,(x'))\ES(A,(x),U(t)).

(4)

v« ef(x)
Res3 is the set difference of the I/O edges that are expected no matter which following state is taken (flvx' ef(x)ES(A,(x),A,(x'))) and the I/O edges that have been observed (ES(A,(x),u(t))). Each rising or falling edge that must occur when the estimated actual state is left but has not been observed is part of Res3. The expected behavior is represented by the intersection of each possible following behavior. It is also possible to give a less strict formulation of the expected behavior by using the union operation instead of the intersection:

Res4(x,u(t)) = U ES(A,(x),A,(x'))\ES(A,(x),u(t)).

(5) Vx' ef(x) Since Res3 c;;Res4, the result of Res4 is usually less restrictive than Res3, it contains more elements. Especially if production systems with many I/Os are consid ered, the residuals can help to get a relatively small set of I/Os that could be related to the fault. A maintenance operator can ,-'"'-/0, _ 0-7 observed \!" --.. results from ES(A(X,).u(I» X\ /0,_0/0,_' Fig. 10. Example for a missed behavior.

Example. The situation in

then check the possibly faulty sensors or the related actuators. Actuators are related to a sensor if their activation or deactivation can have an influence on the sensor value.

Reduction of the residual sets

Once a fault has been detected. the four residuals are available to localize controller I/Os that are related with the fault. It has been explained that the residuals analyze two different generic fault symptoms: missed and unexpected behavior. In general. it is not a priori possible to tell if a fault was caused by missed or by unexpected behavior. Hence. all residuals must be calculated which results in a set of possible fault candidates. The system operator has to check each of the I/Os according to the priority explained before. This operation can be supported by the infor mation if the considered I/O is suspected to have shown a missed or an unexpected behavior. In order to give a more precise estimation of which I/O is possibly affected, it is possible to take further I/O vectors into account that follow the vector that led to fault detection. Of course. this approach is only possible if the system does not immediately have to be stopped after fault detection like for fault tolerant systems. The reduction of the fault candidate set is based on the following heuristic:

After a fault has been detected it is possible to delete a fault candidate from the candidate list if it shows a following behavior that can reasonably be considered as fault-free.

To decide if the following behavior of a possibly affected I/O can reasonably be considered as fault-free. a special state estimation algorithm is proposed. The idea of this algorithm is to mask out the possibly faulty part of the I/O vector and to perform state estimation only based on the 'healthy' I/Os. This state estimation results in a set of states that represent possible descriptions for the non-faulty part of the current system state. If a possibly affected I/O has the same value in each of these states and in the currently observed I/O vector, it is assumed that this can be considered as the result of a fault-free I/O behavior. The I/O showed a behavior that is completely explicable with fault free system states which have been determined without considering the I/O itself. Hence this I/O can be deleted from the fault candidate list.

The candidate set of all I/O edges that are reported by the residuals is given by

Resu(x.u(t)) = ES(l(x).u(t)) u U ES(l(x).l(x'»).

(6) VI( eN') It is obvious that Res., = Res1 URes2URes3 URes4. In order to get the controller I/Os that have edges in the candidate set the following function is applied:

JOUst = (JOdOOi_1 v JOj_O) E Resu}. (7)
This JOUstcontains the part of the system I/O vector that is possibly influenced by the fault based on the result of the residuals. The possibly faulty I/Os from JOList can create the following edges:

CandEdges = {JOj_1,JOi_0 E EIJOj E 10List}. (8)
Since the state estimation algorithm will work on the basis of the non-affected part of the I/O vector. this part has to be calculated by an I/O vector projection: Definition 9 (I/O vector projection). The I/O vector projection of an I/O vector u to a list of I/Os (JOUst) that is to be masked out is defined as uri] if JOi¢JOUst, JOPIOLisr(u)[zl = { * if JO; E JOList vi = 1, ....m (with m denoting the number of controller I/Os in vector u).

~OPIOLjSr(U)[l] 10PIOList(U) = : { JOPIOList(u)[m]
applies the function to the entire vector. JOPIOLisrtU) contains for each I/O that has to be masked out (i.e. an I/O in the JOList) the do not care symbol *. For each I/O that has to be considered, JOPIOList(U) has the according I/O value in vector u.

Algorithm 2. State estimation and candidate set reduction.

Require JOUst. CandEdges. X {Initialization} 1:

X t _ 1 <-X {State estimation} 2:
Wait for the next observed I/O vector denoted as u(t) 3:

Xt <-(x EXt_I IJOPIOList(l(x)) =JOPIOList(U(t))} 4: Xt<-XtU{x'EX 13x EX t_1 :x'Ef(x)AlOP,oList(l(x')) = JOPIOList(U(t))} 5: XReach <-(x(i+n) E XI3(x(i)..... x(i+n») :
x(i) E X t /\ vi ~j ~ i +n : (x(j+ 1) Ef(x(j))/\ES(l(x(J)).

l(x(j+ 1))) end if

With these sets it is possible to perform Algorithm 2 instead of Algorithm 1 in order to estimate the current system state. Since the algorithm starts with considering the first I/O vector after fault detection, the former state estimate is empty and must be reinitialized. Each NDAAO state must be considered as possible estimate (line 1). When the next I/O vector is observed it is checked which NDAAO state of the former estimate is still con form with the 'healthy' part of the new I/O vector (line 3). Possibly affected [/Os do not contribute any information to determine the new estimate X t • After this, NOMO states are determined which can be reached starting in the former estimate X t -1 by considering the I/O vector projection of the new I/O vector. Each following state of the former estimate (x' E!(x)lx EX t -1 ) that has the same 'healthy' part of the [/0 vector as output as it has been observed is added to the new estimate (line 4). Since the state estimation should only be restricted with the non-affected part of the I/O vector, the set XReach is calculated in line 5. It contains each NOMO state that can be reached starting in the current estimation X t if during the state trajetory only possibly affected I/O edges occur (CandEdges). This set is added to the current estimate. This ensures that the state estimation is not blocked by I/Os that are possibly related to the fault. In line 7, the set Xt-1 for the analysis of the next I/O vector is prepared.

After the state estimation part is finished, the candidate reduction starts. For each I/O in the /oUst it is checked if it has the same value in each state of the estimate and in the currently observed I/O vector. If this is true, the I/O gets removed from the JOUst since the currently observed I/O value is explicable with each estimated NOMO state. Since the estimated states were determined only using the 'healthy' part of the I/O vector, it is assumed that it is highly probable that the considered I/O is not faulty as it fits well in each fault-free interpretation of the current system state. After the JOUst has been updated, the CandEdges set must also be recalculated. Then, the calculated sets are inter preted. If the /oUst is empty, the fault detection was probably a false alert since the current observation can perfectly be explained by the estimated fault-free model states. Most often, in this case IXt! = 1 also holds. Hence, the normal observation can go on (Algorithm 1). If at least one I/O is left in the /oUst and X t is not empty, the algorithm loops back and waits for the next I/O vector. If it was not possible to determine a state estimation, the algorithm cannot proceed and stops. The operator must check the remaining I/Os. Note that the number of candidates in JOUst is typically much smaller than the number of I/Os in the system: Before the algorithm is started, JOUst is calculated according to Eq. ( 7) as the union of l/Os contained in the residuals. Processing further I/O vectors, JOUst can only be reduced in the algorithm as can be seen from line 8 to 12.

The algorithm only calculates a new state estimation and an updated candidate set if a new I/O vector is observed. The time to recalculate these two sets depends thus on the time between the observation of the new I/O vectors and cannot be calculated a priori. Since the algorithm has to compare each value of an [/0 in an observed vector with the according I/O in the model states, the complexity of the algorithm with respect to the number of I/Os is linear.

case study

System description

In order to demonstrate the relevance of the proposed method for existing closed-loop DES, a case study is presented in this section. The considered system is a laboratory facility at the Institute of Automatic Control, University of Kaiserslautem. The purpose of the system depicted in Fig. 11 is to treat work pieces that are stored in the feeder (left most station) with three machine tools. The system is controlled using a Siemens S300 PLC (Programmable Logic Controller) with 14 inputs and 15 outputs. The controller inputs and the corresponding sensors can be seen in Fig. 11. Inputs that are written in italic are connected to sensors with a specific technology that delivers a logical 0 if they detect something and a logical 1 if they do not detecting anything. All sensors except of those for work piece positions have this technology. The outputs are given in Table 1. If they are set to 1, the according actuator is activated. In order to make the I/O names easier to read they are not labeled 10; Horizontal milling machine milling motor 03.5 Horizontal milling machine conveyor on but with 11.2 etc. to indicate the second input that belongs to the first machine and 02.4 for the fourth output of the second station.

A production cycle (or system evolution in terms of Section 3) consists of the treatment of two work pieces that are stored in the feeder. At the beginning of a system evolution the feeder pushes the first work piece to the conveyor. Then it is transported to the drilling station (station 1). In the drilling station the work piece is stopped and two holes are drilled. After this, the work piece is transported to the vertical milling machine (station 2) and the next piece is taken from the feeder and transported to station 1. The first work piece gets treated by the vertical milling machine (each of its three milling tools is applied to the piece) and the second work piece by the drilling station. After the treatment at the vertical milling machine has finished the according work piece is transported to the horizontal milling station (station 3). The second work piece moves from the drilling station to the vertical milling machine. When the work piece has been treated in the horizontal milling machine, it is stored in the storage at the right side in Fig. 11. This process continues until both work pieces have been treated by each of the three stations.

Data collection for identification and online diagnosis

As explained in Section 2, the proposed method works on the basis of signals exchanged between controller and plant (see also Fig. 3). The data collection method that was chosen for the case study works with a standard communication processor that is part of the PLC. When controlling the system evolution, the PLC performs the following three steps cyclically: reading the inputs, executing the program in order to determine the new output setting and finally writing the outputs. The transfer of the I/O vector to the diagnosis computer takes place at the end of the second step (program execution). At the same time when the newly determined outputs are transfered to the output card of the PLC, the new I/O vector is sent to a standard PC using the communication processor via a UDP-connection. Since the I/O vectors are captured at the end of the PLC cycle and are thus sampled, it is possible to have different I/Os changing their value between two I/O vectors. This method can be applied to each system with a controller that is able to send data via an Ethernet data link which is a common feature of almost all industrial controllers today.

For the application of the proposed diagnosis method it is necessary to have a model that describes the fault-free behavior of the considered closed-loop system. For the case study system such a model has been identified using the method outlined in Section 3. The data base for the identification consists of 60 fault-free system evolutions (60 treatments of two work pieces) that have been observed using the data link described above. The evolution of the according observed languages (sets of I/O vector sequences of varying length) are shown in Fig. 4 in Section 3.3. It has been explained that with k=2 a model can be identified which is not expected to deliver false alerts during online diagnosis. The identi fied automaton has 121 states and 162 transitions. Fault-free and faulty system behavior can be distinguished since behavior that has not been seen during the model learning phase is considered as a fault symptom. Parts ofthe identified NDMO can be seen in Figs. 12 and 13. Since the state output ),(x) is too large to be depicted (I/O vector with 29 I/Os), only the evolution set that is obtained when comparing the I/O vectors of two connected states is shown.

Fault localization examples

The first example is a fault at the sensor connected with controller input 13.3. During normal operation, this controller input changes its value from 1 to °when the head of the horizontal milling station reaches the bottom position. The fault that has been introduced artificially is a short circuit causing the sensor to switch its signal from 1 to 0. The fault has been produced when the first work piece was treated by the second machine (vertical milling). At the same time, the second work piece gets drilled in the first station. In this situation no work piece is in front of the horizontal milling station and its milling head is in its home position at the top. Hence, none of the inputs belonging to station 3 should change their value. The short circuit has been introduced when milling with the second of the three tools in station 2 started and the vertical milling head left the top position. This situation is depicted in Fig. 12. After the observation of 12.2_1 indicating that the milling head just left its home position, the evaluator algorithm delivers state XIS as single state estimate. Hence, no fault is detected. Then the falling edge 13.3_0 is observed that occurs due to the fault. The evaluator algorithm does not find a following state of XIS that can be reached by this edge. Since X t is empty, a fault is detected due to Eq. ( 1).

Consequently the residuals are applied with x= XIS as the former unique estimate:

Res1(x,u(t» = {13.3_0J, Res2(x,u(t» = {13.3_0}, Res3(x,u(t» = {},
Res4(x,u(t» = {J1.5_1,J2.3_0,Ol.3_1,Ol.4_1,Ol.5_1,02.2_1,02.3_0}.

Res1 and Res2 both result in the same unexpected edge that is caused by the fault. Res3 results in an empty set and Res4 shows a union of the legal following behavior of state XIS. In this case, Res3 and Res4 do not contribute to fault localization. Res1 and Res2 provide the "right" candidate for fault localization.

The second example is a fault at the sensor connected to 12.4. This sensor detects if a work piece is at the left most position of station 2 (12.4= 1). The artificially introduced fault prevents the sensor from switching back to °when the work piece has left the left most position. The fault was introduced when the first work piece gets transported from station 2 to station 3 and the second work piece from station 1 to station 2. Fig. 13 shows two parts of the identified NDMO. The lower part represents the described situation. The evaluation algorithm (Algorithm 1) determines X49 as the actual estimate when the second work piece arrives at position 12.4 and produces the rising edge 12.4_1. Now, the fault prevents 12.4 from switching again to °when the work piece leaves the position due to the start of the conveyor (02.5_1). The fault is detected when the rising edge of 13.5 and its correspond ing output changes are observed before 12.4 changed back to 0. The residuals are calculated: In this case, Res1 and Res2 do not contribute to fault localization since they show the observed behavior that is not responsible for the fault. The controller input related to the blocked sensor is localized by both Res3 and Res4 since it is a missed behavior.

A fault that leads to exactly the same symptoms is a defect at the motor of the conveyor in front of the vertical milling machine (02.5).lfthe conveyor motor does not start due to the motor fault, the symptoms are comparable: The sensor connected to /2.4 does not change its value since the conveyor does not transport the work piece away from the left most position of station 2. Hence, if Res3 or Res4 report a missing change in value at a sensor, it is reasonable to not only check the according sensor but also the actuators that have a direct influence to this sensor.

In the second example, Res3 and Res4 contained the controller input that was connected to the faulty sensor or the sensor that did not change its value due to a fault at a related actuator. Additionally, the residuals contain controller outputs. These controller outputs are set if a controller input leads to a fulfilled logical condition in the control program and are usually not directly related to the fault. Nevertheless, they can be useful to decide if the fault candidates reported by Res1 and Res2 or the candidates of Res3 and Res4 should be analyzed at first. If controller inputs appear in Res1 or Res2, they are supposed to be unexpected. Indeed, the fault-free system model did not expect their change in value at the time of its observation. Nevertheless, if a change in value of an input together with change in values of controller outputs is observed, it is possible that this was not completely unexpected: The controller program was in an internal state where it waited for the change in value of the input and thus triggered some outputs upon its observation.

Given a situation like in the second example, this can help to decide that the observed change in value of 13.5 was not completely unexpected by the controller since it was in an internal state that allowed setting corresponding outputs. Hence, the candidates deliv ered by Res3 and Res4 should be checked first which delivers the right fault candidate.

Apart from using this strategy to decide which candidate should be analyzed first, the candidate set reduction algorithm presented in Section 4.4 can be applied. As explained before, it must be possible to run the system even if a fault was detected. This is a scenario that is often possible in industrial system like e.g. the system in [START_REF] Supavatanakul | Diagnosis of timed automata: Theory and application to the damadics actuator benchmark problem[END_REF]. As an example to illustrate Algorithm 2 we consider the fault of the second example (12.4 does not switch back to 0). Based on the residuals calculated before, the combined set Resu can be determined to Resu (x, u(t)) = {/3.5_1,/2.4_0,03.3_1,03.4_1,03.5_0,O1.5_0} according to Eq. ( 6).

Algorithm 2 is started to reduce the candidate set and to estimate the current system state using the next observed I/O vector following fault detection. In line 1, the complete NDMO state space is assigned to the former estimate. The I/O vector that is observed after fault detection leads to the evolution set {12.5_1, 02.3_1,02.4_1, 02.5_0} when compared to the former I/O vector. The I/O vector projection of the current vector and the state outputs of the former estimate (Iine3) leads to the following estimated states: X17 (upper part of Fig. B), XS1, X93, X113 and X118. This estimate remains unchanged in line 4. XReach representing the states that can be reached starting in one of the states of the current estimate and only producing edges of possibly affected I/Os is determined to state X113. Starting from XSI it is possible to reach X113 only producing possibly affected I/Os. Since X113 is already part of the current estimate it does not change in line 6.

Comparing the observed values of the possibly affected l/Os with the values in each ofthe estimated states reveals that none ofthe I/Os can yet be erased from JOUst. Since IX tI > 1 and IJOUstl > 0 the algorithm loops back to line 2.

The next observed I/O vector leads to the edge 13.2_1 when compared with its predecessor. Since 13.2 is not part of the possibly affected l/Os, it does not have the observed value 1 in one of the states from the former estimate. Hence, Xt becomes an empty set in line 3. As state XlOS has the same I/O vector projection as the currently observed I/O vector and XlOS E!(X93) (X93 was part of the former estimate), the current estimate is determined to {XlOS} in line 4. The operation of line 5 results in an empty set. Hence, the current estimate consists only of one state. The candidate set reduction from line 8 to 12 reveals that only the I/Os 12.4 and 01.5 differ in the observed vector and in the estimated state. The other fault candi dates can thus be removed from /oUst since they showed a behavior that can be considered as normal. Since IX tI > 1 and IJOList I> 0 the algorithm loops back to line 2 but has reduced the number of fault candidates. Hence the fault must be searched at the sensor that is connected to 12.4 which corresponds to the artificially introduced fault.

Conclusions and outlook

In this paper a fault localization technique working on the basis of a fault-free nominal system model has been presented. The concept of residuals known from fault diagnosis in continuous systems has has been introduced to solve the fault localization problem in DES. It has been shown how an identified model can be used for fault detection and fault localization. Results from a case study that represents typical characteristics of industrial production systems show that the residual approach is able to deliver a small subset of controller I/Os possibly related to a detected fault. Using a special candidate set reduction algorithm it is possible to reduce this set of fault candidates. The presented approach was developed for monolithic system models. In large applications, automata networks are often used for diagnosis purposes following a decentralized approach. After the development of an appropriate distributed identification approach in Roth, Jean-Jacques, and Litz (2010b), the fault localization approach will be extended for the application in automata networks in future works.

Many faults are characterized by an altered timed behavior. Typical examples are actuator faults often leading to deadlocks. Analyzing the timed behavior for diagnosis purposes necessitates using timed models. Since the determination of timed models is very demanding, current work concentrates on identifying the timed behavior. After the development of appropriate identifica tion algorithms, future works aim at extending the residual approach for timed models.

  Fig. 1. Closed-loop DESconsisting of plant and controller.

Fig. 2 .

 2 Fig. 3. FDI principle for closed loop DES.

  These signals are available by considering the inputs and outputs of the controller. The con troller I/Os are arranged in I/O vectors with several properties summarized in Definition 1. Definition t (Controller I/O vector). Given r different controller inputs II• ... • IT and s different controller outputs 0 1 • __ .• Os. the controller I/O vector U=(J0I •... ,JOm) with m=r+s is given by IOj = Ijvi = 1.. _.r and IO T + i = 0ivi = 1.. _.s. m denotes the dimension of the vector (number of controller I/Os).

:

  Fig. 4. Evolution of the observed word set for the case study.

  5): Definition 6 (Edges). For each controller I/O IO j there exist three edges: 10;_0 to indicate a change from 1 to 0 (falling edge), 1Oi-1 to indicate a change from 0 to 1 (rising edge) and IO j _ B to indicate no change in value: E = {JOj_O U lOi_1 U IO jj'i1 ::0:; i::o:; m}, with m denoting the number of controller I/Os in the system.

  Fig. 5. Example for rising or falling edges of an 1/0.

Fig. 6

 6 Fig.6shows an example for the evolution set resulting from the comparison of two I/O vectors.

Fig. 6 .

 6 Fig. 6. Example fort he evolution set resulting from the two I/O vectors u(1) and u(2).

  Starting from Xl or X3 only X2 can be reached by producing the same evolution set. Hence the state estimation only contains one state (line 2 of Algorithm 1) which makes the fault detection policy declare "OK". The state estimation proceeds with the next observed (~J 10,_0 (~J 10,_1 (~J 10,_1 (:J IOI_O(~J s ]-t--'-..:;...",-.-"--'----'-'---'-'----'-'--.Y, {x,.x,} (x,)

Fig. 8 .

 8 Fig. 8. Example for state estimation and fault detection.

  Fig. 9. Example for an unexpected behavior.

  Fig. lOis considered to iIlustrate Res3 and Res4. The observation of the edge 10 3 _ °leads to fault detection with Xt as estimated actual NOMO state. Applying the residuals results in Res3={/O t_O} and Res4={lOt_OJ02_1}. The same proce dure as explained for Res1 and Res2 should be started: First, the component connected with lOt should be checked. If the system operator does not find a fault at this component, he should proceed with checking 102 which is additionally part of Res4.

  Fig. 11. Lab system.

Fig. 12 .

 12 Fig. 12. Example for a fault at 13.3.

  Res1(x,u(t» = {13.5_1 ,03.3_1 ,03.4_1 ,03.5_0}, Res2(x,u(t» = {13.5_1 ,03.3_1 ,03.4_1 ,03.5_0}.

Fig. 13 .

 13 Fig. 13. Example for a fault at 12.4.


Alexis Benichoux

Prasad Sudhakar

Frédéric Bimbot

Rémi Gribonval

Well-posedness of the permutation problem in sparse filter estimation with ℓ p minimization

Keywords: sparse filter, convolutive blind source separation, permutation ambiguity, ℓ p minimization, Hall's Marriage Theorem, bi-stochastic matrix Filtres parcimonieux, séparation aveugle de sources, mélange convolutif, ambiguïté de permutation, théorème de mariage de Hall, matrice bi-stochastique

Convolutive source separation is often done in two stages: 1) estimation of the mixing filters and 2) estimation of the sources. Traditional approaches suffer from the ambiguities of arbitrary permutations and scaling in each frequency bin of the estimated filters and/or the sources, and they are usually corrected by taking into account some special properties of the filters/sources. This paper focusses on the filter permutation problem in the absence of scaling, investigating the possible use of the temporal sparsity of the filters as a property enabling permutation correction. Theoretical and experimental results highlight the potential as well as the limits of sparsity as an hypothesis to obtain a well-posed permutation problem.

Introduction

Blind source separation and blind source localization are ubiquitous problems in signal processing, with applications ranging from wireless telecommunications to underwater acoustics and sound enhancement.

These problems can be considered as reasonably well understood and solved in simple linear instantaneous settings, where tools such as Independent Component Analysis, as well as techniques exploiting source sparsity, are now mature. However, the convolutive source localization / separation problem remains much more challenging. In particular, without further assumption than statistical independence between sources, the problem is known to be ill-posed because of the so-called frequency permutation (and scaling) problem: at best, one can hope to estimate for each frequency (up to a source and frequency dependent scaling factor) the collection of frequency components of all sources (and of the associated mixing filters); but one cannot match the estimated frequency components from different subbands to globally identify the sources (and mixing filters).

Several practical approaches have been proposed to solve the permutation and scaling problems in practice, by exploiting various properties of either the mixing filters or the sources to match different frequency subbands. While some of these methods may succeed in practice for certain types of sources / filters, there is no known theory guaranteeing the well-posedness of the permutation and scaling problem under appropriate assumptions.

This paper contributes to fill this gap, by providing well-posedness guarantees for the permutation problem under sparsity assumptions on the mixing filters. Sparse filters, associated to impulse responses corresponding to a limited set of echoes, are typically encountered in a number of underwater communication channels [START_REF] Berger | Sparse channel estimation for multicarrier underwater acoustic communication: From subspace methods to compressed sensing[END_REF] or wireless telecommunications scenarios [START_REF] Sharp | Application of sparse signal recovery to pilotassisted channel estimation[END_REF][START_REF] Bajwa | Compressed channel sensing: A new approach to estimating sparse multipath channels[END_REF] which are relevant for blind source localization and separation, and the theoretical results achieved in this paper indicate that this property can potentially be exploited for blind estimation in this context.

Problem formulation and notations

Let x i [t], 1 ≤ i ≤ M be M mixtures of N source signals s j [t], resulting from the convolution with filters a ij [t] of length L such that:

x i [t] = N j=1 (a ij ⋆ s j )[t], 1 ≤ i ≤ M, (1)
where ⋆ denotes convolution. The filter a ij [t] typically models the impulse response between the j th source and the i th sensor. By abuse of notation, Fa ij = {a ij [ω]} 0≤ω<L denotes the discrete Fourier transform of the filter seen as a vector a ij = {a ij [t]} 0≤t<L ∈ C L . Also, the mixing equation (1) can be rewritten as X = A ⋆ S, with A the matrix of filters

A := ({a ij [t]} 0≤t<L ) 1≤i≤M, 1≤j≤N ,
X the observation matrix and S the source matrix. In this context, blind filter estimation refers to the problem of obtaining estimates of the filters A from the mixtures X, without any explicit knowledge about the sources S. Mixing filters estimation is relevant for several purposes such as deconvolution, source localization, etc. [START_REF] Benesty | Springer Handbook of Speech Processing[END_REF]. It also has a relationship with the problem of Multiple-Input-Multiple-Output (MIMO) system identification in communications engineering [START_REF] Chen | Frequency domain blind mimo system identification based on second and higher order statistics[END_REF].

Frequency domain filter estimation

Estimating the mixing parameters is made easier when all filters are instantaneous, that is to say of length L = 1, as the convolution product in ( 1) is replaced by the usual product. However, things get complicated in the general setting of convolutive mixtures.

A common approach for filter estimation then relies on the transformation of the mixing model in Eq. ( 1) into the time-frequency domain, converting a single convolutive filter estimation problem into several complex instantaneous filter estimation problems. Using standard techniques for instantaneous mixing parameter estimation [START_REF]Handbook of Blind Source Separation, Independent Component Analysis and Applications[END_REF], complex mixing filter coefficients

Ã[ω] = {ã ij [ω]} 1≤i≤M, 1≤j≤N
are estimated for each frequency bin 0 ≤ ω < L.

Permutation and scaling ambiguities

Without further assumption on either the filters a ij [t] or the sources s j [t], one can at best hope to find an estimation à = (ã ij ) where for every frequency ω we have ãij

[ω] = λ j [ω]a iσω (j) [ω], (2)
with λ j [ω] a scaling ambiguity and σ ω ∈ S N a permutation ambiguity, where S N is the set of permutations of the integers between one and N . Several methods [START_REF] Pedersen | A survey of convolutive blind source separation methods[END_REF] attempt to solve for these ambiguities by exploiting properties of either the sources S or the filters A [8, 9, 10].

Exploiting sparsity to solve the permutation ambiguity

The focus of this article is the use of the sparsity of A in the time domain to find σ 0 . . . σ L-1 ∈ S N , assuming the scaling is solved, i.e., λ j [ω] = 1. Assuming that A is sparse means that each filter a ij has few nonzero coefficients, as measured by the ℓ 0 pseudo-norm

a ij 0 := ♯{0 ≤ t < L, a ij [t] = 0} = t |a ij [t]| 0 .
The approach considered in this article is to to seek permutations σ 0 , . . . σ L-1 yielding the sparsest estimated time-domain matrix of filters A = ( a ij ) where a ij [ω] := ãi σω (j) [ω]. Besides the ℓ 0 pseudo-norm A 0 := ij a ij 0 , the following ℓ p quasi-norms will be used to quantify the sparsity of A:

A p p := ij a ij p p = ijt | a ij [t]| p , 0 < p ≤ 1.

Inria

Main results

Our main result (Theorem 2) is a theoretical guarantee that when the filter length L is prime, k-sparse filters (i.e., such that a ij 0 ≤ k) uniquely minimize the ℓ 0 norm of A (up to a global permutation) if k L ≤ α(N ), where N is the number of sources. To reach this bound we exploit uncertainty principles as well as the bistochastic structure of the problem through an apparently new quantitative result on bistochastic matrices (Lemma 2).

Structure of the paper

The main theorems are stated in Section 2, and the main ingredients of their proofs are described in Section 3. In Section 4 we discuss the strength of the assumptions used in the theorems, and how much these could be relaxed. In Sec. 5, a naive combinatorial ℓ p minimization algorithm is proposed to resolve filter permutations and used for Monte-Carlo simulations. We conclude with a discussion of the potential, as well as the limits, of sparsity as a hypothesis to solve permutation problems, in connection with the theoretical and empirical results. All proofs are gathered in the appendix.

Theoretical guarantees

Given an M × N filter matrix A, made of filters of length L, and an L-tuple (σ 0 , . . . σ L-1 ) ∈ S N of permutations, we let à be the matrix obtained from A by applying the permutations in the frequency domain, as in (2), without scaling (λ j [ω] = 1).

The effect of the permutations is said to coincide with that of a global permutation π ∈ S N of the columns of A if ãij = a iπ(j) , ∀i, j, or equivalently in the frequency domain:

ãij [ω] := a iσω (j) [ω] = a iπ(j) [ω], 0 ≤ ω < L, ∀i, j.
This is denoted A ≡ Ã. First, we show that for filters with disjoint time-domain supports, permutations cannot decrease the ℓ p norm, 0 ≤ p ≤ 1:

Theorem 1 Let Γ ij ⊂ {0, . . . , L -1} be the time-domain support of a ij . Suppose that for all i and j 1 = j 2 we have

Γ i,j1 ∩ Γ i,j2 = ∅. (3)
Then, for 0 ≤ p ≤ 1, we have à p ≥ A p .

Note that filters with disjoint supports need not be very sparse: M filters of length L can have disjoint supports provided that max j a ij 0 ≤ L/M . Yet, disjointness of filter supports is a strong assumption, and Theorem 1 only indicates that frequency permutations cannot decrease the ℓ p norm. Thus, the minimum value of the ℓ p norm might not be uniquely achieved (up to a global permutation). In our main result, we consider k-sparse filters of prime length, and p = 0:

Theorem 2 Let A be an M × N matrix of filters of prime length L. Assume that max ij a ij 0 ≤ k, (4)
where

k L ≤ α(N ) := 2 N (N +2) if N is even, 2 (N +1) 2 if N is odd. (5)
Then, up to a global permutation, A uniquely minimises the ℓ 0 pseudo-norm among all possible frequency permutations.

3 Main elements of the proof of Theorem 2

The proof of Theorem 2 relies on a measure of the "amount" of incurred permutation, on uncertainty principles, and on combinatorial arguments related to bi-stochastic matrices, involving Hall's Marriage Theorem.

Measures of the amount of incurred permutations

To measure the "amount" of incurred permutation, one can count the number of frequency bands where a non-trivial permutation is incurred, with respect to the best matching reference global permutation π, i.e., min π ♯{ω, σ ω = π}. However, this generally yields the maximum count L -1.

An alternative is to count the "size" of the incurred permutations, given a reference global permutation π, as the maximum number of frequencies where each estimated filter actually differs from the (globally permuted) original filters, yielding:

∆( Ã, A|π) := max i,j F(ã ij -a iπ(j) ) 0 (6) ∆( Ã, A) := min π∈SN ∆( Ã, A|π). (7)
Note that ∆( Ã, A) = 0 iff à ≡ A.

Exploitation of an uncertainty principle

With this notation, we have the following Lemma:

Lemma 1 Assume that à ≡ A, that L is a prime integer, and that (4) holds with 2k + ∆ ≤ L.

Then à 0 > A 0 and ãij 0 ≥ a ij 0 , ∀i, j. The latter inequality is strict when ãij = a ij . For a general L (not necessarily prime), the same conclusions hold when the assumption (8) is replaced with

2k • ∆ < L. (9)
The skilled reader will rightly sense the role of uncertainty principles [START_REF] David | Uncertainty principles and signal recovery[END_REF][START_REF] Elad | A generalized uncertainty principle and sparse representation in pairs of bases[END_REF][START_REF] Tao | An uncertainty principle for cyclic groups of prime order[END_REF], Theorem 1] in the above lemma.

Inria

Combinatorial arguments

Using Lemma 1 with prime L, a simple combinatorial argument can be used to obtain a weakened version of Theorem 2, with the more conservative constant α ′ (N ) := 1/2N !: by the pigeonhole principle, for any L-tuple of frequency permutations among N sources, at least L/N ! permutations are identical; as a result, ∆( Ã, A) is universally bounded from above by L -L/N !; hence if k ≤ L/2N ! we obtain 2k + ∆ ≤ L and we can conclude thanks to Lemma 1.

The proof of Theorem 2 with the constant α(N ) exploits a stronger universal upper bound ∆( Ã, A) ≤ L(1 -2α(N )), obtained through an apparently new quantitative application of Hall's Marriage Theorem [START_REF] Hall | On representatives of subsets[END_REF] to bi-stochastic matrices.

Definition 1 (Bi-stochastic matrix) An N ×N matrix B is called bi-stochastic if all its entries are non-negative, and the sum of the entries over each row as well as the sum of the entries over each column is one.

Lemma 2 Let B be an N × N bi-stochastic matrix: there exists a permutation matrix P such that all the entries of B on the support of P exceed the threshold

2α(N ) = 4 N (N +2) if N is even, 4 (N +1) 2 if N is odd. (10)
Corollary 1 Let σ 0 , . . . , σ L-1 ∈ S N be L permutations. There exists a global permutation π such that

C jπ(j) = ♯{ℓ : σ ℓ (j) = π(j)} ≥ 2Lα(N ), ∀1 ≤ j ≤ N.

Discussion

The reader may have noticed that Theorem 2, while dropping the disjoint support assumption from Theorem 1, introduces new restrictions: the assumption that L is prime, and the restriction to p = 0 compared to 0 ≤ p ≤ 1 in Theorem 1. How important are these restrictions ? Could they be relaxed while exploiting sparsity together with the disjoint support assumption ? This is discussed in this section.

Extending Theorem 2 to non-prime filter length L?

As indicated by Lemma 3 below, for even L ≥ 4, there exists sparse matrices of filters that are the sparsest but not unique (even up to a global permutation) solution of the considered problem: certain frequency permutations provide an equally sparse but not equivalent solution.

Lemma 3 For any integer k such that 2k divides L, there exists a matrix of k-sparse filters A and a set of L/2k frequency permutations resulting in à ≡ A such that for all 0 ≤ p ≤ ∞: à p = A p , and

ãij p = a ij p , ∀i, j. (11)
We have 2k

• ∆( Ã, A) = L.
The fact that the filter matrices A and à satisfy 2k •∆( Ã, A) = L shows the sharpness of Lemma 1 for the case when L is even: the strict inequality in ( 9) cannot be improved.

Specializing Lemma 3 to k = 1 for even L ≥ 4 yields ideally 1-sparse filters a ij and a set of L/2 frequency permutations such that: ãij are 1-sparse; Ã is not equivalent to A and cannot be discriminated from it by any ℓ p norm.

Stronger guarantees with disjoint supports and sparsity ?

Could one get improved results by combining the disjoint support assumptions from Theorem 1 and the sparsity assumption from Theorem 2 ? For even L ≥ 4, Lemma 4 below indicates the existence of sparse matrices of filters with disjoint supports that are the sparsest but not unique (even up to a global permutation) solution of the considered problem: certain frequency permutations of "size" ∆ = L/2k provide an equally good but not equivalent solution.

Lemma 4 For any integers k ′ < k ≤ L/2 such that 2k ′ divides L, there exist a matrix of k-sparse filters A with disjoint supports (3), and a set of L/2k ′ frequency permutations resulting in à ≡ A, such that for all 0 ≤ p ≤ ∞:

à p = A p and ãij p = a ij p , ∀i, j. (12)
We have

2k ′ • ∆( Ã, A) = L.
Specializing Lemma 4 to k ′ = 1 and k = 2 for even L ≥ 4 yields 2-sparse filters a ij and a set of L/2 frequency permutations such that: ãij are 2-sparse; Ã is not equivalent to A and cannot be discriminated from it by any ℓ p norm. This shows that even by adding the disjoint support assumption, for even L ≥ 4, there is little margin to improve Lemma 1: at best, one can hope to replace the strict inequality in ( 9) with a large one. Can this actually be done ? This is partially answered by the following results:

Lemma 5 Assume that à ≡ A, that (4) holds with 2k • ∆( Ã, A) = L ( 13
)
and that the filters in A have disjoint supports [START_REF] Bajwa | Compressed channel sensing: A new approach to estimating sparse multipath channels[END_REF]. Then, either à 0 > A 0 , or each row of à is obtained by permuting pairs of distinct filters a ij , a ij ′ from the corresponding row of A such that a ija ij ′ is proportional to a modulated and translated Dirac comb with 2k spikes.

For filter matrices with a single row, since à ≡ A means that the filters ã1j are permuted versions of a 1j , we obtain Corollary 2 Consider A with a single row (M = 1). Assume that à ≡ A, that (4) holds with

2k • ∆( Ã, A) = L ( 14
)
and that the filters in A have disjoint supports [START_REF] Bajwa | Compressed channel sensing: A new approach to estimating sparse multipath channels[END_REF]. Then à 0 > A 0 .

Inria

Excessive pessimism?

The counter-examples built in Lemmata 3-4, which are associated to Dirac combs, are highly structured. They provide worst case well-posedness bounds, but existing probabilistic versions of uncertainty principles (see, e.g., the nice survey [START_REF] Tropp | On the linear independence of spikes and sines[END_REF]) lead us to conjecture that if the sparse filters in A are drawn at random (e.g. from Bernoulli-Gaussian distribution), the uniqueness guarantee of Theorem 2 will hold except with small probability O(L -β ), provided that k < c(β)L/ log L, for large L. This is left to further theoretical investigation.

Numerical experiments

The results achieved so far are theoretical well-posedness guarantee, but do not quite provide algorithms to compute the potentially unique (up to global permutation) solution of the frequency permutation problem. We conclude this paper with the description of a relatively naive optimization algorithm, an empirical assessment of its performance with Monte-Carlo simulations, and a discussion of how this compares with the theoretical uniqueness guarantees achieved above.

Proposed combinatorial algorithm

Given a "permuted" matrix Ã, one wishes to find a set of frequency permutations yielding a new matrix A with minimum ℓ p norm. The proposed algorithm starts from A 0 = Ã. Given A n , a candidate matrix A n+1,π can be obtained by applying a permutation π at frequency ω n ≡ n [mod L]. Testing each possible permutation π and retaining the one π n which minimises A n+1,π p yields the next iterate A n+1 := A n+1,πn . The procedure is repeated until the ℓ p norm A n ceases to change. Since there is a finite number of permutations to try, the stopping criterion is met after sufficiently many iterations.

Choice of the ℓ p criterion

In theory, it could happen that the stopping criterion is only met after a combinatorially large number of iterations. However, the algorithm stops much sooner in practice. In fact, if we were to use the ℓ 0 norm, the algorithm would typically stop after just one iteration, because the ℓ 0 norm attains its maximum value M × N × L for most frequency permutations except a few very special ones. For this reason, we chose to test the algorithm using ℓ p norms p > 0, which are not as "locally constant" as the ℓ 0 norm. To our surprise, the experiments below will show that the best performance is not achieved for small p, but rather for p = 2ǫ with small ǫ > 0. For p = 0 and for p ≥ 2, the algorithm indeed completely fails.

Monte-Carlo simulations

For various filter length L, sparsity levels k and dimensions M , N , random sparse filter matrices A made of independent random k-sparse filters were generated. Each filter was drawn by choosing: a) a support of size k uniformly at random; b) i.i.d. Gaussian coefficients on this support. For each configuration (L, k, M, N ), 200 such random matrices A were drawn. For each A, independent random frequency permutations were applied to obtain Ã. The algorithm was then applied to obtain A. The performance was measured using the SNR between A and the best permutation of A. 2) associated with the well-posedness guarantee (using an ℓ 0 criterion) of Theorem 2. Surprisingly, one can observe that the success rate increases when 0 < p < 2 is increased. The maximum success rate is achieved when p = 2ǫ with small ǫ > 0.

Beyond the well-posedness regime suggested by the theory (i.e., to the right of the vertical dashed line) the algorithm can succeed, but at a rate that rapidly decreases when the relative sparsity k/L increases. In the regime where the problem is proved to be well-posed, the proposed algorithm is often successful but can still fail to perfectly recover the filters, especially -and surprisingly-for small values of k. This phenomenon is strongly marked for p < 1 and essentially disappears for p > 1. It remains an open question to determine the respective roles of the ℓ p criterion and of the naive greedy optimization algorithm in this limited performance when the problem is well-posed with respect to the ℓ 0 norm.

Role of the filter length L

Figure 3 shows the results for different L values with p = 1.9, M = N = 2. One can see that the average performance does not seem to depend on whether L is prime or not. As L increases, the performance for "small" k/L slightly increases, but the success rate degrades for "large" k/L close to α(2).

Role of the number of channels M

Figure 4 shows the results for increasing numbers of channels M , with a filter length L = 512, N = 2 sources, p = 1.9. One can observe that the success rate substantially increases when M is increased from M = 1 to M = 2, and slightly increases as M further increases. Although the worst-case well-posedness guarantees are the same, the algorithm seems to benefit from added filter diversity across channels.

Role of the number of sources N

Figure 5 shows the success rate as a function of the relative sparsity k/L, for N ∈ {2, 3, 4}, with L = 31, M = 5 with p = 1.9. The well-posedness limits k/L ≤ α(N ) associated to Theorem 2 are indicated with vertical dashed lines. The empirical curves confirm that the algorithm can still succeeed beyond the worstcase well-posedness guarantees, but with a rapidly decreasing rate of success. When the well-posedness guarantees hold, the algorithm can fail, but its rate of success is high when the relative sparsity is sufficiently small compared to the bound provided by Theorem 2.

Computation time

The algorithm evaluates the ℓ p norm of the N ! permutations of the sources for each of the L frequencies. Hence, the complexity of each sweep through the set of all frequencies is O(N !M N L 2 log 2 L). This is rather expensive because the computational cost grows in factorial with the number of sources and in square with the filter length, but it is tractable for small problem sizes and very efficient compared to the brute force approach that would require O((N !) L-1 M N L log 2 L) operations to test all (N !) L-1 possible permutations up to a global permutation.

Figure 6 shows the average computation time over 200 trials for various filter length. The red dashed line corresponds to its prediction using the theoretical cost estimation as C × L 2 log 2 L with C ≈ 40 nanoseconds.

Conclusions

It is now well known that a sufficient sparsity assumption can be used to make under-determined linear inverse problems well-posed: without the sparsity assumption, the problem admits an affine set of solutions, which intersects at only one point with the set of sparse vectors. Besides this well-posedness property, a key factor that has lead to the large deployment of sparse models and methods in various fields of science is the fact that a convex relaxation of the NP-hard ℓ 0 minimization problem can be guaranteed to find this unique solution under certain sparsity assumptions. The availability of efficient convex solvers then really makes the problem tractable.

The problem considered in this paper is not a linear inverse problem. Even though it is a simplification of the original permutation and scaling problem arising from signal processing, it remains a priori a much harder problem than linear inverse problems in terms of the structure of the solution set: each solution comes with a herd of solutions that are equivalent up to a global permutation.

The fact that we managed to obtain well-posedness results in this context is encouraging, but this is at best the beginning of the story: even if the solution is unique, how do we efficiently compute it? Can one hope to extend these results to the original permutation and scaling problem? Why does the proposed naive algorithm perform better for p > 1? Answers to these questions are likely to have an impact in fields such as blind source separation with sparse multipath channels.

A Proof of Theorem 1

First, notice that for each frequency ω and channel i, permutations preserve the equality j a ij [ω] = j ãij [ω]. Thus, the same holds in the time-domain : j a ij = j ãij . By the disjoint supports hypothesis and the quasi-triangle inequality for ℓ p quasi-norms we have

j a ij p p = j a ij p p = j ãij p p ≤ j ãij p p . (15)
We conclude by summing over all channels i.

B Proof of Lemma 2 and Corollary 1

Let us consider the matching count matrix C with entries Theorem 3 (Birkhoff -Von Neumann Theorem, [START_REF] Birkhoff | Three observations on linear algebra[END_REF][START_REF] Neumann | A certain zero-sum two-person game equivalent to the optimal assignment problem[END_REF]) Every bi-stochastic matrix is in the convex hull of permutation matrices.

C jn := ♯{0 ≤ ℓ < L : σ ℓ (j) = n}, 1 ≤ j, n ≤ N.
Theorem 4 (Carathéodory Theorem [START_REF] Berger | Géométrie[END_REF]) Let X be a non empty subset of an affine space of dimension n ≥ 1. Then every element of the convex hull of X is a convex combination of p elements of X, with p ≤ n + 1.

The set of bi-stochastic matrices is an affine subspace of R N 2 . It is defined by 2N equations, but these equations are linearly dependent since the sum of the sum over rows is the sum of the sum over columns. Hence its affine dimension is n ≤ N 2 -(2N -1) = (N -1) 2 , and we conclude from Carathéodory's theorem that every bi-stochastic matrix is a convex combination of at most (N -1) 2 + 1 permutation matrices. One of the coefficients of this combination must therefore exceed 1/(1 + (N -1) 2 ), and this leads to a version of Lemma 2 with 2α ′′ (N ) = 1 1+(N -1) 2 , as claimed. Yet, this bound is suboptimal. The optimal bound in Lemma 2 follows from Hall's Marriage Theorem, which by the way is also a key ingredient in the proof of the Birkhoff-Von Neumann theorem.

Inria

Theorem 5 (Hall's Marriage Theorem [START_REF] Hall | On representatives of subsets[END_REF][START_REF] Oxley | Matroid theory[END_REF] ) Let (A j ) j∈J be a family of subsets of a set finite S. There exists a bijection π : J → S such that π(j) ∈ A j for all J if, and only if, for all

E ⊂ J ♯ ∪ j∈E A j ≥ ♯E
The bijection π is often referred to as a transversal for S.

Proof (Lemma 2) For shortness of notation we write α for α(N ). Define the sets J = S = 1, N , and A j := {n : B jn ≥ 2α}, j ∈ J, and consider the property

P k : ∀E ⊂ J, ♯E ≤ k ⇒ ♯ ∪ j∈E A j ≥ ♯E.
We wish to prove that P k holds true for all 1 ≤ k ≤ N : then, by Hall's Marriage Theorem, there exists a bijection π : j → π(j) such that π(j) ∈ A j for all j, yielding in turn the permutation matrix P with ones at the entries (j, π(j)). We proceed by contradiction: assume that P N does not hold true. Since P 1 holds true, without loss of generality, for some 1 ≤ k 0 < N :

♯ ∪ 1≤j≤k0 A j ≥ k 0 , and ♯ ∪ 1≤k≤k0+1 A j ≤ k 0 .
Hence, without loss of generality:

∪ 1≤k≤k0 A j = 1 k 0 ⊃ A k0+1 .
It follows that for n > k 0 and j ≤ k 0 + 1, we have n / ∈ A j , hence B jn < 2α. Now we use the bi-stochasticity of B ( j B jn = n B jn = 1, B jn ≥ 0) to obtain

k 0 ≥ n≤k0 j≤k0+1 B jn = j≤k0+1 n≤k0 B jn = j≤k0+1 1 - n>k0 B jn > j≤k0+1 (1 -(N -k 0 )2α) = (k 0 + 1)(1 -(N -k 0 )2α) = k 0 + 1 -(k 0 + 1)(N -k 0 )2α .
This implies 2α > 1/(k 0 + 1)(Nk 0 ). However, this yields a contradiction, since a simple functional study shows that

max 1≤k0<N 1 (k 0 + 1)(N -k 0 ) = 2α.
Equipped with Lemma 2, we can now prove Corollary 1.

Proof (Corollary 1) Since C = L • B where B is bi-stochastic, there is a permutation π such that C jπ(j) ≥ 2Lα(N ).

We conclude this section by showing the sharpness of Corollary 1 through the construction of permutations that reach the bound. Consider N an integer, and k 0 := N/2 (N even) or k 0 := (N -1)/2 (N odd). Let L be a multiple of (k 0 + 1)(Nk 0 ). Consider the L × N matrix: where: a) the left part, of size L × (k 0 + 1), is filled with the column vectors i ∈ R L/(k0+1) made of constant entries equal to the integer 1 ≤ i ≤ k 0 and the vector U ∈ R L/(k0+1) made of the vertical concatenation of the Nk 0 column vectors j ∈ R L/(k0+1)(N -k0) with constant entries k 0 + 1 ≤ j ≤ N ; b) the rows of the the right part, of size L × (Nk 0 -1), include exactly once each integer 1 ≤ ℓ ≤ N which does not already appear in the corresponding row of the left part. By construction, the L rows of the matrix Σ are associated to L permutations σ ℓ . We now show that, for any global permutation π, there is at least one column 1 ≤ j ≤ k 0 + 1 such that

Σ :=          1 U k 0 . . . 2 × . . . × 2
♯{ℓ : σ ℓ (j) = π(j)} ≤ L/(k 0 + 1)(N -k 0 ) = L2α(N ).
Applying again the pigeonhole principle yields: among the k 0 + 1 indices j to consider, at least one, j ⋆ , must be mapped to an integer π(j ⋆ ) ≥ k 0 + 1. By construction, the columns of Σ are such that column j ⋆ contains at most (in fact: exactly) L/(k 0 + 1)(Nk 0 ) instances of the value π(j ⋆ ).

C Proof of Theorem 2

We can now conclude the proof of Theorem 2. By Corollary 1, there is a permutation π such that for each j, we have F(ã ija iπ(j) ) 0 ≤ L(1 -2α(N )), hence ∆( Ã, A|π) ≤ L(1 -2α(N )) and finally ∆( Ã, A) ≤ L(1 -2α(N )). Combined with the assumption k ≤ Lα(N ), we obtain 2k+∆ ≤ L, and we conclude thanks to Lemma 1.

D Proof of Lemmata 1 3 4 5

We prove Lemma 1 first, then the statements of Lemma 3 in the following order: 1), 3), 2). We begin by some notations and fact regarding Dirac combs.

D.1 Dirac combs

Let p, q ≥ 1 be two integers and L = pq their product. The unit Dirac comb with q spikes and of step p, denoted x p , is the vector of C L defined by x q [t] = 1/ √ q if t ≡ 0[p], x q [t] = 0 otherwise. Its Fourier transform is the unit Dirac comb with p spikes and of step q: Fx q = x p . For 0 ≤ n < p an integer translation index and 0 ≤ m < q an integer modulation index, one can define the translated and modulated Dirac comb x q,n,m = T n M m x q where T n is the circular shift by n samples, and M m is the frequency modulation (M m u)[t] := u[t] • e 2iπmt/L .

Inria

One can check that the collection {x q,n,m } 0≤n<p,0≤m<q is an orthonormal basis of C L .

D.2 Proof of Lemma 1

Let π 0 be the permutation such that ∆ 0 (A, Ã) = min π∈SN ∆ 0 (A, Ã|π). By abuse of notation we still denote A the matrix obtained by applying π 0 to permute the columns of the original filter matrix. For each channel i and a source index j such that a ij = ãij we obviously have a ij 0 ≤ ãij 0 . Now, since ∆ 0 ≥ 1 we have à ≡ A hence there exists a pair i, j such that ãij = a ij . By the ℓ 0 Dirac-Fourier uncertainty principle [12, Theorem 1], for any vector u ∈ C L we have u 0 Fu 0 ≥ L. Hence, by the hypothesis k < L/(2∆ 0 ) we have

a ij 0 + ãij 0 ≥ ãij -a ij 0 (16) ≥ L/ F(ã ij -a ij ) 0 (17) ≥ L/∆ 0 > 2k (18) ≥ a ij 0 + a ij ′ 0 (19)
where j ′ is an arbitrary source index. Hence for every i, j such that ãij = a ij and any j ′ , ãij 0 > a ij ′ 0 , and we obtain

ãij 0 > max j ′ a ij ′ 0 ≥ a ij 0 .
Overall, we have shown that à 0 > A 0 . When L is prime, a stronger uncertainty principle u 0 + Fu 0 ≥ L + 1 holds [START_REF] Tao | An uncertainty principle for cyclic groups of prime order[END_REF]. Hence, under the assumption 2k + ∆ 0 ≤ L we can replace ( 17)-( 18

) with . . . ≥ L + 1 -F(ã ij -a ij ) 0 ≥ L + 1 -∆ 0 > 2k
to reach the same conclusion.

D.3 Proof of Lemma 3

We shall simply build an example where A = [α, β] is a 1 × 2 matrix of filters. Extensions to A an M × N matrix are trivial by adding mutually distinct sparse columns that are distinct from α and β, and duplicating the first row. We exploit Dirac combs as described in Appendix D. we obtain that, after permuting the Fourier transforms of α and β at the frequencies ω = 2kr, 0 ≤ r < L/2k,

α = b + u β = a + T L/2 u
We choose the vector u to be zero everywhere with two exceptions u

[0] := -a[0], u[ L 2k ] := -b[ L 2k ]. Since T L/2 u = u and a = b, we have {α, β} = { α, β} and à ≡ A. Moreover, ∆ 0 ( Ã, A) = ∆ 1 ( Ã, A) = L/2k.
Lastly, all considered vectors have k entries of equal magnitude, hence α 0 = β 0 = α 0 = β 0 = k, and for any 0 < p ≤ ∞ α p = β p = α p = β p . In particular, Ã p = A p , 0 ≤ p ≤ ∞.

D.4 Proof of Lemma 5

We repeat the construction of the proof of Lemma 3 starting from the Dirac combs a = x k ′ ,0,0 , b = -x k ′ ,L/2k ′ ,0 . Since k ′ < k ≤ L/2, we have ℓ := kk ′ ≤ L/2k ′ hence we can choose an ℓ-sparse vector u which support is outside the support of x 2k ′ and such that T L/2 u and u have disjoint supports. The four vectors {a, b, u, T L/2 u} have mutually disjoint supports, hence α and β have disjoint supports, {α, β} = { α, β} and à ≡ A. Moreover, ∆ 0 ( Ã, A) = ∆ 1 ( Ã, A) = L/2k. Lastly, we have α 0 = β 0 = α 0 = β 0 = k ′ + ℓ = k, and the ℓ p norms of these vectors are also equal, hence à p = A p , 0 ≤ p ≤ ∞.

D.5 Proof of Lemma 4

As in the proof of Lemma 1 we consider A the permuted matrix associated to the optimal permutation π 0 . Using the inequality 2k ≤ L/∆ 1 ≤ L/2∆ 0 instead of 2k < L/∆ 0 we repeat the steps ( 16)- [START_REF] Oxley | Matroid theory[END_REF] to obtain ãij 0 ≥ a ij ′ 0 for any j ∈ E i := {j, a ij = ãij } and any j ′ . As a result ãij 0 ≥ a ij 0 for all i, j. The assumption that à 0 = A 0 implies that ãij 0 = a ij 0 for all i, j.

By assumption, Ã ≡ A hence there are indices i, j such that a ij = ãij . For such i, j, since a ij 0 = ãij 0 , each inequality in ( 16)-( 19) (the inequality L/∆ 0 > 2k being replaced with L/∆ 1 ≥ 2k) must be indeed an equality. This implies that: a ij 0 = ãij 0 = k; 2k divides L and ∆ 1 = L/2k; the nonzero vector b ij := ãija ij must be an equality case of the ℓ 0 uncertainty principle with b ij 0 = 2k and Fb ij 0 = L/2k. As a result [START_REF] Tao | An uncertainty principle for cyclic groups of prime order[END_REF] b ij is a scaled, modulated and translated version of the Dirac comb x 2k made of 2k Diracs spaced every L/2k samples: there exists a scalar γ ij = 0, and two integers 0

≤ n ij < L/2k, 0 ≤ m ij < 2k such that b ij = γ ij • x 2k,nij ,mij .
Moreover since a ij 0 = ãij 0 = k and ãija ij 0 = 2k, the filters ãij and a ij have disjoint supports of size k. Hence, they are the restriction of b ij (resp. of -b ij ) to their respective supports. Now, define E i,n,m := {j ∈ E i , n ij = n, m ij = m}.

Inria

As observed in the proof of Theorem 1, the equality j a ij = j ãij holds, implying j∈Ei b ij = j b ij = 0. Taking inner products with the Dirac comb orthonormal basis x 2k,n,m , 0 ≤ n < L/2k, 0 ≤ m < 2k, yields j∈Ei,n,m

γ ij = 0, (20)
Since γ ij = 0, whenever E i,n,m is not empty it contains at least two distinct indices.

By the disjoint support assumption: for j, j ′ ∈ E i,n,m , j = j ′ , the original filters a ij and a ij ′ have disjoint supports. Moreover, we know that these supports are subsets of the support of x 2k,n,m which is of size 2k, hence

♯E i,n,m • k = j∈Ei,n,m a ij 0 ≤ 2k.
Hence, whenever E i,n,m is not empty, it contains exactly two distinct elements: E i,n,m = {j, j ′ } where j = j ′ .

Further, observe that: a) a ij and a ij ′ have disjoint supports of size k which are subsets of the support of size 2k of x 2k,n,m ; b) a ij and ãij have the same property. As a result, ãij and a ij ′ have the same support, which is disjoint from that of a ij . Similarly, a ij has the same support as ãij ′ . Finally, Eq. ( 20) can be rewritten γ ij + γ ij ′ = 0, and implies b ij + b ij ′ = 0, that is to say ãij + ãij ′ = a ij ′ + a ij . We conclude that ãij = a ij ′ and ãij ′ = a ij .

Figure 1 :

 1 Figure 1: Histogram of SNR between best permutation of A and original A

Figure 1

 1 Figure1shows the histogram of SNR values achieved for L = 31, 1 ≤ k ≤ L, M ∈ {1, 2}, N ∈ {2, 3, 4}, p = 1. It shows that the algorithm either completely succeeds up to machine precision (SNR above 300 dB) or completely fails (SNR of the order of 0 dB). For this reason, in the rest of the experiments the estimation was considered a success when the SNR exceeded 100 dB.

Figure 2 : 9 0Figure 3 :

 293 Figure 2: Filter recovery success as a function of p, 0 ≤ p ≤ 1.9

Figure 4 :Figure 5 :

 45 Figure 4: Filter recovery success as a function of M , for p = 1.9

Figure 6 :

 6 Figure 6: Computation time of the permutation solving algorithm depending on the length L of the filter

  Since j C jn = n C jn = L we have C = L • B where B is bi-stochastic. A weakened version of Lemma 2, with 2α ′′ (N ) = 1 1+(N -1) 2 , can be obtained by combining the Birckhoff -Von Neumann theorem and Carathéodory theorem.

1 .

 1 Define a = x k,0,0 , b = -x k,L/2k,0 . The filters a and b have disjoint support and satisfy a 0 = b 0 = k. Since ab = √ 2 x 2k,0,0 we have a[ω] = b[ω] whenever ω ≡ 0[2k]. Hence, permuting the Fourier transforms of a and b on the L/2k frequencies {ω = 2kr, 0 ≤ r < L/2k} yields ã = b and b = a. Given any u ∈ C L we define perturbations α and β of a and b α := a + u β := b + T L/2 u with T L/2 a circular shift. Noticing that for ω = 2kr (T L/2 u)[ω] = e 2iπ(L/2)ω L u[ω] = e 2iπkr u[ω] = u[ω]

  InriaWell-posedness of the permutation problem in sparse filter estimation with ℓ p minimization 11

			Average success rate with 200 draws for L=131, N=2, M=5
		100				p=0 p=0.1
	Success rate (%)	40 60 80 20	← α(2)			p=0.3 p=0.5 p=0.7 p=0.9 p=1 p=1.1 p=1.3 p=1.5 p=1.7 p=1.9
						p=2.1
		0
		0	0.2	0.4	0.6	0.8	1
				k/L

RR n°7782

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgment

This work was supported by the EU Framework 7 FET-Open project FP7-ICT-225913-SMALL: Sparse Models, Algorithms and Learning for Large-Scale data, and by Agence Nationale de la Recherche (ANR), project ECHANGE (ANR-08-EMER-006).