Marouane Nazih
email: marouane.nazih1@gmail.com

Khalid Minaoui
email: khalid.minaoui@um5.ac.ma

Pierre Comon
email: pierre.comon@grenoble-inp.fr

Using the Proximal Gradient and the Accelerated Proximal Gradient as a Canonical Polyadic tensor decomposition Algorithms in difficult situations

Keywords: CP decomposition, Tensor, Proximal gradient, Accelerated proximal gradient, Optimization

Canonical Polyadic (CP) tensor decomposition is useful in many real-world applications due to its uniqueness, and the ease of interpretation of its factor matrices. This work addresses the problem of calculating the CP decomposition of tensors in difficult cases where the factor matrices in one or all modes are almost collinear -i.e. bottleneck or swamp problems arise. This is done by introducing a constraint on the coherences of the factor matrices that ensures the existence of a best low-rank approximation, which makes it possible to estimate these highly correlated factors. Two new algorithms optimizing the CP decomposition based on proximal methods are proposed. Simulation results are provided and demonstrate the good behaviour of these algorithms, as well as a better compromise between accuracy and convergence speed than other algorithms in the literature.

INTRODUCTION

In a wide range of applications, it is necessary to handle quantities with several indices. These quantities are often referred to as tensors. The definition of a tensor will generally depend on the scientific field in which it is used.

Generally, a tensor is treated as a mathematical object that owns the property of multi-linearity when changing the coordinate system [START_REF] Comon | Tensor decompositionsstate of the art and applications[END_REF]. We consider that a tensor of order N represents a multidimensional array in which every element is addressed via N indices. For example, a scalar is a tensor of order 0, a vector is a tensor of order 1, a matrix is a second-order tensor and a cube of values is a third-order tensor, etc. In this paper, we focus on tensors of order higher than two since they possess properties which are not valid for matrices and vectors.

Among tensor decompositions, we shall be mainly interested in the so-called Canonical Polyadic (CP) decomposition [START_REF] Hitchcock | The expression of a tensor or a polyadic as a sum of products[END_REF], rediscovered forty years and named Parafac [START_REF] Harshman | Foundations of the parafac procedure: Models and conditions for an" explanatory" multimodal factor analysis[END_REF][START_REF] Harshman | Determination and proof of minimum uniqueness conditions for parafac1[END_REF] or Candecomp [START_REF] Carroll | Analysis of individual differences in multidimensional scaling via an n-way generalization of eckart-young decomposition[END_REF]. As pointed out in [START_REF] Kiers | Towards a standardized notation and terminology in multiway analysis[END_REF][START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF], the acronym "CP" decomposition can smartly stand for either "Canonical Polyadic" or "CanDecomp/Parafac", and we shall follow this terminology. The CP decomposition has been used in various fields, such as Chemometrics [START_REF] Bro | Parafac. tutorial and applications[END_REF] [START_REF] Murphy | Fluorescence spectroscopy and multi-way techniques[END_REF], Telecommunications [START_REF] Jutten | Blind separation of sources, part i: An adaptive algorithm based on neuromimetic architecture[END_REF][START_REF] Rouijel | Cp decomposition approach to blind separation for ds-cdma system using a new performance index[END_REF][START_REF] Sidiropoulos | Blind parafac receivers for ds-cdma systems[END_REF][START_REF] Sahnoun | Joint source estimation and localization[END_REF] and in Denoising of Hyperspectral Images [START_REF] Liu | Denoising of hyperspectral images using the parafac model and statistical performance analysis[END_REF]. The most interesting feature of the CP decomposition is its essential uniqueness for orders strictly larger than two [START_REF]Three-way arrays: rank and uniqueness of trilinear decompositions, with application to arithmetic complexity and statistics[END_REF] [START_REF] Stegeman | On kruskals uniqueness condition for the candecomp/parafac decomposition[END_REF], which permits parameter identification.

Another decomposition of interest is Tucker3 [START_REF] Tucker | Some mathematical notes on three-mode factor analysis[END_REF], and in particular its implementation as a High-Order Singular Value decomposition (HOSVD), where changes of bases are orthogonal. Both HOSVD and CP decomposition reduce to the usual Singular Value decomposition (SVD) in the case of matrices [START_REF] Golub | Cf van loan, matrix computations[END_REF], which are second order tensors. HOSVD is mainly used for the purposes of compression since it is not unique.

There are many situations where CP decomposition algorithms may suffer from the absence or slowness of convergence, and which can be due to degeneracies of the tensor. These situations have been well classified by Richard Harshman [START_REF] Kruskal | How 3-mfa data can cause degenerate parafac solutions, among other relationships[END_REF] in the following three cases:

• Bottleneck: A bottleneck arises when two or more factors in one of the modes are almost collinear [START_REF] Rajih | Enhanced line search: A novel method to accelerate parafac[END_REF].

• Swamp: The swamp phenomenon occurs when all modes have at least two quasi-collinear factors [START_REF] Rajih | Enhanced line search: A novel method to accelerate parafac[END_REF] [21] [START_REF] Rayens | Two-factor degeneracies and a stabilization of parafac[END_REF], which can be seen as a general case of bottleneck.

• CP-degeneracies: CP-degeneracies can be considered as special case of swamps phenomenon, where some of the factors diverge and, at the same time, tend to cancel each other, leading to a better quality of the fit progresses [START_REF] Harshman | The problem and nature of degenerate solutions or decompositions of 3-way arrays[END_REF] [START_REF] Paatero | Construction and analysis of degenerate parafac models[END_REF].

There are many algorithms for calculating the CP decomposition. Due to its simplicity of implementation, the most popular in the literature is the Alternating Least Squares (ALS) originally proposed in [START_REF] Carroll | Analysis of individual differences in multidimensional scaling via an n-way generalization of eckart-young decomposition[END_REF], which is an iterative optimization process that alternately estimates the factor matrices by updating each matrix individually while keeping the others fixed. In this way, the system to be solved is then turned into simple least square (LS) sub-problems. The ALS algorithm is known to converge towards a local minimum under mild conditions [START_REF] Li | Some convergence results on the regularized alternating least-squares method for tensor decomposition[END_REF]. However, the ALS algorithm is highly sensitive to initialization, and its convergence to the global minimum can sometimes be slow, if ever met. In addition, the convergence of the algorithm may, in some cases, fall in swamps [START_REF] Mitchell | Slowly converging parafac sequences: swamps and two-factor degeneracies[END_REF],

where the convergence rate is very low and the error between two consecutive iterations does not decrease.

Several variants of the ALS algorithm have been proposed in the literature in order to reduce the slow convergence of the ALS algorithm. A proper solution has been proposed in [START_REF] Sanchez | Tensorial resolution: a direct trilinear decomposition[END_REF][START_REF] Leurgans | A decomposition for three-way arrays[END_REF][START_REF] Sidiropoulos | Blind parafac receivers for ds-cdma systems[END_REF]. This solution is also known as direct trilinear decomposition [START_REF] Sanchez | Tensorial resolution: a direct trilinear decomposition[END_REF]; it consists in obtaining a first estimate of the factor matrices by constructing a Generalized Eigenvalue (GEVD) problem from two tensor slices. However, such an initialization requires that the two factor matrices are of full rank, and that the third one does not contain zero elements. In [START_REF] Lathauwer | A link between the canonical decomposition in multilinear algebra and simultaneous matrix diagonalization[END_REF], the estimation of the factor matrices of the CP decomposition is linked to the simultaneous matrix diagonalization problem. Another way to improve the convergence speed of the ALS algorithm is to resort to a TUCKER3 compression [START_REF] Andersson | Improving the speed of multi-way algorithms:: Part i. tucker3[END_REF][START_REF] Bro | Improving the speed of multiway algorithms: Part ii: Compression[END_REF]; this is useful at initialization when the dimensions of the tensor are large. In [START_REF] Sidiropoulos | Blind parafac receivers for ds-cdma systems[END_REF], an algorithm that accelerates ALS convergence is proposed, and applies the TUCKER3 compression method followed by an initialization based on the proper analysis; this has now become a usual practice to reduce the computational burden [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF]. The convergence of the ALS algorithm was also improved by the Enhanced Line Search (ELS) method [START_REF] Rajih | Enhanced line search : A novel method to accelerate Parafac[END_REF], which has proven its usefulness when the tensor is affected by degenerative factors (Factor degeneracies), or to decrease its sensitivity to initialization [START_REF] Zarzoso | Robust independent component analysis[END_REF].

The above algorithms improve the speed of the ALS algorithm but remain inadequate to overcome difficult cases when the factor matrices in one or all modes are highly collinear, i.e, when the problem of swamp or bottleneck arises.

Recent works [START_REF] Sahnoun | Joint source estimation and localization[END_REF][START_REF] Farias | Coherence constrained alternating least squares[END_REF] have demonstrated that the introduction of coherence constraints avoids this issue. The proposition in [START_REF] Farias | Coherence constrained alternating least squares[END_REF] is a direct modification of the ALS based on the Dykstra projection algorithm on all correlation matrices, while in proposal [START_REF] Sahnoun | Joint source estimation and localization[END_REF], which consists in simultaneously estimating factor matrices, this method has proven to be one of the most useful methods to overcome the difficulty where the estimated factors are highly collinear.

In this paper, we propose two algorithms to improve both accuracy and convergence speed of the CP decomposition in difficult situations, where swamp and bottleneck problems arise. These algorithms are based on proximal methods [START_REF] Parikh | Proximal algorithms[END_REF] [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF], which are now reliable and powerful optimization tools, leading to a variety of proximal algorithms, such as the proximal point algorithm, the proximal gradient algorithm, the accelerated proximal gradient algorithm, and several others including linearization and/or splitting. We shall be particularly interested in the proximal gradient algorithm and its accelerated version, since they fulfill the assumptions of the CP decomposition problem.

The rest of the paper is organized as follows. The next section presents notations and some properties of third-order tensors. In section 3, we define the coherence constraint ensuring the existence of a best low-rank approximation and then we introduce the proximal methods used in our context. In Section 4 we introduce our new optimization algorithms. In Section 5 we deal with analysis of the simulation results and finally Section 6 concludes the paper.

NOTATIONS AND PRELIMINARIES

Notations and definitions

Let us begin by introducing some key notations and definitions that will be used in this document. Tensors are denoted by calligraphic letters, e.g., T , matrices are denoted by boldface capital letters, e.g., M, vectors are denoted by boldface lowercase letters, e.g., a and Scalars are denoted by lowercase letters, e.g., a. In addition, the p th column of matrix A is denoted by a p , the p th element of a vector a is denoted by a p , the entry of a matrix A in position (i, j) is denoted by A ij and the entry of a tensor T in position (i, j, k) is denoted by

T ijk .
T = a ⊗ b ⊗ c (1)
In Equation (1) above, ⊗ represents the tensor outer product, so that entry (i, j, k) of tensor T is defined by the product

T ijk = a i b j c k .
A tensor T ∈ C I×J×K is said to be rank-1 (also referred to as a decomposable tensor [START_REF] Catalisano | Ranks of tensors, secant varieties of Segre varieties and fat points[END_REF][START_REF] Silva | Tensor rank and the ill-posedness of the best low-rank approximation problem[END_REF][START_REF] Comon | Tensors: a brief introduction[END_REF]) if each of its elements can be represented as :

T ijk = a i b j c k ,
in other words, if it can be expressed as the outer product of three vectors, which will be denoted in a compact form as in [START_REF] Comon | Tensor decompositionsstate of the art and applications[END_REF].

110 Definition 2. The scalar product between two tensors with the same size, X , Y ∈ C I×J×K , is defined as:

X , Y = I i=1 J j=1 K k=1 X * i,j,k Y i,j,k (2)
where * stands for the complex conjugation.

Definition 3. The Frobenius norm . F of a tensor T ∈ C I×J×K is derived from the scalar tensor product:

T F = T , T = ( i,j,k | T ijk | 2 ) (3)
Consequently, the quadratic distance between two tensors X and Y of the same size I × J × K can be determined by the quantity:

X -Y 2 F ( 4
)
Definition 4. Let T ∈ C I×J×K be a tensor, then vec{T } ∈ C IJK×1 represents the column vector defined by :

vec{T } i+(j-1)I+(k-1)IJ = T ijk (5)

PRELIMINARIES

115

A tensor of order N is a mathematical entity defined on a product between N linear spaces, and once the bases of these spaces are fixed, then the tensor may be represented by a N -way array of coordinates [START_REF] Comon | Tensors: a brief introduction[END_REF]. For the sake of simplicity, we use the term tensor in a restricted sense, i.e. as a three-dimensional array of complex numbers (i.e. N = 3), but the generalization to N th order tensors, N ≥ 3, is straightforward. Let us consider a tensor T of order 3 with size I × J × K, its CP-decomposition is defined as follows: where coefficients λ r can always be chosen to be real positive, and decomposable tensors D(r) to have unit norm.

T = R r=1 λ r D(r), (6)
= + +... + I J K a 1 b 1 c 1 a 2 b 2 c 2 a R b R c R λ 1 λ 2 λ R
From Definition (1) of decomposable tensors, another writing of the CP decomposition makes it more explicit :

T = R r=1 λ r (a r ⊗ b r ⊗ c r ), (7)
where λ = [λ 1 , λ 2 , ..., λ R ] is a vector containing the scaling factors λ r and the three matrices

A= [a 1 , a 2 , ..., a R ] ∈ C I×R , B= [b 1 , b 2 , ..., b R ] ∈ C J×R and C= [c 1 , c 2 , ..., c R ] ∈ C K×R represent
the factor matrices. When the number R of terms in ( 7) is minimal, it is called the rank of tensor T , and decomposition [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF] is called the Canonical Polyadic Decomposition (CPD) of T [START_REF] Comon | Tensor decompositionsstate of the art and applications[END_REF]. Note that (7) can be rewritten in terms of tensor entries as:

T ijk = R r=1 λ r a ir b ir c ir . (8)
The CP decomposition is visualized for third-order tensors in Figure 1.

CP-DECOMPOSITION AND PROXIMAL OPERATOR

CP-DECOMPOSITION

Low rank

Even if the tensor of interest is of low-rank, it is often necessary to look for an approximation of low rank R of the observed tensor because of the presence of noise. In the latter case, the observed tensor is indeed generally of generic rank, strictly larger than R [START_REF] Comon | Tensors: a brief introduction[END_REF]. In the low-rank approximation problem [START_REF] Rouijel | Cp decomposition approach to blind separation for ds-cdma system using a new performance index[END_REF],

the goal is to minimize an objective function Υ of the form:

Υ(A, B, C; λ) = X -X 2 F = X - R r=1 λ r (a r ⊗ b r ⊗ c r ) 2 F (9)
Alternatively, the minimization of ( 9) can be written using vectorization defined in ( 5) with x = vec{X } as:

min x Υ( x) = min x x -x 2 F = min A,B,C,λ x - R r=1 λ r (a r b r c r ) 2 F ( 10
)
where symbol represents the Kronecker product [START_REF] Kolda | Tensor decompositions and applications[END_REF][START_REF] Comon | Tensors: a brief introduction[END_REF].

Coherence

Let's introduce the concept of coherence that we will use later, when talking about the conditioning of the problem. The notion of coherence has received different names in the literature: the mutual incoherence of two dictionaries [START_REF] Donoho | Optimally sparse representation in general (nonorthogonal) dictionaries via 1 minimization[END_REF], the mutual coherence of two dictionaries [START_REF] Candes | Sparsity and incoherence in compressive sampling[END_REF], the coherence of a subspace projection [START_REF] Candès | The power of convex relaxation: Near-optimal matrix completion[END_REF], etc. The version here follows that of [START_REF] Gribonval | Sparse representations in unions of bases[END_REF], which has been used particularly as a measure of the ability of algorithms such as basis pursuit and matching pursuit to accurately identify the correct representation of a sparse signal.

Mathematically, let H be a Hilbert space endowed with the scalar product a, b = a H b, and A ⊆ H be a finite set of elements a i with unit norm. The coherence of A is defined as the maximum absolute value of the cross-correlations between the columns of A:

µ(A) = max i =j |a H i a j | (11)
It is obvious that 0 ≤ µ(A) ≤ 1, that µ(A) = 0 if and only if a 1 , ..., a R are orthonormal, and that µ(A) = 1 if and only if A contains at least two collinear vectors, i.e. ∃i = j, λ = 0 : a i = λa j . It should be noted that the L ∞ norm may be bounded by the L 2p norms for large p. This result will allow to bound the coherence by a differentiable quantity:

µ(A) ≤ µ p (A) def = i =j |a H i a j | 2p 1 2p (12)

Conditioning of the problem

One of the interesting properties of tensors of order higher than two, N > 2, is that their CP decomposition is often unique, which is not the case of matrix decompositions [START_REF] Harshman | Foundations of the parafac procedure: Models and conditions for an" explanatory" multimodal factor analysis[END_REF] [39] [START_REF] Sidiropoulos | Blind parafac receivers for ds-cdma systems[END_REF] (the decomposition of a matrix into a sum of rankone matrices also exists, but it is not unique, unless some strong constraints are imposed, such as orthogonality or non-negativity).

Uniqueness means that there is only one set of rank-1 tensors whose sum exactly equals tensor T [START_REF] Domanov | On the uniqueness of the canonical polyadic decomposition of third-order tensors-part i: Basic results and uniqueness of one factor matrix[END_REF]. Note that permutation indeterminacy is inherent in a sum, and that scaling indeterminacy is inherent in the construction of rank-1 tensors [START_REF] Comon | Tensors: a brief introduction[END_REF]. When the decomposition is unique in that sense, factor matrices are sometimes said to be essentially unique. From the above definition of the CP decomposition, it is clear that decomposition [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF] is insensitive to:

• Permutation of the rank-1 terms, which refers to the permutation indeterminacy.

• Scaling of vectors a r , b r and c r , provided the produt of the scaling factors is equal to 1, which refers to the scaling indeterminacy.

In numerical algorithms, it is useful to fix indeterminacies. For instance, columns of factor matrices can be normalized and their norm stored in scaling factor λ [START_REF] Kolda | Tensor decompositions and applications[END_REF]. Another approach described in [START_REF] Rouijel | Cp decomposition approach to blind separation for ds-cdma system using a new performance index[END_REF] and used in our present paper consists in calculating the optimal value of the scaling factor λ to properly control the conditioning of the problem. For this purpose, and for a given matrices A, B and C, the optimal value λ o minimizing the error Υ is determined by cancelling the gradient of (9) w.r.t. λ, which then results to the linear system:

G λ o = f ( 13
)
where G is the Gram matrix of size R×R defined by: G pq = (a p b p c p ) H (a q b q c q ) and f is the R-dimensional vector defined by: f r = Σ ijk T ijk A ir B jr C kr .

Note that entries of G can preferably be obtained by G pq = a H p a q b H p b q c H p c q . Equation ( 13) indicates that coherence plays a central role in the conditioning of the problem. And we can see that scalar products do not appear individually but via their products. Such a statement has profound consequences, particularly on the existence and uniqueness of the solution to problem (9).

Existence

It has been demonstrated in [START_REF] Lim | Blind multilinear identification[END_REF] that if

µ(A) µ(B) µ(C) ≤ 1 R -1 (14)
then the best rank-R approximation exists and that the infimum of ( 9) is attained. The reason is that the error (9) becomes coercive once ( 14) is satisfied, and then must reach its minimum since it is continuous.

Constraint [START_REF] Liu | Denoising of hyperspectral images using the parafac model and statistical performance analysis[END_REF] contains max operators, which are not differentiable. This difficulty can be circumvented by using L 2p norms defined in [START_REF] Sidiropoulos | Blind parafac receivers for ds-cdma systems[END_REF], and imposing differentiable constraints C p (x):

C p (x) def = 1 -R + 1 µ p (A) µ p (B) µ p (C) 0 (15)

Uniqueness

There is a sufficient condition for ensuring the uniqueness of the CP-decomposition [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF] in which the coherences are involved [START_REF] Lim | Blind multilinear identification[END_REF]. However, the use of the following condition [START_REF] Lathauwer | A link between the canonical decomposition in multilinear algebra and simultaneous matrix diagonalization[END_REF] is less restrictive

R ≤ K & R(R -1) ≤ I(I -1)J(J -1) (16)

PROBLEM FORMULATION

As pointed out in [START_REF] Jutten | Blind separation of sources, part i: An adaptive algorithm based on neuromimetic architecture[END_REF], for the sake of simplicity, columns of factor matrices are gathered in a single vector x = vec{[A T , B T , C T ]}. The objective to be minimized consists of two terms: one related to 'data fidelity', which is defined in ( 9) by the cost function Υ, and another one linked to a priori information on model parameters, named 'constraint', which is presently defined in [START_REF]Three-way arrays: rank and uniqueness of trilinear decompositions, with application to arithmetic complexity and statistics[END_REF] and which will be represented by G. One way of reformulating this constrained problem is to consider the unconstrained minimization of:

F(x) = Υ(x) dataf idelity + G(x) constraint ( 17
)
where the function G acts as a barrier in the interior point methods [46] [47];

for example G may be defined as the logarithmic term -ln(C p (x)) weighed by a barrier parameter η > 0. But in order to be able to fully compare our proposed methods with the one proposed in [START_REF] Sahnoun | Joint source estimation and localization[END_REF] that is based on gradient descent algorithm, we have chosen to keep the same objective function, namely:

F(x) = Υ(x) + η exp(-γC p (x)) ( 18
)
where γ is a parameter that controls the sharpness of the penalty C p (x) and η is a penalty weight, which decreases through iterations.

PROXIMAL MAPPING

Proximal mapping has a simple definition yet has long been a powerful tool in optimization, leading to a wide range of algorithms, such as the proximalpoint algorithm, the proximal-gradient algorithm, and many other algorithms involving linearization and/or splitting.

Given a function G, the proximal mapping (or proximal operator) [START_REF] Parikh | Proximal algorithms[END_REF] maps an input point x to the minimizer of G restricted to small proximity to x. The definition of the proximal operator is recalled hereafter.

Definition

Let G ∈ Γ 0 (R N )1 , for every x ∈ R N , the minimization problem (with pa-

rameter µ > 0) minimize y∈R N G(y) + 1 2µ x -y 2 2 ( 19
)
admits a unique solution, which is denoted by prox µG (x).

The operator prox µG : R N -→ R N thus defined is the prox-operator of G, and the parameter µ controls the stretch to which the proximal operator maps points towards the minimum of G, with larger values of µ associated with mapped points near the minimum, and smaller values giving a smaller movement towards the minimum [START_REF] Parikh | Proximal algorithms[END_REF].

In other words, when we evaluate the prox-operator of G, we are attempting to reduce the value of G by penalizing this reduction to not deviate too much from x (The parameter µ is like a "step-size" that controls how much we are penalized for moving away from x).

The proximal mapping gets really handy for composite problem when the minimization has the form of Υ + G with Υ convex and differentiable and G convex with "simple" proximal mapping in the sense that G does not require to be differentiable and that its proximal operator can be evaluated efficiently, i.e. its prox-operator admits a closed form [START_REF] Parikh | Proximal algorithms[END_REF]. This general composite problem can be solved with one of the proximal methods mentioned above, but those that reply to the assumptions of our specific problem [START_REF] Golub | Cf van loan, matrix computations[END_REF], and which are the main proposals in this paper, are the proximal gradient method and the accelerated proximal gradient method. A natural strategy of those methods is firstly to reduce the value of Υ by using unconstrained iterative optimization methods such as descent methods or Newton's method, followed by the reduction of the value of G by applying the prox-operator of G (using the same step-size) and repeat until convergence to a minimizer (under some further conditions). This strategy yields the following iteration:

x (k+1) = prox µ (k) G (x (k) + µ (k) d (k) ) (20)
This last strategy describes the general principle of the proximal gradient algorithm. The accelerated proximal gradient algorithm performs an extrapolation at each iteration, by taking into account information from current and previous iterations.

PROPOSED OPTIMIZATION METHOD

In this article, we propose two methods to overcome swamp and bottleneck problems based on the proximal gradient (PG) and the accelerated proximal gradient (APG) in order to solve the optimization problem defined in [START_REF] Golub | Cf van loan, matrix computations[END_REF].

The general principle of the proposed approaches is detailed in the following paragraphs and summarized in Algorithm 1 and Algorithm 2. There are essentially of two basic steps:

• A gradient step associated with the data fidelity term (denoted by the function Υ).

• A proximal step related to the coherence constraint term (denoted by the function G).

Gradient step

This step improves the approximate solution, by acting only on the data fidelity -regardless of the constraint. This step also involves two other stages.

(i) The first one consists in calculating the descent direction d (k) of Υ, yielding the steepest decrease. In this paper, we propose to use gradient descent method to keep the minimization as simple as possible, and despite its slow convergence, especially for large data set [START_REF] Lee | Applying gradient projection and conjugate gradient to the optimum operation of reservoirs 1[END_REF], these drawbacks will be overcome by the second step of the proximal algorithm. The direction is defined as follows:

d (k) = -∇Υ(x (k) ) ( 21
)
where gradient expressions required to determine the direction of descent d (k) are of the form:

∂Υ ∂A = 2AM A -2N A (22) with M A pq def = Σ jk λ p B jp C kp C * kq B * jq λ * q N A ip def = Σ jk T ijk B * jp C * kp λ * p ( 23
)
Expressions for gradients w.r.t B and C are similar.

(ii) The second stage concerns the determination of the step-size ρ (k) along the chosen direction d (k) . Among many different methods of searching for a good step-size, Backtracking, is widely used. It depends on two parameters α and β, with 0 < α < 0.5 and 0 < β < 1. The idea is to start with a sufficiently large step-size ρ (k) (e.g. ρ = 1) at the beginning, and then reduce it as ρ ← ρ * β, until the following Armijo condition [START_REF] Zhang | Global convergence of a modified fletcher-reeves conjugate gradient method with armijo-type line search[END_REF] is verified:

Υ(x (k) + ρ (k) d (k) ) < Υ(x (k) ) (4)
To conclude, the gradient step can be summarized as:

z (k) = x (k) + ρ (k) d (k) . (5)

Proximal step

Since the previous step concerns only the data fidelity term Υ, the proximal step should adjust the search orientation based on the constraint on coherences G. And to do this, we apply the proximal algorithm to the previous point resulting from the previous step of the gradient, i.e. z (k) , as follows:

z (k+1) = prox ρ (k) G (x (k) + ρ (k) d (k) ) = prox ρ (k) G (z (k) ) = arg min x (G(x) + 1 2ρ (k) x -z (k) 2
2 )

H(x) (6)
This step indicates that prox G (z (k) ) is a point that compromises between minimizing G and being near to z (k) . In the general case, the function describing the constraints may not require to be differentiable but it is recommended to be simple in the sense that its proximal operator can be efficiently evaluated.

Now it remains to calculate an approximation of the proximal operator of G.

We proceed in two stages.

Gradient of H. The gradient of H takes the form:

∇H(x) = ∇G(x) + 1 ρ (k) (x -z (k) ) ( 24
)
where the gradient ∇G(X) is calculated as in [START_REF] Sahnoun | Joint source estimation and localization[END_REF]:

∂G ∂A = γ exp(γC p ) L A p A[(A H A) Ω A -I] (25)
where denotes the Hadamard entry-wise product,

L A p def = (Σ p<q |a H p a q | 2p ) -1 2p -1 µ(B, p) -1 µ(C, p) -1 ,
and Ω A pq def = |a H q a p | 2p-2 . The expressions are similar for the components of B and C.

The minimization in the proximal step is stopped as soon as we meet a good candidate that improves z k , meaning that the new iterate z k+1 ensures the minimization of the function G as well as it does not deviate too much from the iterate z k based on data fidelity.

Monitoring. the Accelerated Proximal Gradient (APG) algorithm first extrapolates a point y k by a combination of the current point and the previous point as:

y k = x k + t k-1 t k (z k -x k ) + t k-1 -1 t k (x k -x k-1 )
and then solves a proximal operator problem. However, for a bad extrapolation y k , F(x k ) may be arbitrarily greater than F(x k+1 ); this may occur because APG is not a monotonous algorithm [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. For this purpose, we introduce a monitor that ensures the descent property F(x (k+1) ) < F(x (k) ); it is defined as follows:

x k+1 =      z k+1 if F(z k+1 ) < F(x k ) x k otherwise. (26)

RESULTS AND DISCUSSION

In this section, we test the performance of the proposed algorithms on three scenarios: (i) in the first one, the limit point lies in the admissible region, i.e. when the constraint is not active (C p ≥ 0) at convergence, (ii) the second scenario is a concrete example in which it is necessary to deal with the swamp problem, i.e. where the factor matrices are highly correlated in all three modes and the constraint is active (C p ≤ 0) at convergence, and (iii) the last scenario addresses the bottleneck problem, where the factors of the first mode are chosen to be co-linear.

To see the interest of our algorithms, we compare them to two other CPD algorithms: i) their counterpart without constraint and without calculating the optimal value of λ defined in [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF], ii) the constrained gradient descent algorithm [START_REF] Sahnoun | Joint source estimation and localization[END_REF].

The performances are evaluated according to three criteria: the error between the calculated and actual factor matrices, the best sum of the "congruences" [START_REF] Tomasi | A comparison of algorithms for fitting the parafac model[END_REF], and the execution speed. The best congruence sum requires finding Algorithm 1: Proximal Gradient algorithm (PG) to minimize [START_REF] Golub | Cf van loan, matrix computations[END_REF] Initialize (A 0 , B 0 , C 0 ) to matrices with unit-norm columns ; calculation of the optimal scaling factor λ * using ( 13) such as:

G 0 λ * = f 0 ;
for k ≥ 1 and subject to a stopping criterion, do 1. Gradient Step (a) Compute the descent direction d (k) as the gradient according to (23) w.r.t. x k :

d (k) = -∇Υ(x k ) (b)
Calculate a step-size ρ k using the backtracking method such:

Υ(x k + ρ k d (k) ) < Υ(x k ) (c) Update z k = y k + ρ k d (k)

Proximal

Step (a) Compute the approximate proximal operator of G at z k using (24) such as:

x (k+1) = prox ρ k G (z k ) 3.
Initialize (A 0 , B 0 , C 0 ) to matrices with unit-norm columns, and the interpolation parameters to : t 0 = 0, t 1 = 1 ; calculation of the optimal scaling factor λ * using ( 13) such as:

G 0 λ * = f 0 ;
for k ≥ 1 and subject to a stopping criterion, do

y k = x k + t k-1 t k (z k -x k ) + t k-1 -1 t k (x k -x k-1 ); t k = 4t 2 k + 1 + 1 2 1. Gradient Step
z k = y k + ρ k d (k)
2. Proximal Step (a) Compute the approximate proximal operator of G at z k using (24) such as:

z (k+1) = prox ρ k G (z k ) (b) Monitoring x k+1 =      z k+1 if F(z k+1 ) < F(x k )
x k otherwise.

3. Extract the three blocks of X k+1 : A k+1 , B k+1 and C k+1 4. Normalize the columns of A k+1 , B k+1 and C k+1 5. calculation of the optimal scaling factor λ * using ( 13) such as: Gλ * = f using [START_REF] Sahnoun | Joint source estimation and localization[END_REF] end for the best permutation σ among the factor matrix columns; it is defined as :

max σ = R r=1 |a H r a σ(r) | a r a σ(r) |b H r b σ(r) | b r b σ(r) |c H r c σ(r) | c r c σ(r) (27)
In order to obtain comparable results and to visualize the behaviour of each algorithm, the starting points for all methods -namely the unconstrained algorithm, denoted by "Algo Unconstrained", constrained algorithm resolved by gradient descent, denoted by "Algo Constrained Gradient", proximal gradient algorithm denoted by ("Algo PG") and the accelerated proximal gradient, denoted by ("Algo APG") -are initialized using the same initial points and share the same stopping criteria:

• The tolerance on the Frobenius norm of the gradient divided by the number of entries in the gradient is set at 10 -8 .

• The maximum number of iteration is set to 10 3 .

In all experiments, The computations are run in Matlab on a computer with Intel i5 CPU (2.6GHz) and 10GB memory running 64bit Mac OS. And the results are obtained from 100 Monte Carlo runs for all scenarios. At each iteration and for each SNR value, a new noise realization is drawn. Note also that we fix the same heuristic choices ρ = 10 and γ = 5 as in [START_REF] Sahnoun | Joint source estimation and localization[END_REF].

Simulation 1

In this experience, we generate a random CP model of size 4 × 3 × 6 with rank 3 and with coherences µ(A) = µ(B) = µ(C) = 0.99, which implies that

µ(A)µ(B)µ(C) = 0.97 is larger than 1 R-1 = 1 2 .
This configuration describes the case where the constraint is active (C p ≤ 0) at convergence. η is varied through iterations. More specifically in this first experiment, η is initialized to 1, and is divided by 10 when Υ(x) is reduced by less than 10 -4 . Figure 2 illustrates the estimation errors of the matrix as a function of SNR. From these results obtained with 100 different initial points, we can see that the unconstrained algorithm produces poor results compared to other algorithms. On the other hand, we can also observe that the two proposed proximal algorithms are more accurate than the constrained gradient algorithm.

Figure 3 represents the best sum of congruences as a function of SNR. We can visualize that with the same initializations, the results are much more promising, especially at a high SNR, which confirms the accuracy of the proposed proximal algorithms.

In order to show a deeper difference between the Algorithms, we will examine the convergence speed of the tensor reconstruction according to the number of iterations, and we will also explore the CPU time to highlight the run time of each algorithm to achieve an accuracy of 10 -6 . Figure 4 indicates that the in terms of number of iterations, followed by proximal gradient algorithm, then by constrained gradient algorithm and finally the unconstrained algorithm which requires a lot of iterations to achieve an accuracy of 10 -6 . This convergence speed result is also clearly presented in Table 1, which indicates the machine time required for each algorithm to achieve an accuracy of 10 -6 , and where it is clearly observed that accelerated proximal gradient yields superior results to deal with the swamp phenomenon.

Simulation 2

In this second experience, we generate a random CP model of size Alternating Least Square (ALS). In this experiment, η is initialized to 0.1, and is divided by 100 when Υ(x) is reduced by less than 10 -4 .

Figure 5 reports the Reconstruction error (9) as a function of the number of iterations. It can hence be observed that the accelerated proximal gradient algorithm converges faster than other algorithms in terms of number of iterations, followed by proximal gradient algorithm, then by the by constrained gradient algorithm, and finally by the unconstrained algorithm. A deeper inspection also reveals that the accelerated proximal algorithm, the proximal gradient algorithm and the constrained gradient algorithm yield 99% of correct estimations whereas the unconstrained algorithm yields 97%.

Simulation 3

In this last experience, we generate a random CP model of size 4 × 3 × 6 with rank 3 and with coherences µ(A) = 0.99, µ(B) = 0.6, µ(C) = 0.6. This configuration addresses the problem of bottleneck, where we chose factors of the first mode (i.e, A) to be almost co-linear. In this experiment, η is initialized to 0.1, and is divided by 100 when Υ(x) is reduced by less than 10 -4 . Figure 6 illustrates the estimation errors of the matrix as a function of SNR. And from these results obtained with 100 different initial points, we can still observe that the unconstrained algorithm produces mediocre results compared to other algorithms. On the other hand, the proximal gradient algorithm produces results similar to those of the constrained gradient algorithm in higher SNR values, par-350 ticularly for the factors of the correlated matrix A. This may be explained by the fact that in higher SNRs, the proximal gradient algorithm is sensitive to the starting points; however the accelerated proximal algorithm remains insensitive to starting points and provides more accuracy than other algorithms. This convergence speed result is also clearly presented in Table 2, which indicates the CPU time required for each algorithm to reach an accuracy of 10 -7 , and it should be noted that the proximal gradient algorithm is a slightly faster compared to the constrained gradient algorithm, while the accelerated proximal gradient algorithm remains the fastest to achieve an accuracy of 10 -7 .

Conclusions

We dient (APG), with a simple and effective monitoring strategy capable of calculating the minimal CP decomposition of three-way arrays. We performed a complete comparison based on computer experiments, which proved the good behaviour of both algorithms in terms of accuracy and convergence speed, even in the presence of bad conditioning, compared to other iterative algorithms available in the literature.

Definition 1 .

 1 The outer product of two vectors a ∈ C I and b ∈ C J defines a matrix M ∈ C I×J M = a ⊗ b = ab T Similarly, the outer product of three vectors a ∈ C I , b ∈ C J and c ∈ C K produces a third order decomposable tensor T ∈ C I×J×K :

Figure 1 :

 1 Figure 1: Visualization of the CP decomposition for a third-order tensor

C k+1 4 .Algorithm 2 :

 42 Extract the three blocks of X k+1 : A k+1 , B k+1 and Normalize the columns of A k+1 , B k+1 and C k+1 5. calculation of the optimal scaling factor λ * using (13) such as: Gλ * = f end for Accelerated Proximal Gradient algorithm (APG) to minimize

( a )

 a Compute the descent direction d (k) as the gradient according to (23) w.r.t. y k : d (k) = -∇Υ(y k ) (b) Calculate a step-size ρ k using the backtracking method such: Υ(y k + ρ k d (k) ) < Υ(y k ) (c) Update

Figure 2 :

 2 Figure 2: Matrix estimation errors, with a random tensor of size 4 × 3 × 6 and rank 3.

Figure 3 :

 3 Figure 3: Congruence versus SNR results, µ(A) = µ(B) = µ(C) = 0.99 up to a precision of 10 -6 and results with 100 random initializations at each SNR.

( a )

 a Up to 1000 iteration (b) Up to a precision of 10 -6

Figure 4 :

 4 Figure 4: Reconstruction error (9) as a function of the number of iterations. For a tensor of size 4 × 3 × 6 and rank 3.

Figure 5 :

 5 Figure 5: Reconstruction error (9) as a function of the number of iterations. For a tensor of size 4 × 3 × 6 and rank 3.

Figure 6 :

 6 Figure 6: Matrix estimation errors, with a random tensor of size 4 × 3 × 6 and rank 3.

Figure 8

 8 Figure8indicates that the accelerated proximal gradient algorithm con-

Figure 7 :

 7 Figure 7: Congruence versus SNR results, µ(A) = 0.99, µ(B) = 0.6, µ(C) = 0.6 up to a precision of 10 -8 and results with 100 random initializations at each SNR.

Figure 8 :

 8 Figure 8: Reconstruction error (9) as a function of the number of iterations. For a tensor of size 4 × 3 × 6 and rank 3.

Table 1 :

 1 CPU time (in seconds) versus SNR results, µ(A) = µ(B) = µ(C) = 0.99 up to a precision of 10 -6 and results with 100 random initializations at each SNR.

	4 × 3 × 6

  have described two methods to overcome swamp and bottleneck prob-

				SNR
	Algorithms	0	10	20	30	40
	Unconstrained	1.85	1.66	1.19	0.95	0.62
	Constrained Gradient 1.38	1.16	0.79	0.57	0.35
	PG	1.37	1.15	0.77	0.52	0.36
	APG	1.35 1.12 0.74 0.47 0.30

Table 2 :

 2 CPU time (in seconds) versus SNR results, µ(A) = 0.99 and µ(B) = µ(C) = 0.6 up to a precision of 10 -7 and results with 100 random initializations at each SNR.

is the class of lower semicontinuous convex functions from R N to ]-∞; +∞] such that dom G = ∅.

with rank 3 and with coherences µ(A) = µ(B) = µ(C) = 0.6, which implies that

This configuration describes the case where the constraint is not active (C p ≥ 0) at convergence. The starting points for all methods are randomly generated, followed by 5 iterations of the