Ludivine Nus

Sebastian Miron

David Brie

An ADMM-based algorithm with minimum dispersion constraint for on-line blind unmixing of hyperspectral images

Keywords: Hyperspectral imaging, Pushbroom acquisition system, On-line unmixing, Alternating Direction Method of Multipliers, Minimum dispersion constraint

Pushbroom imaging systems are emerging techniques for real-time acquisition of hyperspectral images. These systems are frequently used in industrial applications to control and sort products on-the-fly. In this paper, the on-line hyperspectral image blind unmixing is addressed. We propose a new on-line method based on Alternating Direction Method of Multipliers (ADMM) approach, adapted to pushbroom imaging systems. Because of the generally ill-posed nature of the unmixing problem, we impose a minimum endmembers dispersion constraint to regularize the solution; this constraint can be interpreted as a convex relaxation of the minimum volume constraint and therefore, presents interesting optimization properties. The proposed algorithm presents faster convergence rate and lower computational complexity compared to the algorithms based on multiplicative update rules. Experimental results on synthetic and real datasets, and comparison to state-ofthe-art algorithms, demonstrate the effectiveness of our method in terms of rapidity and accuracy.

Introduction

Hyperspectral imaging, which combines the power of digital imaging and spectroscopy, is an effective tool in a wide range of applications (see e.g., [START_REF] Van Der Meer | Multi-and hyperspectral geologic remote sensing: a review[END_REF][START_REF] Qin | Line-scan hyperspectral imaging techniques for food safety and quality applications[END_REF][START_REF] Khan | Modern trends in hyperspectral image analysis: a review[END_REF]). Each pixel in a hyperspectral image provides local spectral information about a scene of interest across a large number of contiguous bands. Because of the limited spatial resolution of the sensor, mixed pixels (pixels containing the contributions of several components) are often encountered in hyperspectral data. Thus, hyperspectral unmixing is an important technique for hyperspectral data interpretation, as it allows to decompose a mixed pixel into a collection of spectral signatures (also called endmembers) and their relative proportions (also called abundances).

This paper addresses the problem of sequential (or on-line) unmixing of hyperspectral Near InfraRed (NIR) images acquired by a pushbroom imager [START_REF] Lu | Medical hyperspectral imaging: a review[END_REF], by means of Non-negative Matrix Factorization (NMF)-like approaches. This problem is encountered e.g., in real-time industrial systems, for product quality control applications. In particular, the approach proposed in this paper is part of a project that aims at predicting and classifying in realtime, the rendering quality of pieces of wood in an industrial process. This is to be done by on-line unmixing method of hyperspectral images of wood surface, which presents several advantages over the off-line methods: i) it is well-adapted to real-time data processing for on-line industrial acquisition systems; ii) it allows to alleviate computational burden and reduce memory requirements for big hyperspectral data cubes; iii) it permits to track the spatial/time variability of the endmembers in a hyperspectral imaging application.

Hyperspectral image unmixing

Hyperspectral unmixing methods are based on a mixture model describing how the endmembers are combined in the acquired image. Depending on the application and/or on the data, this model can be linear or non-linear [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: geometrical, statistical, and sparse regression-based approaches[END_REF]. In this article, we focus on the Linear Mixing Model (LMM) [START_REF] Keshava | Spectral unmixing[END_REF]; it is the most commonly used model in hyperspectral unmixing because it is simple to understand and represents a good approximation of the physical reality in many applications. The LMM model can be generally expressed in the following form: , where R + denotes the set of non-negative real numbers. In hyperspectral imaging, the P columns of X represent the data samples (pixels) recorded at L wavelengths. S is a matrix containing on its columns the R normalized endmembers and A is a matrix containing on its rows the abundances for the recorded samples. Under the non-negativity and the sum-to-one constraint i.e., S ≥ 0, A ≥ 0 and A T 1 = 1 (where 1 is an all-ones vector), the LMM admits a geometric interpretation. Figure 1 illustrates this geometrical representation for a mixture of R = 3 endmembers: all the pixels (red points) belongs to the (R -1)-dimensional simplex (blue lines) whose vertices are the endmembers s 1 , ..., s R (gray points). When the sum-to-one constraint does not hold, the observations belong to the positive cone defined by the endmembers. There are many methods in the literature for hyperspectral unmixing. Most of them are based on the pure pixel assumption, i.e., the existence of pixels containing a single source, see e.g., [START_REF] Winter | An algorithm for fast autonomous spectral endmember determination in hyperspectral data[END_REF][START_REF] Nascimento | Vertex Component Analysis: a fast algorithm to unmix hyperspectral data[END_REF]. However, the pure pixel assumption is a strong requirement that does not hold in general for highly mixed data, and finding the endmembers is therefore a more challenging task. For this case, another family of methods has been developed, based on the volume minimization idea introduced by Craig in 1994 [START_REF] Craig | Minimum-volume transforms for remotely sensed data[END_REF]. These approaches consist in minimizing the volume of the simplex containing the data by "creating" virtual endmembers such as in [START_REF] Berman | ICE: a statistical approach to identifying endmembers in hyperspectral images[END_REF][START_REF] Li | Minimum volume simplex analysis: a fast algorithm to unmix hyperspectral data[END_REF]. However the performance and computational efficiency of these methods are often limited because of complicated simplex volume calculations, sensitivity to initialization and lack of rigorous performance analysis. A detailed analysis can be found in [START_REF] Bioucas-Dias | Hyperspectral unmixing overview: geometrical, statistical, and sparse regression-based approaches[END_REF].

X ≈ SA, (1)
To overcome the up-mentionned problems, Non-negative Matrix Factorization (NMF) [START_REF] Lee | Learning the parts of objects by Non-negative Matrix Factorization[END_REF] has been applied to hyperspectral data unmixing. For a non-negative matrix X, the NMF consists in estimating two matrices, S ≥ 0 and A ≥ 0, satisfying (1). In general, the NMF is not unique and therefore, to reduce the size of the set of admissible solutions, it is necessary to add regularization terms. The most effective constraint approaches are volume-regularized NMF such as [START_REF] Miao | Endmember extraction from highly mixed data using minimum volume constrained Non-negative Matrix Factorization[END_REF][START_REF] Yu | Minimum distance constrained Non-negative Matrix Factorization for the endmember extraction of hyperspectral images[END_REF][START_REF] Schachtner | Minimum determinant constraint for Non-negative Matrix Factorization[END_REF][START_REF] Zhou | Minimum-volumeconstrained Non-negative Matrix Factorization: enhanced ability of learning parts[END_REF][START_REF] Fu | Robust volume minimization-based matrix factorization for remote sensing and document clustering[END_REF][START_REF] Ang | Volume regularized Non-negative Matrix Factorizations[END_REF][START_REF] Fu | Non-negative matrix factorization for signal and data analytics: identifiability, algorithms, and applications[END_REF][START_REF] Leplat | Minimum-volume rank-deficient Nonnegative Matrix Factorizations[END_REF], which can be considered as state-of-the-art methods in blind hyperspectral unmixing.

Pushbroom acquisition scheme

In pushbroom imaging systems, hyperspectral data cubes are acquired slice by slice, sequentially in time. Each slice is an image, characterized by a spectral dimension and a spatial dimension (one line of the scene). Figure 2 illustrates the test bench that we consider for the real-time characterization of the wood material. Hyperspectral images are acquired line by line, as the pieces of wood are carried under the imager, via a conveyor. The stream of spectral-spatial data arrays is then stacked to form the hyperspectral data cube. For each acquisition time k (k = 1, ..., K), the new slice (represented by a dotted line in Figure 3) is a matrix of dimensions L × P , where L denotes the spectral dimension (wavelengths) and P , the across track spatial dimension (one line of the scene). The goal of the on-line blind hyperspectral unmixing is to produce realtime estimates of the endmember and abundance matrices for each new incoming slice, at the pushbroom system acquisition rate.

On-line NMF methods

On-line NMF algorithms sequentially update the endmembers and abundances as the data size increases. In that respect, they are perfectly adapted to the processing of hyperspectral data streaming as they allow to maintain a low and controlled computational complexity. These algorithms can be gathered into two main categories, depending on the considered assumptions on the endmembers. In [START_REF] Mairal | On-line learning for matrix factorization and sparse coding[END_REF][START_REF] Wang | Efficient document clustering via Online Non-negative Matrix Factorizations[END_REF][START_REF] Guan | On-line Non-negative Matrix Factorization with robust stochastic approximation[END_REF][START_REF] Zhao | On-line Non-negative Matrix Factorization with outliers[END_REF], the endmembers do not vary from one sample to another, while in [START_REF] Bucak | Incremental subspace learning via Non-negative Matrix Factorization[END_REF][START_REF] Lefèvre | On-line algorithms for Non-negative Matrix Factorization with the Itakura-Saito divergence[END_REF][START_REF] Zhou | On-line blind source separation using Incremental Non-negative Matrix Factorization with volume constraint[END_REF][START_REF] Wang | On-line learning parts-based representation via incremental orthogonal projective Non-negative Matrix Factorization[END_REF][START_REF] Wu | Visual tracking via On-line Non-negative Matrix Factorization[END_REF][START_REF] Nus | On-line blind unmixing for hyperspectral pushbroom imaging systems[END_REF], the endmembers may evolve between successive samples. In particular, the Incremental NMF (INMF) [START_REF] Bucak | Incremental subspace learning via Non-negative Matrix Factorization[END_REF] considers that the endmembers evolve slowly between two consecutive acquisitions; this is now the most widely used assumption adopted in on-line NMF algorithms. In the context of hyperspectral unmixing, it allows to account for the so-called spectral variability. However, as in the off-line case, the uniqueness of the solution is not guaranteed, which led to the development of on-line regularized NMF such as [START_REF] Zhou | On-line blind source separation using Incremental Non-negative Matrix Factorization with volume constraint[END_REF][START_REF] Wang | On-line learning parts-based representation via incremental orthogonal projective Non-negative Matrix Factorization[END_REF][START_REF] Wu | Visual tracking via On-line Non-negative Matrix Factorization[END_REF]. Recently, the on-line Minimum Volume Constraint-NMF algorithm (OMVC-NMF) was introduced in [START_REF] Nus | On-line blind unmixing for hyperspectral pushbroom imaging systems[END_REF]. It is a straightforward adaptation of [START_REF] Zhou | On-line blind source separation using Incremental Non-negative Matrix Factorization with volume constraint[END_REF] specially designed for pushbroom hyperspectral imaging system. To the best knowledge of the present authors, OMVC-NMF is the only unmixing algorithm in the literature, adapted to on-line processing of hyperspectral images. Thus, it will serve as benchmark for the on-line volume regularized NMF method proposed in this paper. Most of the NMF algorithms mentioned above are based on multiplicative update rules which are known to be highly sensitive to initialization and often suffer from slow convergence rate [START_REF] Lin | Projected gradient methods for Non-negative Matrix Factorization[END_REF]. These limitations are restrictive in particular for the considered application; indeed, the on-line unmixing algorithms must be fast in order to fulfill the real-time industrial production constraints. An interesting alternative to multiplicative updates is represented by the optimization methods based on Alternating Direction Method of Multipliers (ADMM) in [START_REF] Zhao | On-line Non-negative Matrix Factorization with outliers[END_REF]; ADMM proved its superiority over multiplicative updates with respect to both reconstruction accuracy and convergence rate [START_REF] Xu | An alternating direction algorithm for matrix completion with non-negative factors[END_REF][START_REF] Sun | Alternating direction method of multipliers for Non-negative Matrix Factorization with the beta-divergence[END_REF][START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF].

Main contributions

In this paper, we introduce a new algorithm for on-line blind unmixing of hyperspectral images, specially devised for pushbroom acquisition systems. The original contributions of this algorithm compared to our previous work [START_REF] Nus | On-line blind unmixing for hyperspectral pushbroom imaging systems[END_REF] and to state-of-the-art on-line approaches (see e.g., [START_REF] Zhao | On-line Non-negative Matrix Factorization with outliers[END_REF]) are:

1. The addition of the Minimum Dispersion Constraint (MDC) [START_REF] Berman | ICE: a statistical approach to identifying endmembers in hyperspectral images[END_REF], to regularize the problem. MDC can be interpreted as a convex relaxation of the minimum volume simplex constraint. Thanks to its convexity, this constraint offers interesting optimization properties and enables explicit updates of the parameters. 2. The integration of a tracking capability to the algorithm which allows to model dynamic content changes. 3. The use of ADMM-based optimization in the context of on-line volume regularized NMF.

The remainder of this paper is organized as follows: Section 2 is dedicated to the formulation of the on-line NMF problem for the acquisition scheme of a pushbroom imager. Section 3 provides a theoretical study of different volume regularizers. In particular, we highlight the differences between the minimum dispersion and the determinant-based constraint. Section 4 presents the derivation of the proposed on-line ADMM approach; convergence and computational complexity are also discussed. Finally, in Section 5, we give extensive experimental results on both synthetic and real data. Comparison to off-line (or batch) algorithms as well as to an on-line matrix factorization method, are also provided. We conclude the paper in Section 6.

2. On-line hyperspectral unmixing for pushbroom acquisition scheme

Data model

The principle of the proposed on-line method is to alternatively update the endmember and abundance matrices estimated at time instant k when a new sample (slice) arrives at time instant k + 1. One way to handle the problem is to unfold the hyperspectral image as shown in Figure 4, where X(1) = X (1) is the first slice of the hyperspectral image and X(k) is the k th slice. The entire dataset at time instant k +1, i.e., X (k+1) , can be represented as the concatenation of the first k samples with the new incoming sample i.e., X (k+1) = X (k) X(k+1) . Similarly, we define S (k+1) = S (k) S(k+1) and

A (k+1) = A (k) Ã(k+1) .
X [START_REF] Van Der Meer | Multi-and hyperspectral geologic remote sensing: a review[END_REF] ... X(k) X(k+1) L kP P P P Then, the on-line NMF model is given by: X(k+1) ≈ S(k+1) Ã(k+1) .

X (k)
(2)

Cost function

The simplest way to fit the NMF model to the data is to minimize the least square distance between the data and the model. Let J (k) denote the cost function corresponding to the first k samples:

J (k) S (k) , A (k) = 1 2 k ℓ=1 X(ℓ) - S(ℓ) Ã(ℓ) 2 F .
When the k + 1 th sample, X(k+1) arrives, the corresponding cost function can be decomposed as follows:

J (k+1) S (k+1) , A (k+1) = 1 2 k ℓ=1 X(ℓ) - S(ℓ) Ã(ℓ) 2 F + 1 2 X(k+1) - S(k+1) Ã(k+1) 2 F .
(3) Without further assumptions, (3) is just a set of k independent least squares problems and in this case, the on-line setup has no particular interest. However, a natural assumption is that the endmembers vary only slightly between consecutive samples i.e., S(k+1) ≈ S(k) , ∀k. Thus, the cost function (3) can be expressed as:

J (k+1) S(k+1) , Ã(k+1) A (k) = J (k) S(k+1) A (k) + J (k+1) S(k+1) , Ã (k+1)
, with

J (k) S(k+1) A (k) = 1 2 k ℓ=1 X(ℓ) - S(k+1) Ã(ℓ) 2 F , J (k+1) S(k+1) , Ã(k+1) = 1 2 X(k+1) - S(k+1) Ã(k+1) 2 F .
In order to add some tracking capability to the algorithm, a weighting coefficient α (0 ≤ α ≤ 1) is incorporated into the cost function as:

J (k+1) S(k+1) , Ã(k+1) = αJ (k) + (1 -α) J (k+1) . (4)
A version of the cost function ( 4) is used by the INMF algorithm proposed in [START_REF] Bucak | Incremental subspace learning via Non-negative Matrix Factorization[END_REF]. Nevertheless, without additional constraint, the solution of (4) is not unique. The uniqueness of the NMF relies on the sparsity of the underlying latent variables. In particular, if either S and/or A has only non-zero entries, the NMF factorization is not unique. To reduce the size of the set of admissible solutions, we propose to impose on the endmembers matrix S(k+1) , a geometric constraint which forces the simplex bounded by the endmembers to circumscribe the data as closely as possible. For now, we denote this constraint by Vol S(k+1) and integrate it into the cost function as follows:

J (k+1) Vol S(k+1) , Ã(k+1) = αJ (k) + (1 -α) J (k+1) + µ Vol S(k+1) , (5)
where J Vol becomes the volume regularized criterion and µ ≥ 0 controls the trade-off between the data fitting term and the volume regularizer. The choice of Vol S(k+1) is discussed in the next section. Thus, the approach proposed in this paper aims at solving the following optimization problem: minimize

S(k+1) ≥0, Ã(k+1) ≥0 J (k+1) Vol S(k+1) , Ã(k+1) , (6)
for a particular choice of the volume regularizer term.

Volume regularization

The classical measure for the minimum volume constraint is the determinant, i.e. Vol S(k+1) = det S(k+1) [START_REF] Strang | Linear algebra and its applications[END_REF]. However, this determinant is defined only if S(k+1) is a square matrix. In the case where S(k+1) is a tall matrix, a pertinent choice for Vol S(k+1) is [START_REF] Schachtner | Minimum determinant constraint for Non-negative Matrix Factorization[END_REF][START_REF] Zhou | Minimum-volumeconstrained Non-negative Matrix Factorization: enhanced ability of learning parts[END_REF]:

Vol S(k+1) = det S(k+1) T S(k+1) , (7)
or its variants presented in [START_REF] Miao | Endmember extraction from highly mixed data using minimum volume constrained Non-negative Matrix Factorization[END_REF][START_REF] Chan | A convex analysisbased minimum-volume enclosing simplex algorithm for hyperspectral unmixing[END_REF]. The expression [START_REF] Winter | An algorithm for fast autonomous spectral endmember determination in hyperspectral data[END_REF] is mathematically justified because it can be interpreted as the square volume of the simplex bounded by the endmembers [37, Theorem 7].

Identifiability

A lot of work has been done in the last decade to understand the identifiability of the NMF model. The results in [START_REF] Donoho | When does Non-negative Matrix Factorization give a correct decomposition into parts?[END_REF][START_REF] Laurberg | Theorems on positive data: on the uniqueness of NMF[END_REF][START_REF] Huang | Non-negative Matrix Factorization revisited: uniqueness and algorithm for symmetric decomposition[END_REF] have shown that if the matrices S(k+1) and Ã(k+1) both contain null elements, according to a certain pattern, then the model is unique. However, in hyperspectral imaging, the endmembers S(k+1) are likely to be strictly positive and thus, the model is not indentifiable. This is the case where volume minimization approaches can be successfully used. By applying the recent results of [START_REF] Fu | Robust volume minimization-based matrix factorization for remote sensing and document clustering[END_REF] to our on-line mixture model, the following sufficient identifiability conditions can be formulated:

Sufficient identifiability conditions: if ∀ k, rank S(k) = rank Ã(k) =
R and Ã(k) is sufficiently scattered (see [START_REF] Fu | Robust volume minimization-based matrix factorization for remote sensing and document clustering[END_REF] for the exact definition of sufficiently scattered), then the model ( 2) is identifiable under the minimum volume constraint. This condition stipulates that there is no rank loss over the different slices of the hyperspectral image.

Volume regularizers

In this section, to simplify the notations, the upper indices indicating the slice number are omitted. The use of the minimum volume constraint in [START_REF] Winter | An algorithm for fast autonomous spectral endmember determination in hyperspectral data[END_REF] makes the S subproblem non-convex and therefore, more difficult to solve. To tackle the difficulties due to determinant minimization, others volume regularizer surrogates can be considered. An alternative formulation for Vol S , proposed in [START_REF] Bioucas-Dias | A variable splitting augmented Lagrangian approach to linear spectral unmixing[END_REF] is:

log det ST S . (8)
The choice of the logarithm of the determinant rather than the determinant itself is mainly motivated by algorithmic reasons, since it drastically simplifies the update rule for S. To avoid strong negative values of (8) when S becomes close to singularity, i.e. det ST S → 0, a modified version was proposed in [START_REF] Fu | Robust volume minimization-based matrix factorization for remote sensing and document clustering[END_REF][START_REF] Leplat | Minimum-volume rank-deficient Nonnegative Matrix Factorizations[END_REF]:

log det ST S + ǫI , (9)
where ǫ > 0 is a specified small value and I is the identity matrix. While ( 8) is a concave function, it appears that (9) has a stationary point in 0 and is convex in a neighborhood of 0. Indeed, let {λ 2 i , i = 1, ..., R} be the ordered eigenvalues of ST S, that is the λ i 's are the singular values of S, then:

log det ST S + ǫI = R i=1 log λ 2 i + ǫ = R i=1 f (λ i ), with f (λ i ) = log λ 2 i + ǫ , ∂f (λ i ) ∂λ i = 2λ i λ 2 i + ǫ .
Clearly, lim

λ i →0 ∂f (λ i ) ∂λ i = 0. In addition ∂ 2 f (λ i ) ∂λ 2 i = 2ǫ-2λ 2 i (λ 2 i +ǫ) 2 ≥ 0 if λ 2 i ≤ ǫ.
In fact, in the neighborhood of λ i = 0, we can write:

log det ST S + ǫI = R log(ǫ) + R i=1 log λ 2 i ǫ + 1 ≈ R log(ǫ) + 1 ǫ R i=1 λ 2 i = R log(ǫ) + 1 ǫ trace ST S . (10)
In other words, in the neighborhood of λ i = 0, log det ST S + ǫI essentially behaves as trace ST S , the extent of the neighborhood being controlled by the value of ǫ. This has an important consequence: while log det ST S will favor rank deficiency as the regularization parameter µ increases [START_REF] Fazel | Log-det heuristic for matrix rank minimization with applications to Hankel and Euclidean distance matrices[END_REF], the use of log det ST S + ǫI will preserve the full column rank of the solution S even for large values of the regularization parameter.

Finally, a convex surrogate for Vol S was proposed in [START_REF] Berman | ICE: a statistical approach to identifying endmembers in hyperspectral images[END_REF] as:

trace SP ST , (11)
where P = I-1 R 11 T (1 is an all-ones column vector of size R×1). Expression (11) can be interpreted as a measure of the dispersion of the endmembers around their centroid. In the following, we will denote the constraint log det ST S + ǫI by MVC (for Minimum Volume Constraint), and trace SP ST by MDC (for Minimum Dispersion Constraint).

In order to illustrate the difference between these two regularization functions, consider a matrix S composed of two unit column vectors:

S = 1 cos(β) 0 sin(β) ,
with β an angle varying between 0 et π 2 . If β = 0, the vectors are collinear (the volume is zero); a contrario, if β = π 2 , the vectors become orthogonal (the volume is maximal and equal to 1). We plotted in Figure 5 the two functions: MVC with varying ǫ and MDC, for increasing values of β; the curves are normalized to set their minimum value to 0 and their maximum value to 1. For small values of ǫ, MVC promotes rank-deficient solutions but ǫ should no be chosen too small (e.g. 10 -9 ) to avoid a bad conditioning of ST S + ǫI. For larger ǫ values (e.g. 0.1), rank deficiency is no longer promoted and both MDC and MVC preserve the full column rank of S. This rank preserving property is essential to ensure that the sufficient identifiability condition of section 3.1 is satisfied.

OMDC-ADMM algorithm

In the following, we only consider MDC (rather than MVC) because of its convexity which favors full column rank solutions and allows an efficient implementation of on-line ADMM. The resulting algorithm will be termed as On-line MDC-ADMM (OMDC-ADMM).

Algorithm derivation

Let us define:

J (k+1) Vol S(k+1) , Ã(k+1) = αJ (k) + (1 -α) J (k+1) + µ trace S(k+1) P S(k+1) T . ( 12
)
By introducing the auxiliary variables Ũ and Ṽ, problem ( 6) is equivalent to: minimize

S(k+1) , Ã(k+1) , Ṽ(k+1) , Ũ (k+1)
J (k+1) Vol S(k+1) , Ã(k+1) + I R + Ṽ(k+1) +I R + Ũ(k+1) , subject to S(k+1) = Ũ(k+1) and Ã(k+1) = Ṽ(k+1) , (13)
where I R + is the indicator function of R + ensuring the non-negativity of the endmembers and abundances. For algorithmic convenience, we use the scaled form of ADMM [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] in which the linear and quadratic terms are combined in the augmented Lagrangian and the dual variables are scaled. The augmented Lagrangian L for the problem (13) is given by:

L Ã(k+1) , S(k+1) , Ṽ(k+1) , Ũ(k+1) , Π(k+1) , Λ(k+1) = J (k+1) Vol S(k+1) , Ã(k+1) + ρ 2 Ã(k+1) - Ṽ(k+1) + Π(k+1) 2 F - ρ 2 Π(k+1) 2 F + ρ 2 S(k+1) - Ũ(k+1) + Λ(k+1) 2 F - ρ 2 Λ(k+1) 2 F + I R + Ṽ(k+1) + I R + Ũ(k+1) . ( 14
)
The parameter ρ > 0 controls the convergence rate of the method. Λ(k+1) and Π(k+1) are the scaled versions of the dual variables corresponding to the equality constraints S(k+1) = Ũ(k+1) and Ã(k+1) = Ṽ(k+1) , respectively. ADMM optimization alternately minimizes the augmented Lagrangian [START_REF] Yu | Minimum distance constrained Non-negative Matrix Factorization for the endmember extraction of hyperspectral images[END_REF] with respect to Ã(k+1) , Ṽ(k+1) and S(k+1) , Ũ(k+1) and then, updates the dual variables Π(k+1) and Λ(k+1) . The updates at data slice k + 1 can be expressed as:

Ã(k+1) = (1 -α) S(k) T S(k) + ρI -1 (1 -α) S(k) T X(k+1) + ρ Ṽ(k) - Π(k) , (15)
Ṽ(k+1) = max 0, Ã(k+1) + Π(k) , (16)
Π(k+1) = Π(k) + Ã(k+1) - Ṽ(k+1) , (17)
S(k+1) = N (k+1) + ρ Ũ(k) - Λ(k) M (k+1) + 2µP + ρI -1 , (18)
Ũ(k+1) = max 0, S(k+1) + Λ(k) , (19)
Λ(k+1) = Λ(k) + S(k+1) - Ũ(k+1) , (20)
where

N (k+1) = α k l=1 X(l) Ã(l) T + (1 -α) X(k+1) Ã(k+1) T and M (k+1) = α k l=1 Ã(l) Ã(l) T + (1 -α) Ã(k+1) Ã(k+1) T .
Following [START_REF] Bucak | Incremental subspace learning via Non-negative Matrix Factorization[END_REF], under the assumption S(k+1) ≈ S(k) , we can write

N (k) ≈ k l=1 X(l) Ã(l) T and M (k) ≈ k l=1 Ã(l) Ã(l) T
. Thus, N (k+1) and M (k+1) can be expressed as:

N (k+1) = αN (k) + (1 -α) X(k+1) Ã(k+1) T , (21)
M (k+1) = αM (k) + (1 -α) Ã(k+1) Ã(k+1)T . ( 22
)
Algorithm 1 summarizes the proposed OMDC-ADMM algorithm. It includes two main loops: the outer loop produces estimates of all parameters at each new slice. These estimates are iteratively refined in the inner loop using a fixed number of iterations N iter . For notation simplification, the indices k + 1 in the updates rules are omitted. Outputs: A; S; while New sample k + 1 available do X = X(k+1) ;

t = 1; while t < N iter do à = (1 -α) ST S + ρI -1 (1 -α) ST X + ρ Ṽ -Π ; Ṽ = max 0, à + Π ; Π ← Π + à -Ṽ; Ñ = αN + (1 -α) X ÃT ; M = αM + (1 -α) à ÃT ; S = Ñ + ρ Ũ -Λ MS + 2µP + ρI

Convergence

We provide in this section a partial result of convergence for the proposed algorithm: we show that any stationary point generated by a sequence of iterations satisfies the Karush-Kuhn-Tucker (KKT) conditions [START_REF] Boyd | Convex optimization[END_REF]. To simplify the notation, we gather all the variables to be estimated in W = Ã, S, Ṽ, Ũ, Π, Λ . Following [START_REF] Boyd | Convex optimization[END_REF][START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF] and after some basic algebraic manipulations, it can be shown that a point W is a KKT point for the problem (13) if:

(1 -α) ST S Ã -(1 -α) ST X + Π = 0, (23a) S M -Ñ + 2µ SP + Λ = 0, (23b) Ã -Ṽ = 0, (23c) S -Ũ = 0, (23d) Π ≤ 0 ≤ Ṽ, Π ⊙ Ṽ = 0, (23e) Λ ≤ 0 ≤ Ũ, Λ ⊙ Ũ = 0, (23f)
where ⊙ represents the Hadamard product. Let W t := Ãt , St , Ṽt , Ũt , Πt , Λt be the estimate of W at iteration t of Algorithm 1. From ( 15), the following relation can be written between Ãt et Ãt+1 :

(1 -α) ST t St + ρI Ãt+1 -Ãt = -(1 -α) ST t St Ãt -(1 -α) ST t X + ρ Ãt -Ṽt + ρ Πt . ( 24
)
Suppose that the algorithm reaches a stationary point i.e., W t+1 = W t = W ⋆ ; this implies Ãt+1 = Ãt = V t = Ã⋆ . By replacing it in [START_REF] Zhao | On-line Non-negative Matrix Factorization with outliers[END_REF], we obtain:

(1 -α) S⋆ T S⋆ Ã⋆ -(1 -α) S⋆ T X + ρ Π⋆ = 0.
Using similar rationale for the other parameters of W, it can be shown that the first four equalities from the KKT conditions (23) are satisfied for every limit point W ⋆ := Ã⋆ , S⋆ , Ṽ⋆ , Ũ⋆ , Π⋆ , Λ⋆ . To prove (23e), we can write: max 0, Ã⋆ + Π⋆ = Ṽ⋆ . If Ã⋆ = Ṽ⋆ = 0, then max 0, Π⋆ = 0 which leads to Π⋆ < 0. If Ã⋆ = Ṽ⋆ > 0, then Π⋆ = 0. The same kind of rationale also applies to (23f). Thus, we have shown that for the problem (13), any stationary point W ⋆ given by Algorithm 1 satisfies the KKT conditions.

The working assumption to establish this result is S(k+1) ≈ S(k) , meaning that the proposed convergence result is valid only when the steady state is reached; it does not provide any insights into the transient behavior. This point will be addressed through numerical simulations in sections 5.3 and 5.4.

Computational complexity

In this section we evaluate the computational complexity of the proposed algorithm; this is done by taking into account only the matrix multiplication operations, that dominate the algorithm complexity. By considering K slices of size (L × P ), N iter iterations and a decomposition rank R, the computational complexity for OMDC-ADMM algorithm is of the order of 2KN iter (RP L + (L + P )R 2 + R 3 ). For comparison, we developed a batch version of OMDC-ADMM, that we called BMDC-ADMM (for Batch MDC-ADMM ), which processes the entire hyperspectral data cube at once; this algorithm has complexity of 2N iter (RP KL + (L + P K)R 2 + R 3 ). Moreover, the OMVC-NMF algorithm presented in [START_REF] Nus | On-line blind unmixing for hyperspectral pushbroom imaging systems[END_REF], based on multiplicative update rules with minimum volume constraint [START_REF] Craig | Minimum-volume transforms for remotely sensed data[END_REF], has a computational complexity of the order of KN iter (2RP L + (L + P )R + (6L + 2P )R 2 ). Note that the complexity of OMDC-ADMM is comparable to that of BMDC-ADMM and OMVC-NMF. However, as we show in the sequel, OMDC-ADMM requires fewer iterations to converge compared to the two other methods, and thus significantly reduces the computational cost.

Experimental results

In order to evaluate the performance of OMDC-ADMM for hyperspectral unmixing, we conducted several experiments on simulated and real hyperspectral images. For these experiments, we used Matlab (R2016a) on a 2.7 GHz Macbook Pro with 4-core processor and 16 GB of RAM. These experiments have the following objectives:

1. Illustrating the rank preserving properties of MDC and MVC. 2. Showing the efficiency of OMDC-ADMM by examining its convergence speed and comparing it to that of the OMVC-NMF algorithm [START_REF] Nus | On-line blind unmixing for hyperspectral pushbroom imaging systems[END_REF].

The sensitivity of the convergence speed of the proposed method to parameters α and ρ is also studied.

Comparing OMDC-ADMM to its batch counterpart BMDC-ADMM to

assess the advantages of on-line processing of hyperspectral images. 4. Studying the spectral tracking capability of OMDC-ADMM.

Validating the unmixing performance of our algorithm on real hyper-

spectral images of wood.

Performance criteria

Three performance criteria were used for these experiments: the residual error, Spectral Angle Distance (SAD) and Root Mean Square Error (RMSE).

The residual error was calculated for each slice k as follows:

1 2 X -Ŝ Â 2 F , (25)
where Ŝ and  are the estimated endmembers and abundances, respectively.

The SAD measures the similarity between the original endmembers S and the estimated endmembers Ŝ. For a given slice k, it is written as follows:

1 R R r=1 cos -1 sT r ŝr ||s r || 2 || ŝr || 2 , ( 26
)
where R is the number of endmembers, ||.|| 2 , the ℓ 2 -norm, S = [s 1 , ..., sR ] and Ŝ = ŝ1 , ..., ŝR . A low value of SAD indicates that the estimate is close to the ground truth.

The RMSE measures the error between the original abundances A and the estimated abundances  and is computed as:

1 R R r=1 1 KP KP d=1 (a rd -ârd ) 2 , ( 27
)
where a rd and ârd are the ground truth and the estimated abundance respectively of the r th endmember at pixel d (d = 1, ..., KP ).

Note that expressions [START_REF] Lefèvre | On-line algorithms for Non-negative Matrix Factorization with the Itakura-Saito divergence[END_REF] and ( 27) are average values. For some experiments, the SAD and the RMSE values are calculated independently for each source r = 1, ..., R.

Rank preserving properties of MVC and MDC

We performed a numerical experiment to illustrate the rank preserving properties of both volume regularizers, MVC [START_REF] Craig | Minimum-volume transforms for remotely sensed data[END_REF] and MDC [START_REF] Li | Minimum volume simplex analysis: a fast algorithm to unmix hyperspectral data[END_REF]. We used, for the simulations, the endmembers shown in Figure 8(d). None of the three endmembers has any zero value, which results in a non-unique NMF problem. In other words, without additional constraints the model is non-identifiable. The three abundance maps are matrices of size 36 × 36, randomly drawn from a continuous uniform distribution on the interval [0, 1]. By doing so, we ensure (with high probability) that the sufficiently scattered and full column rank conditions are fulfilled (see subsection 3.1). For µ ranging from 0.0001 to 0.003, the quantities trace ŜP ŜT and log det ŜT Ŝ + ǫI were computed, where Ŝ is the endmember matrix estimated by either OMVC-NMF or OMDC-ADMM at the last slice k of the image. For MVC, we used four different values of ǫ: 10 -9 , 10 -4 , 10 -3 and 10 -2 . The normalized curves (maximum value equal to 1, minimum value equal to 0) are shown in Figure 6. In Figure 6, several behaviors can be observed. For MVC with ǫ = 10 -9 , the plot presents a staircase shape and, as the value of ǫ increases, the steps tend to disappear. Actually, MVC with low values of ǫ is favoring rank deficiency as µ increases; this is reflected on the curves by abrupt staircase shape variations. When ǫ increases, rank deficiency is no longer promoted. In Figure 7, we compared the true and estimated simplexes for different values of µ, for MVC (ǫ = 10 -9 , 10 -3 , and 10 -2 ) and MDC. For ǫ = 10 -9 , as µ increases, the volume of the estimated simplex tends to zero, meaning that the endmember matrix is rank deficient. In the limit case (µ → +∞), all endmembers become collinear and the simplex is transformed into a point. On the other hand, as ǫ increases, the rank of the decomposition is preserved for a larger interval of µ. For example, for ǫ = 10 -2 and µ = 0.03, although the estimated simplex is included in the true simplex, the rank is preserved which avoids numerical instabilities; this shows the interest of the ǫ parameter in the context of on-line blind unmixing. Similarly, for MDC (Figure 7(d)), the simplex decreases progressively, while preserving the rank of the decomposition. We used for both OMVC-NMF and OMDC-ADMM the same coefficient α = 0.99 and regularization parameter µ = 0.003. The parameter ρ for OMDC-ADMM algorithm was set to 0.001 and ǫ for MVC to 0.001. To compare the convergence speeds of the two algorithms, the residual error vs. time sample (from 1 to 40) was evaluated for different values of N iter . The results are shown in Figure 9. Note that a logarithmic scale was used for the vertical axis. In Figure 9, we observe that for all the considered values of N iter , OMDC-ADMM has a faster convergence rate than OMVC-NMF. The latter reaches asymptotically an error close to OMDC-ADMM for approximately N iter = 1000 iterations. The residual error for OMVC-NMF and OMDC-ADMM as functions of number of samples (for the data represented in Figure 8)

In order to quantify the estimation quality of the endmembers and of the abundances for different values of N iter , we computed the SAD and the RMSE performance criteria. The SAD was calculated using the endmembers estimated at the last slice, which is coherent in the case of stationary sources. The obtained values are given in Tables 1 and2. One can see that 100 iterations are sufficient for OMDC-ADMM algorithm to converge to a relevant solution; in fact, beyond 100 iterations, the estimation accuracy of OMDC-ADMM no longer improves. OMDC-ADMM processes the entire image in about 0.3 s. OMVC-NMF algorithm requires at least 1000 iterations to estimate correctly the endmembers and the abundances, for a processing time of about 1.3 s. This slow convergence rate can be attributed to the use of multiplicative update rules in OMVC-NMF. These results show that there is a real interest in terms of rapidity and accuracy in using OMDC-ADMM for on-line blind unmixing. The values of RMSE for OMVC-NMF and OMDC-ADMM for different number of iterations (for the unmixing of the data represented in Figure 8)

OMVC-NMF OMDC-ADMM

Another interest of the sequential algorithm for hyperspectral unmixing is the processing of large hyperspectral datasets. To illustrate this aspect, we compared the performance of OMDC-ADMM algorithm with its batch counterpart BMDC-ADMM. We used the same synthetic hyperspectral image (corrupted by a low-level noise) presented in Figure 8. Both algorithms used the same ρ = 0.001; after several trials, the minimum dispersion constraint was set to µ = 0.003 for the on-line version and to µ = 0.5 for the batch version, for optimum performance. The value of α was set to 0.99. The BMDC-ADMM method was applied to the unfolded version of the hyperspectral image, of size 119 × 1600, while the on-line version processed the image slice by slice. To study the convergence speeds of both algorithms, the SAD was evaluated for different values of N iter , ranging from 10 to 2000. Note that for a given value of N iter , the SAD for OMDC-ADMM was calculated for each slice of the image and then averaged, in order to make a fair comparison between the two algorithms. The results obtained by averaging over 20 trials are illustrated in Figure 10. One can see that our on-line algorithm converges much faster than its batch counterpart. Indeed, 100 iterations are enough for OMDC-ADMM to yield accurate estimates, while BMDC-ADMM requires at least 800 iterations. This has strong consequences on the computational cost. Consider the computational complexity of the on-line and batch versions established in section 4. We also examined the influence of parameters α and ρ on the convergence speed of OMDC-ADMM, using the same synthetic dataset in the noise-free case. We first varied α from 0.99 to 0.8 while fixing the values of ρ = 0.001 and µ = 0.003. The values of SAD as functions of N iter , for different values of α, are illustrated in Figure 11. As α increases, fewer iterations are required to converge to the correct solution; this can be explained by the fact that α and µ parameters are strongly linked. Actually, if α changes, the value of µ must also be modified, in order to obtain optimum performance.

Next, we varied ρ from 0.01 to 0.0003 while keeping the value of α = 0.99, µ = 0.003 and N iter = 100. Figure 12 shows that, when ρ decreases, the asymptotic SAD decreases, while the convergence speed does not seem to be affected. As soon as the value of ρ becomes too low (e.g. 0.00003), we notice that the SAD increases, a sign that the algorithm becomes unstable. For the proper functioning of OMDC-ADMM, there is a range of suitable values for the parameter ρ; for this experiment, a value of ρ between 0.001 and 0.00005 ensures the stability of the algorithm. The last aspect studied by numerical simulations is the sensitivity of the algorithm to initial conditions. We used 20 different random initializations for OMDC-ADMM on the synthetic data corrupted by a low-level noise. The SAD and the RMSE were computed for each initialization, using N iter = 100 and µ = 0, 0.0001, 0.0005, 0.003 and 0.01. The results are shown in Figure 13 as boxplots. When µ is too small, different initializations are likely to produce different estimates of the endmembers. When the value of µ is adequately chosen, the algorithm is not sensitive to initialization. Finally, when the value of µ is too high (> 0.003), the SAD and the RMSE increase, indicating that the problem is over-regularized.

Tracking the spectral variability

In hyperspectral imaging, the spectral signature of a component may present intra-class variability [START_REF] Zare | Endmember variability in hyperspectral analysis: addressing spectral variability during spectral unmixing[END_REF][START_REF] Halimi | Unsupervised unmixing of hyperspectral images accounting for endmember variability[END_REF]. Thus, studying the evolution of the endmembers between successive samples is sometimes crucial in understanding the underlying physical phenomenon. Integrating a tracking capability (via the parameter α) allows the tracking of dynamic spectral changes. To address this point, a hyperspectral image of size 119 × 40 × 80 with nonstationary endmembers was simulated. The abundance maps were identical to those in Figure 8, but the endmembers evolved at each slice, according to the following model: S(k+1) = S(k) + T (k) , where T (k) is a low-level random noise matrix drawn from the standard normal distribution. The simulated data were corrupted by noise with an SNR = 26 dB. The parameters µ and ρ were set to 0.01 and 0.001 respectively. In order to assess the ability of our algorithm to follow the evolution of the endmembers, we varied α from 0.99 to 0.7 and, for each slice k, we computed the SAD. The results are shown in Figure 14. As α decreases, both convergence speed and asymptotic SAD (k → +∞) increase. Note that for very large values of α (e.g. 0.99), an increase of the asymptotic SAD is observed. This indicates that there exists an optimal value of α mitigating at best the transient error and tracking error.

Real data application: on-line analysis of wood surface

In this section, the performance of OMDC-ADMM is tested on two real hyperspectral images of wood, recorded by a pushbroom acquisition system (which is the target application of our method) at CRAN (Centre de Recherche en Automatique de Nancy) laboratory.

Wood is a natural material whose rendering is appreciated in many applications. However, its surface heterogeneity and variability makes it difficult to control the quality of the final product. These technical difficulties generate significant non-quality costs and manual sorting in wood industries. Current industrial technical solutions only detect very pronounced defects, such as knots, using color vision systems. There are other types of defects, which have high colorimetric variability, and thus, are not systematically detected by the existing systems; this is the case of the sapwood (material located between the bark and the heartwood) and of the dark grain (dark color wood fibers), for example. The sapwood is generally undesirable in many applications due to its porosity, while dark grain-type defects are undesirable, particularly in the parquet industry, because they degrade the aesthetic appearance of the product. Thus, the detection of these defects represents an important economical issue, in many wood industry branches. A relevant technique to overcome the limitations of color vision systems is NIR hyperspectral imaging, which is the target application of this work. For the next experiment, we considered the wood samples shown in Figure 15, measuring 100 cm × 15 cm. The sample 1 had a sapwood-type defect in the top left corner and in the bottom right corner. The sample 2 had a dark grain area in the middle, over its entire length. These defects were hardly distinguishable by the colorimetric systems and were manually labeled by experts. The two pieces of wood were scanned by a NIR pushbroom imager working in the spectral range of 900-1600 nm, resulting in hyperspectral images of dimension 224 × 264 × 1253 (for sample 1) and 224 × 200 × 1101 (for sample 2), where 224 represents the number of wavelengths and the other two numbers, the (spatial × time) dimensions. The acquired images were processed sequentially by OMDC-ADMM, slice by slice. The parameters of the algorithm were set as follows: µ = 3, α = 0.99, ρ = 0.001 and N iter = 15. The number of endmembers to be extracted was set to R = 3, after several preliminary tests.

The three abundance maps generated by OMDC-ADMM for the sample 1 and for the sample 2 are presented in Figure 16 and 17, respectively. On the first abundance map in Figures 16 and17, we observe that the areas corresponding to the sapwood and to the dark grain are much more contrasted than the rest of the image. These are very promising results, as, for the considered samples, the color of the sapwood/dark grain was very close to the color of heartwood, and therefore, the wood planks were considered as compliant with industrial quality standards by color vision systems. For the considered datasets, the processing time by OMDC-ADMM was of 3 s. The processing time can be largely reduced by low-level implementation of our algorithm directly on the production line to satisfy the real-time industrial constraint of 1 linear m of wood /s. To further assess the performance of the proposed OMDC-ADMM algorithm, we also applied it to a real dataset, called Jasper Ridge (http://lesun.weebly.com/hyperspectral-data-set.html) with available ground truth. This image is commonly used in hyperspectral unmixing community for benchmarking purposes [START_REF] Zhu | Hyperspectral unmixing: ground truth labeling, datasets, benchmark performances and survey[END_REF]. The results are given in the AppendixA.

Conclusions

We proposed a new algorithm, OMDC-ADMM, specially designed for the on-line unmixing of pushbroom hyperspectral images. Tests on simulated data have shown that this new algorithm outperforms the state-of-the-art methods based on multiplicative update rules, in terms of convergence speed and estimation accuracy. Moreover, these experiments revealed that OMDC-ADMM makes it possible to track the spectral variability of the endmembers over time and significantly reduces the processing time compared to its offline counterpart, which is a crucial feature for real-time data processing. From a methodological point of view, we proved the interest of using the minimum dispersion constraint over the minimum volume one, in particular its capacity to regularize the problem and to stabilize the solution. Finally, OMDC-ADMM proved to be an interesting solution for on-line unmixing of pushbroom hyperspectral images, compliant with the real-time constraints of the industrial wood sorting systems.

AppendixA. Test of OMDC-ADMM on Jasper Ridge hyperspectral image

We compared the performance of OMDC-ADMM to its batch version (BMDC-ADMM) and to other state-of-the-art (batch) methods: Vertex Component Analysis (VCA) [START_REF] Nascimento | Vertex Component Analysis: a fast algorithm to unmix hyperspectral data[END_REF], NMF [START_REF] Lee | Learning the parts of objects by Non-negative Matrix Factorization[END_REF], ℓ 1 -NMF [START_REF] Hoyer | Non-negative sparse coding[END_REF] and ℓ 1/2 -NMF [START_REF] Qian | Hyperspectral unmixing via ℓ 1/2 sparsity-constrained Non-negative Matrix Factorization[END_REF]. The results for the last four methods come from simulations conducted in [START_REF] Zhu | Hyperspectral unmixing: ground truth labeling, datasets, benchmark performances and survey[END_REF]. The Jasper Ridge image has a spatial size of 100 × 100 pixels. Each pixel is recorded at 198 wavelengths ranging from 380 nm to 2500 nm. There are four latent endmembers in this data, corresponding to tree, water, soil and road.

In order to respect the pushbroom acquisition scheme that we consider, the image was processed by OMDC-ADMM sequentially, by slice of dimension 198×100, and along the vertical axis. Note that the sequential processing along the horizontal axis was not possible since rank preservation was not guaranteed from one slice to another. The parameters were set as follows: R = 4, µ = 0.05, α = 0.99, ρ = 0.001 and N iter = 200. BMDC-ADMM was applied to the unfolded version of the hyperspectral image, of size 198×10000 with parameters R = 4, ρ = 0.001, µ = 200 and N iter = 2000. For the parameters settings details of the other methods, the reader is referred to [START_REF] Zhu | Hyperspectral unmixing: ground truth labeling, datasets, benchmark performances and survey[END_REF]. In order to evaluate the quality of the estimated endmembers and abundances, the SAD and the RMSE for all methods were computed. Each experiment was repeated 50 times for different random initializations and the average results for the SAD and the RMSE are provided in Table A.3.

It can be seen, from Table A.3, that BMDC-ADMM outperforms the other state-of-the-art methods. Moreover, OMDC-ADMM and BMDC-ADMM yield close estimates for soil, tree, and road. Interestingly, we observed that there is a significant difference between the two methods regarding the SAD of water (0.1186 for OMDC-ADMM vs. 0.3141 for BMDC-ADMM). This difference can be explained by the phenomenon of spectral variability, which can modify locally the spectrum of pure materials. The causes of this variability can be diverse, e.g., the changing illumination conditions during the acquisition, the intrinsic variability of the components or the atmospheric effects. This aspect is illustrated in Figure A.18 which plots the spectra for all estimated endmembers slice by slice for the on-line algorithm. Thanks to its tracking capability, OMDC-ADMM offers the possibility to study the dynamic content changes over time. We note that the spectral signature of water evolves strongly between slices. In Table A.3, the SAD for OMDC-ADMM was calculated using the average over the K time samples of the estimated spectra.

OMDC-ADMM

For the batch algorithm, this spectral variability is not explicitly taken into account, and the estimated endmembers can be interpreted as the average spectra for the entire image. The endmembers and the abundance maps estimated by OMDC-ADMM and BMDC-ADMM, along with the ground truth, are represented in Besides the fact that OMDC-ADMM makes it possible to track the spectral variability, it also requires fewer iterations than the batch version, which significantly reduces the processing time (2.77 s vs. 4.92 s).

  with X ∈ R L×P + , S = [s 1 , ..., s R ] ∈ R L×R + and A ∈ R R×P +

3 Figure 1 :

 31 Figure 1: Geometrical representation of the LMM for R = 3 endmembers

Figure 2 :

 2 Figure 2: Acquisition of hyperspectral images of wood

Figure 3 :

 3 Figure 3: Data structure for pushbroom acquisition

Figure 4 :

 4 Figure 4: Unfolded pushbroom hyperspectral image

Figure 5 :

 5 Figure 5: Illustrations of MVC and MDC penalties as functions of the angle β

Algorithm 1

 1 OMDC-ADMM Inputs: X; P; R; α; µ; ρ; N iter ; Initialization: k = 0; N = zeros(L, R); M = zeros(R, R); S = rand(L, R); Ṽ = zeros(R, P ); Ũ = zeros(L, R); Π = zeros(R, P ); Λ = zeros(L, R); A = [ ]; S = [ ];

- 1 ;Ũ

 1 = max 0, S + Λ ; Λ ← Λ + S -Ũ; t ← t + 1; end while N = Ñ; M = M; A ← A Ã ; S ← S S ; end while

Figure 6 :

 6 Figure 6: The plots of MVC and MDC cost functions for different values of µ

Figure 7 :

 7 Figure 7: Illustration of simulated and true simplexes as functions of µ for MVC and MDC

Figure 8 :

 8 Figure 8: Simulated abundances and endmembers (End)

Figure 9 :

 9 Figure 9: The residual error for OMVC-NMF and OMDC-ADMM as functions of number of samples (for the data represented in Figure 8)

Figure 10 :

 10 Figure 10: SAD for BMDC-ADMM and OMDC-ADMM as functions of N iter (for the unmixing of the data represented in Figure 8)

Figure 11 :

 11 Figure 11: SAD for OMDC-ADMM as functions of N iter for different values of α (for the unmixing of the data represented in Figure 8)

Figure 12 :

 12 Figure 12: SAD for OMDC-ADMM as functions of number of samples for different values of ρ (for the unmixing of the data represented in Figure 8)

3 RMSEFigure 13 :

 313 Figure 13: SAD and RMSE for different initializations of OMDC-ADMM as functions of µ (for the unmixing of the data represented in Figure 8)

Figure 14 :

 14 Figure 14: SAD for OMDC-ADMM as functions of number of samples for different values of α, in the case of non-stationary endmembers

Figure 15 :

 15 Figure 15: Selected wood samples (visible range color image)

Figure 16 :Figure 17 :

 1617 Figure 16: The three abundance maps estimated by OMDC-ADMM (for the wood sample 1 represented in Figure 15(a))

  Figure A.19 and A.20. For OMDC-ADMM, the average values of the endmembers are represented. It can be seen in Figure A.20, that some details are missing on the abundance maps estimated

Figure A. 20 :

 20 Figure A.20: True and estimated abundance maps for OMDC-ADMM and BMDC-ADMM (for the benchmarking image Jasper Ridge)

Table 1 :

 1 The values of SAD for OMVC-NMF and OMDC-ADMM for different number of iterations (for the unmixing of the data represented in Figure8)

	N iter = 50	0.2824	0.0816
	N iter = 100	0.2339	0.0019
	N iter = 150	0.1985	0.0019
	N iter = 1000	0.0028	0.0019
		OMVC-NMF OMDC-ADMM
	N iter = 50	0.2246	0.1312
	N iter = 100	0.2057	0.0029
	N iter = 150	0.1946	0.0028
	N iter = 1000	0.0036	0.0028

Table 2 :

 2

  Table A.3: SAD and RMSE for OMDC-ADMM, BMDC-ADMM, VCA, NMF, ℓ 1 -NMF and ℓ 1/2 -NMF (for the benchmarking image Jasper Ridge)

			BMDC-ADMM	VCA
	Endmembers SAD RMSE	SAD RMSE	SAD RMSE
	Tree	0.0765 0.0579 0.0905 0.0667 0.2565 0.3268
	Water	0.1186 0.0306 0.3141 0.0521 0.2474 0.3151
	Soil	0.1027 0.0993 0.0886 0.0893 0.3584 0.2936
	Road	0.0584 0.0791 0.0473 0.0693 0.5489 0.2829
	Average	0.0891 0.0667 0.1351 0.0694 0.3528 0.3046
		NMF	ℓ 1 -NMF	ℓ 1/2 -NMF
	Endmembers SAD RMSE	SAD RMSE	SAD RMSE
	Tree	0.2130 0.1402 0.0680 0.0636 0.0409 0.0707
	Water	0.2001 0.1106 0.3815 0.0660 0.1682 0.1031
	Soil	0.1569 0.2557 0.0898 0.2463 0.0506 0.2679
	Road	0.3522 0.2450 0.4118 0.2344 0.3670 0.2737
	Average	0.2305 0.1879 0.2378 0.1526 0.1567 0.1789

We benefit from the support of the ANR-OPTIFIN (Agence Nationale de la Recherche-OPTimisation des FINitions) project.

by BMDC-ADMM; this is not the case for the maps generated by OMDC-ADMM.


Abir Mayoufi

Stephane Victor

Rachid Malti

Mohamed Aoun

Mohamed Aoun Dynamic Harmonic

Dynamic Harmonic Regression approach for global terrestrial temperature system identification

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Environnement

INFORMATIONS COMPLEMENTAIRES

Si nécessaire

RESUME ETENDU

Deux pages maximum, en français ou en anglais

The current and future evolution of the global climate raises many questions about how it works. In this context, a terrestrial climate model was proposed by De Larminat [START_REF] De Larminat | Climate Change: Identification and Projections[END_REF] which identifies the dynamic relationship between solar activity, volcanic activity, human activity and the averaged measures of temperature. On one hand the total solar irradiance (TSI, or solar activity) has an 11 year cycle, on the other hand, the aerosol emissions (volcanic activity) and the greenhouse gas emissions (anthropogenic activity) do not have this cyclic component.

Trying to enrich knowledge of this model, this study aims to find all the cyclic components of the TSI signal and the output temperature signal.

Terrestrial climate system

Consider a MISO terrestrial climate model [START_REF] De Larminat | Climate Change: Identification and Projections[END_REF] such as the one given in Figure 1.

Fig. 1: Terrestrial climate model

The input signals are composed of the solar, volcanic and human activities. These inputs can be natural (the variations of radiation intensity (solar activity) and aerosol emissions (volcanic activity)) or anthropogenic through greenhouse gas emissions (carbon dioxide).

The climate system reacts to these forcing factors with temperature variations.

Figure 2 shows a set of input/output reconstructions considered as the connection of two recordings: satellite recording with calibrated proxies records from National Oceanic and Atmospheric Administration (NOAA).

A. Human activity (CO2 concentration)

Carbon dioxide is probably the most well-known greenhouse gas. The signal considered as human activity is considered as the connection of two records, measurement of the centre of Mauna Loa and proxy measurement, to obtain a single piece of data set.

B.

Natural forcing a. Solar activity (Total Solar Irradiance) The intensity and movement of sunspots presents an 11 year cycle, from these sunspots a reconstruction of the Total Solar Irradiance (TSI) was introduced by the TIM instrument (Total Irradiance Monitor), another proxy reconstruction was determined based on isotopic proxies. These two reconstructions were combined to obtain a single piece of data set.

b. Volcanic activity (Aerosol Optical depth)

Aerosols change the climate in several ways. The first effect is observed through the diffusion and absorption of solar radiation, which modifies the irradiative balance of the planet. The Aerosols Optical Depth (AOD) more specifically describes how aerosols affect the passage of light through the atmosphere, for a given wavelength (in this manuscript the wavelength is 550 nm).The signal of volcanic activity is considered as the connection of two recordings: satellite recording with calibrated sulphate records in Antarctic ice cores.

C. Temperature

These forcing factors are considered as the input signals of the terrestrial climate system and temperature, considered as the connection of proxy records and HadCRUT records, as the output.

Dynamic Harmonic Regression model

A Dynamic Harmonic Regression (DHR) model helps to reveal the presence of dynamically, non-stationary, quasi-cyclical components. This Unobserved Component (UC) model takes the following form: where is the signal to be analyzed, is the trend, is a cyclical or quasi-cyclical component and is an irregular where is the number of peaks determined from the Auto-regressive spectrum and and are stochastic time variable parameters (TVP). The trend component can be considered as a stochastic time variable parameter also that could be incorporated into the cyclical component as a zero-frequency term ( with set equal to zero). The DHR model is an extension of the classical harmonic regression model with constant parameter (Fourier series) where the gain and phase of the harmonic components vary according to changes in the parameters and .

The complete DHR estimation algorithm consists of the following steps: first, the dominant harmonic frequencies are extracted from the AR spectrum of the observed signal. Then, the unknown parameters are estimated from a least squares optimization function. Once these parameters have been estimated, the trend and harmonic components varying over time can be obtained using the Kalman filter approach [START_REF]Dynamic harmonic regression -Young[END_REF].

Applying this approach to the temperature output signal, a spectral analysis is carried out by using the auto-regressive spectrum (see Figure 3). The associated spectrum is characterized by peaks, which respectively correspond to 21.96, 6.08, 4.78, 4.23, 3.60, 3.17, 2.86, 2.58, and 2.36 years/cycle. This study consists to isolate the fundamental frequency characterized by the peaks at 21.95 by applying the DHR method to extract the trend and cyclical component (see Figure 4). A clearly defined peak can be noticed at 21.96 years/cycle in the estimated spectrum of the DHR model for the output signal (see Figure 3), which apparently is present in the TSI input. Therefore, a cyclical component at 21.96 years/cycle is extracted on both temperature signal and TSI signal (see Figures 4/5). For future works, the resulting reconstructed signals can be used along the other exogenous inputs for classical parameter system identification purposes. The DHR method provides an approach that helps analyzing non-stationary time-series and allows identifying all unobserved components. More specifically, a cyclic component has been characterized and identified so that the temperature output becomes freed of this specific component.

  NOM ET PRENOM DES AUTEURS MAYOUFI Abir/ VICTOR Stéphane/ MALTI Rachid/ AOUN Mohamed ETABLISSEMENT / LABORATOIRE IMS-Bordeaux, France /MACS-Gabès, Tunisie TITRE DE LA PRESENTATIONDynamic Harmonic Regression approach for global terrestrial temperature system identification BLOC DISCIPLINAIREIndiquer la discipline à laquelle vos travaux appartiennent --Maintenance -Sureté de fonctionnement -Supervision Autre (préciser ci-dessous) au sens de sa diffusion dans la société)Autre (préciser ci-dessous)

Fig. 2

 2 Fig. 2 Input/output reconstruction

Fig. 3

 3 Fig.3 pseudo-spectrum 'fitted' (red) vs. empirical spectrum (blue) Fig.4 Data, trend and cyclical component for the temperature output signal

  )

Fig. 5

 5 Fig.5 Data, trend and cyclical component