Online Canonical Polyadic Decomposition: Application of Fluorescence Tensors with Nonnegative Orthogonality and Sparse Constraint

Abstract

The canonical polyadic decomposition (CPD) is now widely used in signal processing to decompose multi-way arrays. In many applications, it is important to add constraints to quickly converge on an optimal solution. In contrast to classical CPD, we then focus on online CPD. In this context, the number of relevant factors is usually unknown and can vary with time. We propose two algorithms to compute the online CPD based on sparse dictionary learning. We also introduce an application example in environmental sciences and evaluate the performances of the proposed approaches in this context on real data.

1. Introduction

In many fields such as psychometric [1], data mining [2], neuroscience [3], chemometric [4], telecommunications [5], computer vision [6], and biomedical image processing [7], use canonical polyadic decomposition (CPD) also known as PARAllel FActor analysis (PARAFAC).

The data from these fields can be put in a multidimensional data array (or tensor). Then, the CPD permits to decompose this array into factors (which is composed of the number of components or rank depending to the application) that can be interpreted by the user.

CPD algorithms can be summarized in three mains methods:

• Direct methods. These methods are based on algebraic computation such as Direct TriLinear Decomposition (DTLD) [8], SEmi-algebraic framework for approximate CP decompositions via SImultaneous matrix diagonalization (SECSI) [9], and DIrect AlGorithm for canonical polyadic decomposition (DIAG) [10].

• Alternating methods. These methods are based on least square methods, and we can cite alternating least square (ALS) [1], hierarchical alternating least square (HALS) [11, 12, 13], and alternating direction method of multipliers (ADMM) [14].

• Descent methods. These methods are based on traditional gradient descent. We can cite some algorithms using stochastic gradient descent [15, 16, 17] or Nadam optimizer [18].

By considering these algorithms, we can say that the rank of the decomposition, namely the relevant number of factors, is known, but we find in practice that this is not the case. Indeed, this value is unknown so we make an estimate on the order of magnitude while a bad choice can have an impact on the factorial estimates [19]. Traditionally, two opposing approaches are used to address this situation.

• Rank estimation. Rank estimation consists in estimating the appropriate CPD rank before the decomposition. Among these methods, we can cite the CORe CONsistency DIAgnostic (CORCONDIA) [20], split half validation [21], or AutoTen [22]. Though these approaches do not always allow to clearly decide between several possible rank values.

• Overfactoring. In this method, an overestimated value is chosen, i.e., one higher than the current rank and the CPD algorithms are designed to produce additional factors with zero contribution [19]. Thus, in contrast, overfactoring is a posterior rank estimation method because the appropriate rank value is inferred from the CPD.

At this stage, we have talked about offline CPD. However, in this chapter, we are interested in online CPD. In online CPD, the data increases with the time, and the decomposition must follow this augmentation [23].

In the context of fluorescence data and environmental sciences, the online CPD can be summarized as follow:

• 3-way array is used. This tensor is built by concatenation of a new matrix on the last mode. At each time we got a tensor call sub-tensor as shown in Figure 1.

• In practice, some components or the rank of the decomposition can change at interval time from one sub-tensor to another. This phenomenon can be seen as appear and/or disappear of component. We talk about the rank variation but these variations are unknown.

• The goal of online CPD in this context is to update at each new time interval the CPD factors using the factors estimated previously without performing the CPD of the whole.

Figure 1.

Example of a sub-tensor.

In the literature, some online CPD algorithms exist. For example, in [23], the authors introduced two adaptive algorithms: simultaneous diagonalization tracking (SDT) which tracks the SVD of the unfolded tensor and recursive least squares tracking (RLST). In [24], grid-based tensor factorization algorithm (GridTF) is introduced for the large tensor. Also, in [25] authors proposed OnlineCP based on ALS algorithms. In [26] authors introduced an online Tucker decomposition with a fixed rank. This list is not exhaustive, and more online algorithms can be found in [27].

For the rest of the chapter, we propose an approach based on sparse dictionary learning to compute the CPD with sparsity, nonnegativity, and orthogonality constraint. These constraints permit to handling unknown rank variation over time and the better convergence. Considering the most general case, we make no assumptions about these variations. The goal of the proposed approach is to combine dictionary learning with LASSO [28] in a simple but appropriate way for nonnegative online CPD. We have selected the factors of the two fixed modes from two dictionaries which are learned and updated throughout the online process. The appropriate number of vectors is selected from a sparsity constraint on the two atoms matrices (i.e., the CPD rank). We derived from this solution three different algorithms (one for offline CPD and the rest for online CPD). They differ in the way that the dictionaries are updated from one sub-tensor to another. Next, we propose an experiment semi-controlled to obtain real fluorescence data online with rank variations. We provide out an evaluation od our approaches with the state of the art.

Notations are introduced while recalling some basic definitions and the problem of online CPD in the particular context of fluorescence spectroscopy and environmental sciences in Section 2. In Section 3, the new proposed approach is described in the form of three algorithms for computing the nonnegative CPD of a third-order tensor. Two of these algorithms were already presented in a previous conference paper [18] and this journal paper [27]. In Section 4, an experimental design is proposed to evaluate the approaches. Results and comparisons with reference approach are provided and discussed. Section 5 concludes.

2. Problem formulation

We denote tensors and sub-tensors with letter $\mathcal{T}$ and are of size $\left(I\times J\times K\right)$ meaning that $\mathcal{T}$ gather K matrix of size $\left(I\times J\right)$ on its last mode.

2.1 Offline CPD: Example of a fluorescence tensor

In environmental sciences, especially for fluorescence data tensors the third-order tensor is most used. The CPD of this tensor can be used in order to characterize fluorescent dissolved organic matters (FDOM) in natural water samples [21, 27, 29, 30, 31] (see Figure 2).

Figure 2.

Example of the rank 4 CPD of the fluorescence tensor [18, 27].

Fluorescence data tensors are a set of 2D signals called emission and excitation matrices (EEMs) of fluorescence measured from a set of samples. Each sample is then a mixture of an unknown number of fluorescent components (fluorophores), and each part of an EEM corresponds to the fluorescence intensity of one sample at a given couple of excitation and emission wavelengths. At low concentrations, the nonlinear model based on the Beer–Lambert law can be linearized, and it then follows the CPD model [32]. To recover each individual emission and excitation spectra of the fluorophores present in the different samples along with their respective contributions can be found using the CPD. In some applications like in our case, the decomposition factors have physical meaning and are known to be nonnegative. This constraint allows better convergence especially when the columns of the factors are collinear. Nonnegativity constraint can be imposed on the factors during optimization by projecting the values of each factor matrix on ${\mathbb{R}}^{+}$ [33]. This method of applying nonnegativity constraint sometimes shows its limits in certain cases due to the change of scale. In [34], the authors propose to solve this problem by involving two factors, $ϵ$ and $\theta \in \left\{-1;1\right\}$, which multiply each factors. Nonnegativity can also be imposed by choosing an optimization method [14, 17, 35] or even using an exponential variable change [31]. In the literature, the projection method is mostly used such as in approach [14, 17, 18, 33, 36, 37]. In signal or data processing, the nonnegative CPD is usually used as a mathematical model to fit a data tensor $\mathcal{T}$ of size (I, J, K):

$$\forall i,j,k,{\mathcal{T}}_{i,j,k}\simeq {\widehat{\mathcal{T}}}_{i,j,k}\left(\mathbf{A},\mathbf{B},\mathbf{C}\right)=\sum _{r=1}^R{\mathbf{A}}_{\mathit{ir}}{\mathbf{B}}_{\mathit{jr}}{\mathbf{C}}_{\mathit{kr}}$$E1

The matrices $\mathbf{A}\in {\mathbb{R}}_{+}^{I\times R},\mathbf{B}\in {\mathbb{R}}_{+}^{J\times R}$ and $\mathbf{C}\in {\mathbb{R}}_{+}^{K\times R}$ are named the factor matrices. R is the CPD rank. Each column of the factor matrices A, B, and C, defines the CPD factors. We can distinguish two values of R if the CPD factors have a physical meaning:

• The define value of the tensor $\mathcal{T}.$

• The maximal value of R for which all the factors have a physical meaning. We name this value the physical rank of $\mathcal{T}.$ It is usually much smaller than the tensor rank.

In the considered application, these matrices are usually full column rank so that the Kruskal condition [38] is fulfilled and guarantees that the decomposition is unique up to trivial scaling and permutation indeterminacy.

The euclidia distance F is the most used to the reconstruction error. With nonnegative constraint, the offline CPD problem is thus given by

$$\begin{array}{l}min\left\{F\left(\mathbf{A},\mathbf{B},\mathbf{C}\right)=\parallel \mathcal{T}-\widehat{\mathcal{T}}(\mathbf{A},\mathbf{B},\mathbf{C})\parallel {}_F^2\right\}\\ s.t.\mathbf{A}\ge 0,\mathbf{B}\ge 0,\mathbf{C}\ge 0\end{array}$$E2

2.2 Fluorescence tensor in online CPD

In the online case, tensor $\mathcal{T}$ can be seen as

• A composition of sub-tensors ${\mathcal{T}}_n$ of size $\left(I,J,K_n\right)$ acquired at various time intervals.

• Each sub-tensor can have zero matrix in common between last and the present sub-tensor. We call this case partition with nooverlapping.

• In opposite case, we talk about overlapping case. As explained on Figure 3, two consecutive sub-tensors share some matrices, i.e, in this case, ${\mathcal{T}}_n$ also contains the last $\eta$ EEMs of ${\mathcal{T}}_{n-1}.$ Sub-tensor overlapping should help the factor estimations but then the online decomposition becomes slower because we find ourselves with a larger number of sub-tensors.

• Rank variation or fluorophore appearance and disappearance can be explained by sea currents or pollution events, natural degradation…, in a natural marine environment.

Figure 3.

Example of a sub-tensor.

The physical rank of sub-tensor ${\mathcal{T}}_n$ will be denoted ${\tilde{R}}_n$ in this chapter, and we recall that it is unknown. The online CPD problem with rank variation has been clearly described in [18, 27, 39].

In short, the online CPD consists by solving Eq. (2) for each sub-tensor ${\mathcal{T}}_n$ with the nonnegative constraint. We speak of: NonNegative Online Canonical Polyadic Decomposition (NNOnline CPD). The factor matrices estimated at time $t_{n-1}$ with the sub-tensor ${\mathcal{T}}_{n-1}$ are used to update the sub-tensor ${\mathcal{T}}_n$ in order to reduce the computational cost.

3. Algorithms for computing the nonnegative online canonical polyadic decomposition with sparsity and orthogonality constraint

3.1 NNCPD for the first tensor acquired at time $t_0$

The first sub-tensor acquired at time $t_0$ is denoted ${\mathcal{T}}_0.$ At this time, we do not know the value of ${\tilde{R}}_0$ but we suppose that, $\forall n,{\tilde{R}}_n<<R.$ We aim to rewrite the matrices A and B as the product of sparse dictionaries ${\mathbf{D}}^A$ and ${\mathbf{D}}^B$ by atoms ${\mathbf{V}}^A$ and ${\mathbf{V}}^B:$

$$\mathbf{A}={\mathbf{D}}^A{\mathbf{V}}^A\text{and}\mathbf{B}={\mathbf{D}}^B{\mathbf{V}}^B$$E3

with ${\mathbf{D}}^A\left(I,R\right)$ and ${\mathbf{D}}^B\left(J,R\right)$_, respectively, and ${\mathbf{V}}^A\left(R,R\right)$ and ${\mathbf{V}}^B\left(R,R\right).$ We expect that ${\mathbf{D}}^A$ and ${\mathbf{D}}^B$ contain $\tilde{R}$ true components on their columns. The factors without physical meaning are $R-\tilde{R}.$ While ${\mathbf{V}}^A$ and ${\mathbf{V}}^B$ have $R-\tilde{R}$ null columns and that the other columns form generalized permutation matrix. For instance, for $R=3$ and ${\tilde{R}}_0=2$, we could have ideally

$$\mathbf{A}=\stackrel{{\mathbf{D}}^A}{\overbrace{\left(\begin{array}{ccc}9& 7& 2\\ 1& 4& 2\\ 2& 1& 2\end{array}\right)}}.\stackrel{{\mathbf{V}}^A}{\overbrace{\left(\begin{array}{ccc}0.8& 0& 0\\ 0& 0& 0\\ 0& 0.9& 0\end{array}\right)}}=\left(\begin{array}{ccc}7.2& 1.8& 0\\ 0.8& 1.8& 0\\ 1.6& 1.8& 0\end{array}\right)$$

We can easily estimate $\widetilde{R_0}$ by counting the number of non-null columns of A or B. A similar problem has been studied in depth by Cohen and Gillis in [40] but for the case in which the dictionary is known and the atom is a selection matrix. To have zeros columns on the atoms, we thus promote the sparsity constraints on the matrices ${\mathbf{V}}^A$ and ${\mathbf{V}}^B.$ In this purpose, we will aim at minimizing their $L_{1,1}$ norm. Regarding the sparsity constraint, it is generally used for sparse tensors or to ensure the overestimation of the rank, as mentioned above. This constraint is imposed according to the goal we are looking for, and in our case, it is the overestimation of the rank. The constraint can be applied to the matrices factors in the form of a regularization term:

$${min}_x\stackrel{\text{fonction}}{\overbrace{F\left(x\right)}}+\underset{\text{regularization}}{\underbrace{R\left(x\right)}}$$E4

F(x) corresponds here to Eq. (2) and is the attachment part. R(x) is the regularization part on the matrix a, b, and c.

$L_{1,1}$ norm is a classically used to promote the sparsity [41]. Indeed, it is more tractable than the $L_{0,0}$ norm. It has been used in similar contexts in [17, 35, 36, 40]. We can thus rewrite our problem as a LASSO regression problem by adding two penalty terms to the cost function (see Appendix A).

For $\mathcal{T}={\mathcal{T}}_0,$ we solve

$$\begin{array}{l}\begin{array}{l}min\left\{F_1\left({\mathbf{D}}^A,{\mathbf{V}}^A,{\mathbf{D}}^B,{\mathbf{V}}^B,\mathbf{C}\right)\right\}s.t.{\mathbf{D}}^A,{\mathbf{V}}^A,{\mathbf{D}}^B,{\mathbf{V}}^B,\mathbf{C}\ge 0\\ \text{where}{F}_1=\frac{1}{2}\parallel \mathcal{T}-\widehat{\mathcal{T}}\left({\mathbf{D}}^A{\mathbf{V}}^A,{\mathbf{D}}^B{\mathbf{V}}^B,\mathbf{C}\right){\parallel }_F^2+\alpha \parallel {\mathbf{V}}^A{\parallel }_{1,1}+\alpha \parallel {\mathbf{V}}^B{\parallel }_{1,1}\end{array}\end{array}$$E5

$\alpha >0$ is a penalty term.

This approach is named sparse nonnegative CPD (SNNCPD). We resort to a stochastic gradient descent (SGD) algorithm called Nadam [42, 43] in order to solve Eq. (5). Details of the Nadam algorithm are given in Appendix B. In order to ensure the nonnegativity of matrix entries, all the element of matrices ${\mathbf{D}}^A,{\mathbf{D}}^B,\mathbf{C},{\mathbf{V}}^A$, and ${\mathbf{V}}^B$ are projected on ${\mathbb{R}}_{+}$ at each iteration [44]. We can thus differentiate the $L_{1,1}$ norm, and the gradients with respect to the different variables are given by

$$\begin{array}{l}\frac{\partial F_1}{\partial {\mathbf{V}}^A}=-{\left({\mathbf{D}}^A\right)}^{\top }\left({\mathbf{T}}_1-{\mathbf{D}}^A{\mathbf{V}}^A{\left({\mathbf{L}}_1\right)}^{\top }\right){\mathbf{L}}_1+\alpha {\mathbf{1}}_{R,R}\\ \frac{\partial F_1}{\partial {\mathbf{V}}^B}=-{\left({\mathbf{D}}^B\right)}^{\top }\left({\mathbf{T}}_2-{\mathbf{D}}^B{\mathbf{V}}^B{\left({\mathbf{L}}_2\right)}^{\top }\right)\left({\mathbf{L}}_2\right)+\alpha {\mathbf{1}}_{R,R}\\ \frac{\partial F_1}{\partial {\mathbf{D}}^A}=-\left({\mathbf{T}}_1-{\mathbf{D}}^A{\mathbf{V}}^A{\left({\mathbf{L}}_1\right)}^{\top }\right){\mathbf{L}}_1{\left({\mathbf{V}}^A\right)}^{\top }\\ \frac{\partial F_1}{\partial {\mathbf{D}}^B}=-\left({\mathbf{T}}_2-{\mathbf{D}}^B{\mathbf{V}}^B{\left({\mathbf{L}}_2\right)}^{\top }\right)\left({\mathbf{L}}_2\right){\left({\mathbf{V}}^B\right)}^{\top }\\ \frac{\partial F_1}{\partial \mathbf{C}}=-\left({\mathbf{T}}_3-\mathbf{C}{\left({\mathbf{L}}_3\right)}^{\top }\right)\left({\mathbf{L}}_3\right)\end{array}$$E6

Matrices ${\mathbf{T}}_1,{\mathbf{T}}_2,{\mathbf{T}}_3$ are obtained by unfolding the tensor $\mathcal{T}$ with respect to the first, second, and third modes, respectively. ${\mathbf{L}}_1=\mathbf{C}\odot \left({\mathbf{D}}^B{\mathbf{V}}^B\right),{\mathbf{L}}_2=\mathbf{C}\odot \left({\mathbf{D}}^A{\mathbf{V}}^A\right)$ and ${\mathbf{L}}_3=\left({\mathbf{D}}^B{\mathbf{V}}^B\right)\odot \left({\mathbf{D}}^A{\mathbf{V}}^A\right).$ The matrices ${\mathbf{V}}^A$ and ${\mathbf{V}}^B$ are initialized as the identity matrix. The matrices ${\mathbf{D}}^A$ and ${\mathbf{D}}^B$ and C are initialized with nonnegative random values. The different steps of SNNCPD are summarized in Algorithm 1.

Algorithm 1. Sparse NonNegative Canonical Polyadic Decomposition (SNNCPD). Initialization step

• Input: $\mathcal{T},$ R overestimated

 while a convergence criterion is not reached do.

 Update the dictionaries ${\mathbf{D}}^A,{\mathbf{D}}^B$, and C using Nadam optimizer.

 Update the atoms ${\mathbf{V}}^A$ and ${\mathbf{V}}^B$ using Nadam optimizer.

end while

• Output: (${\mathbf{D}}_a,{\mathbf{D}}_b,\mathbf{C},{\mathbf{V}}_a,{\mathbf{V}}_b$)

3.2 Online step: NN-CPD for the first tensor acquired at time $t_n$ with $n>0$

In this step, we assume that the factors matrices ${\mathbf{A}}_n$ and ${\mathbf{B}}_n$ are well estimated.

3.2.1 Algorithm 1

In this approach, we allow some modification of the factors estimated at time $t_{n-1}$ thanks to linear combinations as explained below. We first replace all the null columns of ${\mathbf{A}}_{n-1}$ and ${\mathbf{B}}_{n-1}$ by columns of random numbers drawn from the standard normal distribution. We then look for ${\mathbf{A}}_n$ and ${\mathbf{B}}_n$ as

$${\mathbf{A}}_n={\mathbf{U}}^A{\mathbf{A}}_{n-1}{\mathbf{V}}^A\text{and}{\mathbf{B}}_n={\mathbf{U}}^B{\mathbf{B}}_{n-1}{\mathbf{V}}^B$$E7

The matrices ${\mathbf{U}}^A$ and ${\mathbf{U}}^B$ are square with sizes I and J, respectively. Here, ${\mathbf{U}}^A{\mathbf{A}}_{n-1}$ and ${\mathbf{U}}^B{\mathbf{B}}_{n-1}$ are our dictionaries, while ${\mathbf{V}}^A$ and ${\mathbf{V}}^B$ matrices are still sparse atoms. Therefore, the optimization problem becomes

$$\begin{array}{l}min\left\{F_2\left({\mathbf{U}}^A,{\mathbf{V}}^A,{\mathbf{U}}^B,{\mathbf{V}}^B,\mathbf{C}\right)\right\}s.t{\mathbf{U}}^A,{\mathbf{U}}^B,\mathbf{C},{\mathbf{V}}^B,{\mathbf{V}}^B\ge 0\\ \text{with}{F}_2=\frac{1}{2}\parallel \mathcal{T}-\widehat{\mathcal{T}}\left({\mathbf{U}}^A{\mathbf{A}}_{n-1}{\mathbf{V}}^B,{\mathbf{U}}^B{\mathbf{B}}_{n-1}{\mathbf{V}}^B,\mathbf{C}\right){\parallel }_F^2\\ +\alpha \parallel {\mathbf{V}}^A{\parallel }_1+\alpha \parallel {\mathbf{V}}^B{\parallel }_1\end{array}$$E8

We call this algorithm online sparse and nonnegative CPD (OSNCPD). The different steps are.

summarized in Algorithm 2. Most details can be found in [18, 27].

Algorithm 2. Online Sparse and NonNegative CPD

• STEP 1: Initialization phase

Input: ${\mathcal{T}}_0,$ R overestimated.

Solve Eq. (5) with $\mathcal{T}={\mathcal{T}}_0$.

Compute matrices ${\mathbf{A}}_0$ and ${\mathbf{B}}_0$ using Eq. (3).

Output: ${\mathbf{A}}_0,{\mathbf{B}}_0$ and ${\mathbf{C}}_0$

• STEP 2: Online phase at time $t_n$ with $n>0$

Input: ${\mathcal{T}}_n,$ and ${\mathbf{A}}_{n-1},{\mathbf{B}}_{n-1}$.

Fill null columns of ${\mathbf{A}}_{n-1},{\mathbf{B}}_{n-1}$ by random numbers and Solve Eq. (8).

Update matrices ${\mathbf{A}}_n$ and ${\mathbf{B}}_n$ using Eq. (7).

Output: (${\mathbf{A}}_n,{\mathbf{B}}_n,{\mathbf{C}}_n$ and ${\tilde{R}}_n$)

• STEP 3: Return to Step 2 with n = n + 1

The advantage of this approach is therefore its flexibility. The inconvenient is that the matrices to optimize are larger which make the algorithm more expensive in time. We call this algorithm online sparse and nonnegative CPD (OSNCPD). The different stages are summarized in the Algorithm 2.

3.2.2 Algorithm 2

By taking the Eq. (7) as well as the cost function 8, the optimization of the algorithm becomes costly because of the dictionaries ${\mathbf{U}}^A$ and ${\mathbf{U}}^B.$ This is caused by the fact that no constraint is applied to these matrices to restrict the research space. In this approach, we will apply an orthogonality constraint on these dictionaries. The orthogonality constraint on the factors is also useful in certain applications [26, 45]. It is possible to apply the orthogonality constraint to the factors by using the Eq. 4. If, for example, A must be orthogonal then: ${\mathbf{I}}_A={\mathbf{A}}^t\mathbf{A}.$ We therefore seek to minimize

$$\underset{\mathbf{A}\mathbf{,}\mathbf{B}\mathbf{,}\mathbf{C}}{min}{‖\mathcal{T}-\widehat{\mathcal{T}}\left(\mathbf{A},\mathbf{B},\mathbf{C}\right)‖}_F^2+{‖{\mathbf{A}}^T\mathbf{A}-{\mathbf{I}}_A‖}_F+{‖{\mathbf{B}}^T\mathbf{B}-{\mathbf{I}}_B‖}_F+{‖{\mathbf{C}}^T\mathbf{C}-{\mathbf{I}}_C‖}_F$$E9

To impose orthogonality constraint on the factors matrices for $\mathcal{T}={\mathcal{T}}_N,$ the cost function becomes

$$\begin{array}{l}min\left\{E_2\left({\mathbf{U}}_n^A,{\mathbf{V}}_n^A,{\mathbf{U}}_n^B,{\mathbf{V}}_n^B,\mathbf{C}\right)\right\}\\ \mathrm{o}\grave{\mathrm{u}}{E}_2=\frac{1}{2}\parallel \mathcal{T}-\widehat{\mathcal{T}}\left({\mathbf{U}}_n^A{\mathbf{A}}_{n-1}{\mathbf{V}}_n^B,{\mathbf{U}}_n^B{\mathbf{B}}_{n-1}{\mathbf{V}}^B,\mathbf{C}\right){\parallel }_F^2+\alpha \parallel {\mathbf{V}}_n^A{\parallel }_{1,1}+\alpha \parallel {\mathbf{V}}_n^B{\parallel }_{1,1}+\\ \parallel {\mathbf{I}}_A-{\left({\mathbf{U}}_n^A\right)}^T{\mathbf{U}}_n^A{\parallel }_F^2+\parallel {\mathbf{I}}_B-{\left({\mathbf{U}}_n^B\right)}^T{\mathbf{U}}_n^B{\parallel }_F^2\end{array}$$E10

The gradient of the cost Eq. (10) is given by

$$\begin{array}{l}\frac{\partial E_2}{\partial {\mathbf{V}}_n^A}=-{\left({\mathbf{U}}_n^A{\mathbf{A}}_{n-1}\right)}^{\top }\left({\mathbf{T}}^A-{\mathbf{U}}_n^A{\mathbf{A}}_{n-1}{\mathbf{V}}_n^A{\left({\mathbf{Z}}^A\right)}^{\top }\right){\mathbf{Z}}^A+\alpha {\mathbf{1}}_R\\ \frac{\partial E_2}{\partial {\mathbf{V}}_n^B}=-{\left({\mathbf{U}}_n^B{\mathbf{B}}_{n-1}\right)}^{\top }\left({\mathbf{T}}^B-{\mathbf{U}}_n^B{\mathbf{B}}_{n-1}{\mathbf{V}}_n^B{\left({\mathbf{Z}}^B\right)}^{\top }\right){\mathbf{Z}}^B+\alpha {\mathbf{1}}_R\\ \frac{\partial E_2}{\partial {\mathbf{U}}_n^A}=-\left({\mathbf{T}}^A-{\mathbf{U}}_n^A{\mathbf{A}}_{n-1}{\mathbf{V}}_n^B{\left({\mathbf{Z}}^A\right)}^{\top }\right){\mathbf{Z}}^A{\left({\mathbf{V}}_n^A\right)}^{\top }{\left({\mathbf{A}}_{n-1}\right)}^{\top }-\left(2{\mathbf{U}}_n^A\left({\mathbf{I}}_A^{\top }-{\left({\mathbf{U}}_n^A\right)}^{\top }{\mathbf{U}}_n^A\right)\right)+\\ \left(2{\mathbf{U}}_n^A\left({\mathbf{I}}_A{\mathbf{U}}_n^A-{\left({\mathbf{U}}_n^A\right)}^{\top }{\mathbf{U}}_n^A\right)\right)\end{array}$$

$$\begin{array}{c}\frac{\partial E_2}{\partial {\mathbf{U}}_n^B}=-\left({\mathbf{T}}^B-{\mathbf{U}}_n^B{\mathbf{B}}_{n-1}{\mathbf{V}}^B{\left({\mathbf{Z}}^B\right)}^{\top }\right){\mathbf{Z}}^B{\left({\mathbf{V}}_n^B\right)}^{\top }{\left({\mathbf{B}}_{n-1}\right)}^{\top }-\left(2{\mathbf{U}}_n^B\left({\mathbf{I}}_B^{\top }-{\left({\mathbf{U}}_n^B\right)}^{\top }{\mathbf{U}}_n^B\right)\right)+\\ \left(2{\mathbf{U}}_n^B\left({\mathbf{I}}_B{\mathbf{U}}_n^B-{\left({\mathbf{U}}_n^B\right)}^{\top }{\mathbf{U}}_n^B\right)\right)\\ \frac{\partial E_2}{\partial \mathbf{C}}=-\left({\mathbf{T}}^C-\mathbf{C}{\left({\mathbf{Z}}^C\right)}^{\top }\right){\mathbf{Z}}^C\end{array}$$E11

The advantage of this approach is its speed of convergence because we reduce the research space for our dictionaries. For example, if ${\mathbf{A}}_n$ is similar to ${\mathbf{A}}_{n-1},$ the matrix ${\mathbf{U}}_n^A$ is a diagonal matrix. We call this algorithm orthogonality online sparse CPD (OOSCPD). The different steps are similar with the Algorithm 2.

4. Experiment for fluorescence real data acquisition

We plan an experiment to obtain real data acquisition. Four well-known fluorophores are injected on a reservoir at different time intervals under quasi-real conditions. At all, we get 50 matrices or EEMs os size [21, 36] online with a spectrofluorimeter Hitachi F7000. The corresponding fluorescence tensor was partitioned into four successive sub-tensors (${\mathcal{T}}_0\dots {\mathcal{T}}_3$). Two kinds of partitions are considered. In the first partition (50% overlapping), each sub-tensor contains 20 EEMs, and two consecutive sub-tensors have ten EEMs in common. In the second partition (no overlapping, see Figure 3), sub-tensors have no EEMs in common. The initialization tensor ${\mathcal{T}}_0$ contains 20 EEMs. The other sub-tensors contain ten EEMs. During the acquisition, an extra fluorophore appears in all the EEMs due to the experimental device. A preliminary study showed that this extra fluorophore can be represented by a single additional factor in the CPD model. Thus, the rank of the sub-tensors can be 3, 4, or 5. We recall that SNNCPD is used for sub-tensors ${\mathcal{T}}_0$ in initialization phase.

4.1 Results and discussions

OSNCPD has already been compared with others online algorithms in [18, 27]. We compare the values of the physical ranks estimated by the two online approaches (OSNCPD, OOSCPD) with NN-CPD algorithm proposed in [31] and the actual values. Results are reported in Tables 1 and 2 for the 50% overlapping case and no overlapping case, respectively. The penalty coefficient term ($\alpha$) used in our two algorithms and in NN-CPD is also given. We can see that NN-CPD overestimated the rank in the last sub-tensors in both cases and this means that NN-CPD creates duplicate factors. In contrast, our approaches (OSNCPD and OOSCPD) find the correct rank.

Symbols	Definition
$\mathcal{T},{\mathcal{T}}_n$	Tensor, sub-tensor at time $t_n$
A, ${\mathbf{A}}^{⊺}$, a, a	Matrix, transposed matrix, column vector, scalar
${\mathbf{1}}_R$	A matrix (R, R) of 1
${\mathbf{I}}_A$	Identity matrix with size of matrix A
$\mathrm{ℝ}$	Set of real numbers
${‖.‖}_F,\parallel .{\parallel }_{1,1}$	Frobenius norm, $L_{1,1}$ norm
$\mathbf{A}\ge 0$	Means that all the elements of matrix A are nonnegative
${\tilde{R}}_n,$ $R_n$	Physical rank of the sub-tensor ${\mathcal{T}}_n,$ CPD rank of the sub-tensor
$\odot$	Khatri-Rao product

Table 1.

Main notations used in the paper.

Sub-tensor	0	1	2	3	$\alpha$
True rank	3	4	5	5	—
NN-CPD	3	4	7	10	0.5
OSNCPD & OOSCPD	3	4	5	5	1

Table 2.

Rank estimation in the overlapping case.

When the null columns of $\widehat{\mathbf{A}}$ and $\widehat{\mathbf{B}}$ and the corresponding columns of $\widehat{\mathbf{C}}$ are removed, we can compare the factor matrices estimated by our algorithms with the true factors by means of normalized root mean squared errors:

$$E_A=\frac{\parallel \mathbf{A}-\widehat{\mathbf{A}}{\parallel }_F}{\parallel \mathbf{A}{\parallel }_F},E_B=\frac{\parallel \mathbf{B}-\widehat{\mathbf{B}}{\parallel }_F}{\parallel \mathbf{B}{\parallel }_F}\text{and}{E}_C=\frac{\parallel \mathbf{C}-\widehat{\mathbf{C}}{\parallel }_F}{\parallel \mathbf{C}{\parallel }_F}$$E12

Results are reported in Tables 3 and 4 for the 50% overlapping case and no overlapping case, respectively Table 5.

Sub-tensor	0	1	2	3	$\alpha$
True rank	3	4	5	5	—
NN-CPD	3	4	5	10	0.1
OSNCPD & OOSCPD	3	4	5	5	0.5

Table 3.

Rank estimation in the no overlapping case.

Algorithm 1 (OSNCPD)					Algorithm 2 (OOSCPD)
Sub-tensor	0	1	2	3	0	1	2	3
Mean $E_A$	0.17	0.19	0.19	0.2	0.17	0.18	0.2	0.16
Mean $E_B$	0.14	0.15	0.15	0.16	0.14	0.16	0.14	0.14
Mean $E_C$	0.25	0.25	0.16	0.20	0.25	0.21	0.16	0.18

Table 4.

Mean estimations errors in the overlapping case.

Algorithm 1 (OSNCPD)					Algorithm 2 (OOSCPD)
Sub-tensor	0	1	2	3	0	1	2	3
Mean $E_A$	0.17	0.19	0.20	0.20	0.17	0.18	0.21	0.21
Mean $E_B$	0.14	0.15	0.16	0.19	0.14	0.16	0.19	0.23
Mean $E_C$	0.25	0.32	0.31	0.32	0.25	0.31	0.32	0.25

Table 5.

Mean estimations errors in the no overlapping case.

In addition, to give some physical meaning to these error terms that may seem high, we have plotted in Figure 4 the factors obtained from our approach in the sub-tensor ${\mathcal{T}}_4$ along with the true factors. We observe a good agreement between the true and estimated factors.

Figure 4.

Emission and excitation and spectra of the fluorophores for ${\mathcal{T}}_4.$ Top: Emission spectra, bottom: Excitation spectra. Red dots: estimated spectra, black lines: true spectra.

In the overlapping case, both algorithms give similar results for sub-tensors, and this can be explained that the algorithms use past information. Also, these are close to those obtained from the initialization sub-tensor ${\mathcal{T}}_0$ meaning that the online phase has been performed correctly.

In the no overlapping case OSNCPD and OOSCPD give good results although no sub-tensors are share any information. The error in matrix C seems to be high but in the no overlapping case, because the algorithms do not have the a priori information.

5. Conclusions

We have introduced two algorithms (OSNCPD and OOSCPD) for the online nonnegative canonical polyadic decomposition of sub-tensors. We have also introduced an offline nonnegative algorithm with sparsity constraint. Our algorithms are based on dictionary learning and can deal with unknown rank variations. OOSCPD incorporates orthogonality constraint and can be seen as an extension of OSNNCPD which incorporates nonnegativity and sparsity constraint for a better convergence.

These algorithms are presented in the particular case of the nonnegative online CPD of third-order fluorescence tensors, but they are not limited to this application field, and they can be easily extended to higher-order tensors. A real online fluorescence spectroscopy experiment was conducted in laboratory to validate our approach and compare with state-of-the-art approaches. The proposed algorithms allow to correctly follow the rank variations in most of the considered situations contrary to reference approaches. Eventually, these encouraging results allow to plan, for example, a monitoring in the natural environment in future work.

Conflict of interest

“The authors declare no conflict of interest.”

Appendix A

It is useful to mention that the CPD can incorporate different constraints depending on the field from which the data comes from. For example, it can be constraints of nonnegativity, orthogonality, or sparsity to name only the most commonly used.

Nonnegativity constraint Nonnegativity constraint is imposed in the CPD when the data to be processed is linked to physical quantities. In this sense, we can cite the domain of imagery [46] where nonnegativity is present because we know that optical signals are positive or zero. This is also the case of fluorescence spectroscopy where the data from the sensor is nonnegative [4, 30, 31, 47]. Nonnegativity constraint can be imposed on the factors during optimization by projecting the values of each factor matrix on ${\mathbb{R}}^{+}$ [33]. This method of imposing nonnegativity sometimes shows its limits in certain cases because of the change of scale. In [34], the authors propose to remedy this problem by involving two factors, $ϵ$ and $\theta \in \left\{-1;1\right\}$ which multiply every factor matrix. Nonnégativity can also be imposed by choosing an optimization method [14, 17, 35] or even using an exponential variable change [48] or square [31].

Sparsity constraint Regarding the constraint of parsimony, it is generally used for parsimonious tensors or to ensure the overestimation of the rank, as mentioned above. This constraint is imposed according to the goal we are looking for, and in our case, it is the overestimation of the rank. The constraint can be applied to the matrices factors in the form of a regularization term:

$$\begin{array}{c}{min}_x\underset{\text{function}}{\underbrace{f\left(x\right)}}+\underset{\text{regularization}}{\underbrace{r\left(x\right)}}\end{array}$$E13

F(x) corresponds here to ${‖\mathcal{T}-\widehat{\mathcal{T}}\left(\mathbf{a},\mathbf{b},\mathbf{c}\right)‖}_F^2$ and is the attachment part. R(x) is the regularization part on the matrix a, b, and c.

In the literature, the standard $l_{1,1}$ is the most used as in [17, 35, 41, 46, 49], and it can be defined according to the Eq. (A1) as following:

$$\underset{\mathbf{a}\mathbf{,}\mathbf{b}\mathbf{,}\mathbf{c}}{min}{‖\mathcal{T}-\widehat{\mathcal{T}}\left(\mathbf{a},\mathbf{b},\mathbf{c}\right)‖}_F^2+{‖\mathbf{a}‖}_{1,1}+{‖\mathbf{B}‖}_{1,1}+{‖\mathbf{C}‖}_{1,1}$$E14

Standard $l_{1,1}$ is conventionally used compared to standard $l_{0,0}$ to promote parsimony because it is easier to handle than standard $l_{0,0}.$ In order to minimize the rank of a matrix, we can speak of the nuclear standard which makes it possible to impose a small rank [50, 51, 52] and can be defined as a function of the Eq. (A1) as follows:

$${\mathit{min}}_{\mathbf{a}\mathbf{,}\mathbf{b}\mathbf{,}\mathbf{c}}{‖\mathcal{T}-\widehat{\mathcal{T}}\left(\mathbf{a},\mathbf{b},\mathbf{c}\right)‖}_F^2+{‖\mathbf{a}‖}_{\ast }+{‖\mathbf{B}‖}_{\ast }+{‖\mathbf{C}‖}_{\ast }$$E15

In [53], the authors show that the mixed standard provides a better convex envelope of the rank function than the nuclear standard. It therefore makes it possible to minimize the rank of a matrix. In the literature, the mixed standard is increasingly used as in [54, 55, 56]. It can be defined as a function of the Eq. (A1) as follows:

$$\underset{\mathbf{a}\mathbf{,}\mathbf{b}\mathbf{,}\mathbf{c}}{min}{‖\mathcal{T}-\widehat{\mathcal{T}}\left(\mathbf{a},\mathbf{b},\mathbf{c}\right)‖}_F^2+{‖\mathbf{a}‖}_{2,1}+{‖\mathbf{B}‖}_{2,1}+{‖\mathbf{C}‖}_{2,1}+{‖\mathbf{a}‖}_{1,2}+{‖\mathbf{b}‖}_{1,2}+{‖\mathbf{c}‖}_{1,2}$$E16

Appendix B: Stochastique gradient descent (SGD)

Take a differentiable function f(x). For each iteration k, the update of x in gradient descent is given by

$$\mathbf{x}\left(k+1\right)=\mathbf{x}\left(k\right)-\gamma {\mathbf{f}}^{\text{'}}\left(\mathbf{x}\left(k\right)\right)$$E17

with $\gamma >0,$ the step-size or the learning rate.

B.1 Adaptive moment estimation (ADAM)

Adaptive Moment estimation (ADAM) [43]. Adam optimizer is a variant of SGD. It acts on the gradient component using the exponential mobile average of the gradients and the step-size component by dividing the step-size by the square root of v, the square of exponential average gradients.

Let take a function f(x). For each iteration k, the update of x in gradient descent is given by

$$\begin{array}{l}{\mathbf{g}}_k={\mathbf{f}}^{\prime }\left(\mathbf{x}\left(k\right)\right)\\ {\mathbf{m}}_k={\beta }_1{\mathbf{m}}_{k-1}+\left(1+{\beta }_1\right){\mathbf{g}}_k\\ {\mathbf{v}}_k={\beta }_2{\mathbf{v}}_{k-1}+\left(1+{\beta }_2\right){\mathbf{g}}_k^2\\ {\widehat{\mathbf{m}}}_k=\frac{{\mathbf{m}}_k}{1-{\left({\beta }_1\right)}^k}\\ {\widehat{\mathbf{v}}}_k=\frac{{\mathbf{v}}_k}{1-{\left({\beta }_2\right)}^k}\\ \mathbf{x}\left(k+1\right)=\mathbf{x}\left(k\right)-\frac{\gamma }{\sqrt{{\widehat{\mathbf{v}}}_k+ϵ}}{\widehat{\mathbf{m}}}_k\end{array}$$E18

with $\gamma >0$ the step-size, ${\beta }_1\in \left[0,1\right],{\beta }_2\in \left[0,1\right]$ and m, v is initialized to 0.

B.2 Nesterov accelerated adaptive moment estimation (Nadam)

Nadam optimizer is an extension of the Adam optimization algorithm. The algorithm was described in the 2016 article by Timothy Dozat [42]. For each iteration k of a function f(x), the update of x in gradient descent is given by:

$$\begin{array}{l}\begin{array}{l}{\mathbf{g}}_k={\mathbf{f}}^{\prime }\left(\mathbf{x}\left(k-1\right)\right)\\ {\mathbf{m}}_k={\beta }_1{\mathbf{m}}_{k-1}+\left(1-{\beta }_1\right){\mathbf{g}}_k\\ {\mathbf{n}}_k={\beta }_2{\mathbf{n}}_{k-1}+\left(1-{\beta }_2\right){\mathbf{g}}_k^2\\ {\widehat{\mathbf{m}}}_k=\frac{{\beta }_1{\mathbf{m}}_k}{1-{\left({\beta }_1\right)}^k}+\frac{\left(1-{\beta }_1\right){\mathbf{g}}_k}{1-{\left({\beta }_1\right)}^k}\\ {\widehat{\mathbf{n}}}_k=\frac{{\beta }_2{\mathbf{n}}_k}{1-{\beta }_2}\\ \mathbf{x}\left(k+1\right)=\mathbf{x}\left(k\right)-\frac{\gamma {\widehat{\mathbf{m}}}_k}{\sqrt{{\widehat{\mathbf{n}}}_k+ϵ}}\end{array}\end{array}$$E19

Abbreviations

ALS	alternating least square
ADMM	alternating direction method of multipliers
CPD	canonical polyadique decomposition
NNCPD	nonnegative canonical polyadique decomposition
ELS	enhenced line search
PARAFAC	parallel factor analysis
DOM	dissolve organisme matter
EEM	emission and excitation matrices
pH	potential Hydrogen
NMF	nonnegative matrix factorization
SGD	stochastique gradient descent
Nadam	nesterov accelerated adaptive moment estimation
F	fluorescence
TRP	tryptophane
5S8HQ	5-sulfate-8-hydroxyquinine
RHO	RHOdamine
SVD	singular value decomposition