Franc ¸ois Pacaud

Pierre Carpentier

Jean-Philippe Chancelier

Michel De Lara

Efficacity - † Ensta - ‡ Cermics-Enpc

Stochastic optimal control of a domestic microgrid equipped with solar panel and battery

Microgrids are integrated systems that gather and operate energy production units to satisfy consumers demands. This paper details different mathematical methods to design the Energy Management System (EMS) of domestic microgrids. We consider different stocks coupled together -a battery, a domestic hot water tank -and decentralized energy production with solar panel. The main challenge of the EMS is to ensure, at least cost, that supply matches demand for all time, while considering the inherent uncertainties of such systems. We benchmark two optimization algorithms to manage the EMS, and compare them with a heuristic. The Model Predictive Control (MPC) is a well known algorithm which models the future uncertainties with a deterministic forecast. By contrast, Stochastic Dual Dynamic Programming (SDDP) models the future uncertainties as probability distributions to compute optimal policies. We present a fair comparison of these two algorithms to control microgrid. A comprehensive numerical study shows that i) optimization algorithms achieve significant gains compared to the heuristic, ii) SDDP outperforms MPC by a few percents, with a reasonable computational overhead.

I. INTRODUCTION

A. Context

A microgrid is a local energy network that produces part of its energy and controls its own demand. Such systems are complex to control, because of the different stocks and interconnections. Furthermore, at local scale, electrical demands and weather conditions (heat demand and renewable energy production) are highly variable and hard to predict; their stochastic nature adds uncertainty to the system.

We consider here a domestic microgrid (see Figure 1), equipped with a battery, an electrical hot water tank and a solar panel. We use the battery to store energy when prices are low or when the production of the solar panel is above the electrical demand. The microgrid is connected to an external grid to import electricity when needed. Furthermore, we model the building's envelope to take advantage of the thermal inertia of the building. Hence, the system has four stocks to store energy: a battery, a hot water tank, and two passive stocks being the building's walls and inner rooms. Two kind of uncertainties affect the system. First, the electrical and domestic hot water demands are not known in advance. Second, the production of the solar panel is heavily perturbed because of the varying nebulosity affecting their production.

We aim to compare two classes of algorithms to tackle the uncertainty in microgrid Energy Management Systems (EMS). The renowned Model Predictive Control (MPC) algorithm views the future uncertainties with a deterministic forecast. Then, MPC relies on deterministic optimization algorithms to compute optimal decisions. The contender -Stochastic Dual Dynamic Programming (SDDP) -is an algorithm based on the Dynamic Programming principle. Such algorithm computes offline a set of value functions by backward induction; optimal decisions are computed online as time goes on, using the value functions. We present a balanced comparison of these two algorithms, and highlight the advantages and drawbacks of both methods.

B. Litterature 1) Optimization and energy management systems: Energy Management Systems (EMS) are integrated automated tools used to monitor and control energy systems. The design of EMS for buildings has raised interest in recent years. In [START_REF] Olivares | Trends in microgrid control[END_REF], the authors give an overview concerning the application of optimization methods in designing EMS.

The well-known Model Predictive Control (MPC) [START_REF] Garcia | Model predictive control: theory and practice-a survey[END_REF] has been widely used to control EMS. We refer notably to [START_REF] Oldewurtel | Stochastic model predictive control for energy efficient building climate control[END_REF], [START_REF] Malisani | Pilotage dynamique de l'énergie du bâtiment par commande optimale sous contraintes utilisant la pénalisation intérieure[END_REF], [START_REF] Lamoudi | Distributed Model Predictive Control for energy management in buildings[END_REF] for applications of MPC in buildings. Different solutions are investigated to tackle uncertainties, such as Stochastic MPC [START_REF] Oldewurtel | Stochastic model predictive control for energy efficient building climate control[END_REF] or robust optimization [START_REF] Paridari | Robust scheduling of smart appliances in active apartments with user behavior uncertainty[END_REF].

2) Stochastic Optimization: as we said, at local scale, electrical demand and production are highly variable, especially as microgrids are expected to absorb renewable energies. This leads to pay attention to stochastic optimization approaches. Apart from microgrid management, stochastic optimization has found some applications in energy systems (see [START_REF] Lara | Optimization Methods for the Smart Grid[END_REF] for an overview). Historically, stochastic optimization has been widely applied to hydrovalleys management [START_REF] Pereira | Multi-stage stochastic optimization applied to energy planning[END_REF]. Other applications have arisen recently, such as integration of wind energy and storage [START_REF] Haessig | Dimensionnement et gestion d'un stockage d'énergie pour l'atténuation des incertitudes de production éolienne[END_REF] or insulated microgrids management [START_REF] Heymann | A stochastic continuous time model for microgrid energy management[END_REF].

Stochastic Dynamic Programming (SDP) [START_REF] Bertsekas | Dynamic programming and optimal control[END_REF] is a general method to solve stochastic optimal control problems. In energy applications, a variant of SDP, Stochastic Dual Dynamic Programming (SDDP), has proved its adequacy for large scale applications. SDDP was first described in the seminal paper [START_REF] Pereira | Multi-stage stochastic optimization applied to energy planning[END_REF]. We refer to [START_REF] Shapiro | Analysis of stochastic dual dynamic programming method[END_REF] for a generic description of the algorithm and its application to the management of hydrovalleys. A proof of convergence in the linear case is given in [START_REF] Philpott | On the convergence of stochastic dual dynamic programming and related methods[END_REF], and in the convex case in [START_REF] Girardeau | On the convergence of decomposition methods for multistage stochastic convex programs[END_REF].

With the growing adoption of stochastic optimization methods, new researches aim to compare algorithms such as SDP and SDDP with MPC. We refer to the recent paper [START_REF] Riseth | A comparison of control strategies applied to a pricing problem in retail[END_REF].

C. Structure of the paper

We detail a modelling of the microgrid in Sect. II, then formulate an optimization problem in Sect. III. We outline the different optimization algorithms in Sect. IV. Finally, we provide in Sect. V numerical results concerning the management of the microgrid.

II. ENERGY SYSTEM MODEL

In this section, we depict the physical equations of the energy system model described in Figure 1. These equations write naturally in continuous time t. We model the battery and the hot water tank with stock dynamics, and describe the dynamics of the building's temperatures with an electrical analogy. Such physical model fulfills two purposes: it will be used to assess different control policies; it will be the basis of a discrete time model used to design optimal control policies.

A. Load balance

Based on Figure 1, the load balance equation of the microgrid writes, at each time t:

φ pv (t) + f ne (t) = f b (t) + f t (t) + f h (t) + d el (t) . (1)
We now comment the different terms. In the left hand side of Equation ( 1), the load produced consists of

• the production of the solar panel φ pv (t),

• the importation from the network f ne (t), supposed nonnegative (we do not export electricity to the network).

In the right hand side of Equation ( 1), the electrical demand is the sum of

• the power exchanged with the battery f b (t),

• the power injected in the electrical heater f t (t),

• the power injected in the electrical hot water tank f h (t),

• the inflexible demands (lightning, cooking...), aggregated in a single demand d el (t).

B. Energy storage

We consider a lithium-ion battery, whose state of charge at time t is denoted by b(t). The state of charge is bounded:

b ≤ b(t) ≤ b . (2)
Usually, we set b = 30% × b so as to avoid empty state of charge, which proves to be stressful for the battery. The battery dynamics is given by the differential equation

db dt = ρ c (f b (t)) + - 1 ρ d (f b (t)) -, (3)
with ρ c and ρ d being the charge and discharge efficiency and f b (t) denoting the power exchange with the battery. We use the convention f + = max(0, f ) and f -= max(0, -f ).

As we cannot withdraw an infinite power from the battery at time t, we bound the power exchanged with the battery:

-f b ≤ f b (t) ≤ f b . (4)

C. Electrical hot water tank

We use a simple linear model for the electrical hot water tank dynamics. At time t, we denote by T h (t) the temperature inside the hot water tank. We suppose that this temperature is homogeneous, that is, that no stratification occurs inside the tank.

At time t, we define the energy h(t) stored inside the tank as the difference between the tank's temperature T h (t) and a reference temperature

T ref h(t) = ρV h c p T h (t) -T ref , (5)
where V h is the tank's volume, c p the calorific capacity of water and ρ the density of water. The energy h(t) is bounded:

0 ≤ h(t) ≤ h . (6)
The enthalpy balance equation writes

dh dt = β h f h (t) -d hw (t) , (7)
where

• f h (t)
is the electrical power used to heat the tank, satisfying

0 ≤ f h (t) ≤ f h , (8)
• d hw (t) is the domestic hot water demand,

• β h is a conversion yield. A more accurate representation would model the stratification inside the hot water tank. However, this would greatly increase the number of states in the system, rendering the numerical resolution more cumbersome. We refer to [START_REF] Schütz | A comparison of thermal energy storage models for building energy system optimization[END_REF] and [START_REF] Beeker | Discrete-time optimal control of electric hot water tank[END_REF] for discussions about the impact of the tank's modeling on the performance of the control algorithms.

D. Thermal envelope

We model the evolution of the temperatures inside the building with an electrical analogy: we view temperatures as voltages, walls as capacitors, and thermal flows as currents. A model with 6 resistances and 2 capacitors (R6C2) proves to be accurate to describe small buildings [START_REF] Berthou | Development of building models for load curve forecast and design of energy optimization and load shedding strategies[END_REF]. The model takes into account two temperatures:

• the wall's temperature θ w (t),

• the inner temperature θ i (t). Their evolution is governed by the two following differential equations

c m dθ w dt = θ i (t) -θ w (t) R i + R s Exchange Indoor/Wall + θ e (t) -θ w (t) R m + R e Exchange Outdoor/Wall + γf t (t) Heater + R i R i + R s Φ int (t) Radiation through windows + R e R e + R m Φ ext (t) Radiation through wall , (9a) c i dθ i dt = θ w (t) -θ i (t) R i + R s Exchange Indoor/Wall + θ e (t) -θ i (t) R v Ventilation + θ e (t) -θ i (t) R f Windows + (1 -γ)f t (t) Heater + R s R i + R s Φ int (t) Radiation through windows , (9b)
where we denote

• the power injected in the heater by f t (t),

• the external temperature by θ e (t),

• the radiation on the wall by Φ ext (t),

• the radiation through the windows by Φ int (t). The time-varying quantities θ e (t), Φ int (t) and Φ ext (t) are exogeneous. We denote by R i , R s , R m , R e , R v , R f the different resistances of the R6C2 model, and by c i , c m the capacities of the inner rooms and the walls. We denote by γ the proportion of heating dissipated in the wall through conduction, and by (1 -γ) the proportion of heating dissipated in the inner room through convection.

E. Continuous time state equation

We denote by x = (b, h, θ w , θ i ) the state, u = (f b , f t , f h ) the control, and w = (d el , d hw , φ pv ) the uncertainties. The continuous state equation writes

ẋ = F (t, x, u, w) , (10)
where the function F is defined by Equations ( 3)-( 7)-( 9).

III. OPTIMIZATION PROBLEM STATEMENT

Now that we have described the physical model, we turn to the formulation of a decision problem. We aim to compute optimal decisions that minimize a daily operational cost, by solving a stochastic optimization problem.

A. Decisions are taken at discrete times

The EMS takes decisions every 15 minutes to control the system. Thus, we have to provide decisions in discrete time.

We set ∆ = 15mn, and we consider an horizon T 0 = 24h. We adopt the following convention for discrete processes:

for t ∈ {0, 1, • • • , T = T0
∆ }, we set x t = x(t∆). That is, x t denotes the value of the variable x at the beginning of the interval [t∆, (t + 1)∆[. Otherwise stated, we will denote by [t, t + 1[ the continuous time interval [t∆, (t + 1)∆[.

B. Modeling uncertainties as random variables

Because of their unpredictable nature, we cannot anticipate the realizations of the electrical and the thermal demands. A similar reasoning applies to the production of the solar panel. We choose to model these quantities as random variables (over a sample space Ω). We adopt the following convention: a random variable will be denoted by an uppercase bold letter Z and its realization will be denoted in lowercase z = Z (ω). For each t = 1, . . . , T , we define the uncertainty vector

W t = (D el t , D th t , Φ pv t ) , (11)
modeled as a random variable. The uncertainty W t takes value in the set W t = R 3 .

C. Modeling controls as random variables

As decisions depend on the previous uncertainties, the control is a random variable. We recall that, at each discrete time t, the EMS takes three decisions:

• how much energy to charge/discharge the battery F b t , • how much energy to store in the electrical hot water tank F h t , • how much energy to inject in the electrical heater F t t . We write the decision vector (random variable)

U t = (F b t , F h t , F t t ) , (12)
taking values in U t = R 3 . Then, between two discrete time indexes t and t + 1, the EMS imports an energy F ne t+1 from the external network. The EMS must fulfill the load balance equation ( 1) whatever the demand D el t+1 and the production of the solar panel Φ pv t+1 , unknown at time t. Hence F ne t+1 is a recourse decision taken at time t + 1. The load balance equation (1) now writes as

F ne t+1 = F b t + F t t + F h t + D el t+1 -Φ pv t+1 P -a.s. , (13)
where Pa.s. indicates that the constraint is fulfilled in the almost sure sense. Later, we will aggregate the solar panel production Φ pv t+1 with the demands D el t+1 in Equation ( 13), as these two quantities appear only by their sum.

The need of a recourse variable is a consequence of stochasticity in the supply-demand equation. The choice of the recourse variable depends on the modeling. Here, we choose the recourse F ne t+1 to be provided by the external network, that is, the external network mitigates the uncertainties in the system.

D. States and dynamics

The state becomes also a random variable

X t = (B t , H t , θ w t , θ i t ) . (14)
It gathers the stocks in the battery B t and in the electrical hot water tank H t , plus the two temperatures of the thermal envelope (θ i t , θ w t ). Thus, the state vector X t takes values in X t = R 4 .

The discrete dynamics writes

x t+1 = f t (x t , u t , w t+1 ) , (15)
where f t corresponds to the discretization of the continuous dynamics (10) using a forward Euler scheme, that is,

x t+1 = x t + ∆ × F (t, x t , u t , w t+1 ).
By doing so, we suppose that the control u t and the uncertainty w t+1 are constant over the interval [t, t + ∆[. The dynamics [START_REF] Riseth | A comparison of control strategies applied to a pricing problem in retail[END_REF] rewrites as an almost-sure constraint:

X t+1 = f t X t , U t , W t+1 P -a.s. . (16)
We suppose that we start from a given position x 0 , thus adding a new initial constraint: X 0 = x 0 .

E. Non-anticipativity constraints

The future realizations of uncertainties are unpredictable. Thus, decisions are functions of previous history only, that is, the information collected between time 0 and time t. Such a constraint is encoded as an algebraic constraint, using the tools of Probability theory [START_REF] Kallenberg | Foundations of Modern Probability[END_REF]. The so-called non-anticipativity constraint writes

σ(U t ) ⊂ F t , (17)
where σ(U t ) is the σ-algebra generated by U t and

F t = σ(W 1 , • • • , W t )
the σ-algebra associated to the previous history (W 1 , . . . , W t ). If Constraint [START_REF] Beeker | Discrete-time optimal control of electric hot water tank[END_REF] holds true, the Doob lemma [START_REF] Kallenberg | Foundations of Modern Probability[END_REF] ensures that there exists a function π t such that

U t = π t (X 0 , W 1 , . . . , W t ) . (18)
This is how we turn an (abstract) algebraic constraint into a more practical functional constraint. The function π t is an example of policy.

F. Bounds constraints

By Equations ( 2) and ( 6), the stocks in the battery B t and in the tank H t are bounded. At time t, the control

F b t must ensure that the next state B t+1 is admissible, that is, b ≤ B t+1 ≤ b by Equation (2), which rewrites, b ≤ B t + ∆ ρ c (F b t ) + + 1 ρ d (F b t ) -≤ b . (19)
Thus, the constraints on F b t depends on the stock B t . The same reasoning applies for the tank power F h t . Furthermore, we set bound constraints on controls, that is,

-f b ≤ F b t ≤ f b , 0 ≤ F h t ≤ f h t , 0 ≤ F t t ≤ f t t . ( 20
) Finally the load-balance equation ( 13) also acts as a constraint on the controls. We gather all these constraints in an admissible set on control U t depending on the current state X t :

U t ∈ U ad t (X t ) P -a.s. . (21)

G. Objective

At time t, the operational cost L t : X t × U t × W t+1 → R aggregates two different costs:

L t (x t , u t , w t+1 ) = p e t × f ne t+1 + p d t × max(0, θ i t -θ i t ) . (22)
First, we pay a price p e t to import electricity from the network between time t and t + 1. Hence, electricity cost is equal to p e t × F ne t+1 . Second, if the indoor temperature is below a given threshold, we penalize the induced discomfort with a cost p d t × max(0, θ i t -θ i t ), where p d t is a virtual price of discomfort. The cost L t is a convex piecewise linear function, which will prove important for the SDDP algorithm.

We add a final cost K : X T → R to ensure that stocks are non empty at final time

T K(x T ) = κ × max(0, x 0 -x T ) , ( 23
)
where κ is a positive penalization coefficient calibrated by trials and errors.

As decisions U t and states X t are random, the costs L t (X t , U t , W t+1 ) become also random variables. We choose to minimize the expected value of the daily operational cost, yielding the criterion

E T -1 t=0 L t (X t , U t , W t+1 ) + K(X T ) , (24)
yielding an expected value of a convex piecewise linear cost.

H. Stochastic optimal control formulation

Finally, the EMS problem translates to a generic Stochastic Optimal Control (SOC) problem

min X,U E T -1 t=0 L t (X t , U t , W t+1 ) + K(X T ) , (25a)
X 0 = x 0 , (25b) X t+1 = f t X t , U t , W t+1 P -a.s. , (25c)
U t ∈ U ad t (X t ) P -a.s. , (25d) σ(U t ) ⊂ F t . (25e)
Problem [START_REF] Widén | A high-resolution stochastic model of domestic activity patterns and electricity demand[END_REF] states that we want to minimize the expected value of the costs while satisfying the dynamics, the control bounds and the non-anticipativity constraints.

IV. RESOLUTION METHODS

The exact resolution of Problem ( 25) is out of reach in general. We propose two different algorithms that provide policies π t :

X 0 × W 1 × • • • × W t → U t
that map available information x 0 , w 1 , . . . , w t at time t to a decision u t .

A. Model Predictive Control (MPC)

MPC is a classical algorithm commonly used to handle uncertainties in energy systems. At time t, it considers a deterministic forecast (w t+1 , . . . , w T ) of the future uncertainties (W t+1 , . . . , W T ) and solves the deterministic problem

min (ut,••• ,u T -1 ) T -1 j=t L j (x j , u j , w j+1 ) + K(x T ) , (26a)
x j+1 = f j x j , u j , w j+1 , (26b)

u j ∈ U ad j (x j ) . (26c)
At time t, we solve Problem [START_REF] Lloyd | Least squares quantization in PCM[END_REF], retrieve the optimal decisions (u t , . . . , u T -1 ) and only keep the first decision u t to control the system between time t and t + 1. Then, we restart the procedure at time t + 1.

As Problem ( 26) is linear and the number of time steps remains not too large, we are able to solve it exactly for every t.

B. Stochastic Dual Dynamic Programming (SDDP)

1) Dynamic Programming and Bellman principle: the Dynamic Programming method [START_REF] Bellman | Dynamic Programming[END_REF] looks for solutions of Problem (25) as state-feedbacks π t : X t → U t . Dynamic Programming computes a serie of value functions backward in time by setting V T (x T ) = K(x T ) and solving the recursive equations

V t (x t ) = min u∈U ad t (xt) Wt+1 L t (x t , u, w t+1 )+ V t+1 f t (x t , u, w t+1 ) µ of t+1 (dw t+1 ) , (27)
where µ of t+1 is a (offline) probability distribution on W t+1 . Once these functions are computed, we compute a decision at time t as a state-feedback:

π t (x t ) ∈ arg min u∈U ad t (xt) Wt+1 L t (x t , u, w t+1 )+ V t+1 f t (x t , u, w t+1 ) µ on t+1 (dw t+1 ) , (28)
where µ on t+1 is an online probability distribution on W t+1 . This method proves to be optimal when the uncertainties W 1 , . . . , W T are stagewise independent and when µ on t = µ of t is the probability distribution of W t in (27).

2) Description of Stochastic Dual Dynamic Programming: Dynamic Programming suffers from the well-known curse of dimensionality [START_REF] Bertsekas | Dynamic programming and optimal control[END_REF]: its numerical resolution fails for state dimension typically greater than 4 when value functions are computed on a discretized grid. As the state X t in §III-D has 4 dimensions, SDP would be too slow to solve numerically Problem [START_REF] Widén | A high-resolution stochastic model of domestic activity patterns and electricity demand[END_REF]. The Stochastic Dual Dynamic Programming (SDDP) can bypass the curse of dimensionality by approximating value functions by polyhedral functions. It is optimal for solving Problem [START_REF] Widén | A high-resolution stochastic model of domestic activity patterns and electricity demand[END_REF] when uncertainties are stagewise independent, costs L t and K are convex and dynamics f t are linear [START_REF] Girardeau | On the convergence of decomposition methods for multistage stochastic convex programs[END_REF].

SDDP provides an outer approximation V k t of the value function V t in (27) with a set of supporting hyperplanes

{(λ j t , β j t )} j=1,••• ,k by V t (x t ) = min θt∈R θ t , (29a)
λ j t , x t + β j t ≤ θ t , ∀j = 1, • • • , k . (29b)
Each iteration k of SDDP encompasses two passes.

• During the forward pass, we draw a scenario x 0 , . . . , w k T of uncertainties, and compute a state trajectory x k t t=0•••T along this scenario. Starting from position x 0 , we compute x k t+1 in an iterative fashion: i) we compute the optimal control at time t using the available V k t+1 function

u k t ∈ arg min u∈U ad t (xt) Wt+1 L t (x k t , u, w t+1 )+ V k t+1 f t (x k t ,
: V k+1 t = max{V k t , λ k+1 t , . + β k+1 t }.
Otherwise stated, SDDP only explores the state space around "interesting" trajectories (those computed during the forward passes) and refines the value functions only in the corresponding space regions (backward passes).

3) Obtaining online controls with SDDP: in order to compute implementable decisions, we use the following procedure.

• Approximated value functions V t are computed with the SDDP algorithm (see §IV-B2). These computations are done offline. • The approximated value functions V t are then used to compute online a decision at any time t for any state x t . More precisely, we compute the SDDP policy π sddp t by

π sddp t (x t ) ∈ arg min u∈U ad t (xt) Wt+1 L t (x t , u, w t+1 )+ V t+1 f t (x t , u, w t+1 ) µ on t+1 (dw t+1 ) , (
32) which corresponds to replacing the value function V t+1 in Equation ( 28) with its approximation V t+1 . The decision π sddp t (x t ) is used to control the system between time t and t + 1. Then, we resolve Problem (32) at time t + 1.

To solve numerically problems (30)-( 31)-(32) at time t, we will consider distributions with finite support w

p s L t (x t , u, w s t+1 )+ V t+1 f t (x t , u, w s t+1 ) . ( 33
)
V. NUMERICAL RESULTS A. Case study 1) Settings: we aim to solve the stochastic optimization problem [START_REF] Widén | A high-resolution stochastic model of domestic activity patterns and electricity demand[END_REF] over one day, with 96 time steps. The battery's size is 3 kWh, and the hot water tank has a capacity of 120 l. We suppose that the house has a surface A p = 20 m 2 of solar panel at disposal, oriented south, and with a yield of 15%. We penalize the recourse variable F ne t+1 in [START_REF] Baetens | Modelling uncertainty in district energy simulations by stochastic residential occupant behaviour[END_REF] with on-peak and off-peak tariff, corresponding to Électricité de France's (EDF) individual tariffs. The building's thermal envelope corresponds to the French RT2012 specifications [START_REF] Officiel | Arrêté du 28 décembre 2012 relatif aux caractéristiques thermiques et aux exigences de performance énergétique des bâtiments nouveaux[END_REF]. Meteorological data comes from Meteo France measurements corresponding to the year 2015.

2) Demands scenarios: we have scenarios of electrical and domestic hot water demands at 15 minutes time steps, obtained with StRoBe [START_REF] Baetens | Modelling uncertainty in district energy simulations by stochastic residential occupant behaviour[END_REF]. Figure 2 displays 100 scenarios of electrical and hot water demands. We observe that the shape of these scenarios is consistent: demands are almost null during night, and there are peaks around midday and 8 pm. Peaks in hot water demands corresponds to showers. We aggregate the production of the solar panel Φ pv and the electrical demands D el in a single variable D el to consider only two uncertainties (D el t , D hw t ). 3) Out of sample assessment of strategies: to obtain a fair comparison between SDDP and MPC, we use an out-ofsample validation. We generate 2,000 scenarios of electrical and hot water demands, and we split these scenarios in two distinct parts: the first N opt scenarios are called optimization scenarios, and the remaining N sim scenarios are called assessment scenarios. We made the choise N opt = N sim = 1, 000.

First, during the offline phase, we use the optimization scenarios to build models for the uncertainties, under the mathematical form required by each algorithm (see Sect. IV). Second, during the online phase, we use the assessment scenarios to compare the strategies produced by these algorithms. At time t during the assessment, the algorithms cannot use the future values of the assessment scenarios, but can take advantage of the observed values up to t to update their statistical models of future uncertainties.

B. Numerical implementation

1) Implementing the algorithms: we implement MPC and SDDP in Julia 0.6, using JuMP [START_REF] Dunning | JuMP: A modeling language for mathematical optimization[END_REF] as a modeler, StochDynamicProgramming.jl as a SDDP solver, and Gurobi 7.02 [START_REF] Inc | Gurobi Optimizer Reference Manual[END_REF] as a LP solver. All computations run on a Core i7 2.5 GHz processor, with 16Go RAM.

2) MPC procedure: Electrical and thermal demands are naturally correlated in time [START_REF] Widén | A high-resolution stochastic model of domestic activity patterns and electricity demand[END_REF]. To take into account such a dependence across the different time-steps, we chose to model the process W 1 , . . . , W T with an auto-regressive (AR) process.

a) Building offline an AR model for MPC: we fit an AR(1) model upon the optimization scenarios (we do not consider higher order lag for the sake of simplicity). For i ∈ {el, hw}, the AR model writes

d i t+1 = α i t d i t + β i t + ε i t , (34a)
where the non-stationary coefficients (α i t , β i t ) are, for all time t, solutions of the least-square problem

(α i t , β i t ) = arg min a,b Nopt s=1 d i,s t+1 -ad i,s t -b 2 2 . ( 34b
)
The points

(d i,1 t , . . . d i,Nopt t
) correspond to the optimization scenarios. The AR residuals (ε el t , ε hw t ) are a white noise process.

b) Updating the forecast online: once the AR model is calibrated, we use it to update the forecast during assessment (see §IV-A). The update procedure is threefold: i) we observe the demands w t = (d el t , d hw t ) between time t -1 and t, ii) we update the forecast w t+1 at time t + 1 with the AR model

w t+1 = d el t+1 , d hw t+1 = α el t d el t +β el t , α hw t d hw t +β hw t ,
iii) we set the forecast between time t + 2 and T by using the mean values of the optimization scenarios:

w τ = 1 N opt N opt i=1 w i τ ∀τ = t + 2, • • • , T .
Once the forecast (w t+1 , . . . , w T ) is available, it is fed into the MPC algorithm that solves Problem [START_REF] Lloyd | Least squares quantization in PCM[END_REF].

3) SDDP procedure: even if electrical and thermal demands are naturally correlated in time [START_REF] Widén | A high-resolution stochastic model of domestic activity patterns and electricity demand[END_REF], the SDDP algorithm only relies upon marginal distributions.

a) Building offline probability distributions for SDDP: rather than fitting an AR model like done with MPC, we use the optimization scenarios to build marginal probability distributions µ of t that will feed the SDDP procedure in (30)-(31).

We cannot directly consider the discrete empirical marginal probability derived from all N opt scenarios, because the support size would be too large for SDDP. This is why we use optimal quantization to map the N opt optimization scenarios to S representative points. We use a Loyd-Max quantization scheme [START_REF] Lloyd | Least squares quantization in PCM[END_REF] to obtain a discrete probability distribution: at each time t, we use the N opt optimization scenarios to build a partition Ξ = (Ξ b) Computing value functions offline: then, we use these probability distributions as an input to compute a set of value functions with the procedure described in §IV-B2.

c) Using the value functions online: once the value functions have been computed by SDDP, we are able to compute online decisions with Equation (33) 1 . SDDP, on the contrary of MPC, does not update the online probability distribution µ on t during assessment to consider the information brought by the previous observations.

4) Heuristic procedure: we choose to compare the MPC and SDDP algorithms with a basic decision rule. This heuristic is as follows: the battery is charged whenever the solar production Φ pv is available, and discharged to fulfill the demand if there remains enough energy in the battery; the tank is charged (F h t > 0) if the tank energy H t is lower than H 0 , the heater F t t is switched on when the temperature is below the setpoint θ i t and switched off whenever the temperature is above the setpoint plus a given margin.

C. Results

1) Assessing on different meteorological conditions: we assess the algorithms on three different days, with different meteorological conditions (see Table I). Therefore, we use three distinct sets of N sim assessment scenarios of demands, one for each typical day.

Date

Temp. ( The production of the solar panel varies accordingly.

2) Comparing the algorithms performances: during assessment, we use MPC (see [START_REF] Lloyd | Least squares quantization in PCM[END_REF]) and SDDP (see (32)) strategies to compute online decisions along N sim assessment scenarios. Then, we compare the average electricity bill obtained with these two strategies and with the heuristic. The assessment results are given in Table II: means and standard deviation σ are computed by Monte Carlo with the N sim assessment scenarios; the notation ± corresponds to the interval ±1.96 σ √ N sim , which is a 95% confidence interval. We observe that MPC and SDDP exhibit close performance, and make better than the heuristic. If we consider mean electricity bills, SDDP achieves better savings than MPC In addition, SDDP achieves better savings than MPC for the vast majority of scenarios. Indeed, if we compare the difference between the electricity bills scenario by scenario, we observe that SDDP is better than MPC for about 93% of the scenarios. This can be seen on Figure 3 that displays the histogram of the absolute gap savings between SDDP and MPC during Summer day. The distribution of the gap exhibits a heavy tail that favors SDDP on extreme scenarios. Similar analyses hold for Winter and Spring day. Thus, we claim that SDDP outperforms MPC for the electricity savings. Concerning the performance on thermal comfort, temperature trajectories are above the temperature setpoints specified in §III-G for both MPC and SDDP.

In term of numerical performance, it takes less than a minute to compute a set of cuts as in §IV-B2 with SDDP on a particular day. Then, the online computation of a single decision takes 1.5 ms on average, compared to 0.5 ms for MPC. Indeed, MPC is favored by the linearity of the optimization Problem [START_REF] Lloyd | Least squares quantization in PCM[END_REF], whereas, for SDDP, the higher the quantization size S, the slowest is the resolution of Problem (32), but the more information the online probability distribution µ on t carries. 3) Analyzing the trajectories: we analyze now the trajectories of stocks in assessment, during Summer day. The heating is off, and the production of the solar panel is nominal at midday.

Figure 4 displays the state of charge of the battery along a subset of assessment scenarios, for SDDP and MPC. We observe that SDDP charges earlier the battery at its maximum. On the contrary MPC charges the battery later, and does not use the full potential of the battery. The two algorithms discharge the battery to fulfill the evening demands. We notice that each trajectory exhibits a single cycle of charge/discharge, thus decreasing battery's aging.

Figure 5 displays the charge of the domestic hot water tank along the same subset of assessment scenarios. We observe a similar behavior as for the battery trajectories: SDDP uses more the electrical hot water tank to store the excess of PV energy, and the level of the tank is greater at the end of the day than in MPC. This analysis suggests that SDDP makes a better use of storage capacities than MPC.

VI. CONCLUSION

We have presented a domestic microgrid energy system, and compared different optimization algorithms to control the stocks with an Energy Management System.

The results show that optimization based strategies outperform the proposed heuristic procedure in term of money savings. Furthermore, SDDP outperforms MPC during Winter and Summer day -achieving up to 35% costs savings -and displays similar performance as MPC during Spring day. Even if SDDP and MPC exhibit close performance, a comparison scenario by scenario shows that SDDP beats MPC most of the time (more than 90% of scenarios during Summer day). Thus, we claim that SDDP is better than MPC to manage uncertainties in such a microgrid, although MPC gives also good performance. SDDP also makes a better use of storage capacities.

Our study can be extended in different directions. First, we could mix SDDP and MPC to recover the benefits of these two algorithms. Indeed, SDDP is designed to handles the uncertainties' variability but fails to capture the time correlation, whereas MPC ignores the uncertainties' variability, but considers time correlation by means of a multistage optimization problem. Second, we are currently investigating the optimization of larger scale microgrids -with different interconnected buildings -by decomposition methods.

Figure 1 .

 1 Figure 1. Electrical microgrid

Figure 2 .

 2 Figure 2. Electrical (left) and domestic hot water (right) demand scenarios.

Figure 3 .

 3 Figure 3. Absolute gap savings between MPC and SDDP during Summer day

Figure 4 .Figure 5 .

 45 Figure 4. Battery charge trajectories for SDDP and MPC during Summer day

  , • • • , Ξ S ), where Ξ is the solution of the optimal quantization problem . Then, we set for all time t = 0, • • • , T the discrete offline distributions µ of t = S s=1 p s δ w s t , where δ w s t is the Dirac measure at point w s t and p s = card(Ξ s )/N opt is the associated probability weight. We have chosen S = 20 to have enough precision.

	where w s t = of Ξ	1 card(Ξs)	w i t ∈Ξs w i t is the so-called centroid
				S
			min Ξ	s=1	t ∈Ξs w i	w i t -w s t	2 2	(35)

s

  Winter day, the heating is maximal, whereas it is medium during Spring day and null during Summer day.

			• C) PV Production (kWh)
	Winter Day	February, 19th	3.3	8.4
	Spring Day	April, 1st	10.1	14.8
	Summer Day	May, 31st	14.1	23.3
			Table I
	DIFFERENT METEOROLOGICAL CONDITIONS
	These three different days corresponds to different heating
	needs. During

  Summer day and Winter day, but SDDP and MPC display similar performances during Spring day.

		SDDP	MPC	Heuristic
	Offline time	50 s	0 s	0 s
	Online time	1.5 ms	0.5 ms	0.005 ms
	Electricity bill (e)
	Winter day	4.38 ± 0.02 4.59 ± 0.02 5.55 ± 0.02
	Spring day	1.46 ± 0.01 1.45 ± 0.01 2.83 ± 0.01
	Summer day	0.10 ± 0.01 0.18 ± 0.01 0.33 ± 0.02
		Table II
	COMPARISON OF MPC, SDDP AND HEURISTIC STRATEGIES
	during

1

In practice, the quantization size of µ on t is bigger than those of µ of t , to have a greater accuracy online


Jesús Angulo
email: jesus.angulo@mines-paristech.fr

Santiago Velasco-Forero
email: santiago.velasco@mines-paristech.fr

Non-Negative Sparse Mathematical Morphology

Keywords: mathematical morphology, sparse modeling, non-negative matrix factorization, sparse NMF, non-linear filtering, hyperspectral image processing

Sparse modeling involves constructing a succinct representation of initial data as a linear combination of a few typical atoms of a dictionary. This paper deals with the use of sparse representations to introduce new nonlinear image filters which efficiently approximate morphological operators. Reasons why non-negative matrix factorization (NMF) is a dimensional reduction (i.e., dictionary learning) paradigm particularly adapted to the nature of morphological processing are given. In particular, Sparse-NMF representations are studied and used to introduce first approximations to binary dilations/erosions and then to openings/closings. The idea behind consists of processing exclusively the image dictionary and then, the result of processing each image is approximated by multiplying the processed dictionary by the coefficient weights of the current image. These operators are then extended to gray-scale images and their interest for feature detection is illustrated. The practical relevance of our approach is considered for two applications on multivariate image processing. The first case deals with multispectral texture modeling using Boolean random set theory; the second case with multi-scale decomposition of hyperspectral images and its interest in spectral-spatial pixel classification.

Introduction

Modern image processing techniques should deal with images datasets of huge volume, great variety (e.g., optical, radar and hyperspectral sensors, sequence of images often multi-scale, etc.) and different veracity (i.e., incomplete images, or different uncertainty in the acquisition procedure). This paper addresses the problem of simultaneous treatment of a set of images via mathematical morphology. Our approach is motivated by sound developments in sparse signal representation (or coding), which suggest that the linear relationships among 1 high-resolution signals can be accurately recovered from their low-dimensional projections (also called dictionary) [START_REF] Donoho | Compressed sensing[END_REF]. Sparse coding and dictionary learning, where data is assumed to be well represented as a linear combination of a few elements from a dictionary, is an active research topic which leads to state-of-the-art results in image processing applications, such as image denoising, inpainting, demosaicking or compression [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF][START_REF] Mairal | Sparse representation for color image restoration[END_REF]Yu et al., 2011). For a detailed self-contained overview on sparse modeling for image and vision processing the reader is referred to [START_REF] Mairal | Sparse Modeling for Image and Vision Processing[END_REF]. Inspired by this paradigm of parsimony representation of images, the aim of this work is to explore how image sparse representations can be useful to efficiently calculate approximations to morphological operators.

Mathematical morphology [START_REF] Serra | Image Analysis and Mathematical Morphology. Vol I, and Image Analysis and Mathematical Morphology[END_REF][START_REF] Heijmans | Morphological Image Operators[END_REF] is a nonlinear image processing methodology based on the application of lattice theory to spatial structures. Morphological filters and transformations are useful for various image processing tasks [START_REF] Soille | Morphological Image Analysis[END_REF], such as denoising, contrast enhancement, multi-scale decomposition, feature extraction and object segmentation. In addition, morphological operators are defined using very intuitive geometrical notions which allows us the perceptual development and interpretation of complex algorithms by combination of various operators.

Notation. Let E be a domain of points, which is considered here as a finite digital space of the pixels of the image, i.e., E ⊂ Z 2 such that N = |E| is the number of pixels. Image intensities are numerical values, which ranges in a closed subset T of R = R ∪ {-∞, +∞}; for example, for an image of discrete L values, it can be assumed T = {t 1 , t 2 , • • • , t L }. A binary image X is modeled as a subset of E, i.e., X ∈ P(E); a grey-level image f (p i ), where p i ∈ E are the pixel coordinates, is a numerical function E → T , i.e., f ∈ F(E, T ). In this paper, we are interested in operators ψ as a map transforming an image into an image. There are thus operators on binary images, i.e., maps P(E) → P(E); or on grey-level images, i.e., maps F(E, T ) → F(E, T ).

Motivation and outline of the approach. In order to illustrate our motivation, let us consider the diagrams depicted in Fig. 1. In many practical situations, a collection F = {f 1 (p i ), • • • , f M (p i )} of M binary or grey-level images (each image having N pixels) should be analyzed by applying the same morphological operator ψ (or a series of operators) to each image of F, i.e., ψ (F) = {ψ (f j )} 1≤j≤M , Fig. 1(a). If one considers that the content of the various images is relatively similar, we can expect that the initial collection can be efficiently projected into a R dimensionality reduced image space, using typically sparse modeling techniques, i.e., F ≈ ΦH, in order to learn the dictionary composed of atom images Φ = {ϕ 1 (p i ), • • • , ϕ R (p i )} and the corresponding coding or weights H, Fig. 1(b). In this context, we would like to apply the morphological operator ψ (or an equivalent operator) to the reduced set of images of the projective space Φ, i.e., ψ (Φ) = {ψ (ϕ k )} 1≤k≤R , in such a way that the original processed images are approximately obtained by projecting back to the initial space, i.e., ψ (F) ≈ ψ (Φ) H, Fig. 1(c). Note that if the number of operators to be applied increases, e.g., Φ A , Φ B , Φ C , etc., then such paradigm becomes much more efficient since R ≪ M , Fig. 1(d). At this point, we should remark that our approach is different from the standard use of sparse modeling in image processing [START_REF] Mairal | Sparse Modeling for Image and Vision Processing[END_REF], since in the classical approaches the dictionary is not processed. In fact, sparse representation itself leads to a solution for the tasks of denoising, regularization, coding, etc. Typical examples of image families which can be fit in our framework are: i) a collection of binary shapes, ii) a database of registered grey scale images, for instance, faces, cells (from biomicroscopy applications) or galaxies (in astronomy), iii) a set of patches of a large image, iv) a time series of registered images, v) the spectral bands of a multi/hyper-spectral image, etc. Some of these examples are used to illustrate the relevance of our approach. The rationale behind this methodology is that the intrinsic dimension1 of the image collection is lower than M . Usually the subspace representation involves deriving a set of basis components using linear techniques like Principal Component Analysis (PCA), Independent Component Analysis (ICA), or sparse dictionary learning techniques like K-SVD [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF]. The projection coefficients for the linear combinations in the above methods can be either positive or negative, and such linear combinations generally involve complex cancellations between positive and negative numbers. Therefore, these representations lack the intuitive meaning of "adding parts to form a whole". This property is particularly problematic in the case of mathematical morphology since the basic binary operator, the dilation of a set, is defined as the operator which commutes with the union of parts of the set. This principle is widely discuss in the paper. In practice, that means that the paradigm discussed above in Fig. 1 cannot be used with sign unconstrained sparse modeling methods. Non-negative matrix factorization (NMF) [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF] imposes the non-negativity constraints in learning basis images: the pixels values of resulting images as well as the coefficients for the reconstruction are all non-negative. This ensures that NMF is a procedure for learning a parts-based representation [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF]. In addition, by incorporating sparse constraint to NMF, it will be constructed a succinct representation of the image data as a combination of a few typical patterns (few atoms of the dictionary) learned from the data itself. Hence, our approach of sparse mathematical morphology is founded on sparse NMF modeling.

NMF state-of-the-art. The state-of-the on NMF is nowadays vast. Besides the pioneering work [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF], several NMF variants have been developed by introducing additional constraints and properties to the original NMF. In particular, being relevant for our purpose, we can mention, on the one hand, the local NMF [START_REF] Li | Learning spatially localized parts-based representations[END_REF], which learns spatially localized, parts-based subspaces for images (i.e., visual patterns); and on the other hand, the sparse NMF [START_REF] Hoyer | Non-negative Matrix Factorization with Sparseness Constraints[END_REF], which incorporates explicit sparseness constraints. These powerful and already classical algorithms are reviewed in next section. We note also that in hyperspectral imaging, there is a common task called linear unmixing (or factorization into physical space) which consists in detecting the spectra of the pure materials and estimating their relative abundances. Linear unmixing can be seen as an equivalent problem to NMF [START_REF] Esser | A convex model for nonnegative matrix factorization and dimensionality reduction on physical space[END_REF].

Formulation of PCA-like optimization with nonnegative and sparse constraints have also considered by several works, see for instance [START_REF] Zass | Nonnegative sparse PCA[END_REF]. Motivated by a combination of NMF and kernel theory (emmbedding into polynomial feature space), various nonlinear variants of NMF have been proposed [START_REF] Zafeiriou | Nonlinear non-negative component analysis algorithms[END_REF]. Similarly, algorithms have been also introduced for nonnegative ICA [START_REF] Plumbley | Algorithms for nonnegative independent component analysis[END_REF]. NMF has been extended to factorization of non vector spaces, typically for data sampled from a low-dimensional manifold embedded in a high-dimensional space, using affinity graph matrix non-negative decomposition [START_REF] Cai | Graph regularized nonnegative matrix factorization for data representation[END_REF] and manifold regularized NMF [START_REF] Guan | Manifold regularized discriminative nonnegative matrix factorization with fast gradient descent[END_REF]. NMF has been also extended to tensor factorizations [START_REF] Cichocki | Nonnegative matrix and tensor factorizations[END_REF], and has considered for applications to exploratory multi-way data analysis and blind source separation. In general, NMF is a NP-hard problem [START_REF] Vavasis | On the complexity of nonnegative matrix factorization[END_REF]. However, starting from the notion of separable NMF, introduced by [START_REF] Donoho | When does non-negative matrix factorization give a correct decomposition into parts?[END_REF], which in geometric terms involves that the conical hull of a small subsets of the columns of the data matrix contain the rest of the columns, an emerging line of work in NMF, based on computational geometry approaches, is leading to efficient solutions. It can be then solved by linear programming [START_REF] Arora | Computing a nonnegative matrix factorizationprovably[END_REF], conical hull algorithms [START_REF] Kumar | Fast conical hull algorithms for near-separable nonnegative matrix factorization[END_REF] and random projections [START_REF] Damle | Random projections for non-negative matrix factorization[END_REF]. When the data to be factorized is binary or Boolean (i.e., only "0" or "1" values), the NMF problem is called discrete basis problem [START_REF] Miettinen | The discrete basis problem[END_REF]. This discrete basis problem which is equivalent to the binary k-median problem is NP-hard and some greedy algorithms [START_REF] Miettinen | Sparse Boolean matrix factorizations[END_REF] have been proposed to deal with. Other binary matrix factorization approaches are straightforward related to classical NMF [START_REF] Zhang | Binary matrix factorization with applications[END_REF] either by a thresholding algorithm (smooth approximation to Heaviside function) or by binary-constraint penalty function. A binary version of the sparse NMF can be also mentioned [START_REF] Yuan | Binary Sparse Nonnegative Matrix Factorization[END_REF], based on windowed image parts with Haar-like box functions. Extension of NMF to binary data factorization has been also considered by means of a Bernoulli distribution of binary observations [START_REF] Schachtner | A nonnegative blind source separation model for binary test data[END_REF], or more recently using a logistic PCA approach [START_REF] Tomé | A logistic non-negative matrix factorization approach to binary data sets[END_REF].

Paper organisation. The methodology described here was initiated in [START_REF] Angulo | Sparse Mathematical Morphology using Non-Negative Matrix Factorization[END_REF] for binary erosion/dilation. In the present work, we develop also sparse approximations to binary opening/closing as well as the full counterpart of sparse approximation to morphological operators for gray-scale images. The applications considered in the last part of the paper are also new with respect to our previous work [START_REF] Angulo | Sparse Mathematical Morphology using Non-Negative Matrix Factorization[END_REF].

The rest of paper is structured as follows. Section 2 reviews the formal definition of NMF and various algorithms proposed in the state-of-the-art, including the sparse variant introduced in [START_REF] Hoyer | Non-negative Matrix Factorization with Sparseness Constraints[END_REF], which is the one used in the our study. The use of NMF representations for implementing sparse pseudo-morphological binary operators is introduced and theoretically justified in Section 3. The section discusses also why sparse NMF is appropriate in morphological image processing. Section 4 generalizes the framework of morphological sparse processing to gray-scale images using two alternatives. The first one is based on level set decomposition of numerical functions. The second approach, more useful in practice, uses a straightforward representation of sparse NMF basis from the gray-scale images. Applications are studied in Section 5. The first case-study deals with the sparse processing and modeling of a multivariate Boolean texture. The second application focusses on morphological sparse processing of hyperspectral images. Conclusions and perspectives close the paper in Section 6.

NMF and Sparse Variants

Definition of NMF on vector space

Let us assume that our data consists of M vectors of N non-negative scalar variables. Denoting the column vector v j , j = 1, • • • , M , the matrix of data is obtained as

V = (v 1 , • • • , v M )
(each v j is the j-th column of V), with |v j | = N . If we analyze M images of N pixels, these images can be stored in linearized form, so that each image will be a column vector of the matrix.

Given the non-negative matrix V ∈ R N ×M , V i,j ≥ 0, NMF is defined as a linear nonnegative approximate data decomposition into the two matrices

W ∈ R N ×R and H ∈ R R×M such that V ≈ WH, s.t. W i,k , H k,j ≥ 0, (1)
where usually R ≪ M , which therefore involves a dimensionality reduction. Each of the R columns of W contains a basis vector w k and each row of H contains the coefficient vector (weights) h j corresponding to vector v j :

v j = ∑ R k=1 w k H k,j = Wh j .
Using the modern terminology from sparse coding theory, the matrix W contains the dictionary and H the encoding.

A theoretical study of the properties of NMF representation has been achieved in [START_REF] Donoho | When does non-negative matrix factorization give a correct decomposition into parts?[END_REF]) using geometric notions. Hence, NMF is interpreted as the problem of finding a simplicial cone which contains the data points in the positive orthant, or in other words, NMF is a conical coordinate transformation. It is interesting to note that NMF is somehow equivalent to Kernel k-means clustering and Laplacian-based spectral clustering [START_REF] Ding | On the equivalence of nonnegative matrix factorization and spectral clustering[END_REF].

Basic NMF algorithms

The factorization V ≈ WH is not necessarily unique, and the optimal choice of matrices W and H depends on the cost function that minimizes the reconstruction error. The most widely used is the Euclidean distance, i.e., minimize

∥V -WH∥ 2 2 = ∑ i,j (V i,j -(WH) i,j ) 2 ,
with respect to W and H, and subject to the constraints W, H > 0. Although the minimization problem is convex in W and H separately, it is not convex in both simultaneously. In [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF] a multiplicative good performance algorithm to solve (1) is proposed. They proved that the cost function is nonincreasing at the iteration and the algorithm converges at least to a local optimal solution. More precisely, the update rules for both matrices are:

H k,j ← H k,j (W T V) k,j (W T WH) k,j ; W i,k ← W i,k (VH T ) i,k (WHH T ) k,j .
Another useful cost function, also considered in [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF], is the Kullback-Leibler (KL) divergence, which leads also quite simple multiplicative update rules.

In [START_REF] Li | Learning spatially localized parts-based representations[END_REF], a variant of KL divergence NMF was proposed, which is named Local NMF (LNMF), aiming at learning spatially localized components (by minimizing the number of basis R to represent V and by maximizing the energy of each retained components) as well as imposing that different bases should be as orthogonal as possible (in order to minimize redundancy between the different bases). The multiplicative update rules for LNMF are given by

H k,j ← √ H k,j ∑ i V i,j W i,k (WH) i,k ; W i,k ← W i,k ∑ j V i,j H k,j (WH) i,j ∑ j H k,j ; W i,k ← W i,k ∑ i W i,k
.

As we have discussed in the introduction, there are many other alternative NMF algorithms, but we focus on the following technique which provide a sparse variant.

NMF with sparseness constraints

A very powerful framework to add an explicit degree of sparseness in the basis vectors W and/or the coefficients H was introduced in [START_REF] Hoyer | Non-negative Matrix Factorization with Sparseness Constraints[END_REF]. First of all, the sparseness measure σ of a vector v ∈ R N ×1 used in [START_REF] Hoyer | Non-negative Matrix Factorization with Sparseness Constraints[END_REF]) is based on the relationship between the L 1 norm and the L 2 norm:

σ(v) = √ N -∥v∥ 1 /∥v∥ 2 √ N -1 .
This function is maximal at one if and only if v contains only a single non-zero component, and takes a value of zeros if and only if all components are equal (up to signs). Then, matrix W and H are solved by the problem (1) under additional constraints σ(w k ) = S w and σ(h j ) = S h , where S w and S h are respectively the desired sparseness of W and H. The algorithm introduced in [START_REF] Hoyer | Non-negative Matrix Factorization with Sparseness Constraints[END_REF]) is a projected gradient descent algorithm (additive update rule), which takes a step in the direction of the negative gradient

W i,k ← W i,k -µ W (W i,k H k,j -V i,j ) H T j,k ; H k,j ← H k,j -µ H W T k,i (W i,k H k,j -V i,j ) ;
and subsequently projects each column of W and each row of H onto the constraint space.

The most sophisticated step is therefore how to find, for a given vector v, the closest nonnegative vector u with a given L 1 norm and a given L 2 norm. The algorithm works as follows [START_REF] Hoyer | Non-negative Matrix Factorization with Sparseness Constraints[END_REF]. The vector is projected onto the hyperplane of L 1 . Next, within this space, one projects to the closest point on the joint constraint hypersphere (intersection of the sum and the L 2 constraints), by moving radially outward for the center of the sphere (the center is the point where all the components have equal values). If the components of this projection point are not completely non-negative, the negative values must be fixed at zero and a new point found in similar way under the same constraints. Sparseness is controlled explicitly with a pair of parameters that is easily interpreted; in addition, the number of required iterations grows very slowly with the dimensionality of the problem. In fact, for all the empirical tests considered in this paper, we have used the MATLAB code for performing NMF and its various extensions (LNMF, sparse NMF) provided by [START_REF] Hoyer | Non-negative Matrix Factorization with Sparseness Constraints[END_REF].

Besides the sparseness parameters (S w , S h ), a crucial parameter to be chosen in any NMF algorithm is the value of R, that is, the number of basis of projective reduced space.

Any dimensionality reduction technique, such as PCA, requires also to fix the number of components. In PCA, the components are ranked according to the second-order statistical importance of the components and each one has associated a value of the represented variance; whereas in NMF the selection R can be evaluated only a posteriori, by evaluating the error of reconstruction. This issue is not considered in the paper and we empirically fix R.

A few properties of NMF

We conclude this review on NMF algorithms by discussing some well-known propeorties which are useful for the sequel.

• Boundedness property (for any NMF variant) [START_REF] Zhang | Binary matrix factorization with applications[END_REF]

: We say that V is bounded if 0 ≤ V i,j ≤ 1.
If V is bounded, then the factor matrices W and H are also bounded, i.e., 0 ≤ W i,k ≤ 1 and 0 ≤ H k,j ≤ 1.

• Indeterminacies (for any NMF variant) [START_REF] Theis | First results on uniqueness of sparse non-negative matrix factorization[END_REF][START_REF] Huang | Non-negative matrix factorization revisited: Uniqueness and algorithm for symmetric decomposition[END_REF]: Positivity and sparseness are invariant under permutation and scaling of columns of W (and correspondingly of the rows of H), i.e., V = WH = (WP -1 L -1 )(LPH), where P is a permutation matrix and L a positive scaling matrix.

• Uniqueness (only for Sparse-NMF) [START_REF] Theis | First results on uniqueness of sparse non-negative matrix factorization[END_REF][START_REF] Huang | Non-negative matrix factorization revisited: Uniqueness and algorithm for symmetric decomposition[END_REF]: Under sparsity constraints, projection step proposed in [START_REF] Hoyer | Non-negative Matrix Factorization with Sparseness Constraints[END_REF]) has a unique solution, which is found. Hence, under non-degenerate data, Sparse-NMF a.s. produces a unique factorization. This is another fundamental reason why we think that Sparse-NMF is an excellent choice for our learning our non-negative dictionaries.

3 Sparse approximation to binary morphological operators

Let X = {X 1 , • • • , X M }
be a collection of M binary shapes, called a family of shapes, i.e., X j ∈ P(E). For each shape X j , let x j (i) : I → {0, 1}, with i ∈ I = {1, 2, • • • , N } and N = |E|, be its indicator vector :

∀X j ∈ P(E), ∀p i ∈ E, x j (i) = { 1 if p i ∈ X j 0 if p i ∈ X c j (2)
Then the family of shapes X has associated a data matrix V ∈ {0, 1} N ×M , where each indicator vector corresponds to one of its columns, i.e., V i,j = x j (i).

Sparse NMF approximations of binary sets

After computing the NMF representation on data V, for a given constant R > 0, an approximation to V is obtained. More precisely, if we denote by ϕ k (p i ) : E → R + the basis images associated to the basis matrix W, i.e., ϕ k (p i ) = W i,k , the following image is obtained as

a X j (p i ) = R ∑ k=1 ϕ k (p i )H k,j (3)
It is obvious that without any additional constrains, function a X (p i ) is not a binary image. By the boundedness property discussed above, we have 0 ≤ a X (p i ) ≤ 1. Hence, a thresholding operation at value α is required to impose a binary approximate set X j to each initial shape X j , i.e.,

X j N M F ----→ X j : p i ∈ X j if a X (p i ) > α (4)
We propose to fix, for all the examples of the paper, the threshold value to α = 0.45, in order to favor the reconstruction of X j against its complement. Let us consider a practical example from a binary image collection X , composed of M = 100 images from the Fish shape database (N = 400 × 200), see in Fig. 2 a few examples of images. Fig. 3 depicts the corresponding basis images for various NMF algorithms: we have fixed R = 10 for all the cases (factor 10 of dimensionality reduction). We observe that standard NMF produces a partial part-based representation, which also includes almost complete objects for outlier shapes (basis 2-upper-center and 5-center-center). As expected, LNMF produces more local decompositions, however the orthogonality constraints involves also an atomization of some common parts. A similar problem arises for Sparse-NMF when S w ̸ = 0 (constraint of sparsity in basis matrix W). When the sparsity constraint is limited to the coding S h , with a typical value around 0.6, the obtained dictionary of shapes is less local, but in exchange, this constraint involves that each binary shapes is reconstructed using a limited number of atoms. The various groups of fish shapes are therefore better approximated by the latter case than using the other NMF algorithms. The comparison of Fig. 4 illustrates qualitatively the better performance of Sparse-NMF (S w = 0, S h = 0.6) with respect to the others. We have also included in Fig. 3

Sparse max-approximation to binary dilation and erosion

Dilation and Erosion. The two fundamental morphological operators are the dilation and the erosion, which are defined respectively as the operators which commute with the union and the intersection. Given a structuring element B ⊆ E, i.e., a set defined at the origin which introduces the shape/size of the operator, the dilation of a binary image X by B and the erosion of a binary image X by B are defined respectively by [START_REF] Serra | Image Analysis and Mathematical Morphology. Vol I, and Image Analysis and Mathematical Morphology[END_REF][START_REF] Heijmans | Morphological Image Operators[END_REF][START_REF] Bloch | Mathematical Morphology[END_REF]:

δ B (X) = ∪ {B(p i )| p i ∈ X} , (5)
and

ε B (X) = {p i ∈ E| B(p i ) ⊆ X} , ( 6
)
where B(p i ) is the structuring element centered at pixel p i . In the case of numerical functions f ∈ F (E, T ), which are considered in detail in next section, the dilation of a grey-level image is defined by [START_REF] Heijmans | Morphological Image Operators[END_REF][START_REF] Soille | Morphological Image Analysis[END_REF][START_REF] Bloch | Mathematical Morphology[END_REF]:

δ B (f )(p i ) = { f (p m )| f (p m ) = sup [f (p n )] , p n ∈ B(p i ) } = sup pn∈ B(p i ) {f (p n )} , ( 7
)
and the dual grey-level erosion is given by [START_REF] Heijmans | Morphological Image Operators[END_REF][START_REF] Soille | Morphological Image Analysis[END_REF][START_REF] Bloch | Mathematical Morphology[END_REF]) The indicator function of set X, denoted X : E → {0, 1}, is defined by

ε B (f )(p i ) = {f (p m )| f (p m ) = inf [f (p n )] , p n ∈ B(p i )} = inf pn∈B(p i ) {f (p n )} , ( 8
∀X ∈ P(E), ∀p i ∈ E, X (p i ) = { 1 if p i ∈ X 0 if p i ∈ X c (9)
Obviously, we have X c = 1 -X . Given two sets X, Y ∈ P(E), one has the two following basic properties of indicator function: For a function f : E → T , the thresholded set at value t ∈ T is a mapping from F(E, T ) to P(E) given by [START_REF] Serra | Image Analysis and Mathematical Morphology. Vol I, and Image Analysis and Mathematical Morphology[END_REF]:

X∩Y = min{ X , Y } = X • Y ; X∪Y = max{ X , Y } = X + Y -X • Y = min{1, X + Y }.
ϖ t (f ) = {p i ∈ E| f (p i ) ≥ t} . ( 10
)
Using these transformations it is obvious that the binary dilation (6) can be computed using the numerical operator (7), i.e.,

δ B (X) = ϖ 1 (δ B ( X )(p i )) . ( 11
)
By the fundamental property of dilation [START_REF] Serra | Image Analysis and Mathematical Morphology. Vol I, and Image Analysis and Mathematical Morphology[END_REF], given a set defined as the union of a family of sets, i.e., X = ∪ k∈K X k , its corresponding binary dilation by B is

δ B (X) = δ B (∪ k∈K X k ) = ∪ k∈K δ B (X k ). ( 12
)
From ( 11) and the property of the indicator function for the union of sets, it is easy to see that ( 12) can be rewritten using the dilation of functions as

δ B (X) = ϖ 1 ( δ B ( min { 1, ∑ k∈K X k (p i ) })) = ϖ 1 ( min { 1, ∑ k∈K δ B ( X k (p i )) })
, and finally it is obtained that

δ B (X) = ϖ 1 ( ∑ k∈K δ B ( X k (p i ))
)

. ( 13
)
Hence, the justification for using NMF-based part decompositions in sparse mathematical morphology arises from equations ( 12) and ( 13).

Coming back to the NMF reconstruction, expressions (3) and ( 4), we can now write

X j ≈ X j = ϖ α ( R ∑ k=1 ϕ k (p i )H k,j ) . ( 14
)
Hence, we introduce the following nonlinear operator, named sparse max-approximation to binary dilation:

D B (X j ) = ϖ α ( R ∑ k=1 δ B (ϕ k ) (p i )H k,j ) . ( 15
)
Note that by the nonnegativity of H k,j , we have

δ B (ϕ k (p i )H k,j ) = δ B (ϕ k ) (p i )H k,j . We can say that if X j ≈ X j then δ B (X j ) ≈ D B (X j ),
or in other words, the degree of approximation to the dilation depends on the degree of approximation to the original set. However neither the increasiness nor the extensitivity of D B (X j ) w.r.t. X j can be guaranteed and consequently, operator ( 15) is not a morphological dilation in an algebraic sense [START_REF] Heijmans | The algebraic basis of mathematical morphology I. Dilations and erosions[END_REF].

In conclusion, in order to approximate the dilation by B of any of the M sets X j , we only need to calculate the dilation of the R basis images. In addition, if sparsity is imposed to H, that involves that only a limited number of dilated atoms are required for each X j .

Dual sparse max-approximation to binary erosion. One of the most interesting properties of mathematical morphology is the duality by the complement of pairs of operators, and in particular the duality between the dilation and the erosion. Hence, the binary erosion of set X by structuring element B can be defined as the dual operator to the dilation: ε B (X) = (δ B (X c )) c , where the complement set is X c = X = E \ X. Using this property, we propose to define the sparse max-approximation to binary erosion as

E B (X j ) = ϖ α ( R ∑ k=1 [ δ B ( [ϕ k ] ) (p i ) ] H k,j ) , ( 16
)
where the complement basis images are defined by

[ϕ k (p i )] = max(W i,k )-ϕ k (p i )+min(W i,k ).
The results of sparse max-approximations to D B (X j ) and E B (X j ) for three examples of the Fish shapes, compared to the exact binary dilation and erosion, are given in Fig. 5. We have compared in particular the sparse max-approximation for the standard NMF and for the Sparse-NMF.

Sparse max-approximation to binary opening and closing

Composition of dilation and erosion produces two other operators, called opening and closing, which are the fundamental bricks for morphological filtering. More precisely, opening of a binary image X by B and the dual closing of a binary image X by B are defined as [START_REF] Serra | Image Analysis and Mathematical Morphology. Vol I, and Image Analysis and Mathematical Morphology[END_REF][START_REF] Heijmans | Morphological Image Operators[END_REF]:

γ B (X) = δ B (ε B (X)) = ∪ {B(p i )| p i ∈ E and B(p i ) ⊆ X} , ( 17
) φ B (X) = ε B (δ B (X)) = (γ B (X c )) c . ( 18
)
Similarly for the case of numerical functions f ∈ F(E, T ), opening and closing of a grey-level image by B and its dual grey-level closing are respectively given by [START_REF] Serra | Image Analysis and Mathematical Morphology. Vol I, and Image Analysis and Mathematical Morphology[END_REF][START_REF] Heijmans | Morphological Image Operators[END_REF]:

γ B (f )(p i ) = δ B (ε B (f )) (p i ) = sup pm∈ B(pn) inf pn∈B(p i ) {f (p m )} , φ B (f )(p i ) = ε B (δ B (f )) (p i ) = inf pm∈B(pn) sup pn∈ B(p i ) {f (p m )} .
In order to apply the non-negative decomposition underlaying the NMF coding, we first introduce sufficient conditions for a compatibility of opening with set decomposition, which are based on the geometric formulation of binary opening provided by (17

). Let X = ∪ k∈K X k such that • either ∂X i ∩ ∂X j = ∅, ∀i, j ∈ K, • or X i ⊆ X j or X j ⊆ X i , ∀i, j ∈ K then γ B (X) = γ B ( ∪ k∈K X k ) = ∪ k∈K γ B (X k ). ( 19
)
Hence, opening commutes with the union only in the cases where the subsets are either totally disjoint or totally contained. We call these cases as separable subsets.

Under the assumption of separability of the atoms of the dictionary {ϕ k } 1≤k≤R , by combining ( 17) and ( 19) (equivalence only valid on the separable case), i.e.,

γ B (X) = ∪ k∈K ∪ {B(p i )| p i ∈ E and B(p i ) ⊆ X k } = ∪ { B(p i )| p i ∈ E and B(p i ) ⊆ ∪ k∈K X k } ,
and similarly to the dilation and erosion, we introduce the sparse max-approximation to binary opening and to binary closing respectively as

G B (X j ) = ϖ α ( R ∑ k=1 γ B (ϕ k ) (x i )H k,j ) , ( 20
)
F B (X j ) = ϖ α ( R ∑ k=1 [ γ B ( [ϕ k ] ) (p i ) ] H k,j ) . ( 21
)
As discussed below, in Section 4, even in the case when the separability is not always fulfilled, operators G B (X j ) and F B (X j ) are interesting to approximate the effect of the opening γ B (X j ) and the closing φ B (X j ).

Consistency and noise robustness of sparse morphological operators

Let us discuss the two following important properties which have been empirically observed.

From our viewpoint, they prove the pertinence of the sparse max-approximation to binary dilation based on Sparse-NMF representations.

Consistency. Behind this notion of consistency between Sparse-NMF and morphological binary dilation, we mean the fact that the Sparse-NMF basis from an image collection X should be stable to dilation of X , i.e., δ B (X ) = {δ B (X 1 ), • • • , δ B (X M )}. Fig. 6 illustrates the consistency: it is compared in (a) the Sparse-NMF basis from the original image data set X to the Sparse-NMF basis obtained using the same parameters from two dilated image data set δ B (X ), in (b) for B a square of 5 × 5 and in (c) for a square 7 × 7. We have noted that, for a symmetric structuring element B of size relatively small with respect to the size of the objects, i.e.,

|X j \ δ B (X j )|/|X j | ≪ 1, ∀j, the consistency involves that ∑ i ∥ϕ ⊕B k - δ B (ϕ k )∥(p i )/ ∑ i ϕ k (p i ) ≈ 0, ∀k, where X Sparse-N M F ---------→ {{ϕ k } 1≤k≤R ; H} ⇒ δ B (X ) Sparse-N M F ---------→ { {ϕ ⊕B k } 1≤k≤R ; H } .
Using the semi-group property of multi-scale dilation δ nB by homothetic convex structuring elements, i.e., δ nB •δ nB = δ (n+m)B , a consequence of such consistency is depicted in the example of Fig. 6(d): the sparse max-approximation to dilation of size 7 from X is close the sparse max-approximation to dilation of size 5 from a Sparse-NMF representation on δ 2B (X ), i.e., D 5B (δ 2B (X j )) ≈ D 7B (X j ). More generally, given a convex symmetric structuring element B, if one has |X j \ δ mB (X j )|/|X j | ≪ 1 for a scale m, then we have

D (n+m)B (X j ) ≈ ϖ α ( R ∑ k=1 δ nB ( ϕ ⊕mB k ) (p i )H k,j
) . Robustness against noise.The perturbation associated to X is now related to the presence of noise; in particular, since we are dealing with binary images, salt-and-pepper noise is considered. Accordingly, a given percentage p of pixels is corrupted in our dataset to obtain X noise-p . Concerning the noise, it is well known that morphological operators are very sensitive to noise and even a small amount of salt-and-pepper provides a strong perturbation of dilation/erosion. As it is illustrated in Fig. 7, Sparse-NMF dictionary learning produces atoms or basis which are quite robust against this kind of noise. If we denote by ϕ noise-p k the basis obtained from noisy data X ′ p , we have observed that for p up to 20% -30%, we get

∑ i ∥ϕ noise-p k -ϕ k ∥(p i )/ ∑ i ϕ k (p i ) ≈ 0,

Sparse approximation to numerical morphological operators

Approaches introduced in Section 3 are extended here to gray-levels images. In this case, we deal with families of discrete gray-level images, i.e.,

F = {f 1 (p i ), • • • , f M (p i )}, with f j (p i ) ∈ F(E, T ), T = {t 1 , t 2 , • • • , t L } with (t l+1 -t l ) = ∆t.
Fig. 8 provides some examples of ORL face database [START_REF] Guo | Face recognition by support vector machines[END_REF] which is used to illustrate and compare our techniques. We consider in particular two alternative paradigms: i) each gray-scale function f j is represented as a stack of upper level sets and Sparse-NMF processing is applied on the upper level sets, followed by image recomposition from processed upper level sets; ii) straightforward Sparse-NMF representation and processing of gray-scale images

Sparse-NMF processing of upper level sets

The thresholded set of f j at each t l , i.e., X t l j = ϖ t l (f j ), is called the upper level-set at t l of f j . The set of upper level sets constitutes a family of decreasing sets [START_REF] Serra | Image Analysis and Mathematical Morphology. Vol I, and Image Analysis and Mathematical Morphology[END_REF]:

t λ ≥ t µ ⇒ X t λ ⊆ X tµ and X t λ = ∩{X tµ , µ < λ}.
Moreover, any (semi-continuous) image f j can be viewed as an unique stack of its upper level sets, which leads to the following reconstruction property [START_REF] Serra | Image Analysis and Mathematical Morphology. Vol I, and Image Analysis and Mathematical Morphology[END_REF]:

f j (p i ) = sup{t l | p i ∈ X t l j }, t l ∈ T .
We prefer here to consider the alternative reconstruction using a numerical sum of the indicator function of upper level sets [START_REF] Wendt | Stack Filters[END_REF][START_REF] Ch | Bounded variation in posets, with applications in morphological image processing[END_REF]:

f j (p i ) = ∆t L ∑ l=1 X t l j (p i ). ( 22
)
It is well known in mathematical morphology that any binary increasing operator, such as the dilation and erosion, can be generalized to gray-level images by applying the binary operator to each cross-section, and then by reconstructing the corresponding gray-level image [START_REF] Serra | Image Analysis and Mathematical Morphology. Vol I, and Image Analysis and Mathematical Morphology[END_REF][START_REF] Ch | Bounded variation in posets, with applications in morphological image processing[END_REF], i.e.,

δ B (f j )(p i ) = ∆t L ∑ l=1 δ B ( X t l j ) (p i ), ( 23
) ε B (f j )(p i ) = ∆t L ∑ l=1 ε B ( X t l j ) (p i ). (24)
Consider now that each image of the initial gray-level family F of M images is decomposed into its L upper level set. Hence, we have

F = {f 1 , • • • , f M } → X = {X t 1 1 , X t 2 1 • • • , X t L 1 , • • • , X t L M -1 , X t 1 M , • • • , X t L M },
where X is a family of M ′ = M ×L binary images. Therefore, we can use NMF algorithms, for a given dimension R, to approximate each set X t l j and then approximate the corresponding function f j (p i ). Thus, using the results of the previous Section, we are able now to introduce the following definition for the sparse max-approximation to grey-level dilation and erosion given respectively by:

D B (f j )(p i ) = ∆t L ∑ l=1 D B ( X t l j ) (p i ), (25)
and

E B (f j )(p i ) = ∆t L ∑ l=1 E B ( X t l j ) (p i ). (26)
with

D B ( X t l j ) = ϖ α ( R ∑ k=1 δ (ϕ k ) (p i )H k,j+l ) , E B ( X t l j ) = ϖ α ( R ∑ k=1 [ δ B ( [ϕ k ] ) (p i ) ] H k,j+l
) .

Similarly, a sparse max-approximation to gray-scale opening and closing can be achieved by means of the expressions (20) and ( 21) for respectively G B (X j ) and F B (X j ) applied to the upper level sets of f j .

Behavior of sparse max-approximation to gray-level dilation and erosion based on upper level set decomposition is shown in Fig. 9. Note that, for this example, the number of initial images M = 20, and each of them has been quantized in L = 10 gray-levels, i.e., M ′ = 20 × 10 = 200 dimensions and them for Sparse-NMF the number of atoms is fixed to R = 75. Hence a reduction factor which is not significant. Some of the corresponding atoms are given in Fig. 9(a). As expected the quality of the sparse max-approximation to dilation and to erosion depends on the quality of the initial Sparse-NMF reconstruction of the image. For instance, the first face (man with glasses) is not well approximated with the learned NMF basis and hence, its approximated dilation and erosion are also unsatisfactory. On the contrary, in the case of the last image (woman), the results are more relevant. Nevertheless, at this point, we can conclude that the sparse approximation using upper level set stack, which is theoretically consistent with the binary framework discussed in Section 3 is not very useful in practice.

Sparse-NMF representation and processing of gray-scale images

As an alternative to the previous formulation, we propose a straightforward use of Sparse-NMF representation from the image gray-scale image collection

F = {f 1 (p i ), • • • , f M (p i )}.
That is, each column of the data matrix V corresponds to an image, i.e., V i,j = f j (p i ). Then, after Sparse-NMF dictionary learning of dimension R, each image is approximated as

f j (p i ) Sparse-N M F ---------→ f j (p i ) = R ∑ k=1 ϕ k (p i )H k,j .
According to our principle, for any image f j of the family F and given a structuring element B, we can introduce the following operators.

The sparse max-approximation to gray-level dilation and to gray-level erosion are respectively given by

D B (f j )(x i ) = R ∑ k=1 δ B (ϕ k ) (x i )H k,j , ( 27
) E B (f j )(x i ) = R ∑ k=1 [ δ B ( [ϕ k ] ) (x i ) ] H k,j . ( 28
)
In addition, the sparse max-approximation to gray-level opening and gray-level closing are respectively defined as

G B (f j )(x i ) = R ∑ k=1 γ B (ϕ k ) (x i )H k,j , ( 29
)
F B (f j )(x i ) = R ∑ k=1 [ γ B ( [ϕ k ] ) (x i ) ] H k,j . ( 30
)
We must point out again that these approximate nonlinear operators do not satisfy the standard properties of grey-level dilation and erosion.

Similarly to the binary case, let us discuss the conditions under which these max-approximations to dilation and opening are more relevant. First, let us rewrite (27) as supremum of cylinders [START_REF] Heijmans | Morphological Image Operators[END_REF][START_REF] Bloch | Mathematical Morphology[END_REF]:

δ B (f ) = ∨ { C B(p i ), f (p i ) : p i ∈ E } ,
where C B,t is the cylinder of base B at the origin and height t. From this expression, it is easy to see that dilation commutes with supremum of functions, i.e.,

δ B ( ∨ k∈K f k ) = ∨ { C B(p i ), ∨ k∈K f (p i ) : p i ∈ E } = ∨ { ∨ k∈K C B(p i ), f k (p i ) : p i ∈ E } = ∨ k∈K ∨ { C B(p i ), f k (p i ) : p i ∈ E } = ∨ k∈K δ B (f k )
For the case of the opening ( 29), one has [START_REF] Heijmans | Morphological Image Operators[END_REF][START_REF] Bloch | Mathematical Morphology[END_REF]:

γ B (f ) = ∨ { C B(p i ), t : p i ∈ E, C B(p i ), t ≤ f } .
Thus, a sufficient condition for the commutation of opening with the supremum is the "disjointness" of functions, i.e., ∧ k∈K f k = 0. It is obvious also that for "disjoint functions", one has:

∨ k∈K f k = ∑ k∈K f k .
Therefore, operators ( 27) and (29) will provide a good approximation to dilation and opening in the case where the Sparse-NMF representation is composed of disjoint atoms, i.e., ∧ ϕ k (p i ) ≈ 0, which turns out to be a condition of non-redundancy of the Sparse-NMF basis.

Fig. 10(a) gives the Sparse-NMF representation of the faces, where a collection of M = 80 images has been used in this experiment and dimension of dictionary has been fixed to R = 20 atoms. As we can observe, the effectivity of the max-approximation to dilation and opening is consistent with the approximation to the corresponding original image. Thus, such approximations can be then used to compute evolved morphological operators. We have in particular illustrated in Fig. 10 two operators which can be applied in face image feature extraction. The first one is the morphological gradient, i.e., δ ( B)(f ) -ε ( B)(f ), typically used for contour detection. The second one, the black top-hat, i.e., φ ( B)(f ) -f , extracts the dark structures, which in the present case of faces, correspond just to the eyes and lips. negative representation is compatible with linear physical model of endmembers/abundances. Fig. 13 shows some of the spectral bands of Pavia image in (a) as well as the classical NMF representation compared to the Sparse-NMF counterpart (with S w = 0 and S h = 0.5), where for both cases R = 5. We have assessed the effect of the sparseness parameter S h on the approximation error (in average percentage) to morphological operators: max-approximation to dilation D B (27), to erosion E B (28) and to opening G B (29), for the the M = 103 spectral bands. From this example we clearly observe that with S h = 0.5 the approximation error is much better than without any sparseness. We note also that a low value of S h can eventually produces worst results than unconstrained NMF.

(p i ) (c) Sparse-NMF approximation f j (p i ) (d) Dilated images δ B (f j )(p i ) (e) Sparse max-approximation dilat. D B (f j )(p i ) (f) Opened images B (f j )(p i ) (g) Sparse max-approximation open. G B (f j )(p i ) (h) Black top-hat images φ B (f j )(p i ) -f j (p i ) (i) Sparse max-approximation black top-hat F B (f j )(p i ) -f j (p i ) (j) Morphological gradient images δ B (f j )(x i ) -ε B (f j )(p i ) (k) Sparse max-approximation morpho. gradient D B (f j )(x i ) -E B (f j )(p i )
Mathematical morphology is widely used in remote sensing hyperspectral imaging. One of the most popular applications in the state-of-the-art is the spectral-spatial classification based on the notion of morphological profile [START_REF] Pesaresi | A new approach for the morphological segmentation of high resolution satellite imagery[END_REF], and its various extensions to hyperspectral images (Fauvel et al., 2008;[START_REF] Velasco-Forero | Classification of hyperspectral images by tensor modeling and additive morphological decomposition[END_REF]. Morphological profiles are founded on the granulometry (or opening-based scale-space). We remind that a granulometry [START_REF] Serra | Image Analysis and Mathematical Morphology. Vol I, and Image Analysis and Mathematical Morphology[END_REF] is a family of openings {γ Bn } 0≤n≤L , depending on a (discrete) positive parameter n, representing scale of homothetic structuring element, i.e., B n = nB. Then, given an image f , it produces a multi-scale decomposition of bright structures:

s n (f )(p j ) = f (p j ) -γ Bn (f )(p j ),
or using a differential representation:

r n (f )(p j ) = γ B n-1 (f )(p j ) -γ Bn (f )(p j )
such that:

f (p j ) = γ B L (f )(p j ) + L ∑ n=1 r n (f )(p j ).
Then, this multi-scale decomposition can be used either for tensor-based spatial-spectral dimensionality reduction (Velasco-Forero and Angulo, 2013) (4D tensor = 2D space × 1D spectrum × 1D morphology):

{ {f λ } 1≤λ≤M , {γ B 1 (f λ )} 1≤λ≤M • • • , {γ B L (f λ )} 1≤λ≤M } ;
or for spatial-spectral classification, where the feature vector per pixel p j given by:

(f 1 (p j ), f 2 (p j ), • • • , f M (p j ), γ B 1 (f 1 )(p j ), γ B 1 (f 2 )(p j ), • • • , γ B L (f M )(p j )) ,
is the so-called morphological profile at p j . As it is shown in Fig. 15 for two examples of spectral bands, the max-approximation to opening-based scale-space using Sparse-NMF representation is quite satisfactory and can be naturally used in morphological profile-based spectral spatial classification. Figure 15: Comparison of opening-based scale-space f i (x), s 1 (f i )(x), s 3 (f i )(x), s 5 (f i )(x), s 7 (f i )(x) for exact openings and for Sparse-NMF max-approximations (R = 5, S h = 0.5). In (a) for spectral band f 1 (p j ) and in (b) for f 3 (p j ).

ϕ 1 (p i ) ϕ 2 (p i ) ϕ 3 (p i ) ϕ 4 (p i ) (a) (b) (c) X 1 X 5 (d) δ B (X i ) (e) ∑ 4 k=1 δ B (ϕ k ) (x j )H k,i (f) D B (X i )
• • f 1 (p j ) f 2 (p j ) f 3 (p j ) f 103 (p j ) (a)

Conclusions and perspectives

We have introduced the notion of sparse binary and gray-level max-approximation to morphological operators based on Sparse-NMF representation.

We have linked the role of non-negative decompositions of sets of functions, viewed as vectors in a high dimensional space, to the commutation by supremum of the dilation operator, which is just the counterpart of convolution in max-plus algebra. We have also seen that in the case of the opening operator, additional conditions of separability (or disjointness) are required for an exact compatibility. Empirical examples showed that this assumption can be relaxed since the obtained approximations are satisfactory in practical cases.

We have illustrated the practical interest of our approach for morphological processing of multivariate vector images, namely multispectral/hyperspectral images. Our results are relatively encouraging and they open a new avenue to study how the current paradigm of sparse modeling (based mainly on linear operations) in computer vision can be extended to the nonlinear morphological framework.

As we have discussed, NMF produces linear non-negative decompositions which are well indicated for morphological operators in the case of separable functions, i.e., decomposition by sum is equivalent to decomposition by max.

It is well known that morphological operators are linear in the max-plus algebra. Consequently, the intrinsic linear decomposition in such algebra is the more appropriate one to introduced sparse morphological processing, without separability conditions. Some recent work has addressed the problem of matrix factorization in max-plus algebra, see for instance [START_REF] Hook | Max-plus singular values[END_REF], which is based on some classical results on max-plus spectral theory and max-plus eigenvalues as the asymptotic behavior of standard eigenvalues by means of a nonlinearization of the matrix [START_REF] Gaubert | Generic asymptotics of eigenvalues using min-plus algebra[END_REF][START_REF] Schutter | The QR decomposition and the singular value decomposition in the symmetrized max-plus algebra revisited[END_REF]. Other recent results on matrix factorization over max-times algebra [START_REF] Karaev | Matrix Factorization over Max-Times Algebra for Data Mining[END_REF] can also be used as a representation approach for morphological operators. To explore the interest of decompositions on max-plus and max-times algebras for morphological processing will be the object of future research.

Figure 1 :

 1 Figure 1: Sparse mathematical morphology: motivation and outline.

  Fig.4illustrates qualitatively the better performance of Sparse-NMF (S w = 0, S h = 0.6) with respect to the others. We have also included in Fig.3(a) the 10 first eigenimages obtained by PCA (in red the positive values and in blue the negative values); as expected, the corresponding representation does not fit with a part-based decomposition needed for morphological operators.

Figure 2 :

 2 Figure 2: Examples of original binary images from the fish database, composed of 100 images.

  Figure 3: PCA vs. non-negative representation of binary shapes. A collection of M = 100 shapes has been used in the experiments (see examples in Fig. 2), where the number of reduced dimensions has been fixed to R = 10 (in the examples are given the first 9 basis images).

  Figure 4: Sparse-NMF approximations to binary sets: (a) three original shapes X j ; (b)-(e) Top, reconstructed function a X j and Bottom, approximate set X j .

  Figure 5: Comparison of dilation/erosion (a)/(b) vs. sparse pseudo operators for three examples of the Fish shapes. It is compared in particular the sparse max-approximation to dilation/erosion for the standard NMF (c)/(d) and for the Sparse-NMF (e)/(f), with (S w = 0,S h = 0.6). The structuring element B is a square of 5 × 5 pixels.

Figure 6 :

 6 Figure 6: Consistency of Sparse-NMF representation for binary dilation: (a)-(c) Sparse-NMF basis images (S w = 0,S h = 0.6), M = 100 and R = 10; (d) max-approximation to dilation on original image collection vs. on dilated image collection.

  ∀k, and therefore we obtainD B ( X noise-p j ) ≈ D B (X j ).Fig. 7(d) shows examples of the robustness of sparse max-approximation to dilation for p = 10% and p = 20%. In conclusion, Sparse-NMF representation allows us (pseudo-)morphological processing noisy images.

Figure 7 :

 7 Figure 7: Robustness against noise of Sparse-NMF representation for binary dilation: (a)-(c) Sparse-NMF basis images (S w = 0,S h = 0.6), M = 100 and R = 10; (d) noisy binary images, dilated noisy images and max-approximation to dilation of noisy images.

Figure 8 :

 8 Figure 8: Examples of original numerical images from the ORL face database.

  Figure 9: A collection of M = 20 faces has been used in this experiment, where the number of reduced dimensions for the binary matrix V has been fixed to R = 75 (note that M ′ = 20 × 10 = 200 dimensions). Top, some examples of the Sparse-NMF atoms obtained from the 200 upper level sets. Bottom, four examples of the ORL face database (b) (quantized in L = 10 gray-levels) and (c) corresponding approximated images using Sparse-NMF. Comparison of dilation/erosion (d)/(f) vs. sparse max-approximation to dilation/erosion for the Sparse-NMF (e)/(g). The structuring element B is a square of size 3 × 3 pixels.

Figure 10 :

 10 Figure 10: Sparse-NMF representation and processing of gray-scale images F: A collection of M = 80 faces has been used in this experiment and dimension of dictionary has been fixed to R = 20. Top, the obtained Sparse-NMF atoms. Bottom, four examples of the ORL face database (b) and (c) corresponding approximated image using Sparse-NMF. Comparison of dilation, opening, black top-hat and gradient (d)/(f)/(h)/(j) vs. sparse max-approximation to this operators. The structuring element B is a square of size 5 × 5 pixels for the dilation and opening, size 9 × 9 for the black top-hat.

Figure 12 :

 12 Figure 12: Dimensionality reduction of multivariate Boolean set X from Fig. 11: (a) Sparse-NMF basis with R = 4; (b) three spectra used to simulate multivariate set X; (c) coding vectors H associated to basis of (a); (d) two dilated components of X; (e) approximated dilation using the Sparse-NMF basis; (d) thresholded images.

Figure 14 :

 14 Figure 13: Pavia hypersectral images (crop of a part of the whole image): (a) some of the spectral bands, (b) classical NMF representation, (c) Sparse-NMF representation (S w = 0, S h = 0.5). In both cases, left, basis {ϕ k (p j )} 1≤k≤5 and right, coding H.

Roughly, the intrinsic dimension of an image f is defined as the number of "independent" parameters needed to represent f .

Applications

We consider in this Section to potential applications of Sparse-NMF morphological processing to problems arising from multi/hyperspectral image processing.

Sparse processing of multivariate Boolean textures

The motivation of this study is to deal with random image models on multivalued images, and more precisely with the case of multivariate Boolean random set model [START_REF] Jeulin | Multivariate Random Image Models[END_REF]. This approach is an extension of the classical Boolean random model [START_REF] Matheron | Random Sets and Integral Geometry[END_REF]. Multivariate texture images which fit with this model appear in different fields: multi/hyperspectral images, energy dispersive spectrometry, electron energy-loss spectrometry, etc.

Let us consider a Poisson point process in a set E ⊂ R n with intensity (average number of point per unit of volume) θ. In each point x k ∈ E of the process, an independent realization of a multivariate random compact set , 255], θ|E| = 30, and where A i are random disks is given in Fig. 11. In practical situations, this kind of binary textures can be for instance obtained after segmentation of multi/hyperspectral images.

In the framework of the theory of random sets, these structures and their models are fully characterized by a functional called Choquet capacity [START_REF] Matheron | Random Sets and Integral Geometry[END_REF]. In the present model, the multivariate Choquet capacity T X (K), with

, is defined as [START_REF] Jeulin | Multivariate Random Image Models[END_REF]:

µ n being the Lebesgue measure in R n and where X ⊕ K is the dilation of set X ∈ P(E) by structuring element K, i.e.,

where K = {-p i : p i ∈ K}. Hence, an experimental estimation of T X (K) can be obtained by morphological dilations on realizations of components of K. If a single grain is used for each Poisson point, the components are independent Boolean random sets [START_REF] Matheron | Random Sets and Integral Geometry[END_REF][START_REF] Jeulin | Multivariate Random Image Models[END_REF]; otherwise some correlations between components are present. As a fundamental property, this model is stable by the union; such that any number of components (even correlated) of a Boolean multivariate random set is still a Boolean

random set. Hence a notion of dimensionally reduction by supremum of primary schemes is compatible with the model. In our current example of Fig. 11, the A i random disks have been simulated using three vector prototypes, see Fig. 12(b) which consequently involves correlations between the components. Using our theory of Sparse-NMF representation and max-approximation to dilation we note that, on the one hand, the Sparse-NMF basis {ϕ k (p j )} gives just a dimensionality reduction of the initial multivariate Boolean model. In other words, thresholding the four components of Fig. 12(a) produces a Boolean model of dimension 4 instead of dimension 20 of the original model. We observe in the coding vectors H associated to basis, Fig. 12(c), that the four components are consistent with the decomposition of the original spectra Fig. 12(d). On the other hand, approximated dilations to the initial components of the model can be naturally obtained using expression (15), Fig. 12(d)-(f).

In summary, in order to model a multispectral Boolean set X, we need to compute the characteristic curves related to Choquet capacity T X (K), for different K, and these computations can be efficiently approximated in a consistent theoretical way with our Sparse-NMF dilations.

Sparse processing of hyperspectral images

High-spatial resolution hyperspectral images are used nowadays in remote sensing and other application domains. They constitute high dimensional datasets: for instance, the Pavia image, used in this case-study, has dimensions N = 340 × 610 and M = 103. As we have discussed in the introduction, NMF is widely used in hyperspectral imaging since the non-