Yann Busnel
email: yann.busnel@univ-nantes.fr

Nothing can compare with a population, besides agents

Keywords: Population Protocols, Multi-Agents System, Model Equivalence

Leveraging the resemblances between two areas explored so far independently enables to provide a theoretical framework for distributed systems where global behaviors emerge from a set of local interactions. The contribution of this paper arise from the observation that population protocols and multi-agent systems (MAS) bear many resemblances. Particularly, some subclasses of MAS seem to fit the same computational power than population protocols. Population protocols provide theoretical foundations for mobile tiny device networks. On the other hand, from long-standing research study in distributed artificial intelligence, MAS forms an interesting model for society and owns a broad spectrum of application field, from simple reactive system to social sciences. Linking the both model should offers several extremely interesting outcomes.

Introduction

With the recent expansion of distributed system the two last decades, a boundless set of models has been developed to capture diversity and richness of these systems. Among this set, some models, while proposed in different context, catch the same range of computability. Some recent studies deal with the comparison and equivalence between models [START_REF] Basu | Systematic correct construction of self-stabilizing systems: a case study[END_REF][START_REF] Bertier | On gossip and populations[END_REF][START_REF] Marzougui | A new formalism for modeling a multi agent systems: Agent petri nets[END_REF][START_REF] Perathoner | Influence of different system abstractions on the performance analysis of distributed real-time systems[END_REF][START_REF] Poslad | Specifying protocols for multi-agent systems interaction[END_REF]. In this paper, we investigate the connections between two famous models: Population Protocols and Multi-Agent Systems (MAS).

Population protocols [START_REF] Angluin | Computation in networks of passively mobile finite-state sensors[END_REF] provide theoretical foundations for distributed systems in which global behavior emerges from a set of simple interactions between their agents. Originally developed in the context of mobile tiny devices, typically sensors, in this model, agents are considered anonymous, and therefore, undistinguishable. Many variants of population protocols have been proposed [START_REF] Angluin | Fast computation by population protocols with a leader[END_REF][START_REF] Angluin | The computational power of population protocols[END_REF][START_REF] Angluin | Self-stabilizing population protocols[END_REF][START_REF] Aspnes | An introduction to population protocols[END_REF][START_REF] Delporte-Gallet | When birds die: Making population protocols fault-tolerant[END_REF]. Among them, community protocols [START_REF] Guerraoui | Names trump malice: Tiny mobile agents can tolerate byzantine failures[END_REF] augment the original model by assigning agents a unique identifier and letting nodes remember a limited number of other identifiers. This not only significantly increases the computation power of the system but also provides a way to tolerate a bounded number of byzantine failures. In the sequel, the class of population protocols and variants will be referred as population protocols, and the original model as basic population protocol.

More specifically, the population protocol model consists in a finite space of agent's states, a finite set of inputs, a finite set of outputs and a transition function. The set of possible node's interactions is represented by a graph. When two agents are sufficiently close for a sufficiently long time, they interact by exchanging their local information, and update their state according to the transition function. For instance, if agents are small devices embedded on animals, an interaction takes place each time two animals are in the same radio range. The interaction patterns, orchestrated by a scheduler, are considered as unpredictable. Yet, the scheduler is assumed to be fair i.e., it ensures that any reachable global system state can be reached infinitely often. In the absence of global knowledge, agents cannot usually verify that the protocol has terminated, therefore the model considers convergence (of the distributed output) rather than termination.

On the other hand, largely inspired by the social behavior analysis of insect colony, Multi-Agent System model aims at studying the concept of collective and/or distributed intelligence [START_REF] Ferber | Multi-Agent Systems[END_REF]. The former model and its several extensions should be viewed as a crossroads in artificial intelligence (AI), distributed AI, distributed systems, software engineering, and smart objects.

As proposed in [START_REF] Ferber | Multi-Agent Systems[END_REF] and enhanced in [START_REF] Sycara | The many faces of agents[END_REF], a MAS is made up of an environment, wherein a set of active objects, so called agents, interact with passive objects. A set of relation links objects between them, according to their activity, as a set of operations defines what is possible to apply to objects. Finally, some operators are in charge to reflect the operation mapping on the environment (usually denoted "universe laws"). It exists a specific case where the environment is empty and its does not exists passive object. This kind of system is called purely communicating MAS, as the set of relations defines a network between agents, for which actions are reduced to communication. This specific structure is commonly used in distributed AI, mostly for collaboration between units design for problem solving of expert systems.

Usually, MAS are studied owing to their noteworthy ability for self-organization, and are sometimes surprisingly able to reproduce certain form of complex social systems. MAS could be split in two main classes: Cognitive and Reactive agents. These classes come from the design of the active agents, according to their algorithmic ability. More precisely, the cognitive class is commonly used for distributed AI as the reactive one permits to study virtual life [START_REF] Gechter | A reactive agent-based problem-solving model: Application to localization and tracking[END_REF]. In this paper, we mostly consider the latter one that shares several similarities with the Population Protocol model. Recent interest in the probabilistic convergence of MAS raises [START_REF] Martinez | A convergence result for multiagent systems subject to noise[END_REF][START_REF] Tang | Convergence of a class of multi-agent systems in probabilistic framework[END_REF]. A motivation of our work relies on the fact that the latter convergence has been extensively studied in the Population Protocol model [START_REF] Angluin | Fast computation by population protocols with a leader[END_REF][START_REF] Angluin | The computational power of population protocols[END_REF], and results should be simply bring into some MAS models.

As in [START_REF] Bertier | On gossip and populations[END_REF] in which authors bridge the gap between population protocols and gossip-based protocols, the aim of this paper is then to correlate the two afore-mentioned model: Population Protocols and Multi (reactive) Agents Systems, and more specifically, the ones made up of eco-agents. They both rely on finitestate agents. Both aim at achieving an emerging global behavior from a set of local interactions in a fully decentralized manner. The main contribution of this paper is to acknowledge these similarities and leverage them in both contexts.

This paper is organized as follow. In Section 2, we present a classification of Multi-Agent Systems. Second, Section 3 provides some background on basic population and community protocols. Then, we prove their equivalence in Section 4 and discuss some opportunities of leveraging it, before concluding in Section 5.

Different Kinds of Agent

The Kenetic [START_REF] Ferber | Multi-Agent Systems[END_REF] claims to be the science and technics of artificial organizations. Using MAS as a designing tool, this field of research possesses the global control of all the experimental system parameters as an advantage. Yet, if the designer is able to tune any possible parameters, she has no access to any determinism aspect, regarding the system evolution. Quite the contrary, as changing configuration in MAS are subject to chaotic phenomena. This implies that any variation in initial parameter, or any adding of insignificant random variable, should lead the same system to exactly opposite final state. In fact, these tiny changes are dramatically magnified through interactions between agents, avoiding some precise system state prediction.

All the Kenetic problems are located at a meeting point of the notion of agent and of society. In other words, it is at a crossroads in the relationship between individual behavior and globally observed phenomena. In this context, cooperation, conflict, collaboration and coordination of actions make sense. Let us restrain the possible field by defining what we mean by agent and multi-agent system.

MAS Model

MAS

are not yet subject to any law. Any work in this model is located in the center of a duality: Agent and Organization. Any organization comes from agent's interactions, but the local behavior of these agents should be modified in return from the constraints imposed by the created organization structures. Except predefined organizations by designers, task coordination, assignment and distribution are resulted from agents themselves. Some emerging properties appear, without any initial programming. It is necessary to introduce two schools of thought according to the model appliance, mainly relying on architectural differences; these differences concern the algorithmic nature of used agents. One school of thought is specialized in distributed AI, so called cognitive school, as the other one, denoted reactive school, is more interested in virtual life. The one considered in this paper corresponds to this reactive one.

Cognitive agents Cognitive MAS design aims to reach communication and cooperation in classical expert systems. In this case, the MAS will be made up of a small number of "intelligent" agents. Each one owns a basic knowledge that includes a set of mandatory information and technical know-how in order to succeed in its task, and to manage interactions with its environment and other agents. They are sometimes called intended, i.e., they own some goals, and explicit plans permitting to reach them. In the context of cognition, mostly planned, MAS leads to some group of individuals, govern by predefined social rules (for instance, in case of conflict, agents should be forced to negotiate). Researches in this field often rely on sociological studies about organization and small groups. In this paper, we do not consider further this kind of MAS.

Reactive Agents

The main class considered in this paper relies on the fact that a globally intelligent system should not required individual intelligence of its agents. Simple reaction mechanism to stimuli should resolve complex problem, without short-term plan, or explanation of aims and objectives. These MAS permit to raise some self-organization, as the famous example of ant's colony. In the latter, as all ants are undistinguishable, and without any global authority, coordination raised from agents' action in order to develop and save the colony. This collective entity is able to resolve complex issues as construct the anthill, find food, take care of eggs and larva, etc. [START_REF] Drogoul | Manta: New experimental results on the emergence of (artificial) ant societies[END_REF]. Emerging organization is principally related to reactive MAS, as they are not characterized by predefined structures (which is seldom, if ever, the case of cognitive MAS).

Formalism

Multi-Agent System

We propose to model a MAS as proposed in [START_REF] Ferber | Multi-Agent Systems[END_REF].

Definition 1. A Multi-Agent System (MAS) is a system made up of:

-An environment E, i.e., generally a metric space; -A set of objects O. These objects are usually located, i.e., at each time, it is possible for any object to associate it with a position in E. Most of these objects are passive, i.e., they can be sensed, created, removed or modified by agents; -A set of agents A such that A ⊆ O, which represent active entities of the system; -A set of relations R that links objects (and thus agents) between them.

-A set of operations Op that lets agents in A to perceive, produce, consume transform and manipulate objects in O. -Operators responsible for applying these operations and the world reaction to this modification. We denote these reactions Universe Laws.

It exists a specific case of MAS where A = O (i.e., there is no passive objects) and where E = ∅. This kind of system is call purely communicating MAS, as the set of relations R defines a network between agents, for which actions are reduced to communication. This specific structure is commonly used in distributed AI, mostly for collaboration between units and designed for problem solving of expert system. This organization emulates the modus operandi of a social structure as administration for instance.

Computational power of specific MAS relies on the potential of the environment [START_REF] Weyns | Environment as a first class abstraction in multiagent systems[END_REF]. The environment gives structure to the MAS, by providing a kind of shared memory, or a physical/spatial structure. In this study, we mainly focus on virtual environment using discrete representation. The environment could also abstract the communication structure, which is composed by an infrastructure for message passing, stigmergy or implicit communication. In cognitive MAS, it should also provide some social structure, defining roles, groups and communities of agents. We introduced then the notion of autonomous environment. An autonomous environment is defined as an environment with extended functionality compared with a simple spatial structure. The simplest way to model the environment is to provide only a metric space in order to estimate the ability of communication of any two agents, according to their distance between their respective locations. Any enhanced environment will be denoted as autonomous in the following.

Agent The core entity of MAS is the agent itself. Let us introduce a generic definition: Definition 2. An agent is a virtual or physical entity that:

is able to act in an environment; -is able to communicate directly toward other agents; -is directed by a set of trends (as individual objectives or satisfaction function, or even survival); -owns proper resources; -is able to perceive (in a limited manner) its environment; -owns a partial representation of this environment (possibly empty); -owns abilities and should offer services; -can be able to reproduce; -has a behavior striving for satisfying its aims, by tacking into account available resources and skills, according to its perceptions, representations and communications.

As introduced in Section 2.1, agents should be cognitive or reactive. However, this splitting is sometimes to simplistic. We need to enhance it according to two axes, which are summarized in Table 1: -Agent's comportments are usually leaded by their teleonomic behavior1 in contrast with the reflex infers by perception. The trends that govern agents should come either from the environment, or explicitly expressed in agents.

We mention it as a reflex behavior in the first case, and a teleonomic one in the second one.

❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤

Eco Problem Solving

The Eco Problem Solving (EPS) is a decentralized approach of problem solving, using emergence of interacting reactive agents [START_REF] Drogoul | Eco-problem-solving model: results of the n-puzzle[END_REF]. A stable state is then sought and considered as the solution of the problem. Any eco-agent has only one purpose: to be satisfied. Its environment is composed by the other neighboring agents. The latter agents should be in an acquaintance or in a dependency network, interacting using perception of actions. Two specific perceptions of eco-agent exist: being aggressed and being hampered. Three different actions can then be chosen: reach satisfaction (carrying the action that reach its objective), aggress another agent, or take flight. It can be represented as a finite state automaton, with four-intern state: Satisfied (S), Seeking satisfaction (SS -usually the entry state), Seeking to take flight (ST ) and Taking flight (T F ). A transition relation is then set between these states.

The actions made by an agent in its acquaintance or its dependency network follow the rules hereafter:

if an agent is satisfied, it informs his dependencies that they should be satisfied also; -if an agent cannot be satisfied yet, it seeks the hampering agents among its acquaintance and aggresses them; -if an agent seeks to take flight, it seeks the hampering agents among its acquaintance and aggresses them.

EPS can be applied in very different context [START_REF] Drogoul | Eco-problem-solving model: results of the n-puzzle[END_REF], from task scheduling (using "task agents" and "resources agents"), to production line optimization, including simulation of emergent structural turbulence stability in fluid dynamics. Depending on the applicative requirement, eco-agent should own an identifier or not. The latter case infers a specific class of EPS so-called Anonymous EPS in the rest of this paper.

However, the scope of possibilities of eco-agent is not restricted to problem solving. It is especially adapted to evolving universe simulation. In fact, in an eco-agent system, an external stimuli is treated as any other entry of the problem. Any agent acts locally in response of a perturbation. It is perfectly illustrated by the Misachieving baby problem [21]. It introduces some malicious external adversaries, which are very powerful and raise the difficulty of eco-resolution.

About Population Protocols

In this section, we briefly present the basic population protocol model and the community protocol variant, which relaxes the assumption on the anonymity of agents.

Basic population protocol

The basic population protocol model, initially introduced in [START_REF] Angluin | Computation in networks of passively mobile finite-state sensors[END_REF], is composed of a collection of agents, interacting pairwise in an order determined by a fair scheduler. Each agent has an input value and is represented by a finite state machine. This agent can only update its state through an interaction. Updates are defined by a transition function that describes the function f computed by the system. At each interaction, the agents compute an output value from their current state and converge eventually to the correct output value, depending to the inputs initially spread to the agents.

More formally, a population protocol is composed of:

a complete interaction graph Λ linking a set of n ≥ 2 agents; -a finite input alphabet Σ; -a finite output alphabet Y ; -a finite set of possible agent's states Q; -an input function ι : Σ → Q mapping inputs to states; -an output function ω : Q → Y mapping states to outputs; -a transition relation δ :

Q × Q → Q × Q on pairs of states.
In the sequel, we call (p, q) → (p ′ , q ′ ) or (p, q, p ′ , q ′ ) a transition if (p, q, p ′ , q ′ ) ∈ δ. A transition can occur between two agents' states only if these two agents have an interaction. The protocol is deterministic if δ is a function (i.e., at most one possible transition for each pair in Q 2 ).

A configuration of the system corresponds to an unordered multi-set, containing states of all agents. We denote C → C ′ the fact that a configuration C ′ can be obtained from C in one step (i.e., with only one transition for one existing interaction). An execution of the protocol is a finite or infinite sequence of population configurations

C 0 , C 1 , C 2 , . . . such that ∀i, C i → C i+1 .
As introduced above, the order of the interactions is unpredictable, and is decided by the scheduler. The scheduler is assumed to be fair, i.e a feasible configuration cannot be endlessly ignored. In other words, if a configuration C appears an infinite number of times during an execution, and there exists a possible step C → C ′ , then C ′ must also appear an infinite number of times in the execution. This ensures that any attainable configuration is eventually reached.

Community protocols

Many variants of the last model exist. In this paper, we focus on the specific extension so-called community protocol [START_REF] Guerraoui | Names trump malice: Tiny mobile agents can tolerate byzantine failures[END_REF], which significantly increases the computational power of the basic population protocol. This model augments the basic population protocol model by assigning unique identifiers to agents. All possible identifiers and a special symbol ⊥ are grouped in an infinite set U . The difference between basic population protocols and community protocols is the definition of the set of states: Q = B × U d where B is the initial definition of the population protocol's set of states collapsed to a memory of d identifiers. As in population protocols, algorithms cannot use any bound on the number of agents and moreover, U is infinite. In order to maintain the population protocol spirit in this extended model, some constraints are added: only existing agent identifiers can be stored in the d slots intended for identifiers of an agent's state and no other structural information about identifiers can be used by algorithms. We consider, for q ∈ Q and id ∈ U , that id ∈ q means that q stores id in one of its d identifier slots. Thus, community protocols have to verify the two following formal constraints:

∀(q 1 , q 2 , q ′ 1 , q ′ 2 ) ∈ δ, id ∈ q ′ 1 ∨ id ∈ q ′ 2 ⇒ id ∈ q 1 ∨ id ∈ q 2 . (1)
Consider π a permutation of U with π(⊥) = ⊥. For q = b, u 1 , u 2 , . . . , u d ∈ Q, let π(q) = b, π(u 1 ), π(u 2 ), . . . , π(u d ) We assume that:

∀(q 1 , q 2 , q ′ 1 , q ′ 2 ) ∈ δ : (π(q 1 ), π(q 2 ), π(q ′ 1 ), π(q ′ 2 )) ∈ δ. (2)
In short, the first assumption ensures that no transition introduce new identifiers and the second one that identifiers can only be stored or compared for equality, but not manipulated in any other way. Any population protocol can be viewed as a community protocol with d = 0.

Finally, a population or community protocol stably computes a function f : Σ + → Y if ∀n ∈ N, ∀σ ∈ Σ n , every fair execution with n agents initialized with the elements of σ, eventually stabilizes to output f (σ). That means that the output value of every agent eventually stabilizes to f (σ).

Computational power

In order to characterize computable functions of population protocols, Angluin et al. study a set of decidable predicate restriction. In fact, the basic population protocol model characterization is precisely identical than the set of semi-linear predicate [START_REF] Angluin | Computation in networks of passively mobile finite-state sensors[END_REF][START_REF] Angluin | Stably computable predicates are semilinear[END_REF]. A set of semi-linear is a subset of N d made up of finite union of linear set such that {

- → b + k 1 - → a 1 + k 2 - → a 2 + • • • + k m -→ a m }, where - → b is the d- dimension reference vector, -
→ a 1 to -→ a m are basis vectors, and k 1 to k m are nonnegative coefficients. A semi-linear predicate on an entry set corresponds to a predicate that is validated exactly on a semi-linear set.

In [START_REF] Angluin | Computation in networks of passively mobile finite-state sensors[END_REF], Angluin et al. prove that basic population protocols are able to compute any predicate define using Presburger's arithmetic, which precisely match with the set of all the semi-linear predicates. Presburger's arithmetic is a firstorder theory of the natural number, in the Peano's arithmetic language without the multiplication operation (i.e., with only the addition operation, eventually ordered, besides zero and successor operator). Moreover, any semi-linear predicate should be compute using a population protocol [START_REF] Angluin | Stably computable predicates are semilinear[END_REF].

Guerraoui and Ruppert provide a proof that their extended model drastically increases the latter power [START_REF] Guerraoui | Names trump malice: Tiny mobile agents can tolerate byzantine failures[END_REF]. Indeed, CPs are equivalent to the class of the symmetric functions in N SP ACE(n log n); i.e., any language decided using a non-deterministic Turing machine using O(S log S) slots, with respect to the following condition. Consider a function f : Σ + → Y . f is symmetric if, for all string y obtained by swapping character from another string x, we have f (x) = f (y). In other words, ∀x ∈ Σ + , ∀π permutation of x,

f (x) = f (π(x)) ⇔ f is symmetric.

Population vs. Reactive Agents

Recently, the probabilistic convergence of MAS has been studied apart [START_REF] Martinez | A convergence result for multiagent systems subject to noise[END_REF][START_REF] Tang | Convergence of a class of multi-agent systems in probabilistic framework[END_REF]. A motivation of our work relies on the fact that the latter convergence has yet been extensively studied in the Population Protocol model [START_REF] Angluin | Fast computation by population protocols with a leader[END_REF][START_REF] Angluin | The computational power of population protocols[END_REF], and these results should be simply bring from one model to another.

In this section, we prove the classification and the relations between all the considered models that are summarized in Figure 1 where R-MAS-AE (resp. R-MAS-NAE) means Reactive MAS with Autonomous Environment (resp. without Autonomous Environment). This provides a refined classification of MAS and population protocols. Let us split the proof in three main theorems.

First of all, the two left horizontal relations of Figure 1 are straightforward by extension of [START_REF] Bertier | On gossip and populations[END_REF], as they are based on the fact that increasing the model with unique agent identifier allows to extend the computational power of the former model. Thus, we obtain that Anonymous-EPS ≺ EPS and PP ≺ CP. On the other hand, the rightest vertical relation between R-MASs is also straightforward as an autonomous environment can obviously simulated a non-autonomous environment, using a restriction on the environment of the R-MAS-AE. The reverse side is impossible as agents in R-MAS-NAE are not able to deterministically compute the unique evolution of their environment. Thus, we trivially obtain that R-MAS-NAE ≺ R-MAS-AE. In the following sections, we provide proofs of remaining relations illustrated on Figure 1.

Anonymous EPS weaker than Basic Population Protocol

We show there that there exists a population protocol P able to simulate any execution of eco-resolution.

Theorem 1. If f is computable by a protocol from Anonymous EPS, then there exists a basic population protocol which can compute f . Proof. Let E an EPS using eco-agents that computes a specific function f . As presented above, a finite-state machine models eco-agents.

Mapping the domain of the transition function: The domain of any relation in R or operation in Op is finite (and corresponds to some Cartesian products of O, the set of objects, with D S the set of possible states of an object in the environment). Moreover, as elements in R∪Op are functions, theirs ranges are also finite by definition. Based on these sets, we define D E (g) a specific subset of the Cartesian product between the domain and the codomain of g ∈ R ∪ Op (i.e., D E (g) contains each ordered pair such that the first entry is in the domain of g and the second entry is the mapped element of this first entry by g). Thus, for any g ∈ R ∪ Op, D E (g) is finite and contains all the possible transitions of object states, based on the knowledge of a surrounding object sub-state in the environment.

Design of the basic population protocol for the purpose of simulation: Consider the following basic population protocol P, represented by the 7-uplet (Λ, Σ, Y, Q, ι, ω, δ). Consider a complete interaction graph Λ. Let the set of agent states in P be identical to the set of object states, i.e., Q = D S . Consider that Σ and Y are the same as the input and output sets of E, if they exist. In this case, ι and ω are the same functions than the ones which respectively associate the input set of E to D S , and D S to the output set of E. Conversely, if no specific input and output sets are defined in E, then Σ = Y = D S and ι ≡ ω corresponds to the identity function. Finally, the transition function δ is defined as follows.

∀g ∈ R ∪ Op, ∀(s l , s r , s ′ l ) ∈ D E (g), ∃(s r , s l , s ′ r ) ∈ D E (g) s.t. (s l , s r , s ′ l , s ′ r ) ∈ δ.
On the environment and the fair scheduler: Communication exchange between agents is an inherent characteristic of EPS. However, in this class, the environment is only used to model the ability for two agents, or an agent and an object, to pairwise interact. Therefore, as the scheduler of P is mandatory fair, any possible interaction eventually happens and makes sufficient condition to obtain a correct execution of E, using the aforementioned P.

Thus, there exists a basic population protocol P, which simulates the considered EPS E and computes the function f . Then, for any function computable by an EPS, there exists a basic population protocol, which stably computes this function.

⊓ ⊔

On the other hand, the reverse of Theorem 1 is not valid. In fact, there exists some functions that are computable by a population protocol, but no anonymous EPS is able to compute the same function. This is mainly due to the global fairness property associated to the Population Protocol model. In some EPS, it is possible to construct a loop-based infinite execution that endlessly avoid a given reachable system state, belying the aforementioned fairness property.

EPS weaker than Community Protocol

Along the same lines, we prove here the following theorem in order to prove that it is possible to simulate any EPS using a community protocol.

Theorem 2. If f is computable by a protocol from EPS, then there exists a community protocol which can compute f . Proof. Let E an EPS using eco-agents that computes a specific function f . Thus, each eco-agent in the system is aware of its unique identifier. Consider the following community protocol C.

Preliminary assumptions on C: We assume that, in the community protocol used on C, agents are uniquely identified, and a unique agent is assigned a specific identifier id L . This agent is considered as the leader by all other agents. In this specific community protocol, we assume that this leader is aware of the size of the system2 (denoted n in the sequel).

Design of the community protocol for the purpose of simulation: The only difference between a basic population protocol and a community protocol consists in the definition of the set of states (i.e., Q = B × U d ) and the two constraints on the state's identifier part (i.e., the part belonging to U d cannot be used freely). Similarly, the only difference between Anonymous-EPS and EPS remains on the use of identifiers.

Then, due to the proof of Theorem 1, the community protocol C has to be designed to simulate a given EPS E. Then, we can still consider the fair scheduler as the environment. In order to do that, the state of an agent in C is the same than for the eco-resolution one, enriched by its own coordinate in the environment.

Dealing with the environment simulation in C: For each interaction, δ is applied on involved agents if and only if these two agents are in adjacent coordinates. Moreover, some specific state required surrounding environment knowledge. For instance, to simulate the Taking Flight state, the corresponding agent must check if there is any space in front of it. In order to verify such property, C uses a synchronization barrier to check it as follows.

Consider an agent id that fall into ST state. In this state, it only waits until it encounters the agent id L . During its next interaction with id L , id sets its state to wait and id L store the current coordinates of id in its own state. For all futher interaction between id L and others agents, id L store the position of this object/agent and the latter also sets its state to wait. Thus, all agents eventually stabilize to the wait state, and the number of coordinated stored on id L eventually converges to n (the number of agents in the population). At this point, all agents have reached the synchronization barrier.

After the barrier, id L enables all agents to resume their own execution by release their wait state as described in [START_REF] Bertier | On gossip and populations[END_REF], and gradually remove stored coordinates after these interactions. The set of these stored coordinated eventually becomes empty and all agents eventually leave the wait state of the last barrier, and are ready for their next action. In the meanwhile, id L waits to first interact with id before starting this releasing phase. This ensure that id will be aware of the position of all the objects in the environment and is then able to check if there is an empty space in front of it.

In conclusion, if E computes the function f , then the community protocol C simulates the behavior of E and also computes the function f . ⊓ ⊔

Again, the opposite of Theorem 2 cannot be verified, due to the same argument about to the fairness property (see Section 4.1).

R-MAS-NAE is equivalent to Community Protocol

Finally, we prove here the following theorem in order to prove the equivalence between community protocols and R-MAS-NAE. Theorem 3. A predicate is computable by a community protocol if and only if it can be computed by a Reactive MAS without Autonomous Environment (R-MAS-NAE).

Proof. The proof of Theorem 3 is directly inferred from the statements of Lemmata 1 and 2, which show respectively both implications of this equivalence. ⊓ ⊔ Consider the first implication of this theorem. Based on the usage of the environment as the scheduler, the following lemma is almost straightforward. Lemma 1. For each f computable by a community protocol, there exists a R-MAS-NAE which compute f . Proof. Let C the community protocol computing f and defined by the 7-tuples (Λ, Σ, Y, Q, ι, ω, δ). Consider the R-MAS-NAE M described below. We have to show that M simulates C. First of all, each agent in Λ is hosted by a specific agent of M. Input, output and state sets are the same in M than in C. A fortiori, both map function ι and ω remain identical in M.

Dealing with the transition function: The M's relations R and operations Op are defined from δ: the universe laws is only composed of the possible transition extracted from δ. Moreover, we consider a purely communication MAS (see Section 2.2) where R is restricted to pairwise communication (i.e., each time agents interact with others, they could only operate a pairwise transformation state that exists in δ). Thus, the sequence of population configurations is valid and represents the community protocol C, as it only stems from the transition function δ using pairwise interactions.

On the fairness assumption: The last assumption to verify is the fairness condition. In R-MAS-NAE, the scheduler is fully defined by the (nonautonomous) environment evolution, which, as we mentioned before, potentially leads to any possible communication according to the current distance between each agent. In that context, every possible finite scheduling has a non-null probability to happen. Thus, every possible transitions between two system configurations C → C ′ has a non-null probability to happen. The fairness assumption is then verified.

We are then able to conclude that M simulates the community protocol C, which computes the function f . ⊓ ⊔ Finally, let now show the opposite of Lemma 1, corresponding to the second part of Theorem 3. Lemma 2. For each f computable by a R-MAS-NAE, there exists a community protocol which computes f . Proof. Let M the given R-MAS-NAE and consider the following community protocol C.

Summary of agent state requirement: Due to space constraints, we do not present the formal mapping between possible agent states and the transition function, which are trivially extended from the proof of Theorem 2.

Dealing with the environment simulation in C: Consider a specific agent in C that store the whole state of the environment of M, which has a finish state (as the set of object is finished). Following the outline of the rendez-vous mechanism presented in the proof of Theorem 2 (and in more details in [START_REF] Bertier | On gossip and populations[END_REF]) with the environment agent as the leader, it becomes simple to simulated M with C. Indeed, as the environment is not autonomous, there is no change in the latter one without direct interaction with an object of M. Thus, each action of any agent of M should infer a synchronization barrier, an update of the environment (on the aforementioned dedicated agent), and finally, a global release that lets the other agent opportunity to also interact with the environment.

In conclusion, if M computes f , then the community protocol C simulates the behavior of M and also computes the function f . ⊓ ⊔

We then have proved all the claims presented in Figure 1.

Conclusion and future works

The main contribution of this paper is to establish a correlation between Population Protocols and Multi (reactive) Agents Systems. Both aim at achieving an emerging global behavior from a set of local interactions in a fully decentralized manner. This parallel between two worlds, explored so far independently, offers several extremely interesting outcomes, acknowledge these similarities and leverage them in both contexts. These results can be leveraged for existing results as well as results to come in both areas. For instance, the intensive research around structured language for MAS (as KQML [START_REF] Finin | Kqml as an agent communication language[END_REF] or FIPA ACL [START_REF] Fipa | Fipa acl message structure specification[END_REF] for instance) should be mapped to the population theory. That should lead to a standardization of high-level communication between agents. From the inverse, the theoretical analysis of Population Protocols models could resolve the issue of a lack of computational power analysis with different MAS models.

Fig. 1 .

 1 Fig. 1. Relationship between Multi-Agent Systems and Population Protocols

Table 1 .

 1 Characterization of agents -Agent's relation with the environment raises the classical problem of subject/object. Is an agent have an explicit and symbolic representation of the surrounding world, on which she is able to reason or this representation is sub-symbolic, i.e., integrated in its sensorimotor intelligence. In the first case, we mention cognitive agent, and in the second, reactive agents.

	Behavior	Environment relation	Cognitive agents	Reactive Agents
	Teleonomic	Intentional Agents	Instinctual Agents
	Reflex		"Modules" Agents	Tropic Agents

"Teleonomy" represents the scientific concept of end in itself (purposefulness and goal-directedness of structures and functions in living organisms).

This assumption is only expected for the synchronization barrier. It can be relaxed using the probabilistic clock phase mechanism proposed for population protocols in[START_REF] Angluin | Fast computation by population protocols with a leader[END_REF].

Acknowledgment

The author would like to warmly thank Philippe Lamarre (INSA Lyon, France) for his valuable discussions at the early steps of this work.


Simon Henrot

Member, IEEE Charles Soussen

Manuel Dossot

Member, IEEE David Brie

Does deblurring improve geometrical hyperspectral unmixing?

Keywords: Hyperspectral unmixing, Minimum-volume simplex, Deconvolution I. INTRODUCTION

In this paper, we consider hyperspectral unmixing problems where the observed images are blurred during the acquisition process, e.g. in microscopy and spectroscopy. We derive a joint observation and mixing model and show how it affects endmember identifiability within the geometrical unmixing framework. An analysis of the model reveals that nonnegative blurring results in a contraction of both the minimum-volume enclosing and maximum-volume enclosed simplex. We demonstrate this contraction property in the case of a spectrally-invariant point-spread function. The benefit of prior deconvolution on the accuracy of the restored sources and abundances is illustrated using simulated and real Raman spectroscopic data.

signatures are mixed in each pixel because of limited spatial resolution. Spectral unmixing usually refers to the estimation of endmembers and their fractional contribution to each pixel, named abundances.

Many unmixing methods in the literature can be categorized as belonging to either a geometrical, a statistical (Bayesian) or a sparse regression framework [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF]. In the former approach, each spectral vector belongs to a simplex whose vertices are the endmembers we seek [START_REF] Craig | Minimum-volume transforms for remotely sensed data[END_REF]. Some algorithms assume that there exists at least one pure pixel (PP) per endmember and exploit this property by identifying the endmembers to the "purest" pixels in the hyperspectral scene. When the pure pixel hypothesis does not hold, the geometrical approach to unmixing then consists in fitting the Minimum Volume Simplex (MVS) to the data using one of many MV algorithms. However, highly mixed data are beyond the reach of geometrical algorithms because spectral signatures are located near the center of the true endmember simplex. One could then resort to statistical or sparse methods to unmix the data set. In this paper, we focus instead on analyzing and alleviating the effects of hyperspectral image blurring, which will be shown to possibly cause highly mixed data.

Hyperspectral images are usually degraded by the observation process, typically including blurring and corruption by noise. Data acquired by spectro-imagers dedicated to geoscience, remote sensing or planetology, result from a physics-based inversion pipeline ensuring that these images are free of blur. However, in other domains, the acquisition process does not include a deblurring step, resulting in blurry hyperspectral images. For instance, microscopy [START_REF] Zhang | Gaussian approximations of fluorescence microscope point-spread function models[END_REF], Raman spectroscopy [START_REF] De Grauw | Axial resolution of confocal raman microscopes: Gaussian beam theory and practice[END_REF] or ground-based astrophysical [START_REF] Bourguignon | Processing MUSE hyperspectral data: Denoising, deconvolution and detection of astrophysical sources[END_REF], [START_REF] Villeneuve | PSF estimation of hyperspectral data acquisition system for ground-based astrophysical observations[END_REF] images are blurred by diffraction through the objective aperture. We consider such hyperspectral images in this paper. Mathematically, each observed image plane is obtained as the twodimensional convolution of the original image and the channel linear shift-invariant blur or Point-Spread Function (PSF). If the size of the PSF is less than that of an image pixel, the blur does not alter the spatial resolution of the data. However, if the PSF support extends beyond a pixel, neighboring pixels are combined in the observed data, resulting in additional mixing.

The latter case arises whenever one aims at obtaining a fine spatial resolution, e.g. in super-resolution reconstruction [START_REF] Akgun | Super-resolution reconstruction of hyperspectral images[END_REF]. In this framework, some scenes that could a priori be geometrically unmixed may lose this property because of the blur. The operation of deconvolution which inverts the observation process may then restore the stack so that low computational complexity geometrical methods qualify to unmix it.

In addition, deconvolution algorithms also provide the additional benefit of restoring images with higher signal-to-noise ratios (SNRs).

The ideal approach to account for these two physical phenomena would be to jointly deconvolve and unmix the data. Let us mention [START_REF] Zhao | Deblurring and sparse unmixing for hyperspectral images[END_REF] where the authors perform joint deconvolution and abundance estimation. They assume that endmembers belong to a large dictionary and the abundance estimation is formulated as a sparse regression problem. However, to the best of our knowledge, no algorithms have been designed to perform joint deconvolution and blind unmixing1 of the data. Nevertheless, many off-the-shelf algorithms are available to perform deconvolution and unmixing separately. We propose to show in this paper that the sequential procedure of (deconvolution/unmixing) yield better results than a (unmixing/deconvolution) strategy.

To this end, we will study how blurring may affect the results of geometrical unmixing methods and assess whether deblurring hyperspectral images before unmixing them yields quantitative improvements.

The remainder of the paper is organized as follows. Section II presents the linear mixing model and recalls its associated geometrical framework ; we also discuss identifiability issues by the geometrical unmixing approach. In section III, we factor in the general observation model and present experimental results on synthetic data to evaluate the gain of a deconvolution step before applying a spectral unmixing algorithm. Section IV focuses on the specific case where the PSF is invariant across the spectral channels.

Finally, we present experimental results on real data in section V and we conclude in section VI.

II. HYPERSPECTRAL GEOMETRICAL UNMIXING FRAMEWORK

In this section, we introduce concepts and notations used throughout the rest of the paper. Most of the results given below are not new and can be found elsewhere in the literature: see [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF] for a thorough review of hyperspectral unmixing and [START_REF] Chan | A convex analysis-based minimum-volume enclosing algorithm for hyperspectral unmixing[END_REF] for an introduction to convex geometry applied to endmember extraction.

A. Linear mixing model

Consider an hyperspectral image measuring radiance on L different spectral bands (channels) and N pixels. We gather the data in a L × N matrix X and use the following notations:

1) x ℓ is the ℓ-th row of x, that is the 2D image at channel ℓ after lexicographical ordering into a row vector of length N ;

2) x n is the n-th column of x, i.e. the L × 1 spectrum of the n-th pixel (also termed spectral vector or pixel vector).

Each spectral vector in the image is a linear combination of an known number P of endmembers {s 1 , . . . , s P }. When unknown, P can be obtained by some model order estimation method such as virtual dimensionality [START_REF] Chang | Estimation of number of spectrally distinct signal sources in hyperspectral imagery[END_REF]. Ignoring noise for now, the linear mixing model (LMM) writes

X = SA
where the pth column of L × P source matrix S indexes endmember s p and the nth column a n of P × N abundance matrix A contains the fractional abundance coefficients for x n :

x n = Sa n = P p=1 A p,n s p . ( 1
)
The LMM is generally based on the following assumptions [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF]:

i) the number of endmembers P is much smaller than the number of bands L, that is P ≪ L;

ii) matrix S is of full column rank, i.e. endmembers {s 1 , . . . , s P } are linearly independent;

iii) Abundance Nonnegativity Constraint (ANC) : A p,n ≥ 0 for all p and n; iv) Abundance Sum Constraint (ASC) : P p=1 A p,n = 1 for all n.

Assumptions (i) and (ii) seem very reasonable in hyperspectral imaging since many bands are collected and the image is made up of a few distinct materials. Assumptions (iii) and (iv) come from the physical interpretation of abundance coefficient A p,n as the fractional spatial area occupied by the pth endmember in the nth pixel.

To model illumination variability due to surface topography [START_REF] Nascimento | Vertex component analysis: a fast algorithm to unmix hyperspectral data[END_REF], equation (1) sometimes writes x n = γ n Sa n where parameter γ n is a scale factor for each pixel in the hyperspectral scene. This phenomenon can be accounted for in model (1) by simply rejecting hypothesis (iv). The geometrical consequences of this choice are detailed below.

B. Dimension reduction and affine projection

From equation (1) and assumptions (i) and (ii), it is immediate that pixel vectors lie in a P -dimensional vector subspace of R L . Suppose the subspace orthonormal basis E = [e 1 , . . . , e P ] has been inferred from the data, e.g. by the HySime algorithm [START_REF] Bioucas-Dias | Hyperspectral subspace identification[END_REF], then the n-th dimension-reduced pixel is given by

xn = E T x n ∈ R P . ( 2
)
From equation [START_REF] Bioucas-Dias | A variable splitting augmented lagrangian approach to linear spectral unmixing[END_REF] and assumption (iii), P -dimensional pixel vectors are enclosed in the convex cone spanned by dimension-reduced endmembers {ŝ 1 , . . . , ŝP } [START_REF] Ifarraguerri | Multispectral and hyperspectral image analysis with convex cones[END_REF] (see figure 1): In order to satisfy hypothesis (iv), the P -dimensional points are then projected onto a (P -1)-dimensional affine plane (see figure 1). The choice of a particular affine plane is arbitrary and may be viewed as imposing the norm of the supporting endmembers.

{x|x = P p=1 θ p s p , θ p 0, p = 1 . . . P }. (3)
Most algorithms use either an orthogonal projection [START_REF] Chan | A convex analysis-based minimum-volume enclosing algorithm for hyperspectral unmixing[END_REF] or the perspective projection [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF] onto the (P -1)-dimensional affine plane : the choice of one method or another is related to the SNR, see [START_REF] Nascimento | Vertex component analysis: a fast algorithm to unmix hyperspectral data[END_REF].

For example, the perspective projection consists in scaling pixel vectors

xn = xn xT n u (4)
to project them onto the hyperplane defined by vT u = 1 for every vector v ∈ R P , where u ∈ R P + is chosen beforehand such that xT n u > 0 for all n. From equations (1), ( 2) and ( 4), the n-th projected pixel vector rewrites

xn = P p=1 α p,n sp (5)
where coefficients α p,n are given by

α p,n = A p,n ŝT p u xT n u = A p,n ŝT p u q A q,n ŝT q u . ( 6
)
It is easy to see that the α p,n 's satisfy the following properties:

1) α p,n ≥ 0 for all p and n, under the ANC;

2) p α p,n = 1 for all n.

It follows that each projected pixel vector xn is a convex combination of the projected endmembers {s p , p = 1 . . . P }. Hence, projected pixel vectors belong to the convex set S:

S = {x|x = P p=1
θ p sp |θ p ≥ 0, p = 1 . . . P,

P p=1 θ p = 1}. (7)
By definition, S is the simplex whose vertices are the projected endmembers {s 1 , . . . ,s P }.

C. Minimum Volume Simplex (MVS)

Most geometrical unmixing algorithms are based on the key idea introduced by Craig [9] that the minimum-volume simplex (MVS) enclosing all the pixel vectors identifies with S and thus the endmembers we seek can be extracted as the vertices of the MVS. The MV constraint allows to regularize the ill-posed nature of the nonnegative source separation problem [START_REF] Moussaoui | Non-negative source separation: range of admissible solutions and conditions for the uniqueness of the solution[END_REF]. Alternatively, one can look for the maximum-volume simplex inscribed in the dataset, provided the image contains pure pixels (PP).

Recovering both simplices can be addressed as a combinatorial problem but is too complex to handle high dimensional data [START_REF] Gritzmann | On the complexity of some basic problems in computational convexity: I. containment problems[END_REF], [START_REF] Zhou | Algorithms for a minimum volume enclosing simplex in three dimensions[END_REF]. Simple and computationally efficient approaches were proposed e.g.

[24], [START_REF] Nascimento | Vertex component analysis: a fast algorithm to unmix hyperspectral data[END_REF] in the MVS case and [START_REF] Craig | Minimum-volume transforms for remotely sensed data[END_REF], [START_REF] Chang | A new growing method for simplex-based endmember extraction algorithm[END_REF], [START_REF] Chan | A convex analysis-based minimum-volume enclosing algorithm for hyperspectral unmixing[END_REF] in the PP case, among many others. We point out that unlike in the PP case, the MVS and the true simplex S whose vertices are the endmembers may be different simplices. Figure 2 illustrates the distinction between the MVS and S in the simple case where P = 3, yielding a 2D affine plane. In what follows, we discuss the problem of endmember identifiability using the MVS approach. Craig's statement can be rewritten in the following way:

Statement 1. the MVS identifies with S.

The first problem with the statement is that there might exist more than one MVS: see figure 2b. This case is not likely to arise in practice if enough data points have been collected so that there exists data points near the vertices of S.

Assuming there exists only one MVS, the following result can be found in the literature [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches[END_REF]:

Sufficient condition 1.
If there is at least one data point on each vertex of S, then statement 1 is true.

This sufficient condition (SC) is known as the pure pixel condition. For P = 3, we make the following statement whose proof is given in appendix A (see figure 2c):

Sufficient condition 2.
If there are at least two data points per facet of S whose distance with the closest vertex is no greater than a third of the corresponding facet's length, then statement 1 is true.

SC 2 is weaker than SC 1. We note that a recent paper establishes necessary and sufficient conditions for the identifiability of endmembers by Craig's criterion in the 3-endmember case, see [START_REF] Lin | On the endmember identifiability of craigs criterion for hyperspectral unmixing: A statistical analysis for three-source case[END_REF]. Formulating such a statement for higher dimensions (that is, for more than three endmembers) is difficult and out of the scope of this paper. Rather, we stress that endmembers are more likely to be accurately identified by a MVS algorithm if there are pixel vectors near the vicinity of the vertices of S in the data. Obviously, if the data is highly mixed, the MVS will be enclosed in S and identifiability conditions will not be satisfied. Such a situation can arise because of two major reasons:

1) the abundance distributions themselves;

2) additional mixing due to the imaging system, as we will demonstrate in the next section.

In the latter case, this undesirable contribution can be removed by image restoration methods, thus improving endmember identifiability conditions.

III. EFFECT OF THE IMAGING PROCESS

A. Blurring of the linear mixture

In this section we now account for the fact that the image is degraded during the acquisition process.

Under the common linear blur assumption, the 2D image y ℓ observed at a given channel ℓ is obtained as the 2D convolution product of the true image and the channel point-spread function H ℓ :

y ℓ = x ℓ H ℓ (8)
where the N × N observation matrix H ℓ is a convolution matrix corresponding to H ℓ . For instance, when the blur is space-invariant, (H ℓ ) T is a block-Toeplitz matrix where each block is Toeplitz [START_REF] Henrot | Fast positive deconvolution of hyperspectral images[END_REF] 2 . Each entry of the observed data matrix Y is given by

Y ℓ,k = N n=1 X ℓ,n H ℓ n,k . (9)
Using equations ( 1) and ( 9), the overall model combining noise, observation blurring and linear mixing of endmembers writes

Y ℓ,k = N n=1 P p=1 S ℓ,p A p,n H ℓ n,k + W ℓ,k (10)
where W is the noise term and model [START_REF] De Grauw | Axial resolution of confocal raman microscopes: Gaussian beam theory and practice[END_REF] assumes that the SNR is high enough for the noise to be additive and i.i.d. Gaussian. We observe that the blurred data do not satisfy the LMM since the mixing coefficients A p,n H ℓ n,k depend on the channel index ℓ. Thus, even in the absence of noise, the y k 's do not necessarily belong to S anymore.

B. Pixel trajectories

Since the abundances vary with the channel index, it is no longer possible to use barycentric coordinates to represent a pixel vector in the simplex formed by the endmembers [START_REF] Honeine | Geometric unmixing of large hyperspectral images: A barycentric coordinate approach[END_REF]. To give some insight into the understanding of the effect of the channel dependent blur on the linear mixture, we propose to use the representation of figure 4. The endmembers correspond to the vertices of the simplex. Note that in this figure the simplices locations are arbitrarily fixed to (0, 0), (0, 1) and (1/2, √ 3/2) resulting in an isoceles triangle. As the mixing coefficients are depending on the channel index, a mixed pixel can no longer be represented as a single point. For each channel, we have different mixing coefficients. We rather represent a channel dependent mixed pixel by a trajectory where each point of the trajectory corresponds to the barycentric coordinates at a particular channel.

Figure 3 displays a few pixel trajectories for a Gaussian PSF whose width linearly increases with the channel index 3 . Interestingly, this representation provides spatial information about the pixel at hand. For instance, a trajectory that goes from a vertex to another may indicate that the pixel is located in the boundary of two regions, each corresponding to a particular endmember; a trajectory pointing towards the center of the simplex corresponds to a highly mixed pixel.

As a consequence, even if applying a linear unmixing algorithm directly to the data is feasible, the resulting endmembers and abundances are literally unmeaningful. We rather propose to perform deblurring prior to the unmixing step and numerical experiments are conducted in the next section to assess the gain of doing so.

To conclude this section, we note that blurring caused by a nonnegative PSF results in a contraction property whose extent increases with the PSF width. We will come back to that point in section IV.

C. Experimental results

To quantitatively assess the effects of deblurring on hyperspectral unmixing, we use synthetic data so that the extracted endmembers may be compared to their true counterparts. First, we select P = 3 true independent endmembers from the United States Geological Survey (USGS) spectral library 4 . Each endmember signature is sampled on 224 spectral bands. We then synthesize a set of P random twodimensional (2D) 30 × 30 abundance maps modeled as a mixture of 2D Gaussian functions with random locations and widths, under the constraint that the abundances for each pixel sum to one. The true data stack is then computed using equation (1). Image degradation is simulated by convolving each channel with a 2D Gaussian PSF whose width varies linearly with the channel index, from one pixel to the entire field of view. The approximation of physical PSFs by 2D Gaussian functions is actually reasonable for numerous optical systems, including fluorescence microscopy [START_REF] Zhang | Gaussian approximations of fluorescence microscope point-spread function models[END_REF], Raman microscopy [START_REF] De Grauw | Axial resolution of confocal raman microscopes: Gaussian beam theory and practice[END_REF] and atmospheric turbulence [START_REF] Zhao | Deblurring and sparse unmixing for hyperspectral images[END_REF]. We then add i.i.d. Gaussian noise at a specific SNR value. We wish to evaluate the effect of deconvolution for different noise levels. Hence, we run three separate experiments with respective SNRs of 5dB, 25dB and 50dB.

The rationale for these values is that at 5dB, the role of the deconvolution algorithm is to both denoise and deblur the data, while at 50 dB its main role will consist in deblurring. A SNR of 25dB is viewed as a more standard situation.

Image deconvolution is achieved using the algorithm described in [START_REF] Henrot | Fast positive deconvolution of hyperspectral images[END_REF] which accounts for spectral and spatial smoothness and non-negativity of multispectral and hyperspectral data. The rationale for this choice is that the priors of the methods are well adapted to the data and its low complexity allows January 30, 2014 DRAFT a computationally efficient simulation. This algorithm requires tuning two regularization parameters, weighting respectively a spatial and spectral smoothness prior, denoted by µ and ν. For each situation, we run the algorithm by setting these parameters on a 10 × 10 logarithmic grid, with values ranging from 10 -4 to 10 2 . Since we only evaluate hyperparameters on a finite grid, deconvolution of the data may not be optimal. However, the deconvolved image is always closer to the true image in the mean square error sense.

For each experiment, spectral unmixing is carried out on both the degraded data set and the deconvolved data set. Endmember extraction is performed using three widespread geometrical-based algorithms :

Vertex Component Analysis (VCA) [START_REF] Nascimento | Vertex component analysis: a fast algorithm to unmix hyperspectral data[END_REF], Minimum-Volume Enclosing Simplex (MVES) [START_REF] Chan | A convex analysis-based minimum-volume enclosing algorithm for hyperspectral unmixing[END_REF] and Simplex Identification via variable Splitting and Augmented Lagrangian (SISAL) [START_REF] Bioucas-Dias | A variable splitting augmented lagrangian approach to linear spectral unmixing[END_REF], using MATLAB code published by the authors. Abundance estimation is automatically carried out with MVES and we use Fully Constrained Least Squares inversion (FCLS) [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF] to compute abundances for VCA and SISAL. While VCA is based on the pure pixel hypothesis, both MVES and SISAL are minimum-volume algorithms, hence the choice of these three methods allows the scrutiny of a wide range of situations. The goal here is to evaluate the increase of performance brought by a deconvolution step for all algorithms, rather than comparing the results of one algorithm to another, which would not be fair as they depend on the situation at hand.

We choose the following performance metrics:

1) the spectral angle distance (SAD) between the true and extracted endmembers:

SAD p = cos -1 s T p ŝp s p 2 2 ŝp 2 2 , p = 1 . . . P. (11)
The total SAD is given by SAD = p SAD p ;

2) The Mean Square Error (MSE) between the true and estimated abundances, given by

MSE = Â -A 2 F . (12)
We then store the unmixing results on the raw data for each point of the hyperparameter grid and select the point which corresponds to the best extraction SAD. This choice is arbitrary : depending on the application, it may be better to select the deconvolution result yielding the smallest abundance estimation MSE. In either case, we point out that both these quantities could provide interesting criteria for hyperparameter selection.

Table I gives the various performances obtained for endmember extraction and abundance estimation by VCA, MVES and SISAL with and without a prior deconvolution step. Deconvolving the image prior to image unmixing improves the results for all algorithms.

IV. CASE OF PSF INVARIANCE

In this section, we focus on the specific case where the PSF is invariant across channels. This case is actually important in practice and arises in many imaging situations, such as Raman microscopy.

A. Bilinear model

In the specific case where the PSF is invariant across channels, the model reduces to the bilinear equation Y = SAH + W where each dimension-reduced pixel is given by

ŷk = N n=1 P p=1 A p,n H n,k ŝp + ŵk . (13)
The k-th projected pixel vector ẙk = ŷk /ŷ T k u is then a noisy bilinear combination of the projected endmembers {s p = ŝp /ŝ T p u, p = 1 . . . P }:

ẙk = N n=1 P p=1 α ′ p,n,k sp + ŵT k u ŷT k u ẘk . (14)
where coefficients α ′ p,n,k are given by

α ′ p,n,k = A p,n H n,k ŝT p u ŷT k u = A p,n H n,k ŝT p u q,m A q,m H m,k ŝT q u . ( 15
)
From equation ( 15), we can see that coefficients α ′ p,n,k satisfy the nonnegativity and sum-to-one properties:

1) α ′ p,n,k ≥ 0 for all p, n and k under the ANC and knowing that PSFs entries are nonnegative; 2)

p,n α ′ p,n,k = 1 for all k = 1 . . . N .

Thus, in the case where the PSF is invariant across channels and in the absence of noise, the observed pixel vectors {ẙ k , k = 1 . . . N } remain in the simplex S (defined in equation ( 7)). However, the additional mixing introduced by the PSF may hinder spectral unmixing in a different way, as shown in section IV-B.

B. MVS contraction

How does the observation process affect the distribution of pixel vectors inside the simplex? The answer to the question obviously depends on the nature of the PSF. Since the entries of H are known to be nonnegative, the blurring process tends to average neighboring pixel intensities. Hence:

1) suppose the data include a pure pixel x n : then xn locates on a vertex of S. If a mixed pixel x m is in the spatial neighborhood of x n , the observed y n will appear mixed in the observed image, which will relocate ẙn towards the inside of S;

2) similarly, if a source s p is not present in a given pixel x n but appears in its neighborhood, the observed pixel vector y n will include a contribution from s p and ẙn will move towards the inside of S.

This phenomenon inclines observed spectral vectors to cluster towards the center of S. As a graphical illustration, consider for a simple multispectral image with 3 × 3 = 9 pixels and 3 spectral bands, each blurred with a 3 × 3 averaging kernel:

H =      1/9 1/9 1/9 1/9 1/9 1/9 1/9 1/9 1/9      . ( 16
)
Suppose our data set satisfies the pure pixel hypothesis, and also includes some mixed pixels, yielding the following 3 × 9 abundance matrix:

A =      1 1/2 0 1/2 1/3 1/6 0 1/6 0 0 0 0 1/2 1/3 1/6 1 2/3 1/2 0 1/2 1 0 1/3 2/3 0 1/6 1/2      . ( 17
)
The resulting simplices for the original and degraded data are given in figure 7 Therefore S * x encloses all pixels {ẙ k , k = 1 . . . N }. Since S * y is the MVS enclosing {ẙ k , k = 1 . . . N }, we necessarily have

V(S * y ) V(S * x ) (19)
(Note however that while the volume of S * y is smaller, S * y is not necessarily included in S * x .) In the case of Gaussian PSFs, Proposition 1 has an interesting consequence for the variant PSF case. Proposition 2. Consider two symmetric separable 2D Gaussian PSFs H 1 and H 2 of respective width satisfying σ 1 > σ 2 . Let {ẙ 1 ,ẙ 2 } denote the respective observed projected pixels and {S * 1 , S * 2 } the respective MVSs. Then the convex hull of ẙ1 is included in the convex hull of ẙ2 and V(S * 1 ) V(S * 2 ).

Proof: The PSFs being separable, H 1 may be written as the 2D convolution product of H 2 and a Gaussian kernel of standard deviation σ 2 1 -σ 2 2 . Since all kernels are nonnegative, Proposition 2 is a straightforward application of Proposition 1.

Getting back to the idea of pixel trajectories outlined in section (III-B), when the PSF width increases with the channel, proposition (2) allows us to deduce that the convex hull of points at a given channel ℓ is included in the convex hull of points at channel ℓ -1. We now resume the discussion of the invariant PSF case.

Hence, we know from section II-C that directly applying a geometrical unmixing algorithm to the observed data may produce incorrect endmembers, and thus the subsequent estimation of abundances will also be biased. In the next section, we study the improvement of unmixing performance after applying a deconvolution algorithm to the data.

C. Experimental results

Here, we consider the effect of the standard deviation of the PSF on our performance metric, using the same synthetic data and experimental protocol as in section III-C. For a fixed SNR value, we consider each kernel size from one pixel (σ = 0.1) to the the entire field of view of the image (σ = 5). For each value of σ, we run the deconvolution algorithm for spatial and spectral regularization parameters on the previously used grid. We display our experimental results on figure 9 and note a few key points:

σ 1 σ 2 σ L s 1 s 2 s 3
Fig. 8. Illustration of the nesting property of convex hulls from channel to channel (variant PSF case).

1) As expected, the performance of both endmember extraction and abundance estimation are better after a deconvolution step, the margin between the two results tending to increase with σ except for a few localized values of σ. The reason why the unmixing algorithms perform slightly better on the raw data for these values certainly owes to the limited tuning of hyperparameters in the deconvolution step.

2) All unmixing results show a decrease in performance, both in terms of endmember extraction (increase of SAD) and abundance estimation (increase of MSE) when σ increases. Owing to the random initializations of the unmixing algorithms, the error curves are not purely monotonic. The decrease in performance seems very intuitive since the severity of the observed data degradation increases with σ.

3) The parameter search on the deconvolution step allows to stabilize the solution. Of course, the simulation makes use of ground truth while computing the best parameter values. In order to perform tuning in real imaging situations, one could use a library of spectral signatures and seek the parameters minimizing the spectral angle distance between the extracted endmembers and the corresponding ground truth endmembers that are a priori known to be present in the image. Let us sum up our experimental results in the invariant PSF case. A prior deconvolution step increases the performance of endmember extraction and abundance estimation for a vast majority of σ values. Of course, the wide variety of situations (abundance maps, endmembers, PSF spread, SNR, . . . ) make general statements difficult to quantify. However, in all our experiments, the unmixing results after deconvolution were always at least comparable and most of the time much better than the results obtained on the raw data.

V. EXPERIMENTAL RESULTS ON REAL DATA

In this section, we evaluate the gain of deconvolution on unmixing real Raman spectroscopy data. Our hyperspectral data set comprises images of size 98 × 131 pixels, each pixel being 100 nm × 100 nm, acquired on 337 bands ranging from 800 cm -1 to 1200 cm -1 . The scene of interest consists in a grain of sodium acetate (CH 3 COONa) covered with sodium carbonate (Na 2 CO 3 ) laying on a silicon layer (Si). Part of the sodium carbonate reacts with water vapor to yield hydrated sodium carbonate. These four chemical compounds are the endmembers we seek. A thorough inspection of the data reveals that the silicon compound contributes to all pixels of the image. The extraction of these endmembers is a challenging problem, since silicon is the only compound for which the pure pixel assumption is fulfilled.

Because of the inherent high mixing of the data, even MV methods are not supposed to produce good results on this data set.

Given the limited spectral range, the PSF is considered to be invariant across channels. It is modeled as a 2D Gaussian function [START_REF] De Grauw | Axial resolution of confocal raman microscopes: Gaussian beam theory and practice[END_REF] with an experimentally measured full-width at half maximum of 300 nm. We apply our deconvolution algorithm to the data by setting regularization parameters through a trial-and-error process to µ = 20 and ν = 5.

VCA and FCLS are carried out on both the raw data matrix and its restored counterpart ; the results are displayed in figure 11. The first extracted endmember corresponds to the silicon layer, which presents a broad band at 910-960 cm -1 due to the 2TO harmonic phonon mode of bulk silicon. The deconvolution step allows to denoise its abundance map and more importantly, uncovers structure that was distributed throughout other abundance maps. The second endmember is sodium acetate, with a peak at 930 cm -1 due to the intense C-C stretching mode of the acetate molecule [START_REF] Wang | Raman spectroscopic studies on single supersaturated droplets of sodium and magnesium acetate[END_REF]. It is well isolated by VCA (both for the raw and restored data), perhaps because the width of the grain make this endmember almost uniquely contribute to some pixels. However, the third and fourth endmembers present mixed spectral contributions from all compounds : both forms of sodium carbonate, hydrated and non hydrated (respective peaks at 1060 and 1080 cm -1 ) as well as silicon and sodium acetate. The deconvolution step again reveals structure in the fourth endmember that was not displayed in the case of raw data.

Further analysis of the spectra also show bands at 840 cm -1 and 1004 cm -1 , which reveals that the data have accidentally been polluted with sodium benzoate (NaC 7 H 5 O 2 ) [START_REF] Veij | Reference database of raman spectra of pharmaceutical excipients[END_REF]. Because this compound is associated with the sodium carbonate phase, experiments with 5 endmembers do not yield improvements, both compounds being mixed in a unique endmember.

Both MVES and SISAL extract endmembers with negative peaks from this data set, rending them physically unmeaningful, so we do not include results from these algorithms. However, we employ another minimum-volume algorithm : the Nonnegative Matrix Factorization (NMF) with volume constraints (NMF-vol) on the inferred simplex [START_REF] Miao | Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization[END_REF]. The algorithm operates on both the raw data (where negative pixels have been clipped to zero since the algorithm imposes a nonnegative data matrix) and restored data. Due to computation times, it was not possible to include this method in our previous simulations on synthetic data; here we set 10000 as the maximum number of iterations. The resulting abundance maps and endmembers are also given in figure 11. The spectral shapes of endmembers extracted from the raw data appear noisier, a problem solved by the deconvolution step. The sodium acetate endmember displays more undesirable contribution from the silicon compound. However, the algorithm is able to separate the sodium carbonate compound (third endmember, characterized by the peak at 1080 cm -1 ) from the hydrated sodium carbonate compound (fourth endmember, peak at 1060 cm -1 ; both are mixed with the sodium benzoate, as expected. The main gain of deconvolution clearly appears on the third endmember, where the silicon contribution is completely suppressed. Another benefit of the restoration step is to reveal structure hidden in raw abundance maps (first and third endmember).

VI. CONCLUSION

In this paper, we have analyzed the effect of a linear observation system on geometrical hyperspectral unmixing. When the observation system varies with the spectral channels, the joint mixing and observation model is no longer bilinear. Blurring by a nonnegative PSF results in the contraction of the minimum volume simplex enclosing the data. We formally prove the contraction property in the case of a spectrally constant PSF and show that its extent increases with the PSF width. Hence, unmixing blurred hyperspectral data is expected to significantly benefit from a deconvolution procedure, which we demonstrate in terms of endmember extraction and abundance estimations using both synthetic and real Raman spectroscopic data. Comparing the benefits and computational costs of a sequential or joint deconvolution and unmixing

approach remains an open problem, especially when one needs to resort to blind deconvolution methods.

Perspectives also include identifying application domains where the issue might have been overlooked, e.g.

VCA MVES SISAL SNR = 5dB

Fig. 1 .

 1 Fig. 1. Illustration in dimension P = 3 of the convex cone spanned by dimension-reduced endmembers {ŝp, p = 1 . . . P } enclosing dimension-reduced pixels xn's. xn is obtained by perspective projection of xn on the affine plane supported by {ŝp, p = 1 . . . P }.

Fig. 2 .

 2 Fig. 2. 2D examples of triangles enclosing the data for P = 3 endmembers. The true unknown simplex S is represented with a straight line while the other candidate triangle is plotted with a dotted line. Data points are represented as hollow circles and the vertices of the simplices are filled circles. (a): Statement (1) is true, i.e. S is the MVS. (b): both triangles are MVS's for this data configuration. (c): the dotted-lined triangle has a smaller volume than S and thus is the MVS.

Fig. 3 .

 3 Fig. 3. Examples of trajectories of a few pixel vectors.

Fig. 4 .

 4 Fig. 4. Selected endmembers.

Fig. 5 .

 5 Fig. 5. Examples of abundance maps generated as spatial mixtures of 2D Gaussian distributions with random means and variances.

Fig. 6 .

 6 Fig. 6. Original data simplex (with no blur and noise) obtained by Principal Component Analysis and orthogonal projection. Pixel vectors are denoted by the symbol '+', endmembers by circles.

Fig. 7 .

 7 Fig. 7. Left : original projected pixel vectors, right : degraded projected pixel vectors. The data cloud is shrinked by the observation process and pure pixels become mixed in the observed data.

Proposition 1 .

 1 . The data cloud clearly suffers a shrinkage inside of S from the observation process. Let us now formally state this contraction property, using the notion of convex hull. The convex hull of the data should not be confused with the MVS which necessarily has P vertices; in contrast, the convex hull has up to N vertices which all are pixels. Let S *x denote the MVS enclosing the true data X, S * y the one enclosing the blurred data Y and V denote the volume. If the PSF is nonnegative, then the convex hull of Y is enclosed in the convex hull of X andV(S * y ) V(S * x )(18)Proof: Assume the PSF coefficients are nonnegative. Then each observed pixel ŷk is a linear nonnegative combination of the original pixels {x n , n = 1 . . . N } (or possibly a subset of X depending on the support of the PSF). Hence ŷk belongs to the conical hull of X for each k. It follows that the convex hull of the projected observed pixel vectors {ẙ k , k = 1 . . . N } is included in the convex hull of {x n , n = 1 . . . N }.

Figure 10

 10 Figure 10 is a geometrical view of the results of figure 9 illustrating the distance between the unknown and recovered sources, both on the raw and restored data for three values of the PSF spread σ. The results clearly show the MVS contraction caused by the blurring process and the reverse effect of the deconvolution step.

Fig. 9 .Fig. 10 .Fig. 11 .

 91011 Fig.9. Performance of unmixing methods as a function of the PSF spread σ. First column : VCA, second : MVES, third : SISAL. For each SNR section, the first row displays the spectral angle divergence between the observed / restored extracted endmembers and their true counterparts (respectively solid and dotted line). The second row displays the mean square error between the observed / restored estimated abundances and their true counterparts.

Fig. 12 .

 12 Fig.[START_REF] Gritzmann | On the complexity of some basic problems in computational convexity: I. containment problems[END_REF]. S is the straight-lined triangle, S ′ is the dotted-line triangle. r = .2 in this illustration.

January 30, 2014 DRAFT

Here, blind unmixing refers to the joint estimation of endmembers and abundances.

The transpose comes from the fact that x ℓ and y ℓ are row vectors.

We restrict our analysis to the case of PSFs whose width increases with the channel index (i.e. wavelength) because this is the relevant case of conventional optical systems.January 30, 2014 DRAFT

Available online from: http://speclab.cr.usgs.gov/spectral-lib.html January 30, 2014 DRAFT

microscopy and ground-based astrophysical imaging, and addressing nonlinear unmixing schemes [START_REF] Honeine | Approches géométriques pour lestimation des fractions dabondance en traitement de données hyperspectales extensions aux modeles de mélange non-linéaires[END_REF]. APPENDIX A SC 2 PROOF Suppose an hyperspectral image contains three endmembers. After dimension reduction and affine projection (see section (II)), pixel vectors are enclosed in a triangle S whose vertices correspond to the endmembers. Suppose further that there are exactly two data points per facet of S and denote these six points by {P 1 . . . P 6 }. Each P i is assumed to be located at a distance r * d from the closest vertex on its facet, where d is the facet length and r ∈ [0, 1/2]. As illustrated in section (II-C), there are two possible triangles fitting {P 1 . . . P 6 }, S and S ′ .

We first perform an affine transformation so that S is transformed into a right isosceles triangle. The respective areas A(S) and A(S ′ ) are scaled by the same number (that is, the determinant of the transform matrix) and the distance ratio r is conserved since the transformation is affine. In the new domain, we have A(S) = 1/2. By first calculating the equations of the facets of S ′ and the coordinates of the vertices of S ′ , it is easy to compute A(S ′ ) as a function of r:

Solving A(S ′ ) > A(S) for r ∈ [0, 1/2] is trivial and yields r < 1/3.

ACKNOWLEDGMENT

This work was supported by the French Agence Nationale de la Recherche, with reference number ANR-09-BLAN-0336-04. The authors are very indebted to the anonymous reviewers for their insightful comments which significantly improved the contents of the paper.