N the set of nonnegative integers.

N n the set of vectors with n-nonnegative integer entries.

N n t the set of vectors with n-nonnegative integer entries bounded by t ∈ N such that

N n t = {α ∈ N n ||α| := n i=1 α i ≤ t}.
R + the set of nonnegative real numbers.

R[x]

the ring of multivariate real polynomials in n variables.

V K V K (I) is the variety associated to I ⊆ R[x], i.e., V K (I) = {x ∈ K n |f (x) = 0 ∀f ∈ I}.

x α the monomial whose degree is |α|, i.e., x α = x α 1 1 . . . x αn n .

p α the coefficient of the monomial x α of the polynomial p(x).

M 0, means that M is semidefinite positive.

⌈•⌉ ⌈x⌉ is the smallest integer not less than x.

•, • the inner product.

• The Frobenius norm.

• 2

The spectral norm.

⊗ the tensor product.

⊠ the Kronecker product.

⊙ the Khatri-Rao product.

deg the degree of a polynomial.

dim the dimension of a subspace.

Ran the range or column space. For instance Ran(M ) is the column space of matrix M . The Canonic Polyadic (CP) tensor decomposition1 has become an attractive mathematical tool these last ten years in various fields, such as blind source separation [Comon 2010], telecommunications [START_REF] De Almeida | [END_REF], chemometrics [Smilde 2005], neuroscience [Becker 2014], sensor array processing [Sahnoun 2015] and data mining [Savas 2008]. The interest in resorting to CP tensor decomposition, compared to more standard matrix-based approaches, lies in the uniqueness of the decomposition.

Let T ∈ K I 1 ו••×I N be a N -order tensor with entries in some field K, the CP decomposition is given by

T = R r=1 λ r (a (r) 1 ⊗ • • • ⊗ a (r) N ), (1.1)
with λ r ∈ K, a (r) j ∈ K I j , 1 ≤ r ≤ R, 1 ≤ j ≤ N , and R is the minimal number of rank-one components, namely, the rank of the tensor [Hitchcock 1927].

Tensors can be represented in a matrix way, which is called unfoldings. Actually, there exists several ways to "matricizing" a tensor. Here, one uses the representation in [Kolda 2009]. For an N -order tensor, one has N different unfoldings, each one under each dimension of the tensor. The n-unfolding of a tensor is given by

T (n) = A n Λ(A N ⊙ • • • ⊙ A n+1 ⊙ A n-1 ⊙ • • • ⊙ A 1 ) T , (1.2)
where

A n = a (1) n • • • a (R) n and Λ = diag λ 1 , • • • , λ R .
Figure 1.1 illustrates a decomposition of a three-way tensor.

The goal of CP decomposition is to find all components λ r and a (r)

i , 1 ≤ i ≤ N, 1 ≤ r ≤ R from tensor T . However, in practice, tensors are corrupted by noise so that one needs to compute an approximate decomposition since the rank is unknown. The approximation problem is stated as

arg min λr, a (r) i =1 T - R r=1 λ r (a (r) 1 ⊗ • • • ⊗ a (r) N ) .
(1.3)

T I 1 Although finding a low rank approximation is an ill-posed problem in general [De Silva 2008, Hillar 2013], some classical algorithms are still useful to compute CP decomposition. The main workhorse for solving CP is Alternating Least Squares (ALS), which is a simple iterative method that updates alternately the factor matrices A n , 1 ≤ n ≤ N in equation 1.2. There are also algorithms based on all-at-once estimation of the matrix factors such as Conjugate Gradient (CG) and Levemberg-Maquardt (LM) methods [Paatero 1997, Comon 2009a]. Some algebraic geometric methods are developed in [Lasserre 2001, Bucero 2014] in order to recover all components, however, there are some limitations with these methods, such as the exponential growth of the number of variables due to convex relaxations which compromises the time of computation for solving the problem. Finally, one cites deflation methods such those in [Cichocki 2009a, Phan 2014], where the main idea is to compute the components of the tensor by applying successive rank-one approximations and subtractions.

I 2 I 3 = + λ 1 • a (1) 1 a (1) 2 a (1) 3 + • • • + λ 2 • a (2) 1 a (2) 2 a (2) 3 λ R • a (R) 1 a (R) 2 a (R) 3
In the following, we present the main motivations regarding the choice of the CP decomposition as my research subject, and also the objectives of the thesis.

Why study the canonical polyadic tensor decomposition?

As mentioned before, the CP decomposition plays a prominent role in a wide number of applications, and it also has advantageous mathematical properties when compared to matrix decompositions. These practical and theoretical points of CP are the main motivations to study this decomposition in my thesis. In the following, we point out the advantage of using tensor approach instead of the matrix one. Next, we draw some example of applications at which the CP decomposition takes place.

Let M be a I 1 × I 2 matrix with entries in some field K. This matrix can be decomposed as a sum of rank-1 matrices as follows

M = R r=1 λ r u r v H r , (1.4)
where u r ∈ K I 1 and v r ∈ K I 2 , ∀r. When rank(M ) = R > 1, there exist infinitely many factors u r , v r , 1 ≤ r ≤ R satisfying the decomposition (1.4). Indeed,

M = λ 1 u 1 v H 1 + λ 2 u 2 v H 2 + • • • + λ R u R v H R = λ 1 (b + u 2 )v H 1 + λ 2 u 2 v H 2 + • • • + λ R u R v H R = λ 1 bv H 1 + u 2 (λ 1 v H 1 + λ 2 v H 2 ) + • • • + λ R u R v H R ,
for b = u 1u 2 ∈ K I 1 . Notice that other choices for b could be made.

In order to ensure a unique decomposition for M , we must impose some conditions. Uniqueness is possible, for instance, if the factors u r and v r , 1 ≤ r ≤ R, are orthogonal to one another (i.e. u i , u j = 0 and v i , v j = 0, i = j), and λ i = λ j , ∀i = j. These constraints define the singular value decomposition (SVD) of M .

The orthogonality of the factors is a strong condition to guarantee the uniqueness of the decomposition of M . In fact, if a problem is modeled as in (1.4) but not with orthogonal factors, the estimation of these factors would be very hard or even impossible. Thus, the CP tensor decomposition arises as a powerful tool that ensures the uniqueness of the factors2 under weaker conditions.

The most known result of uniqueness is a sufficient condition due to Kruskal [Kruskal 1977, Stegeman 2007], which was generalized by other authors to N -th order tensors [Sidiropoulos 2000a]. Let A 1 , A 2 , . . . , A N be the factor matrices of a rank-R complex tensor T . If N n=1 k An ≥ 2R + N -1, then the decompostion of T is essentially unique. k A n is the Kruskal's rank of the unfolding A n . In addition to the Kruskal's condition, necessary conditions to uniqueness were formulated in [START_REF] Liu | [END_REF], De Lathauwer 2006]. Besides, the authors in [Lim 2014] present a condition for uniqueness and existence of low rank tensor approximations based on a coherence measure, giving thus an angular interpretation of the approximate decomposition.

Another motivation for studying the CP tensor decompostion is its applicability in several domains. We described here three examples of applications on digital communications [Lim 2014, Sidiropoulos 2000b[START_REF] De Almeida | [END_REF]. The choice of this domain lies in the fact that part of this thesis was supported by the Department of Teleinformatics Engineering at Federal University of Ceará, so that the interest of developing tensor algorithms keeping in mind this practical domain became central.

Antenna array processing

We show the application in antenna array processing described in [Lim 2014]. Consider a transmission of a wireless communication system composed of n-antenna subarrays, each one with m sensors, and a source in the far field transmitting symbols through a multi-path fading channel. Let b i be the vector defining the spacial position of the sensors for one of the subarrays called reference subarray. The position of the other subarrays with the same sensor structure is defined by the vector ∆ j , 1 ≤ j ≤ n. We set ∆ 1 = 0 since it denotes the reference subarray. Thus, the signal received at discrete time t k , k = 1, 2, . . . , l, on the i-th sensor of subarray j is given by

s ijk = R r=1 σ r (t k ) exp(ψ ijk ), (1.5)
where R is the number of paths, σ r (t k ) is the r-th path of the transmitted signal at time t k , and

ψ ijk = √ -1(ω/c)(b T i d r + ∆ T j d r )
, with ω the pulsation, c the wave celerity and d r the direction of arrival of the r-th path onto the array structure. This problem can be formulated as a CP decomposition. Indeed, set u ir = exp( √ -1(ω/c)(b T i d r )), v jr = exp( √ -1(ω/c)(∆ T j d r )), w kr = σ r (t k )/ σ r , for σ r = [σ r (t 1 ) σ r (t 2 ) . . . σ r (t K )] T , and λ r = σ r . Thus, we have (1.6) where u (r) = [u 1r u 2r . . . u mr ] T , v (r) = [v 1r v 2r . . . u nr ] T , and w (r) = [w 1r w 2r . . . w lr ] T .

s ijk = R r=1 λ r u ir v jr w kr ⇐⇒ S = R r=1 λ r u (r) ⊗ v (r) ⊗ w (r) ∈ C m×n×l ,
Thereby, solving the CP decomposition provides a blind estimation of the directions d r , the signal wave form σ r and the signal λ r power, 1 ≤ r ≤ R.

DS-CDMA communication

The authors in [Sidiropoulos 2000b] use the CP decomposition to estimate some parameters of a DS-CDMA uplink communication between a single antenna composed of I sensors and R users in a multipath fading baseband channel. Let α ir be the gain of the channel between the sensor i and user r, s jr the j-th symbol transmitted by user r, and c kr the k-th chip of the spreading code of the user r. The baseband output of the i-th sensor of the antenna for symbol period j during chip k can be written as r) .

x ijk = R r=1 α ir s jr c kr ⇐⇒ X = R r=1 α (r) ⊗ s (r) ⊗ c (
(1.7)

Hence, the channel gain, the symbols transmitted by the users and the spreading sequence can be estimated by decomposing the received signal X . For more details about interference and synchronization issues of the model see the mentioned article.

Block-fading MIMO channel

Finally, a last example of application in digital communications can be found in [START_REF] De Almeida | [END_REF]. Therein, the authors represent a block-fading MIMO channel as a CP model. Let B be an I ×L matrix denoting the collection of L path gains during I transmission blocks, A R (θ) an M R ×L matrix whose each element consists of the direction of arrival of the l-th path at the receive antenna m r , and a matrix C(τ, φ) that represents a combined space-time channel response depending on the departure angles φ from the transmitted antennas and also on the delay τ of the paths. The dimension of C(τ, φ) is N P × L, where N is the size of a training sequence, and P the ratio between the oversampling factor at the receiver and the symbol rate. Thus, the received signal can be represented as a unfolding matrix of a three-way tensor as follows

X = C(τ, φ)(B ⊙ A R (θ)) T .
(1.8)

By using an algorithm such as ALS, the channel fading, the arrival and departure angles, and the path delay can be estimated from the CP decomposition.

Different objectives and the scope of the thesis

This thesis deals with the CP tensor decomposition problems (1.1) and (B.2) in the real and complex fields. Herein, we intend to contribute with new tools to improve the efficiency of finding solutions of these problems considering orthogonal and non-orthogonal tensors. We also present a chapter describing the reduction of quadratic polynomial systems to best rank-1 tensor approximation problems.

In a first stage, our goal is the computation of an approximate CP decomposition when R = 1. Since we know that the best rank-1 approximation problem is NP-hard, the idea is to propose suboptimal algorithms that gives good rank-1 approximations. For that, we focus on two types of algorithms: finite and iterative ones. In the first case, our idea was to conceive a simple and efficient algorithm compared to state-of-art ones. Therefore, we propose an algorithm called SeROAP, that has been presenting an excellent performance in terms of approximation error for three-way tensors. In the second case, the goal was the conception of an iterative algorithm presenting some advantage over the standard algorithm ALS. The CE algorithm fulfills this role under some conditions.

In a second stage, our objective is to solve the low rank CP decomposition by using an iterative deflation algorithm similar to the hierarchical method described in [Cichocki 2009c]. However, instead of updating alternately the columns of the factor matrices of a tensor, we perform rank-1 approximations. The reason of proceeding with the computation of rank-1 approximations within the iterative deflation procedure is because we noted that the best rank-1 approximation of a tensor can reduce residues very quickly, ensuring a convergence in a few iterations. Again, the difficult of computing the best rank-1 approximation impels us to use finite rank-1 approximations, such as SeROAP, alternatively, which keeps a reasonable computational complexity as residues are reduced with high probability. It is important to mention that before using an iterative deflation algorithm to compute the CP decomposition, we have performed some theoretical studies to ensure the viability of the method.

Another topic described in my thesis is the application of tensors to solve quadratic polynomial systems. The reduction of a multivariate quadratic system into a best rank-1 approximation problem provides a set of new tools to tackle these kind of systems, and also evinces some advantages compared to other standard mathematical tools, such as those in algebraic geometry and Newton-based methods. Moreover, the best rank-1 approximation always exists and it can be computed by standard algorithms. The reduction is made only in the real field. However, we have got extending the approach to a specific complex system but more general than the real one. In the complex case, we propose an algorithm to compute a constrained rank-1 approximation problem.

Finally, we focus on a particular problem of orthogonal CP tensor decomposition, when one of the factor matrices is semi-unitary. The goal of this part of the manuscript is the proposition of a new algorithm based on the combination of the alternating projection method and the deflation procedure, supported by important convergence guarantees.

Structure of the thesis and outline

• Chapter 1: Introduction. The formulation of the CP tensor decomposition. The motivations of studying the CP decomposition: theoretical advantages over matrix tools and a wide number of interesting applications. The description of the topics developed during the thesis. The structure of the manuscript and the publications.

• Chapter 2: Rank-1 tensor approximation. Formulation of the best rank-1 approximation problem. Description of Lasserre's approach to solve this problem. Description of a finite rank-1 approximation algorithm, called SeROAP, that is always at least as good as the state-of-art finite algorithm called ST-HOSVD. Performance of the finites algorithms. Description of an iterative algorithm, called CE, based on a coupled-eigenvalue problem to compute a rank-1 tensor approximation. Performance of the CE algorithm.

• Chapter 3: Iterative deflation. Introduction to the deflation procedure. Description of the iterative deflation algorithm DCPD. Study of the convergence of DCDP under deterministic and probability conditions. Performance of DCDP with finite rank-1 algorithms.

• Chapter 4: Multivariate quadratic systems and the best rank-1 tensor approximation problem. Introduction to quadratic polynomial systems. Reduction from quadratic systems to the best rank-1 approximation problem of three-way real tensors. Discussion on the real tensor approach. Extension to the complex field and reduction to a conjugated partially symmetric rank-1 approximation problem. Description of the CPS algorithm. Examples of real and complex systems with/without perturbations. Performance of the tensor approach for solving real and complex generic systems.

• Chapter 5: Alternating projections on orthogonal CP tensor decomposition. A brief introduction on alternating projection methods and orthogonal tensors. Description of the CAPD algorithm. Theoretical results on the convergence of CAPD. Results on the orthogonal CP decomposition with a column-wise factor matrix.

Publications

CHAPTER 2

Rank-1 tensor approximation

In this chapter we study alternatives to compute rank-1 tensor approximations in the real and complex fields. We start the chapter with the description of the best rank-1 tensor approximation problem, and the presentation of existing methods to tackle it. The three next sections deal with different approaches to solve rank-1 approximation problems: algebraic geometry, finite methods and iterative methods.

In the algebraic geometry approach, we start by adapting the problem to use Lasserre's relaxation method, which consists of satisfying a constraint condition that ensures the existence of a global solution. Then, we verify if a certificate of optimality is ensured for a set of real tensors with entries distributed according to a probability measure, and an example is drawn up. Additionally, we discuss about the limitations of the method to decompose higher-order and higher-dimension tensors. We conclude the section with a discussion about other relaxation methods.

In the part concerning finite methods, we begin describing three algorithms: T-HOSVD, ST-HOSVD and SeROAP. Therein, we show how these methods work and delineate their main features, including advantages and drawbacks when compared to one another. In the sequel, we give a mathematical proof that our proposed method SeROAP performs at least as good as the other ones in terms of approximation error for three-way tensors. Next, we compare the performance of the algorithms by means of simulations for the following scenarios: three-way tensors and higher-order tensors. The goal is to identify in which configurations SeROAP is useful. The section finishes by comparing the estimation of the rank-1 approximation delivered by the three methods with the best rank-1 approximation obtained by means of Lasserre's method, for 2 × 2 × 2 real tensors.

The part about iterative methods deals with two specified methods: the standard ALS algorithm and the proposed CE algorithm. The goal is to propose for three-way tensors that CE can be a competitive algorithm vis-a-vis the standard ALS. First, we detail how our algorithm is formulated, how it works, and what are its mainly features, in particular the improvement of solutions delivered by other rank-1 methods, and convergence. Second, we evaluate the performance of CE and ALS in terms of the approximation error and the average number of iterations to satisfy a stopping criterion, in order to have a meaningful comparison between both methods.

To close this chapter, we boil down to the main points discussed herein and we also point out some directions to future work.

Main contributions of this chapter:

• Adaptation of the rank-1 tensor approximation problem to ensure its global solution can be attained by means of Lasserre's method. We constrain the scalar parameter λ in Problem (B.3) in order to ensure that a specific closed ball be present in the description of the set of constraints, thereby guaranteeing the global minimum of a relaxed version of this problem (see Appendix (A.2)) is attained;

• Finite method to compute rank-1 tensor approximation. We propose an algorithm, called SeROAP, based on successive SVDs and projections, which is less computationally demanding than the T-HOSVD algorithm, and has a competitive complexity compared to ST-HOSVD, for tensors with high dimensions and small order;

• approximation error of SeROAP algorithm. We give an analytic proof that the approximation error obtained by using SeROAP is always at least as good as T-HOSVD and ST-HOSVD, for three-way tensors;

• Limitation of SeROAP algorithm. Heuristics show that SeROAP presents a poorer performance for tensors with order larger than five, when compared to the other algorithms. As a metter of fact, even for 4-th order tensors, ST-HOSVD performs better than SeROAP. However, SeROAP is statistically better than THOSVD with high probability for 4-th order tensors.

• Iterative rank-1 tensor approximation. We propose an iterative method, called CE, based on a coupled-eigenvalue problem that computes rank-1 approximations of three-way tensors. The framework of how construct the method is drawn up and some theoretical results ensuring the improvement and convergence of the objective of the rank-1 approximation problem is also shown;

• Performance of CE method. The method has competitive advantages compared to the standard ALS. We show that the computational complexity of CE is smaller than ALS when one dimension of the tensor is sufficiently larger than the other ones. Moreover, for two types of initializations, simulations show that CE performs better than ALS in terms of approximation errors and average number of iterations in which a stopping criterion is satisfied.

Rank-1 tensor approximation: description and algorithms

Contrary to best rank-r tensor approximations with r > 1, the best rank-1 tensor approximation of any tensor always exists in R or C since the set of rank-1 tensors is in a cone of Segre varieties [START_REF] Abo | [END_REF], Landsberg 2012], which ensures its well-posedness. The complexity class of the best rank-1 approximation problem is NP-hard.

An interest in rank-1 approximations lies, for instance, in deflation of symmetric and orthogonally decomposable tensors [Zhang 2001]. Some applications of the best rank-1 approximation are multi-target tracking [Shi 2013], blind source separation [START_REF] Grellier | [END_REF]], and blind equalization [Grellier 1999]. Recently, the authors in [Yang 2016] apply the rank-1 tensor problem to a class of low rank tensor optimization problems. Another interest lies in their use in iterative deflation technique to compute the CP tensor decomposition, which will be explained in Chapter 3.

Let λ • a 1 ⊗ • • • ⊗ a N be a rank-1 tensor (or decomposable tensor) in K I 1 ×I 2 ×...×I N with unit factors a i , i ∈ {1, 2, . . . , N } and λ ≥ 0. The best rank-1 approximation problem can be stated as

p min = min λ,a i ∈K I i T -λ • a 1 ⊗ • • • ⊗ a N 2 s.t. a i = 1. (2.1)
To compute the global minimum of Problem (B.3), we could resort to algebraic geometric tools such as those described in [Lasserre 2001, Bucero 2014]. However, these techniques are only efficiently applied to small-sized real tensors since they introduce a lot of variables due to relaxations, so that computational complexity and storage requirements become an issue.

Still in the algebraic geometry context, [Nie 2014] describes an algorithm that computes best rank-1 approximations in a reasonable time, but the global minimum is attained only if some restrictive rank conditions are satisfied. Again, the drawback is the slow convergence for tensors with moderate dimensions.

Standard iterative algorithms such those described in Appendix A.1 can be employed to tackle Problem (B.3). However, as mentioned in chapter 1, these algorithms do not guarantee that the global minimum will be attained since they are initialization-dependent algorithms, despite the fact that local convergence is ensured for the standard alternating least squares (ALS) algorithm [Uschmajew 2012]. Another specific iterative algorithm to compute the best rank-1 approximation is introduced by Friedland in [START_REF] Friedland | [END_REF]], which is based on the maximal singular pair of matrices but it does not ensure convergence to the global minimum either.

The best rank-1 tensor approximation can also be computed approximately by finite algorithms, i.e., algorithms terminating within a finite number of steps. The solutions obtained by them can be close to the global minimum of Problem (B.3) and they are generally used for initializing standard iterative algorithms. Moreover, due to their low computational complexity, they are also used in the iterative deflation context. Among existing finite algorithms, we highlight here two algorithms based on the singular value decomposition (SVD). In [De Lathauwer 2000], the authors extend the SVD concept to tensors, and a rank-1 approximation can be obtained by truncating the proposed higher-order singular value decomposition. It is called THOSVD throughout this paper. In [START_REF] Vannieuwenhoven | [END_REF], a new truncated strategy to higher order SVD, called ST-HOSVD, presents a computational complexity smaller than that of THOSVD. Moreover, in terms of approximation error, the authors show that ST-HOSVD is always at least as good as THOSVD for three-way real tensors when at least one of the unfolding matrices of the approximating tensor has a rank equal to 1. Therefore, ST-HOSVD has been an efficient alternative to compute a rank-1 three-way tensor approximation [START_REF] Vannieuwenhoven | [END_REF]]. Other finite rank-1 algorithms are also described in [Grellier 1999[START_REF] Grellier | [END_REF][START_REF] Hackbusch | Tensor spaces and numerical tensor calculus[END_REF]].

For computing rank-1 approximations, we propose a finite algorithm, called SeROAP, that computes a rank-1 tensor approximation by means of a sequence of singular value decompositions of decreasing order tensors followed by a sequence of projections onto Kronecker vectors of increasing size. The goal of the proposed algorithm is to provide an approximation error at most equal to that obtained with ST-HOSVD (and consequently at most equal to that of THOSVD) for three-way tensors, with computational complexity smaller than that of THOSVD, and still competitive with that of ST-HOSVD, at least for tensors of small orders, regardless of their dimensions. These algorithms are used later in an iterative deflation method to compute low rank tensor approximations.

We also propose an iterative algorithm, called CE, based on an alternating coupled eigenvalue problem. We prove that the objective of Problem B.3 convergences to a limit value. Moreover, under some conditions, the computational complex of CE can be lower than that of the standard ALS algorithm.

An algebraic geometric approach

The noteworthiness of many tools on algebraic geometry employed in multivariate polynomial optimization brings about the use of these tools to solve rank-1 CP tensor approximation problems. Particularly, polynomial relaxations [Parrilo 2003, Lasserre 2001, Bucero 2014, Laurent 2009] come up with powerful strategies which turn polynomial-based optimization problems into convex problems with approximative solutions.

Moment relaxation [Lasserre 2001] is an attractive method that approximates real-valued polynomial optimization problems over closed semialgebraic sets as closely as desired to a finite sequence of SDP problems. Frequently, this method obtains exactly the global solution of its unrelaxed problem. We give a brief overview of the moment relaxation technique in Appendix A.2. Some variables employed there are used throughout this section.

The rank-1 problem (B.3) can be reformulated as a polynomial form as follows

p min = min λ∈R + ,a i ∈R I i T (n) -λ • a n (a N ⊠ • • • ⊠ a n+1 ⊠ a n-1 ⊠ • • • ⊠ a 1 ) T 2 s.t. a i 2 = 1, ∀i ∈ {1, 2, . . . , N }, (2.2)
for any unfolding n ∈ {1, 2, . . . , N }. The degree of the cost function is reduced to N + 1 when unit factors a i are allowed for. Indeed,

T (n) -λ • a n (a N ⊠ • • • ⊠ a n+1 ⊠ a n-1 ⊠ • • • ⊠ a 1 ) T 2 = λ 2 -2λ trace(T T (n) a n (a N ⊠ • • • ⊠ a n+1 ⊠ a n-1 ⊠ • • • ⊠ a 1 ) T ) + T 2 .
(2.3) Problem (2.2) has a unique global minimizer in general [Friedland 2014] up to permutation indeterminacies, which accounts for applying the certificate of optimality pointed in Theorem A.2.2, as far as the moment relaxation approach is concerned. However, in Lasserre's relaxation context, we must ensure that the global minimum of relaxed SDP problems is attained, for instance, in the condition of Proposition A.2.1.

In order to satisfy the condition pointed by Proposition A.2.1, we must add on Problem (2.2) the following ball constraint

R 2 -λ 2 + N i=1 a i 2 ≥ 0, (2.4)
where R is the ball radius. The factors a i are unit vectors, therefore only the variable λ need to be bounded to afford a conformed R. Note that λ ≥ 0 and

T (n) -λ • a n (a N ⊠ • • • ⊠ a n+1 ⊠ a n-1 ⊠ • • • ⊠ a 1 ) T 2 = T 2 -λ 2 , if the Frobenius norm is used. Thus, T 2 -λ 2 ≥ 0 =⇒ λ ≤ T . Therefore, λ 2 + N i=1 a i 2 ≤ T 2 + N,
and we can pick any R ≥ T 2 + N to satisfy the inequality (2.4). Now, we construct the closed semialgebraic set K, defined in Appendix A.2, for Problem (2.2). Note that each equality constraint a i 2 = 1 can be written as two inequalities: a i 2 -1 ≥ 0 and a i 2 + 1 ≥ 0. Thus, we have

K = λ ∈ R, a i ∈ R I i λ ≥ 0, a i 2 -1 ≥ 0, -a i 2 + 1 ≥ 0, ∀i ∈ {1, 2, . . . , N }, T 2 + N -(λ 2 + N i=1 a i 2 ) ≥ 0 . (2.5) We set R = T 2 + N. Problem (2.
2) is then reformulated as follows

p min = min K T (n) -λ • a n (a N ⊠ • • • ⊠ a n+1 ⊠ a n-1 ⊠ • • • ⊠ a 1 ) T 2 . (2.6)
The global minimum of its relaxed SDP version, namely p mom t (as described in Appendix A.2), is thus attained with the optimization over K.

As far as Theoreom A.2.2 is concerned, the certificate of optimality is ensured by simulations, as presented next.

Simulations

In spite of the mathematical elegance of the moment relaxation approach, the size of the involved matrices and the number of variables of the relaxed SDP problems impose a serious drawback. As a matter of fact, Problem (2.2) has F = 1 + N i=1 I i variables, whereas the relaxed SDP problem in Lasserre's approach has F +2t 2t variables. This stems from the introduction of a new variable for each monomial considering the lexicographic basis [START_REF] Cox | [END_REF]]. Moreover, the size of the moment matrices is F +t t . These limitations enforce computer simulations of small-sized tensors. Let a sample of 100 real 2 × 2 × 2 tensors whose entries are distributed according to a uniform probability measure in [-1, 1]. For this scenario, d p = 2 and d K = 1, so that t ≥ 2 (vide Appendix A.2). Using the software GloptiPoly c [START_REF] Henrion | [END_REF][START_REF] Henrion | [END_REF] we have observed for all tensors that the certificate of optimality is already ensured for t = 3, i.e., p mom 3 = p min .

We have selected randomly a tensor from the sample, merely as an example, to verify that the solution obtained by solving the SDP problem for t = 3 is the global minimum of Problem (2.2). Let its 1-unfolding T (1) = 0.9483 -0.7513 -0.5911 -0.2967 -0.9169 0.1291 -0.9392 -0.2961 .

In view of Theorem A.2.2, the condition rank(M 2 (y)) = rank(M 1 (y)) does not hold for t = 2. Indeed, rank(M 2 (y)) = 9 and rank(M 1 (y)) = 5. On the other hand, when t = 3, the certificate of optimality is ensured inasmuch as rank(M 2 (y)) = rank(M 1 (y)) = 4. GloptiPoly c extracts the best rank-1 approximation whose a global minimizer is

λ = 1.4791, a 1 = [-0.7701 0.6379] T , a 2 = [0.8944 -0.4473] T and a 3 = [-0.9951 -0.0990] T .
The average time of simulations for computing the best rank-1 approximation with moment relaxation approach for our sample was about 23.4826 seconds, which is overwhelming compared to non algebraic approximate algorithms, such as ALS and SeROAP 1 . As a metter of fact, as mentioned before, the size of the relaxed problem increases exponentially. For 2 × 2 × 2 tensors and t = 2, the relaxed problem has 330 variables and the size of the moment matrix M t (y) is 36 × 36. These sizes rise to 1716 and 120 × 120, respectively, when t = 3.

We are not going to show results for larger tensors because the time for computing best rank-1 approximations with GloptiPoly c becomes exponentially slower. In fact, we have ascertained that even simulations in 3 × 3 × 3 scenarios take more than 30 minutes to compute the global solution.

Other relaxation methods

During our research on applying moment relaxation to the best rank-1 approximation problem, we have found out other researches on this topic that were standing out from other ones, such as those in [Nie 2014] and [Bucero 2014].

The former article proposes semidefinite relaxations algorithms based on SOS, that compute best rank-1 approximations if some rank conditions are satisfied, otherwise a nonlinear optimization method should be applied to improve the solution. In spite of these conditions, the authors claim their methods find best rank-1 approximations most of the time. Results for symmetric and non symmetric real tensors are shown even for moderate-sized ones in a reasonable time. For random tensors with dimensions equal to 15 and order 6, their method finds a best rank-1 approximation in only 4 minutes approximately, which would be unfeasible with Laserre's approach. Nonetheless, the time of simulations becomes a drawback for tensors with dimensions larger than 30, even for three-way tensors.

The latter improves the efficiency of Lasserre's moment approach by using border basis reduction [START_REF] Mourrain | [END_REF]]. The authors take the same examples within [Nie 2014] and show that their method computes best rank-1 approximations. Besides, the method is also applied to find the best rank-2 approximation of a symmetric 3 × 3 × 3 × 3 tensor. The time of simulations is not presented in their results but it is expected to be much faster than Lasserre's.

On account of the auspicious results that had been developing by both aforementioned researches at the time that our research was in progress, we felt ourselves pressured to call off the Lasserre's relaxation approach. As a metter of fact, we have decided to switch to standard tools in linear algebra, as presented in further sections.

Finite rank-1 approximation algorithms

THOSVD and ST-HOSVD are conceived to compute low multilinear rank approximations. The former was proposed in [De Lathauwer 2000] in the context of the generalization of the SVD. It has been playing a prominent role in many applications where compression is desired, such as signal processing [De Lathauwer 2004[START_REF] Muti | [END_REF], computer vision [Vasilescu 2002], and data mining [Savas 2003[START_REF] Savas | [END_REF]]. The latter is an alternative method to compute HOSVD with low computational complexity. Contrary to THOSVD, the approximation error computed with ST-HOSVD depends on the order in which the modes are processed, i.e., different order of modes imply different approximation errors. ST-HOSVD is always at least as good as THSOVD for three-way tensors when at least one multilinear rank is equal to 1. Note that the CP rank-1 tensor approximation is equivalent to the multilinear rank-(1, 1, . . . , 1) approximation, thereby only the rank-1 version of these algorithms are presented in this manuscript.

Our proposed finite method SeROAP computes rank-1 approximations based on a sequence of singular value decompositions followed by a sequence of projections onto Kronecker vectors. The approximation error is at most equal to that of ST-HOSVD for three-way tensors with the same ordering of modes in both algorithms. In this section, we begin describing THOSVD and ST-HOSVD for the rank-1 case and set forth its interest in rank-1 approximations. In a second stage, we draw our attention to SeROAP, describing its formulation and pointing out its advantages.

Truncation of HOSVD

THOSVD is described in the lines of Algorithm 1. Every n-factor is computed by applying the SVD algorithm onto the respective n-mode of T and keeping the dominant left singular vector u n . With all factors, the rank-1 approximation is obtained by computing the tensor product of the N factors, that is U , and performing the contraction of T onto U .

We delineate the following points about THSOVD algorithm.

input : T ∈ K I 1 ×I 2 ו••×I N : input data output: X ∈ K I 1 ×I 2 ו••×I N : rank-1 approximation for n = 1 to N do u n ← first left singular vector of T (n) ; end U ← ⊗ N n=1 u n ; λ ← T , U ; X ← λ • U .
Algorithm 1: THOSVD algorithm -Multilinear rank-(1, 1, . . . , 1).

• For computing the first left singular vector we do not need to compute the complete SVD in every iteration n, but only the dominant singular values. We use Lanczos algorithm running with a number of k steps [Comon 1990] to compute the singular triplet of a m × n matrix, with a complexity O{2kmn}. In practice, taking k equal to 3 or 4 is enough. Thereby, the accumulated complexity computed for all u n in THOSVD is equal to O{2N k N j=1 I j }. The computation of U requires N j=1 I j FLOPS. The contraction to obtain λ (line λ ← T , U ) also needs N j=1 I j operations. To sum up, the total number of operations of THOSVD is of order:

O (2N k + 2) N j=1 I j ;
• Contrary to matrices, the truncated SVD for higher order tensors does not provide a best rank-1 approximation, in spite of keeping a compressed information of T in multilinear sense. Notice also that in every iteration n, the truncated SVD is calculated for modes of size N j=1 I j , which is more computationally expensive when compared with other reduced-order methods (i.e. ST-HOSVD and SeROAP);

Sequentially truncated HOSVD.

The ST-HOSVD algorithm is an alternative strategy that can be used to compute a rank-1 approximation of a tensor [START_REF] Vannieuwenhoven | [END_REF]]. The idea behind this algorithm in this case is to construct the rank-1 tensor approximation with the principal left singular vectors u i , 1 ≤ i ≤ N, as in THOSVD, but computed from a sequence of tensors of smaller and smaller order, which in turn are constructed from the principal right singular vectors. Herein, we omit the effect of singular values, present in the description of ST-HOSVD in [START_REF] Vannieuwenhoven | [END_REF]], because they do not affect the computation of the factors for rank-1 approximations. The approximation error T -X , for an estimated rank-1 tensor X , depends on the order in which the modes are processed.

Let π i , 1 ≤ i ≤ N , be the unfolding operator applied to a tensor T ∈ K I 1 ×I 2 ו••×I N along the i-th mode. The unfolding operator is the procedure of reshaping the entries of a tensor into a matrix form [Kolda 2009]. Herein, T (i) = π i (T ) is a matrix representation of T in the matrix space

K I i ×I 1 I 2 •••I i-1 I i+1 •••I N . Let also π -1
i , 1 ≤ i ≤ N, be its inverse operator applied to a matrix form that recovers the original tensor T , that is T = π -1 i (T (i) ). We assume p = [p 1 p 2 . . . p N ], with p i ∈ {1, 2, . . . , N } and p i = p j for i = j, a vector that defines the ordering of modes processed in ST-HOSVD.

The rank-1 version of ST-HOSVD algorithm is described as follows: the principal left singular vector of the unfolding S (p 1 ) , say u p 1 , is computed as the first factor of the rank-1 approximation. The tensor S = π -1 p 1 (v T p 1 ), where v p 1 is the principal right singular vector of S (p 1 ) , is an N -1-th order tensor, since I p 1 = 1. At the next step, u p 2 is obtained by computing again the principal left singular vector of the unfolding S (p 2 ) of S. The tensor π -1 p 2 (v T p 2 ), where v p 2 is the principal right singular vector of S (p 2 ) , is now an (N -2)-th order tensor (I p 1 = I p 2 = 1). The process continues according to p until all factors u n are computed and π -1 p N (v T p N ) is a scalar number. In order to obtain a suboptimal solution X , we contract T onto U = ⊗ N n=1 u n , viz., λ = T , U , and set the rank-1 approximation as X = λ • U . The method is detailed in Algorithm 2.

input : T ∈ K I 1 ×I 2 ו••×I N : input data, p : processing order output: X ∈ K I 1 ×I 2 ו••×I N : rank-1 approximation S ← T for n ← p 1 , p 2 , . . . , p N do S (n) ← π n (S) u n ← principal left singular vector of S (n) v n ← principal right singular vector of S (n) S ← π -1 n (v T n ) end U ← ⊗ N n=1 u n ; λ ← T , U ; X ← λ • U .
Algorithm 2: ST-HOSVD algorithm -Multilinear rank-(1, 1, . . . , 1).

We have some comments about ST-HOSVD algorithm for rank-1 approximations:

• The ordering vector p impacts on the reduction (or increase) of the computational complexity of ST-HOSVD. For computing the principal singular vectors of S (n) at iteration n, we need O{2k N i=n I p i } operations. Thus, the complexity to compute all factors u n is of order O{2k N -1 n=1 N i=n I p i }. Note that if the dimensions of T were sorted otherwise, the complexity would be different. Moreover, compared to THOSVD, ST-HOSVD has smaller computational complexity. Indeed, as we saw, the main loop of THOSVD has O{2kN N i=1 I i } operations.

• The approximation error of the rank-1 approximation of T also depends on the choice of the ordering vector p. As a metter of fact, for finding the optimal configuration that minimize the approximation error, we should compute the error in N ! different configurations, which is intractable for large order tensors;

• Apart from the conditions brought up before, it is not ensured that ST-HOSVD performs as good as the THOSVD algorithm [START_REF] Vannieuwenhoven | [END_REF]].

Sequential Rank-One Approximation and Projection

The proposed algorithm, called SeROAP, is a competitive finite algorithm that computes a rank-1 tensor approximation. As ST-HOSVD, SeROAP also constructs a sequence of tensors of smaller and smaller order, thereby the approximation error depends on the ordering of modes in the process.

Contrary to ST-HOSVD, SeROAP does not compute the factors at every iteration. Instead, the rank-1 approximation is directly computed after a projection process. The order-N version of SeROAP algorithm goes along the lines depicted in Algorithm 3 (see also [da Silva 2014] for a longer description).

input : T ∈ K I 1 ×I 2 ו••×I N : input data, p : processing order output: X ∈ K I 1 ×I 2 ו••×I N : rank-1 approximation S ← T for n ← p 1 , p 2 , . . . , p N -1 do S (n) ← π n (S) v n ← principal right singular vector of S (n) S ← π -1 n (v T n ) end u ← principal left singular vector of S (p N-1 ) w ← v * p N-1 ⊠ u for n ← p N -2 , p N -3 , . . . , p 1 do X (n) ← S (n) ww H w ← vec(X (n) ) end X ← π -1 p 1 (X (p 1 ) ) Algorithm 3: SeROAP algorithm
The SeROAP algorithm can be described in two different phases: decreasing tensor order and projection. The decreasing order phase (first for loop) is similar to ST-HOSVD algorithm, except that in SeROAP the factors of the rank-1 approximation are still not computed. As a metter of fact, the goal here is to store the unfoldings S (p 1 ) , S (p 2 ) , . . . , S (p N-2 ) to project their rows into Kronecker vectors that will be obtained in the projection phase. After the first loop, a Kronecker vector w is obtained from the last unfolding S (p N-1 ) by computing its rank-1 approximation. For SeROAP, we do not need to take use of the last element p N of the mode ordering.

In the projection phase (second for loop), at the first step, we project the rows of the matrix S (p N-2 ) onto the Kronecker vector w, by performing X (p N-2 ) = S (p N-2 ) ww H . Notice that X (p N-2 ) can be viewed as the unfolding of a rank-1 three-way tensor. Indeed,

X (p N-2 ) =S (p N-2 ) ww H ⇐⇒ X (p N-2 ) = S (p N-2 ) w ⊗ u * ⊗ v p N-1 .
The vector w is updated using the vec operator applied to X (p N-2 ) . Since X (p N-2 ) is the unfolding of a three-way tensor, vec(X (p N-2 ) ) is a Kronecker product of three vectors: vec(X (p N-2 ) ) = v p N-1 ⊠ u * ⊠ S (p N-2 ) w. Thus, at the next iteration, the unfolding S (p N-3 ) is projected now onto a three-way Kronecker vector. This projection generates a matrix X (p N-3 ) that is the unfolding of a 4-th order tensor. By continuing the process, we note that at the end of the projection phase, the vector w is a Kronecker product of N -1 vectors, so that the matrix X (p 1 ) is actually the unfolding of an order-N rank-1 tensor.

We draw some important aspects about SeROAP method:

• The complexity of SeROAP is O{2k N -1 n=1 N i=n I p i } for the decreasing order phase, and O{2 N -1 n=1 N i=n I p i } for the projection phase. Moreover, all matrices S (p 1 ) , S (p 2 ) , . . . , S (p N-2 ) must be stored before the second phase, which gives a total of N -2 n=1 N i=n I p i push operations. Notice that if N is not large, the complexity of SeROAP is smaller than that of THOSVD, and not much larger than that of ST-HOSVD. For instance, for k = 4, N = 3 and I 1 = I 2 = I 3 = 100, we have approximately 24, 8.08, and 10.1 million FLOPS for THOSVD, ST-HOSVD and SeROAP algorithms, respectively. For SeROAP, we need 10 6 floating points in memory to store the unfolding matrix S (p 1 ) .

• As ST-HOSVD algorithm, different ordering vectors p yield different approximation errors;

• The idea behind the computation of the factors u n , for p 1 ≤ n < p N , from the reshaped form of the right singular vectors v n , lies on the fact that the straight line along v n fits the row vectors of S (n) such that the sum of squared distances of the row vectors from the line is minimized. Thus, we try to find approximate vectors w fitting the rows of increasing order unfolded tensors X (n) , and thus obtain a good rank-1 approximation;

• For higher order tensors, the successive projections of SeROAP introduces errors in the rank-1 tensor approximation so that the comparison with THOSVD and ST-HOSVD methods could not be recommended. More details in Section 2.3.2.

Theoretical result for rank-1 approximation methods

For three-way tensors, in SeROAP we only need to compute the principal singular triplet of two matrices: one to reduce the order of the tensor, and the other to construct the Kronecker vector. Yet, only a single projection is performed to obtain the rank-1 approximation. In this case, we present in the following a theoretical result showing that the rank-1 approximation computed with SeROAP is always at least as good as that delivered by ST-HOSVD, for the same ordering of modes.

Theorem 2.3.1 Let T ∈ K I 1 ×I 2 ×I 3 be a three-way tensor. K = C or R. Let also X ST and X Se be the rank-1 approximations delivered by ST-HOSVD and SeROAP algorithms, respectively. Then the approximation error T -X Se ≤ T -X ST holds for any mode ordering defined by p.

Proof: Let p = [p 1 p 2 p 3 ], with p i ∈ {1, 2, 3} and p i = p j for i = j. Let also T (p 1 ) , X ST (p 1 ) and X Se (p 1 ) be the unfolding matrix along mode p 1 of tensors T , X ST and X Se , respectively.

Denoting by u p 1 , u p 2 and u p 3 the factors obtained with ST-HOSVD, the approximation error can be written as

T (p 1 ) -X ST (p 1 ) 2 = T (p 1 ) -λu p 1 (u p 3 ⊠ u p 2 ) T 2 ,
where λ is a positive scalar number for unit factors u n , n ∈ {p 1 , p 2 , p 3 }. Hence,

T (p 1 ) -λu p 1 (u p 3 ⊠ u p 2 ) T 2 = T 2 -λ trace{T H (p 1 ) u p 1 (u p 3 ⊠ u p 2 ) T }- -λ * trace{T (p 1 ) (u * p 3 ⊠ u * p 2 )u T p 1 } + λ * λ.
(2.7)

The scalar λ is actually the contraction of T onto the unit rank-1 tensor u p 1 ⊗ u p 2 ⊗ u p 3 , so that it can be written as

λ = trace{T H (p 1 ) u p 1 (u p 3 ⊠ u p 2 ) T } = u H p 1 T (p 1 ) (u * p 3 ⊠ u * p 2 ).
(2.8)

Plugging (2.8) into equation (2.7), we obtain, after simplifications,

T (p 1 ) -X ST (p 1 ) 2 = T 2 -|λ| 2 . Note that u H p 1 T (p 1 ) = T (p 1 ) 2 v H p 1 for the principal singular triplet (u p 1 , v p 1 , T (p 1 ) 2 ) of T (p 1 )
, where • 2 stands for the spectral norm. Hence

λ = u H p 1 T (p 1 ) (u * p 3 ⊠ u * p 2 ) = T (p 1 ) 2 v H p 1 (u * p 3 ⊠ u * p 2 ) =⇒ |λ| 2 = T (p 1 ) 2 2 |v H p 1 (u * p 3 ⊠ u * p 2 )| 2
, which leads to the following approximation error

T (p 1 ) -X ST (p 1 ) 2 = T 2 -T (p 1 ) 2 2 |v H p 1 (u * p 3 ⊠ u * p 2 )| 2 .
On the other hand for SeROAP, we have X Se (p 1 ) = T (p 1 ) ww H , for w = v * p 2 ⊠ u. Thus,

T (p 1 ) -X Se (p 1 ) 2 = T (p 1 ) -T (p 1 ) ww H 2 = T 2 + T (p 1 ) w 2 w 2 -2 T (p 1 ) w 2 = T 2 -w H T H (p 1 ) T (p 1 ) w.
(2.9)

The vectors u and v p 2 are the principal left and right singular vectors of the unfolding S (p 2 ) , as described in Algorithm 3. The eigenvalue decomposition of T H (p 1 ) T (p 1 ) can be expressed by

T H (p 1 ) T (p 1 ) = T (p 1 ) 2 2 v p 1 v H p 1 + V , (2.10)
where V is a semidefinite positive matrix. Plugging (2.10) into (2.9)

T (p 1 ) -X Se (p 1 ) 2 = T 2 -T (p 1 ) 2 2 w H v p 1 v H p 1 w -c = T 2 -T (p 1 ) 2 2 |v H p 1 w| 2 -c, with c = w H V w ≥ 0.
To complete the proof of the theorem, we just need to show that |v

H p 1 w| 2 ≥ |v H p 1 (u * p 3 ⊠ u * p 2 )| 2 , or equivalently that | w, v p 1 | ≥ | u * p 3 ⊠ u * p 2 , v p 1 |.
This is true, because w is by construction (cf. Algorithm 3) the vector closest to v p 1 among all vectors of the form a ⊠ b where a and b have unit norm.

Corollary 2.3.2 Let X TH be the rank-1 approximation delivered by THOSVD algorithm. Then the approximation error T -X Se ≤ T -X TH holds for any mode ordering defined by p.

Proof: The proof follows directly from Theorem 2.3.1 above and from Theorem 7.2 described in [START_REF] Vannieuwenhoven | [END_REF]], in which the approximation error computed with ST-HOSVD is at most equal to that obtained with THOSVD for three-way tensors, when at least one of the multilinear ranks of the approximating tensor is equal to 1.

Performance of finite rank-1 approximation methods

In this section, we are going to evaluate the performance of SeROAP, THOSVD and ST-HOSVD in terms of approximation error. Firstly, we compare SeROAP and THOSVD for different threeway tensor sizes with entries in the complex field, and considering three different order of modes. The same comparison is made in the following for SeROAP and ST-HOSVD. Secondly, we evaluate all methods under higher order tensor scenarios. Finally, we measure how far the estimated rank-1 approximations are from the best rank-1 approximation, which is obtained from Lasserre's approach.

Scenarios of three-way tensors

To compare the approximation error of the rank-1 approximation methods SeROAP, THOSVD and ST-HOSVD, we have simulated four tensor scenarios: 3 × 4 × 5, 3 × 4 × 20, 3 × 20 × 20 and 20 × 20 × 20, which are labeled as scenarios 1, 2, 3 and 4, respectively. With these scenarios, we have intended to take into account cases with equal dimensions or not. For each of them, a sample of 300 complex tensors with real and imaginary parts uniformly distributed in [-1, 1] was generated thereby ensuring a coherent and meaningful comparison.

As mentioned in previous sections, the approximation error of SeROAP depends on the choice of the order of modes. Therefore, we have computed the solution delivered by SeROAP in three different modes defined by the permutations

p 1 = [1 2 3], p 2 = [2 3 1], and p 3 = [3 1 2].
In order to characterize the comparison between SeROAP and THOSVD methods, we introduce the following metric

∆φ t = T -X TH -T -X Se ,
which is the difference between the approximation errors obtained by both methods. Notice that ∆φ t > 0 implies that SeROAP performs better than THOSVD. As predicted by Proposition 2.3.1, the results show that the approximation error is smaller for SeROAP algorithm. Indeed, it turns out that ∆φ t > 0 for the three order of modes employed in SeROAP.

The behavior of ∆φ t does not change significantly throughout the different permutations. However, with the variation of dimension sizes, we have noted different distribution of ∆φ t values. As a metter of fact, the larger the dimensions, more ∆φ t concentrates around some value. For instance, the scenario with 20 × 20 × 20 tensors, ∆φ t ≈ 0.5. Notice also the distribution of ∆φ t values reduces as the size of dimensions increases (compare scenarios 2, 3 and4). Now, we focus on the comparison of the approximation error obtained with SeROAP and ST-HOSVD. As defined before, we use the following metric that computes the difference between the approximation error evaluated by both algorithms.

∆φ s = T -X ST -T -X Se ,
Again, if ∆φ s > 0, SeROAP algorithm performs better than ST-HOSVD in terms of approximation error. Contrary to THOSVD algorithm, the behavior of the approximation error depends on the order of modes, which can be viewed in our results. Note that for p 1 , the values of ∆φ s in the scenario 3 × 20 × 20 are close to zero, which implies that SeROAP and ST-HOSVD have approximately the same approximation errors. Indeed, for this scenario, 0 ≤ ∆φ s ≤ 0.05 for 97.33 % of tensors.

When the ordering vector used in both algorithms is p 2 or p 3 , SeROAP provides a better performance compared to ST-HOSVD. This can be noted throughout the three last scenarios inasmuch as ∆φ s values are rather far from zero. As expected, results also confirm Theorem 2.3.1 irrespective the employed scenarios.

It is important to mention that our results reinforce the result described in [START_REF] Vannieuwenhoven | [END_REF]] in which ST-HOSVD outpeforms THOSVD in terms of approximation error. Indeed, while ∆φ t varies from 0.1 and 1.1, the values of ∆φ s fluctuates between 0 and 0.4, approximately.

Scenarios of higher order tensors

Theorem 2.3.1 ensures that the error of the rank-1 approximation of SeROAP is smaler than those of THSOVD and ST-HOSVD methods only for three-way tensors. In this section, we evaluate the impact on the approximation errors when the order of the tensors is higher than 3.

As before, we consider a sample of 300 complex tensors with real and imaginary parts uniformly distributed in [-1, 1]. We evaluate two higher-order scenarios: 3× 4 and 3× 5 . We have chosen for both algorithms the ordering p 1 = [1 2 3 4] andp 2 = [1 2 3 4 5 For both scenarios, we note that Theorem 2.3.1 does not hold anymore. Indeed, for the 4-th order tensor scenario, it turns out by our numerical experiments that SeROAP delivers a better rank-1 approximation with high probability than that of THOSVD, but its performance worsens drastically when compared to ST-HOSVD. The approximation errors computed with SeROAP is even worst for 5-th order tensors. Since the successive projections of SeROAP introduces substantial errors in the rank-1 tensor approximation, the comparison is not recommended for higher order tensors.

Comparison to best rank-1 approximation in 2 × 2 × 2 scenario Now, we assess the approximation errors in 2 × 2 × 2 scenario of the previous finite methods by comparing them to the best rank-1 approximation, which is computed using Lasserre's approach. We have generated a sample of 200 real tensors with entries uniformly distributed in [-1, 1]. In simulations, we have used only the permutation vector p = [1 2 3]. In order to have a meaningful comparison, we consider the MSE metric given by

MSE = 1 200 200 n=1 ∆φ [n] m 2 ,
where ∆φ

[n] m = T [n] -X [n] m -T [n] -X [n]⋆ . Here, m ∈ {"T H", "ST ", "Se"} is a tag indicating the method, X [n]
m is the rank-1 approximation of T [n] delivered by the method m, and X [n]⋆ is the best rank-1 approximation of T [n] . Clearly, ∆φ

[n] m ≥ 0, ∀m, n. We also point out the percentage of tensors in which ∆φ

[n] m ≤ 10 -4 in order to measure how close to the best approximation the solutions of the finite methods are. The table bellow depicts these results.

Algorithm

m MSE % ∆φ

[n] m ≤ 10 -4 THOSVD "TH" 0.0230 1.5 % ST-HOSVD "ST" 0.0013 11.0 % SeROAP "Se" 2.4511e-04 26.5 % Table 2.1: Performance of finite rank-1 methods for 2 × 2 × 2 tensors.

As expected, SeROAP exhibits the lowest MSE among all methods, and THOSVD the worst one. Moreover, SeROAP also delivers rank-1 approximations close to the best one more often than ST-HOSVD and THOSVD. While only for 3 tensors the condition ∆φ

[n] m ≤ 10 -4 hold s in THSOVD, a total of 53 tensors have satisfied this condition in SeROAP. ST-HOSVD has got the second better performance.

Iterative rank-1 approximation methods

To compute rank-1 approximations using iterative methods, we can use those ones mentioned in the Introduction. Due to high complexity of Newton and gradient-based methods, ALS algorithm is usually preferred to do this task. As a metter of fact, it is easy to be implemented and ensure satisfactory results in general. Additionally, ALS also presents a general global convergence property for rank-1 approximations [START_REF] Wang | [END_REF]]. Indeed, for a given initialization point, a sequence of iterates always converge to a unique limit point in general.

To compete with ALS, it is necessary to propose some method that combines simple implementability and good convergence features. For this, we propose a new method for three-way tensors, called CE (Coupled-Eigenvalue), that:

• presents less computational complexity than ALS when one dimension is larger enough than the other ones;

• for any initialization, the approximation error always converges to a limit value;

• performs better than ALS in terms of approximation error and convergence rate in most cases according to simulations;

Coupled-eigenvalue rank-1 approximation

This section presents an alternating eigenvalue method for three-way tensors that can improve local solutions obtained from any other rank-1 approximation method (e.g. SeROAP, ST-HOSVD and THOSVD algorithms).

Let T be a tensor in K I 1 ×I 2 ×I 3 and t in be the vectorization of its slice i n , 1 ≤ i n ≤ I n , and n ∈ {1, 2, 3}. Without loss of generality, we assume n = 3. The rank-1 approximation problem B.3 can be stated as

[z opt , x opt , y opt ] = arg min z,x,y Υ(x, y, z) s.t. x = 1, y = 1, (2.11) with Υ(x, y, z) = I 3 i 3 =1 t i 3 -z i 3 (y ⊠ x) 2 and z = [z 1 • • • z I 3 ]. Indeed, T -x ⊗ y ⊗ z 2 = T (3) -z(y ⊠ x) T 2 = I 3 i 3 =1 t i 3 -z i 3 (y ⊠ x) 2 .
Let L 0 be the lagrangian of Problem (2.11) given by

L 0 (x, y, z, η 1 , η 2 ) = I 3 i 3 =1 t i 3 -z i 3 (y ⊠ x) 2 + η 1 ( y 2 -1) + η 2 ( x 2 -1),
where η 1 and η 2 are the Lagrange multipliers.

By computing the stationary points of L 0 , we obtain for each z i 3 the following

∂L 0 ∂z * i 3 = -(y ⊠ x) H t i 3 + z i 3 = 0 =⇒ z i 3 = (y ⊠ x) H t i 3 .
Plugging the value of z i 3 into Problem (2.11), we can rewrite it as the following equivalent maximization problem

[x opt , y opt ] = arg max x,y (y ⊠ x) H M (y ⊠ x) s.t. y 2 = 1, x 2 = 1,
(2.12) where M = I 3 i 3 =1 t i 3 t H i 3 . Now, we decompose M as a sum of Kronecker products. This can be done by reshaping M and applying the SVD decomposition [Van Loan 1993]. Thus, M can be given by

M = R ′ r=1 Q (r) ⊠ P (r) , with matrices P (r) ∈ K I 1 ×I 1 and Q (r) ∈ K I 2 ×I 2 . R ′ is the Kronecker rank of M satisfying R ′ ≤ I 1 I 2 . Substituting M into Problem (2.
12), we have: r) x). This holds because

[x opt , y opt ] = arg max x,y Γ(x, y) s.t. y 2 = 1, x 2 = 1, (2.13) with Γ(x, y) = R ′ r=1 (y H Q (r) y)(x H P (
(y ⊠ x) H R ′ r=1 Q (r) ⊠ P (r) (y ⊠ x) = R ′ r=1 (y ⊠ x) H (Q (r) ⊠ P (r) )(y ⊠ x) = = R ′ r=1 (y H Q (r) ⊠ x H P (r) )(y ⊠ x) = R ′ r=1 (y H Q (r) y) ⊠ (x H P (r) x) = R ′ r=1 (y H Q (r) y)(x H P (r) x).
Now, let L be the Lagrangian of Problem (2.13)given by

L(x, y, η 3 , η 4 ) = -Γ(x, y) + η 3 ( y 2 2 -1) + η 4 ( x 2 2 -1).
By computing the critical points, we obtain a pair of coupled eigenvalue problems r) , and

   y H A (1,1) y • • • y H A (1,I 1 ) y . . . . . . . . . y H A (I 1 ,1) y • • • y H A (I 1 ,I 1 ) y    x = λx (2.14) and    x H B (1,1) x • • • x H B (1,I 2 ) x . . . . . . . . . x H B (I 2 ,1) x • • • x H B (I 2 ,I 2 ) x    y = λy, (2.15) where λ = η 3 = η 4 , A (m,n) = R ′ r=1 P (r) mn Q (
B (k,l) = R ′ r=1 Q (r) kl P (r) , with 1 ≤ m, n ≤ I 1 and 1 ≤ k, l ≤ I 2 .
The coupled-eigenvalue algorithm is presented in Alg. 4. We can initialize the algorithm by computing x 0 and y 0 from the rank-1 approximation solution obtained with SeROAP, THOSVD or any other rank-1 approximation method.

The complexity per iteration of the CE algorithm is dominated by the construction of the matrices in the LHS of (2.14) and (2.15), which is of order

O{min(I 2 1 I 2 2 , I 2 1 I 2 3 , I 2 2 I 2 3 )}. Suppose I 3 is input : X : rank-1 approximation of a tensor T output: X ⋆ : improved rank-1 approximation Compute x 0 from X as x 0 ← X (:, i 2 , i 3 )/ X (:, i 2 , i 3 ) for some i 2 , i 3 ; t ← 0; repeat Set x ← x t in eigenvalue problem (2.15
) and take y t+1 as the eigenvector whose eigenvalue is maximum; Set y ← y t+1 in eigenvalue problem (2.14) and take x t+1 as the eigenvector whose eigenvalue is maximum; t ← t + 1; until some stopping criterion is satisfied;

for i 3 = 1 to I 3 do z ⋆ i3 ← t i3 , y t ⊠ x t ; end X ⋆ ← x t ⊗ y t ⊗ z ⋆ ;
Algorithm 4: CE rank-1 approximation. Above, we chose to start with x, but we could equivalently have started with y.

the largest dimension. If I 3 ≫ I 1 I 2 , then we can take advantage of the CE algorithm in terms of complexity in comparison with the ALS algorithm. Indeed, the complexity per iteration of ALS for rank-1 approximation is of order O(3I 1 I 2 I 3 ), which is higher than that of the CE algorithm in this case. For instance, let I 1 = 5, I 2 = 10 and I 3 = 50. While ALS performs 7500 operations/iteration, the computational complexity of CE is of order 2500. Notice, however, that a properly comparison makes sense if the same initialization is employed in both algorithms.

The following proposition shows that the above algorithm improves (in worst case the solution remains the same) any rank-1 approximation algorithm.

Proposition 2.4.1 Let X be a rank-1 approximation of a three-way tensor T delivered by some method. If X and X ⋆ are the input and output of CE algorithm, then the inequality T -X ⋆ ≤ T -X holds.

Proof: Plugging the expression of A (m,n) into equation (2.14), we obtain, after simplifications, that λ = Γ(x, y), which is the objective function of Problem (2.13). The same result is obtained when the matrix B (k,l) is plugged into equation (2.15). Now ∀t ≥ 1, let λ Γ(x t , y).

Also, the eigenpair (λ

(x)
t+1 , x t+1 ) obtained by setting y = y t+1 in equation (2.14), is solution of the problem λ

(x) t+1 = max x 2 =1
Γ(x, y t+1 ).

Since Γ(x t+1 , y t+1 ) = max

x 2 =1
Γ(x, y t+1 ), it follows in particular that

Γ(x t , y t+1 ) ≤ Γ(x t+1 , y t+1 ), which implies that λ (y) t+1 ≤ λ (x)
t+1 . Similarly, plugging x t+1 into equation (2.15), we can conclude that λ

(x) t+1 ≤ λ (y)
t+2 for the reason that

Γ(x t+1 , y t+1 ) ≤ Γ(x t+1 , y t+2 ).
Hence, the sequence

{Γ t } t∈N = {. . . , λ (y) t , λ (x) t+1 , λ (y) t+1 , λ (x) t+2 , . . .}
is monotonically non-decreasing. The same conclusion would be achieved if we begin by plugging x t into equation (2.15). Now, let X = x 0 ⊗ y 0 ⊗ z 0 be a rank-1 approximation obtained with any other method. Assume x 0 and y 0 are unit vectors, and define λ 0 = Γ(x 0 , y 0 ). By setting x = x 0 in equation (2.15) in the first iteration (a similar operation would be possible for y 0 in equation (2.14)), we clearly have

λ 0 ≤ λ (y) 1 ≤ λ (y)
tmax , where t max is the iteration in which the stopping criterion is satisfied. Since the optimization problems (2.11) and (2.13) are equivalent, z ⋆ i 3 , 1 ≤ i 3 ≤ I 3 , can be obtained by performing the scalar product between vectors t i 3 and y tmax ⊠ x tmax (which is equivalent to contracting tensor T on x tmax and y tmax ). Hence, the tensor

X ⋆ = x tmax ⊗ y tmax ⊗ z ⋆ is a better rank-1 approximation of T than X , implying T -X ⋆ ≤ T -X .
Corollary 2.4.2 For any input X in CE algorithm, the objective of (2.13) always converges to a limit value.

Proof: In the proof of Proposition 2.4.1, we have shown that {Γ t } is monotonically non-decreasing for any input X . Let p ⋆ be the maximum of the objective (2.13). Since the best rank-1 approximation problem always has a solution, then p ⋆ < ∞. But max x Γ(x, y t+1 ) ≤ max x,y Γ(x, y), which implies that {Γ t } t∈N is bounded above by p ⋆ . Since {Γ t } is a real non decreasing sequence bounded above, it converges to a limit Γ ⋆ , Γ ⋆ ≤ p ⋆ .

Performance of iterative methods

The next section deals with the performance of CE algorithm in some simulated scenarios. We compare for a sample of complex tensors our CE method with the standard ALS algorithm using the estimated approximation error. We also determine the average iteration in which both algorithms satisfy a stopping criterion.

Performance CE and ALS methods in general scenarios

In order to compare CE and ALS algorithms, we assume that both algorithms are initialized with the same factors. We assume two types of initializations: random and that delivered by SeROAP method running with the permutation vector p = [1 2 3]. The performance of the algorithms are evaluated in four scenarios: 3 × 4 × 5, 3 × 4 × 20, 3 × 20 × 20 and20 × 20 × 20 complex tensors. As performed in Section 2.3.2, the choice of scenarios intends to account for the variability of tensor dimensions. A sample of 300 tensors with real and imaginary parts distributed according to a uniform measure in [-1, 1] was generated. The stopping criterion of the algorithms is the maximum number of iterations

(it max = 200) or |E it -E it-1 | ≤ 10 -6
, where E it is the approximation error at it-th iteration.

Let

∆φ [n] = T [n] -X [n] ALS -T [n] -X [n]
CE be the difference between the approximation errors of tensor n of the sample obtained with ALS and CE. Notice that if ∆φ [n] ≥ 0, then CE method performs better than ALS for the tensor T [n] .

Tables below present the percentage of tensors in which ∆φ [n] ≥ 0 and the mean iteration in which the stopping criterion is satisfied for two types of initializations. We have rounded the mean iteration to the nearest integer.

Scenario

% ∆φ [n] ≥ 0 mean iteration (ALS) mean iteration (CE) Notice that in the three first scenarios with the SeROAP initialization, CE method performs better than ALS for most of tested tensors. In last scenario, the performance of CE is still better than that of ALS but for less tensors.

For the random initialization, CE carries on with the better performance, but the percentage of ∆φ [n] decrease in all scenarios.

As expected the mean iteration is greater for the scenario with random initialization. Additionally, CE algorithm converges faster than ALS on average in all cases.

Chapter Summary and Directions

We conclude this chapter by presenting the main points observed in the results and also we given some directions for future work.

• Lasserre's approach for computing rank-1 approximations. Despite the fact that the moment approach provides a certificate of optimality to verify whether the estimated solution delivered by the relaxed problem (A.5) is optimum (which is true in general), the method is inefficient even for medium-sized tensors. This is the reason why only 2 × 2 × 2 scenarios were simulated.

As a metter of fact, we confirmed that computation of the best rank-1 approximation for a 3 × 3 × 3 tensor already requires a lot of time (about 30 minutes). In order to mitigate this drawback in future work, we incite the use of the algorithm proposed in [Bucero 2014], which improves Lasserre's method. However, we can not expect a meaningful improvement in terms of computational time for larger tensors, since the exponential increasing of variables due to relaxations is still an inconvenience. Thus, since the monomials of the polynomial describing the objective of the best rank-1 approximation problem are sparse, we expect that optimization can still be done in order to reduce the computational time of the rank-1 estimation.

• SeROAP method. Our proposed method was proved useful for computing rank-1 approximations of three-way tensors. Indeed, we gave a mathematical proof for three-way tensors in which the approximation error is at most equal to those computed with standard methods, such as the truncated HOSVD (THOSVD) and the more recently method called ST-HOSVD. For fourth order tensors, we confirmed by simulations that SeROAP still brings some advantages over THOSVD. However, for higher order tensors, we saw that its performance is hugely degraded. Additionally, SeROAP allowed us to compute a rank-1 approximation closer than the best one in 2 × 2 × 2 scenarios, when compared to the other methods. Unfortunately, we could not see how far the approximation errors delivered by the mentioned finite methods are from the best rank-1 approximation for tensors with larger dimensions, due to limitations of Lasserre's method for computing solutions in a reasonable time. The development of more efficient algorithms allowing us to compute the best rank-1 approximation faster is still a challenge.

• CE algorithm. Although this method has better computational complexity than ALS only in specific conditions (one dimension must be large enough than the other two dimensions), the approximation error and the average number of iterations to convergence obtained with CE is better than those computed with ALS. These results showed that CE is a very competitive iterative method to compute rank-1 approximations for three-way tensors. In the future, we aim to analyze some convergence aspects of CE algorithm. For instance, we still do not have any theoretical result showing that CE algorithm converges to a unique stationary point for any initialization (we have only numerical proofs).

CHAPTER 3

Iterative Deflation

This chapter presents a detailed study of iterative deflation algorithms based on finite rank-1 approximations with the goal of computing low rank tensor approximations or the exact CP decomposition.

First, we introduce the general concept of deflation and present some examples of its usage in some mathematical domains. Specifically, we present the iterative deflation concept on low rank tensor approximations. Next, we perform a detailed description of the so-called DCPD algorithm, whose concept is based on the computation of successive rank-1 tensor approximations and deflation procedures. Contrary to others deflation algorithms, called in the literature as hierarchical algorithms, we use within DCPD the finite rank-1 approximation algorithms described in Chapter 2.

In the sequel, we perform a theoretical study on iterative deflation in order to analyze the convergence of the DCPD algorithm assuming a best rank-1 approximation to update rank-1 components. Some important lemmas, propositions and a conjecture are presented to ensure the convergence of DCPD algorithm under some conditions for the exact CP decomposition.

In the last part of this chapter, we draw some numerical experiments. In a first stage, we evaluate the percentage of successful decompositions of the DCPD algorithm for the three rank-1 approximations algorithms: SeROAP, THOSVD and ST-HOSVD. We compare the performance with the standard ALS and the iterative deflation algorithm HALS. Next, for the same algorithms, we analyze the behavior of residuals and convergence under noisy scenarios, and the average approximation error under rank and order variation. The chapter is closed with some conclusion and perspectives of future works.

Main contributions of this chapter:

• DCPD algorithm. We propose an iterative deflation algorithm to compute low rank tensor approximations. The idea of iterative deflation algorithms, originally called hierarchical algorithms, was already introduced in the literature, but here we replace the original alternating update of factors of the known HALS algorithm by finite rank-1 tensor approximation algorithms, with the goal of reducing residuals in a few iterations. Our algorithm is called DCPD, stands for Deflation Canonical Polyadic Decomposition.

• Theoretical study on the convergence of DCPD algorithm. We show that the norm of residuals is monotonically reduced within the iterative deflation process. We also prove that the DCPD algorithm recovers the exact CP decomposition of a given tensor when residuals do not fall within a cone with an arbitrary small volume. In a second stage, we prove that the iterative deflation method can reduce the norm of the initial residual by a factor smaller than (sin(β)) L-1 (β being the angle of a suitable cone where the residuals can fall in) after L iterations with high probability, when tensors are distributed according to an absolutely continuous probability measure, and the probability function of residuals is continuous on some suitable angular interval. We also present a conjecture stating the existence of probability measures ensuring the convergence of the DCPD algorithm to an exact decomposition with high probability.

• Performance of DCPD algorithm. Our conjecture is reinforced for DCPD algorithm by means of the computation of CP decompositions of tensors with the entries distributed according to a uniform distribution. In general, we also verify by simulations that residuals decrease faster with DCPD + rank-1 approximation algorithm, when compared to the HALS algorithm. We also present a framework of experiments for noisy and noiseless scenarios with rank and order variations.

Deflation on CP decomposition

Generally speaking, deflation is a computational technique applied to a wide range of mathematical problems whose purpose is to remove some mathematical features and hence to deal with less complex problems. As a simple example, consider that some algorithm (Newton, gradient, etc) is employed to compute a root α of a univariate polynomial p(x). This root can be eliminated (defleted) from the polynomial p(x) by performing

q(x) = p(x) x -α ,
where deg(q(x)) = deg(p(x)) -1. Clearly, the roots of q(x) are also roots of the polynomial p(x), but now one root had been removed the division by xα. The same algorithm combined with this deflation procedure can be sucessively applied to obtain all other roots of p(x).

Deflation also takes place in the improvement of conjugate gradient methods, in which a subspace is hidden (deflated) from the algoritm itself in order to improve its convergence speed. See [Nicolaides 1987] for a longer discussion.

A last example is the computation of rank-r approximations of matrices. It is known that the rank-r approximation of an m × n matrix A can be obtained from its r-dominant singular triplets

(σ i , u i , v i ), 1 ≤ i ≤ r, with r ≤ min(m, n). As a metter of fact, r i=1 σ i u i v H i = arg min X,rank{X}=r A -X .
For large-sized matrices, the computation of r-singular triplets all-at-once can be an issue, so that it is suitable to compute part of them. For instance, we can compute one-by-one the singular triplets with deflations by subtraction. Thus, the rank-1 approximation of A, namely σ 1 u 1 v H 1 , can be computed with some algorithm (e.g. Lanczo's [Comon 1990]) and deflated from A, resulting in the matrix Aσ 1 u 1 v H 1 , on which the same algorithm can be used in order to recover the second dominant singular triplet (σ 2 , u 2 , v 2 ). Hence, by repeating r-times both the deflation and the rank-1 approximation, we recover all r-singular triplets, thereby constructing the rank-r approximation of A.

The question that arises is whether a deflation property similar to that applied to matrices can be extended to tensors. In other words, we wonder if the best rank-1 approximation of a tensor T is a rank-1 component of its decomposition. The answer is no in general, and it is shown in [Stegeman 2010]. Therein, the authors show that subtracting a best rank-1 approximation from a tensor can even increase the rank of the resulting tensor. However, if the tensor T is symmetric and orthogonally decomposable, then the deflation procedure can be applied to [Zhang 2001].

Herein, we introduce the concept of iterative deflation, which is defined as the computational iterative technique whose purpose is the extraction of some mathematical features of some problem at any given iteration to be used at the next iteration. Notice that contrary to the preceding definition of deflation, we do not remove the mathematical features, but they are entries of a next problem. This techique can be used to compute the low rank approximation of tensors, as we will present in Section 3.1.1.

In order to understand the iterative deflation procedure, let T = R r=1 X r +E be a decomposition of a tensor T , where X r , 1 ≤ r ≤ R, are rank-1 tensors, and E is a residual tensor. At iteration 1, we compute the best rank-1 approximation of X 1 + E, namely X ′ 1 , and extract the residual

E ′ = X 1 + E -X ′ 1 .
At iteration 2, we add the residual E ′ to tensor X 2 , and compute again the rank-1 approximation but now of the tensor X 2 + E ′ , and extract another residual to be used at iteration 3, and so on. This simple procedure is actually the idea behind of DCDP algorithm, which is explained in details later.

The DCPD algorithm is similar to the nonnegative tensor factorization algorithm, called simple HALS NTF, described in [Cichocki 2009c]. Therein, the algorithm computes the nonnegative factors A r , 1 ≤ r ≤ R of a nonnegative tensor Y, where R is the rank of the approximating tensor. They use an alternating least square procedure for updating each column of each A r , with the negative entries projected to R + . This method can be extend to ordinary tensors by not performing the projection, in which case we call it simply HALS. These authors prefer to use the terminology hierarchical instead of iterative deflation. We do not agree with the former terminology since there is no actually hierarchical relation in their methods. Contrary to HALS, the DCPD algorithm does not use an alternating least square procedure to update rank-1 tensors at every iteration. For this task, we use rank-1 approximation algorithms, such as SeROAP, THOSVD and ST-HOSVD. The goal is to reduce quickly the norm of residuals according as the iteration increases, which not happens with HALS algorithm (see Section 3.3.2). Other tensor deflation algorithms can be found in [Cichocki 2009b, Phan 2015a, Phan 2015b].

In the following, we present a detailed description of DCPD algorithm in the context of low rank tensor approximation.

Description of the DCPD algorithm

The DCPD is an iterative deflation algorithm [da Silva 2015] that computes low rank tensor approximations for real or complex tensors. Contrary to the so-called Hierarchical ALS algorithm proposed in [Cichocki 2009c], rank-1 components are updated by using rank-1 approximation algorithms. The DCPD algorithm goes along the lines depicted in Alg.5.

In the first for loop, the rank-1 tensors X [1, 1], . . . , X [R, 1] are computed by successive rank-1 approximations and subtractions. Since the rank of the tensor does not decrease with subtractions in general [Stegeman 2010], a residual E[R, 1] is then produced. In the iterative process (repeat loop), a new rank-1 component, namely X [1, 2], is computed from the sum of the previous residual and X [1,1]. The tensor Y[1, 2] is updated within the if-else condition, and a new residual

E[1, 2] is produced with the subtraction Y[1, 2] -X [1, 2]
. By applying the same procedure to the other components, we update all R rank-1 tensors, so that another residual E[R, 2] is generated at the end of the second for loop. The second loop continues to execute until some stopping criterion is satisfied, and hence R rank-1 components are computed for T . If rank{T } ≤ R, then an exact decomposition can be recovered.

Study on iterative deflation and best rank-1 tensor approximations

In [da Silva 2015], we proved that the normalized residual ( E[R, l] ) l∈N >0 is a monotonically decreasing sequence when the best rank-1 approximation is assumed within DCPD. In this section, a thorough theoretical study is presented. Based on a geometric approach, we sketch an analysis of the input : T ∈ K I 1 ×I 2 ו••×I N : input data, R: rank parameter. φ: an algorithm computing a rank-1 approximation output:

X r ∈ K I 1 ×I 2 ו••×I N , for r = 1, . . . , R: rank-1 components Y[1, 1] ← T ; for r = 1 to R do X [r, 1] ← φ(Y[r, 1]); if r < R then Y[r + 1, 1] ← Y[r, 1] -X [r, 1]; else E[R, 1] ← Y[R, 1] -X [R, 1]; end end l ← 2; repeat for r = 1 to R do if r > 1 then Y[r, l] ← X [r, l -1] + E[r -1, l]; else Y[1, l] ← X [1, l -1] + E[R, l -1]; end X [r, l] ← φ(Y[r, l]); E[r, l] ← Y[r, l] -X [r, l]; end l ← l + 1; until some stopping criterion is satisfied; foreach r ∈ [1, . . . , R] do X r ← X [r, l]; end
Algorithm 5: DCPD algorithm convergence of the DCDP algorithm, including a conjecture that it converges to an exact decomposition with high probability when tensors within T (R) = {T ∈ T : rank T ≤ R} are distributed according to some absolutely continuous probability measure.

First, let us take a closer look at the 2D geometric interpretation of the DCPD algorithm. Figure 3.1 depicts the [r, l]-iteration for r > 1, so that γ[r, l] is the angle between the tensors E[r -1, l] and , andγ

X [1, l -1]. For r = 1, the residual E[r -1, l] can be just replaced with E[R, l -1] in the figure
[1, l] is then defined from E[R, l -1] and X [1, l -1].
The meaning of angle β and the gray area in the figure will be explained later.

Before stating some theoretical results on the DCPD algorithm, we present a fundamental lemma related to the error in rank-1 approximations of tensors of the form X + E, where X is a rank-1 tensor and E any other tensor, both with entries in some field K.

β γ[r, l] X [r, l -1] Y [ r , l ] E[r -1, l] Figure 3
.1: Visualization of the residual in an n-sphere for some iteration l of DCPD algorithm.

Lemma 3.2.1 Let X be a rank-1 tensor and φ the best rank-1 approximation operator. For any tensor

E, X + E -φ(X + E) ≤ sin(γ) E ,
where γ denotes the angle between E and X .

Proof: Let P X (X + E) be the orthogonal projection of X + E onto span(X ). Because φ(X + E) is a best rank-1 approximation of X +E, P X (X +E) cannot be a strictly better rank-1 approximation than φ(X + E). Thus,

X + E -φ(X + E) ≤ X + E -P X (X + E) .
On the other hand, X +E -P X (X +E) ⊥ X . Hence, we have X +E -P X (X +E) = sin(γ) E by using basic trigonometry. This concludes the proof.

The following results for the DCPD algorithm stems from the previous lemma.

l = R r=1 sin(γ[r, l]), the inequality E[R, l] ≤ c l E[R, l -1] holds.
Proof: By applying R -1 times the result of Corollary 3.2.2, we have

E[R, l] ≤ (sin(γ[R, l]) • • • sin(γ[2, l])) E[1, l] .
From Lemma 3.2.1, we know that

E[1, l] ≤ sin(γ[1, l]) E [R, l -1] . Thus, it follows that E[R, l] ≤ c l E[R, l -1] .
Notice that the same result brought by Proposition (4.4) in [da [START_REF] Da | [END_REF] can be deduced from Corollary 3.2.3 since 0 ≤ c l ≤ 1 for every iteration l, which implies the monotonic decrease of the sequence

{ E[R, l] } l∈N >0 . Corollary 3.2.4 If E[R, l] = E[R, l -1] , then c l = 1. Proof: Because c l ≤ 1 and E[R, l] = E[R, l -1] , one concludes directly from Corollary 3.2.3 that c l = 1.
From Corollary 3.2.4, we note that DCPD might not improve the estimation of the rank-1 components anymore for l ≥ l 0 > 1. And this may occur not only in the presence of noise. As a metter of fact, even for an almost orthogonal case c l ≈ 1, E[R, l] may tend to a stationary non-zero value as l increases. However, the DCPD algorithm converges to an exact decomposition if c l ≤ C, for all l > 1, and some constant C < 1. This will be subsequently detailed by means of a geometric approach.

Figure 3.1 can also be seen as the representation of an n-sphere

of dimension n = I 1 I 2 • • • I N -1 in K n+1 space. β is half the white cone angle defined in [0, π/2].
The direction of the rank-1 tensor X [r, l -1] defines the axis of the white cone and varies with r or l. Thus, the cone axis changes at every iteration [r, l] but the opening of the white cone remains the same. Under a condition on β, we can state an important proposition ensuring the convergence of the DCPD algorithm. Proposition 3.2.5 Let T be a tensor such that rank T ≤ R. An exact decomposition is recovered by the DCPD algorithm if and only if there exists for every (r, l) a half cone of angle β (in white in Fig. 3

.1), 0 ≤ β < π/2, such that β ≥ max l>1 min 1≤r≤R γ[r, l]. Proof: (⇐) For any iteration l > 1, take γ[r 0 , l] = min 1≤r≤R γ[r, l]. Notice that c l ≤ sin (γ[r 0 , l]) from Corollary 3.2.3. By hypothesis, sin (γ[r 0 , l]) ≤ sin(β), which implies that E[R, l] ≤ c l E[R, l -1] ≤ sin(β) E[R, l-1] . Because β is an upper bound for γ[r, l], l > 1, we have sin(β) E[R, 1] ≥ E[R, 2] =⇒ (sin(β)) l-1 E[R, 1] ≥ E[R, l] . Hence, when l → ∞, E[R, l] → 0.
(⇒) Let l 0 be some iteration such that l 0 > 1. Without loss of generality, assume E[R, l] = 0 for l ≥ l 0 (l 0 can be arbitrarily large). Then ( E[R, l] ) l∈N >0 is a strictly monotonically decreasing sequence for 1 < l < l 0 , otherwise the algorithm would converge to a nonzero constant for some iteration smaller than l 0 . Hence, for every (r, l) we can choose β, 0 ≤ β < π/2 such that β ≥ max l>1 min 1≤r≤R γ[r, l], and the proof is complete.

As a conclusion, if for a given iteration l, all tensors E[r, l], 1 ≤ r ≤ R, fall within the gray volume depicted in Fig. 3.1 (the complementary of the white cone), then the sequence E[r, l] does not tend to zero. Even if this gray volume can be made arbitrarily small, it is not of zero measure. In general, the best we can do is to study the convergence of the DCPD algorithm to an exact decomposition under some probabilistic conditions. Lemma 3.2.6 If tensors T are distributed within T (R) according to an absolutely continuous probability measure, then E[r, l] are absolutely continuous random variables.

Proof: Let D = [I 1 • • • I N ]
be a specific size of n-order tensors, and let

T (R) D = {T ∈ K I 1 ו••×I N : T ⊂ T (R) }. Because T (R) ⊃ T (R) D , any tensor T within T (R)
D is also distributed according to an absolutely continuous probability measure. Via the DCDP algorithm, each rank-1 component obtained in successive deflations is also in T (R) D . Hence, since the sum (subtraction) of continuous random variables does not affect the continuity, the residuals E[r, l] are also absolutely continuous random variables. Since the norm is a C 0 function in finite dimension, E[r, l] is also absolutely continuous.

For the next developments, let Z l = E[R, l] and define the following probability for some iteration L > 1:

F L [β] = P Z L ≤ sin(β)Z L-1 ≤ . . . ≤ (sin(β)) L-1 Z 1 . (3.1)
F L [β] can be viewed as the probability that residuals fall within at least one of the R white cones in every iteration l ≤ L.

The following proposition ensures a reduction of Z 1 by a factor smaller than (sin(β)) L-1 after L iterations with high probability, if a condition on the continuity of F L [β] is assumed.

Proposition 3.2.7 Let L be fixed. If ∃β 0 : β 0 ∈ [0, π/2) such that F L [β] is continuous on [β 0 , π/2], then ∀ε : ε ∈ (0, 1], ∃β ∈ [β 0 , π/2) such that F L [β] > 1 -ε. Proof: Since F L [π/2] = 1 and F L [β] is continuous on [β 0 , π/2],
the proof follows directly from the intermediate value theorem.

Although Z l , 1 ≤ l ≤ L, are absolutely continuous random variables and g m (β) = sin m (β) are continuous functions for all m ≥ 0, the continuity of F L [β] in β is not guaranteed due to the dependence of the random variables Z 1 , . . . , Z L (Z l depends on Z l-1 ). For example, for L = 2 and

Z 1 = 2Z 2 with probability 1, it is easy to check that F 2 [β] is not continuous at β = π/2. Indeed, lim β→π/2 -F 2 [β] = 0 whereas F 2 [π/2] = 1.
Now, let β l be the half-angle of the white cone at every iteration l. Contrary to the first definition of β, the angle value is not fixed anymore for all iteration [r, l], but only along the r-iteration. In this case, equation 3.1 becomes

F L [β L ] = P Z L ≤ sin(β L )Z L-1 ≤ L i=L-1 (sin(β i )) Z L-1 ≤ . . . ≤ L i=2 (sin(β i )) Z 1 . (3.2)
The following conjecture claims that there exists absolutely continuous distributions of tensors in T (R) such that the probability F l [β l ] tends to 1 as l → ∞, and at the same time the norm of residuals tends to 0, which is suitable for the convergence of the DCPD algorithm to an exact CP decomposition. Conjecture 3.2.8 There exists absolutely continuous probability measures µ for tensors T within T (R) for which the following holds:

(i). ∀ε : ε ∈ (0, 1], ∀l : l > 1, ∃β l : β l ∈ [0, π/2), such that F l [β l ] > 1 -ε. (ii). ∀l : l > 1, ∃β l : β l ∈ [0, π/2), such that l i=2 sin(β i
) is a strictly monotonically decreasing sequence converging to 0. Conjecture 3.2.8 shows that for some probability distributions of T within T (R) , the residual norm variable Z l → 0 as l → ∞ with high probability. Indeed, subsequent computer simulations support the existence of a uniform probability measure µ for the entries of tensors within T (R) , such that F L [β L ] ≈ 1 and Z L ≈ 0 for large values of L. This reinforces our conjecture.

Estimation of the CP decomposition

In this section, we present some numerical experiments using the iterative deflation algorithm DCPD. The rank-1 components are updated with three finite rank-1 approximation algorithms: SeROAP, ST-HOSVD and THOSVD. We also include in our simulations the standard ALS and the HALS algorithms for the following reasons: the former is known as the "workhorse" algorithm to compute low rank tensor approximations [Kolda 2009]; the latter is the original iterative deflation algorithm extended to complex tensors.

Complexity of the algorithms

Since DCPD is not a finite algorithm, the total complexity is unbounded. Therefore, we have chosen the number of multiplications by iteration as the complexity metric. The complexity is mainly dominated by the rank-1 approximation function φ, which is computed R times. Table 3.1 summarizes the number of operations per iteration of the DCPD algorithm for three finite rank-1 approximation algorithms: SeROAP, ST-HOSVD and THOSVD; and also those of the standard ALS and HALS algorithms.

In Table 3.2, we present the complexity/iteration of the algorithms for all scenarios that will be evaluated in next sections. We assume k = 4. Although in simulations we have defined the ordering vector p = [1 2 . . . N ] for DCPD-STHOSVD and DCPD-SeROAP, we do to need to take it into account here, because the tensors have equal dimensions in all scenarios. To obtain the real value of the complexity per iteration, one should multiply the values on the table by R. From Table 3.2, the order of complexity for the studied scenarios can be establish as ALS ≺ HALS ≺ DCPD-STHOSVD≺DCPD-SeROAP≺DCPD-THOSVD, where ≺ means "less complex than". Notice that DCPD-STHOSVD and DCPD-SeROAP algorithms have approximately the same complexities.

Algorithm complexity/iteration ALS O{3R

N j=1

I j } HALS O{(N + 2)R N j=1 I j } DCPD-THOSVD O{(2N k + 2)R N j=1 I j } DCPD-STHOSVD O{2kR N -1 n=1 N i=n I p i + 2R N j=1 I j } DCPD-SeROAP O{(2k + 2)R N -1 n=1 N i=n I p i }

Percentage of successful decompositions

Table 3.3 presents the percentage of successful decompositions of rank-3 tensors for the algorithms ALS, HALS, DCPD-THOSVD, DCPD-STHOSVD, and DCPD-SeROAP. The ALS and HALS algorithms are randomly initialized (In [Cichocki 2009c], the authors also initialize HALS with ALS, which case will be not considered here for fairness, otherwise DCPD could also be initialized in the same way). Noise is not considered in this case so that the performance is evaluated for the computation of an exact decomposition of 500 complex tensors whose real and imaginary parts are distributed uniformly in [-1, 1]. We consider that a decomposition is succeeded if the residual E ≤ 10 -6 , which is one of the stopping criteria for the algorithms. The algorithms also stop running when a maximum of 1000 iterations is attained.

The main conclusions from the table are the following:

• The iterative deflation algorithm DCPD presents similar performance for all finite rank-1 approximations, mainly for scenarios with larger dimensions;

• HALS presents the worst performance, mainly for the scenarios with smaller dimensions; • The performance of the DCPD algorithm reinforces Conjecture 3.2.84 ;

• The scenario 3 × 3 × 3 presents the poorer performance because the algorithms do not satisfy the residual condition E ≤ 10 -6 before the maximum number of iterations is attained. As a metter of fact, when the rank is of the same order of the dimensions, these algorithms need more iterations to converge. We note in this figure that average residual decreases more quickly for DCPD-SeROAP and DCPD-STHOSVD when compared to the other algorithms. For lower SNRs, most of the algorithms need a few iterations to converge to some stationary value. Notice that HALS takes more iterations to reduce the residual in all scenarios. For DCPD-THOSVD, we see that for an SNR of 0 dB and the iteration greater than 5, the curve lies above the other ones. The reason of that is because for most of tensors the residual E[R, l] does not attain its minimal value before the maximal number of iterations is attained. This fact is more clear explained in the next section.

Residual vs iteration

Convergence vs Iteration

Let χ l be the proportion of 5 × 5 × 5-rank-3 tensors over all 500 ones at which the simulation has finished before l iterations. For the same noisy scenarios evaluated in the previous section, we present in the figure below the metric χ l , 1 ≤ l ≤ 1000 for all aforementioned low rank-1 approximation algorithms. only 40% finished before the same number of iterations for the other algorithms. Notice also that 80% of simulations end up to 500 iterations with ALS, and approximately 65% with HALS and DCPD-SeROAP for the scenario with 0 dB. Among the iterative deflation algorithms, DCPD-SeROAP is the one that presents the more satisfactory performance. Notice that approximately along all scenarios, DCPD-SeROAP is at least as good as the other deflation methods. Apart from the last scenario, the HALS algorithm presents the worst performance. Yet, for the scenario with 0 dB, we note that the algorithms do not converge for some tensors. For instance, less than 60% of simulations converged before the maximum number of iterations for DCPD-THOSVD algorithm, which can explain the reason why the norm of the residual does not decrease in average for the same value obtained with the other algorithms, as described in Figure 3.2.

Residual vs rank

Now, we compare the algorithms for three SNR scenarios by computing µ = E{ E[R, l max ] / T }, where l max is the maximum iteration in which the stopping criterion is satisfied (or the last computed normalized residual). We perform simulations with different ranks for 5 × 5 × 5 complex tensors. The results are shown in Figure 3.4.

In the first scenario (SNR = 40 dB), we note that when the rank increases, the residual µ also increases. In this scenario, HALS presents the higher residual for all rank values. For the middle scenario (SNR = 20 dB), the performance of all algorithms is approximately the same for all ranks, with the residual µ in HALS slightly larger than those of the other algorithms. In the last scenario with lower SNR, the behavior of µ with the rank variation is the opposite of the first scenario. As a metter of fact, the higher is the rank, smaller is the residual µ. The figure also shows that all algorithms present approximately the same performance, except to DCPD-THOSVD, which is slightly worst. As we discussed in Chapter 2, the approximation error of SeROAP algorithm is degraded for higher order tensors. This drawback is reflected into DCPD-SeROAP algorithm, where we clearly note that it does not work for N ≥ 5, manly for high SNR (or noiseless scenario). DCPD-THOSVD also presents a poor performance, manly for the first scenario. On the other hand, DCPD-STHOSVD is robust to order variation, as ALS and HALS are. The latter, however, presents a slight increase of the residual for 6-th order tensors.

Residual vs order

Chapter Summary and Directions

We point out the following conclusions and future directions about iterative deflation:

• DCPD algorithm. we proposed an iterative deflation algorithm, called DCPD, based on finite rank-1 approximation algorithms to compute de CP decomposition and low rank tensor approximations. Instead of using an alternating update of factors for the rank-1 components as it is performed in HALS algorithm, we used the following finite rank-1 approximations: SeROAP, THOSVD and ST-HOSVD. The idea was the fast reduction of residuals in a few iterations, compared to standard methods like ALS and HALS itself. We have cogitated to employ an iterative algorithm to compute the best rank-1 approximation as the update of the rank-1 components in DCPD, in oder to reduce more quickly the residuals. However, the complexity per iteration would be unbounded in that case, so we did not implemented such an algorithm.

• Theoretical study on iterative deflation. We have showed in Corollary 3.2.3 that the norm of residuals E[r, l] is a monotonic decreasing sequence on r and l when the rank-1 components in DCPD algorithm are updated with the best rank-1 approximation. In Proposition 3.2.5 , we have proved that our iterative deflation method can recover an exact CP tensor decomposition if the residuals do not fall within an arbitrary small cone. Another result was the reduction of residuals with high probability when a defined probability function of residuals is continuous in some interval. This was showed in Proposition 3.2.7. Finally, we conjectured that there exist

I 1 = I 2 = • • • I N = 3.
absolute continuous probability measures for the entries of tensors with rank at most R such that the DCDP algorithm, with a best rank-1 approximation as the procedure to update the rank-1 tensors X [r, l], converges to an exact rank R CP decomposition with very high probability (see Conjecture 3.2.8). Hereafter, we intend to work on the proof of our conjecture. Some tools on algebraic geometry and measure theory may be starting points to tackle this problem.

• Performance of iterative deflation. We have evaluated five tensor algorithms on several fronts: complexity per iteration, percentage of successful decompositions, convergence versus iteration, residual versus iteration, residual versus rank, and residual versus order. Based on our experiments, we establish a subjective qualification factor from 1 to 5 for each of the aforementioned performance metrics, denoting the ranking of the algorithms. We set the value 5 for the algorithm considered the one with poorer performance for the studied scenarios, and 1 for the best one. This is summarized in the Table 3.4: Summary of performance factor for the algorithms.

Table 3.4 shows that ALS was the best algorithm whereas HALS the worst one in general for the evaluated scenarios. The DCPD-SeROAP algorithm is considered a robust algorithm in terms of convergence and residual versus iteration, and also in terms of rank variation. However, we saw that it is a very bad algorithm for tackling higher order tensors. Although it was considered the fourth algorithm in terms of percentage of successful decompositions, the performance was very close to ALS (in first position) for higher dimension scenarios (see table 3.2). As HALS, DCPD-THOSVD did not present good results. Finally, DCPD-STHOSVD had a moderate performance.

CHAPTER 4

Multivariate quadratic systems and the best rank-1 tensor approximation

In this chapter we propose a new approach based on rank-1 three-way tensor approximations that extracts solutions of multivariate quadratic polynomial systems in real and complex fields. This approach allows us to make use of standard iterative tensor algorithms [Comon 2009b], which in practice present satisfactory performances in time and convergence.

We start the chapter by presenting a brief introduction on the existing mathematical tools and some applications on polynomial optimization theory, in particular we focus on multivariate quadratic systems. The following sections deal with the formulation, discussion and experiments of our rank-1 tensor approach to tackle this problem.

First, we show the details of reducing a general multivariate quadratic system in the real field into a best rank-1 approximation problem. This is proved along four propositions describing the reduction (or equivalence) between intermediate problems. Our rank-1 tensor approach is discussed in the following section, where the advantages and drawbacks are presented.

Second, we extend our approach to the complex case in which only quadratic monomials of the form x * i x j , 1 ≤ i, j ≤ 2m, where m is the number of equations, are allowed for. In that case, the system is reduced to a conjugated partially symmetric rank-1 tensor approximation. We omit details of the reduction since they follow the same procedures described in the real case. In the sequel, we propose an alternating algorithm to compute the best conjugated partially symmetric rank-1 approximation. This algorithm, called CPS, is based on the stationary equations of a least square problem.

In the parts about examples and performance of the rank-1 tensor approach, we extract solutions of a multivariate quadratic system using ALS for real tensors, and CPS for complex tensors. We present an example of an over-determined real system with and without perturbation, where the solution is obtained by using different initializations methods. We also present an example of computing a solution of an under-determined real system. The last example is a complex over-determined system with a single complex solution that is extracted by CPS algorithm. We also perform a set of Monte Carlo simulations for several scenarios in real and complex case, where we show the computational time and the approximation error of estimations, in order to evince the usefulness of our approach.

Finally, we close this chapter drawing some conclusions and giving some directions to future work.

An introduction to multivariate quadratic systems

Multivariate polynomial optimization (MPO) is still an important topic in mathematics wherein a vast amount of applications takes place. Formally, a MPO problem consists in minimizing a polynomial cost function constrained by a set of polynomial equations and/or inequalities. In the last decades, some mathematical techniques have been applied to solve this problem, such as Gröbner bases, resultants and eigenvalues/eigenvectors of companion matrices [START_REF] David A Cox | Using algebraic geometry[END_REF]], semidefinite relaxations [Parrilo 2003, Lasserre 2001, Bucero 2014], and numerical homotopy [Li 1997, Verschelde 1999].

The complexity of a general polynomial optimization problem is NP-hard, even for quadratic polynomials with affine constraints [Pardalos 1991]. As a metter of fact, any MPO problem can be modified to a quadratic optimization (QO) problem by adding some variables and constraints.

A particular problem in this vein is the multivariate quadratic (MQ) polynomial systems, which are of great interest in various areas. In game theory, the Nash equilibria of a non-cooperative game between two players can be found by solving an MQ system [Lipton 2004]. In cryptography, the security of systems depends on the difficult to solve large quadratic systems in finite fields [START_REF] Courtois | [END_REF], Thomae 2012]. In [START_REF] Lebrun | [END_REF], Lin 2004], the authors present the design of multivariate filter banks modeled by quadratic polynomial systems. A last application lies in multilinear algebra [START_REF] Brachat | [END_REF], where the authors propose an efficient algorithm to decompose a symmetric tensor and show the equivalence between an existence condition of the decomposition and the solution of an MQ system.

Besides the techniques mentioned in the first paragraph, other methods can be employed to solve multivariate quadratique systems. For instance, we cite Newton-based and tensor-based1 algorithms [Dennis Jr 1996[START_REF] Schnabel | [END_REF], low-rank matrix recovery strategy [Davenport 2016], and symbolic computation [Grigoriev 2005].

The goal of this chapter is to present a rank-1 tensor approach to solve general MQ systems in the real field, and also to solve a particular extension in complex fields. In the real case, we show that a general MQ system can be reduced to a best rank-1 three-way tensor approximation problem [Lathauwer 2000[START_REF] Friedland | [END_REF], da Silva 2014], while the particular complex system can be reduced to a conjugated partially symmetric rank-1 three-way approximation. Despite the fact that most rank-1 tensor approximation algorithms do not guarantee to yield the optimum solution of the problem, it is known that standard algorithms, such as alternating least squares (ALS) and gradient-based methods [Comon 2009b] generally give satisfactory solutions in practice. Reducing the MQ problem to a rank-1 approximation can hence be attractive to compute a solution of MQ systems.

To compute a conjugated partially symmetric rank-1 approximation of some complex tensor T , we propose an alternating algorithm, called CPS. Some examples and experiments confirm the efficiency of rank-1 approximation methods to extract solutions of MQ systems.

From quadratic systems to best rank-1 approximations

Let P1 be the general system of quadratic polynomial equations given by P1:

       x T A 1 x + b T 1 x + c 1 = 0 . . . x T A m x + b T m x + c m = 0, where x ∈ R n , A j ∈ R n×n are symmetric matrices, b j ∈ R n and c j ∈ R, 1 ≤ j ≤ m.
Let also I be the ideal in R[x] generated by the set of polynomials of (P1), that is,

I = h 1 (x), h 2 (x), . . . , h m (x) , where h j (x) = x T A j x + b T j x + c j , 1 ≤ j ≤ m.
Thus, the set of solutions of (P1) is defined by the affine variety V (I) := {x ∈ R n | f (x) = 0, ∀f ∈ I}. Also define the optimization problem (P2) as follows P2:

       p ⋆ 2 = min y p(y) s.t. y = 1, where p(y) = m j=1 (y T Q j y) 2 , for y ∈ R n+1 and Q j =   A j b j /2 b T j /2 c j   ∈ R n+1×n+1 .
Define the set of global minimizers of (P2) as S P 2 = {y ∈ R n+1 | p(y) = p ⋆ 2 , y = 1}, and the subset SP 2 ⊆ S P 2 given by SP 2 = S P 2 ∩ (R n × R\{0}). In other words, SP 2 is the set of solutions of Problem (P2) such that y n+1 = 0. Let also N = {z ∈ R n+1 | z = y/y n+1 , ∀y ∈ SP 2 }. Proposition below connects Problems (P1) and (P2).

Proposition 4.2.1 If V (I) = ∅ then V (I) × {1} = N. Proof: By setting y = [x 1] T , it turns out that x T A j x + b T j x + c j = y T Q j y,
∀j ∈ {1, 2, . . . , m}. This shows that the set of solutions of the following system, equivalent to (P1), Since p(y) is by construction the sum of squares of the quadratic polynomials y T Q j y, 1 ≤ j ≤ m, and V (I) = ∅, it follows that p ⋆ 2 = 0. Thus, ∀ȳ ∈ V (I) × {1}, p(ȳ/ ȳ ) = 0, which implies that ȳ/ ȳ ∈ SP 2 and ȳ ∈ N. This proves that V (I) × {1} ⊆ N.

y T Q j y = 0, for 1 ≤ j ≤ m y n+1 = 1 (4.
On the other hand, ∀y ∈ SP 2 , p(y/y n+1 ) = 0, which implies that (y/y n+1 ) T Q j (y/y n+1 ) = 0, 1 ≤ j ≤ m. Thus, y/y n+1 is solution of (4.1), which leads to N ⊆ V (I) × {1}, and the proof is complete. Proposition 4.2.1 says that all solutions of System (P1) can be recovered from Problem (P2) by taking the n first elements of y/y n+1 , ∀y ∈ SP 2 .

Remark 1 When b j = 0, c j = 0, ∀j ∈ {1, 2, . . . , m}, V (I)\{0} can be viewed as a projective variety, since x ∈ V (I) =⇒ αx ∈ V (I) ∀α ∈ R. Therefore, it is enough to solve the system for x = 1. By setting y = x, and Q j = A j , 1 ≤ j ≤ m, all nontrivial solutions of (P1) can be recovered by solving (P2).

Now, let the following optimization problem

P3:

       p ⋆ 3 = max y (y ⊠ y) T M (y ⊠ y) s.t. y = 1, where M ∈ R (n+1) 2 ×(n+1) 2 is a semidefinite positive matrix given by M = λ max I -M Q , where M Q = m j=1 Q j ⊠ Q j ,
I is the identity matrix, and λ max is the largest eigenvalue of M Q .

Proposition 4.2.2 (P2) ⇐⇒ (P3).

Proof: We have

p(y) = m j=1 (y T Q j y) 2 = m j=1 (y T Q j y) ⊠ (y T Q j y) = m j=1 (y ⊠ y) T (Q j ⊠ Q j )(y ⊠ y) = (y ⊠ y) T   m j=1 Q j ⊠ Q j   (y ⊠ y) = (y ⊠ y) T M Q (y ⊠ y).
Hence, it follows that

min y =1 (y ⊠ y) T M Q (y ⊠ y) = max y =1 -(y ⊠ y) T M Q (y ⊠ y) ⇐⇒ max y =1 λ max -(y ⊠ y) T M Q (y ⊠ y) = max y =1 λ max (y ⊠ y) T (y ⊠ y) -(y ⊠ y) T M Q (y ⊠ y) = max y =1 (y ⊠ y) T (λ max I -M Q )(y ⊠ y) = max y =1 (y ⊠ y) T M (y ⊠ y).
Remark 2 Problem (P3) can be viewed as a symmetric best rank-1 approximation of a fourth order tensor. Indeed, it can be written as min

y =1 M -y ⊗ y ⊗ y ⊗ y 2 , where M ∈ R n+1×n+1×n+1×n+1
is a symmetric tensor constructed from M . This problem is difficult to solve and some standard algorithms, such as ALS, could be adapted to deal with it. However, ignoring the symmetry generally leads to faster convergence. For this reason, reducing the system to a standard three-way rank-1 approximation problem is attractive under the algorithmic point of view.

In the following, we establish a link between (P3) and a rank-1 three-way tensor approximation problem. Let (λ k , v k ), 1 ≤ k ≤ K, K = rank M , be the k-th eigenpair of M . Due to the semidefiniteness of M , we can define real vectors

t k = √ λ k v k ∈ R (n+1) 2 . Set T k = Unvec(t k )
. The matrices T k can be viewed as frontal slices of a three-way tensor T ∈ R n+1×n+1×K [Kolda 2009], as illustrated in Figure 4.1.

T = T 1 T 2 T K Figure 4
.1: Frontal slices of a three-way tensor T ∈ R n+1×n+1×K . Now, define the following rank-1 tensor approximation problem with two identical factors.

P4:

       p ⋆ 4 = min y,w T -y ⊗ y ⊗ w . s.t. y = 1.
The following proposition holds.

Proposition 4.2.3 (P3) ⇐⇒ (P4).

Proof: Expand matrix M into its eigencomponents as:

M = K k=1 λ k v k v T k = K k=1 t k t T k
, and write the objective function of (P3) as

(y ⊠ y) T M (y ⊠ y) = K k=1 (y ⊠ y) T t k t T k (y ⊠ y). Set w k = (y ⊠ y) T t k . Problem (P3) can be rewritten as max (y,w)∈C K k=1 w k t T k (y ⊠ y) (4.3)
where

C = {(y, w) ∈ R n × R K , : y = 1, w = [w 1 w 2 . . . w K ] T , w k = (y ⊠ y) T t k , 1 ≤ k ≤ K}. Next, we have the equivalence max (y,w)∈C K k=1 w k t T k (y ⊠ y) ⇐⇒ min (y,w)∈C K k=1 t T k t k -w k t T k (y ⊠ y). Yet, since (y, w) ∈ C, it turns out that t T k t k -w k t T k (y ⊠ y) = t k -w k (y ⊠ y) 2 , 1 ≤ k ≤ K. Therefore, Problem (P3) is equivalent to min (y,w)∈C K k=1 t k -w k (y ⊠ y) 2 (4.4)
or, to minimize the Lagrangian:

Υ(y, w) = K k=1 t k -w k (y ⊠ y) 2 + η( y -1) + K k=1 γ k (w k -(y ⊠ y) T t k ),
where γ k , 1 ≤ k ≤ K, and η are Lagrangian multipliers.

We now need to relax the constraint on w k . Stationary points w.r.t. w k yield ∂Υ(y, w)

∂w k = -2(t k -w k (y ⊠ y)) T (y ⊠ y) + γ k = 0 =⇒ =⇒ (-t k + w k (y ⊠ y)) T (y ⊠ y) + γ k /2 = 0.
Now using all the constraints defined in C (i.e. stationary equations w.r.t. γ k and η), the latter equation leads to γ k = 0, ∀k, which means that the constraints on w k do not need to be imposed. Hence we end up with the simplified Lagrangian:

Υ(y, w) = K k=1 t k -w k (y ⊠ y) 2 + η( y -1). (4.5)
To complete the proof, note that by reshaping the entries of the vectors t k using the Unvec operator, (5) can still be written with the help of matrix slices:

min y =1,w K k=1 T k -w k yy T 2 (4.6)
or in its tensor representation min

y =1,w T -y ⊗ y ⊗ w 2 .
Problem (P4) is a best partially symmetric rank-1 approximation, which cannot be directly tackled by standard tensor approximation algorithms. However, this problem can still be changed in order to drop the constraint of identical factors.

According to [START_REF] Friedland | [END_REF]], the best rank-1 approximation of a real symmetric tensor with respect to a subset of indices (modes) can be chosen symmetric with respect to this subset (see Theorem 1 therein). Thus, if tensor T in Problem (P4) were symmetric with respect to the first two modes, we would not need to impose any constraint on the factors, so that the best rank-1 approximation could be chosen partially symmetric.

In general, T is not symmetric with respect to the first two modes. However, the linear transformation T ′ k = (T k + T T k )/2 on each k-th slice T k of T , ensures that any best partially symmetric rank-1 approximation of T is actually a best rank-1 approximation of a tensor T ′ with slices T ′ k (which is therefore symmetric with respect to the first two modes). In other words, Problem (P4) can be reduced to the problem P5:

       p ⋆ 5 = min x,y,w T ′ -x ⊗ y ⊗ w .
s.t. y = 1. Before proving this, we need the following lemma.

4.2.4 Let f (X, α) = K k=1 A k -α k X 2 and g(X, α) = K k=1 (A k +A T k )/2-α k X 2 with X, A k ∈ R n×n , 1 ≤ k ≤ K, α = [α 1 α 2 . . . α K ] T . If X is symmetric, then arg min X,α f (X, α) = arg min X,α g(X, α). Proof: Let h(X, α) = K k=1 A T k -α k X 2 . Since X is symmetric, it follows that arg min X,α f (X, α) = arg min X,α h(X, α) = arg min X,α f (X, α) + h(X, α).
By expanding the functions,

f (X, α) + h(X, α) = 2 X 2 K k=1 α 2 k -2 K k=1 α k trace{(A k + A T k )X} + 2 K k A k 2 , and g(X, α) = X 2 K k=1 α 2 k - K k=1 α k trace{(A k + A T k )X} + 1 4 K k=1 A k + A T k 2 .
Thus, f (X, α) + h(X, α) = 2g(X, α) + c, for some constant c, and the proof is complete.

Proposition 4.2.5 Let S P 4 and S P 5 be the sets of rank-1 tensors that minimize Problems (P4) and (P5), respectively. Then S P 4 ⊆ S P 5 .

Proof: Since yy T is symmetric, we can apply Lemma 4.2.4 to the cost function of the slice form of Problem (P4) to show that arg min y,w

K k=1 T k -w k yy T 2 = arg min y,w K k=1 T ′ k -w k yy T 2 .
That means that any rank-1 tensor X = y ⊗ y ⊗ w ∈ S P 4 is a partially symmetric rank-1 tensor that minimizes (P5), which is in turn a best rank-1 approximation from Theorem 1 in [START_REF] Friedland | [END_REF]]. Therefore, ∀X ∈ S P 4 , X ∈ S P 5 =⇒ S P 4 ⊆ S P 5 .

Discussion on the real tensor approach

The goal of this section is to point out the advantages and disadvantages of using rank-1 tensor approximation algorithms to solve quadratic polynomials systems in the real field. We start describing the drawbacks, but with the intention of emphasizing some aspects that mitigate part of the deficiencies of the approach.

• Apart from the case described in Remark 1, Problems (P1) and (P2) are not equivalent in general. Thus, solving Problem (P5) does not ensure that the factor y obtained with a rank-1 tensor approximation method satisfies the condition y n+1 = 0. In this case, the standard ALS algorithm can be adapted by setting y n+1 = 1 at every iteration in order to avoid undesirable solutions. However, this adaptation can compromise the convergence of the algorithm.

• Problems (P4) and (P5) are not equivalent either, which means that a non symmetric rank-1 approximation can be a minimizer of (P5). However, routine simulations have shown that ALS delivers (partially) symmetric rank-1 approximations for (partially) symmetric tensors.

• The complexity for constructing tensor T is of order O{mn 4 + n 6 }, which comes mainly from the computation of M and the eigenpairs (λ k , v k ), 1 ≤ k ≤ K.

• Solutions delivered by iterative tensor approximation algorithms have an initialization dependence, which means that such methods do not guarantee that the optimum solution of the best rank-1 approximation problem is attained. However, the solution delivered by Problem (P5) allows us to test whether x satisfies or not System (P1). If it does not, we can run again the algorithms with other initializations.

• Contrary to (numerical) algebraic geometry methods, such as those described in [START_REF] David A Cox | Using algebraic geometry[END_REF], Verschelde 1999], tensor approximation algorithms can find only a solution of (P1). As a metter of fact, the latter should run several times with different initializations in order to find other possible solutions.

On the other hand, we can take advantage of tensor approach for solving quadratic polynomial systems. We emphasize the following points:

• Numerical experiments for generic systems have shown that y n+1 = 0 is always satisfied for the standard ALS method with random initializations.

• The best rank-1 tensor approximation problem is well-posed, which means that Problem (P5) has always a partially symmetric minimizer.

• When System (P1) is inconsistent, the tensor approach gives a solution that fits (P1). Thus, tensor approach can be useful when a system becomes inconsistent due to perturbations.

• Tensor approach can be applied to square, under-determined and over-determined systems, and zero-dimensional and positive-dimensional systems.

• ALS algorithm is easy to implement and converges fast for generic tensors approximations.

Additionally, ALS has also locally convergence properties. When the starting point is in the neighborhood of a global solution of (P5) then the iterates of ALS converge linearly to this solution [Uschmajew 2012].

• Although tensor approximation algorithms can not find all solutions of (P1), they can deal with systems up to 20 variables/equations in a reasonable time, contrary to algebraic geometry methods.

• Newton-based algorithms can not be applied to ill-conditioned systems, which is not a drawback to tensor approximations algorithms.

• Although the construction of T ′ depends on the number of equations m, Problem (5) itself does not depend on m, which is an advantage compared to other methods that take the number of equations into account, such as those in [START_REF] Schnabel | [END_REF], Grigoriev 2005].

Extension to the complex field

The theory developped in Section 4.2 for real tensors can be extended to the following complex system:

P6:

       x H A 1 x + b T 1 x + c T 1 x * + d 1 = 0 . . . x H A m x + b T m x + c T m x * + d m = 0 , where x ∈ C n , A j ∈ C n×n , b j ∈ C n , c j ∈ C n , and d j ∈ C, for j = 1, . . . , m.
Problem (P6) is not a general complex quadratic system since there is no terms of the form

α ij x i x j nor β ij x * i x * j , for α ij , β ij ∈ C, 1 ≤ i, j ≤ n.
As a metter of fact, general quadratic complex systems cannot be formulated as rank-1 tensor approximation problems with the same approach detailed before, unless the systems should be separated in real and imaginary parts, in which case both the number of variables and equations are doubled. Besides, rewrite a complex system by taking separately the real and imaginary parts can be an unnecessary and annoying task, manly if the number of equations/unknowns is large enough. On the other hand, without doubling the number of variables, Problem (P6) can be reduced to a conjugate partially symmetric rank-1 tensor approximation problem, which can be solved by performing an algorithm called CPS (see Section 4.5).

The following lemma shows that each polynomial in (P6) is equivalent to two polynomial equations composed of Hermitian matrices. Lemma 4.4.1 Let y = [x T 1] T ∈ C n+1 . For every j ∈ {1, 2, . . . , m},

x H A j x + b T j x + c T j x * + d j = 0 ⇐⇒ y H Q j y = 0 y H Q j+m y = 0,
where

Q j =   A j + A H j b * j + c j b T j + c H j d j + d * j   and Q j+m = 1 i   A j -A H j c j -b * j b T j -c H j d j -d * j   , with i = √ -1.
Proof: The proof follows directly from the real and imaginary parts of each equation in (P6). Hence, we have

ℜ{x H A j x + b T j x + c T j x * + d j } = 0 ⇐⇒ x H (A j + A H j )x + (b T j + c H j )x + (c T j + b H j )x * + d j + d * j = 0 ⇐⇒ y H Q j y = 0,
and

ℑ{x H A j x + b T j x + c T j x * + d j } = 0 ⇐⇒ 1 i x H (A j -A H j )x + (b T j -c H j )x + (c T j -b H j )x * + d j -d * j = 0 ⇐⇒ y H Q j+m y = 0.

We can define p(y) = 2m

j=1 |y H Q j y| 2 , and write the following optimization problem P7:

       p ⋆ 7 = min y p(y) s.t. y = 1.
(P6) and ( P7) are connected problems as (P1) and ( P2) are for the real case, which means that the solutions of System (P6) are minimizers of (P7). Yet, by defining

M Q = 2m j=1 Q * j ⊠ Q j , we can draw up the equivalent problem P8:        p ⋆ 8 = max y (y * ⊠ y) H M (y * ⊠ y) s.t. y = 1,
where M ∈ C (n+1) 2 ×(n+1) 2 is a semidefinite positive matrix given by M = λ max I -M Q , with λ max the largest eigenvalue of M Q . Thus, the equivalence between (P2) and (P3) is extended to the complex case with (P7) and (P8).

Remark 3 We could define p(y) = 2m j=1 (y H Q j y) 2 since Q j are Hermitian matrices, which would lead to a cost function for (P8) given by (y ⊠ y) H M (y ⊠ y), with M Q = 2m j=1 Q j ⊠ Q j . However, as we will see in Section 4.5, it would imply to conceive a more complex algorithm to tackle complex tensors.

The same reasoning used to construct the partially symmetric rank-1 approximation (P4) from (P3) can be applied here. We omit again the details because they are identical to those of the real case. Thus, (P8) is equivalent to the rank-1 tensor problem given by P9:

       p ⋆ 9 = min y,w T -y ⊗ y * ⊗ w . s.t. y = 1,
with T ∈ C n+1×n+1× K , and K ≤ rank{ M }. Problem (P9) is actually a conjugated partially symmetric rank-1 approximation.

Theorem 1 in [START_REF] Friedland | [END_REF]] cannot be extended to conjugated partially symmetric complex tensors 2 . In other words, we cannot ensure that a best rank-1 approximation of a such tensor can be 2 Even to partially symmetric complex tensors. chosen conjugated partially symmetric. Hence, we cannot derive a best rank-1 approximation problem as (P5), so that we should directly tackle Problem (P9). For this, we propose an algorithm called CPS, based on the stationary equations of (P9). This algorithm is described in the next section.

Conjugated partially symmetric rank-1 approximation

In order to compute the best conjugated partially symmetric rank-1 approximation of a tensor in the complex field, we propose an alternating algorithm, based on the stationary equations of the Lagrangian of the matrix slice form of Problem (P9), which follows

L = µ y 2 -1 + K k=1 T k -w k yy H 2 .
The stationary equations are

                         ∂L ∂µ = 0 =⇒ y 2 = 1. ∂L ∂w * k = 0 =⇒ w k = y H T k y y 4 , ∀k. ∂L ∂y = 0 =⇒ K k=1 w * k T k + w k T H k -2|w k | 2 y 2 I y = µ * y. (4.7)
By using the first stationary equation of (4.7), it follows that

w k = y H T k y and K k=1 w * k T k + w k T H k -2|w k | 2 I y = µ * y. (4.8)
By multiplying both sides of the second equation of (4.8) by y H , we obtain µ = 0. Hence,

K k=1 w * k T k + w k T H k y = γy, (4.9)
where γ = 2 K k=1 |w k | 2 . Notice that the Lagrangian can be rewritten as L = T 2γ, which reveals that to minimize L, we need to maximize γ. The maximization of γ in equation ( 4.9) and the equations w k = y H T k y, 1 ≤ k ≤ K can be solved alternately by the algorithm depicted in Alg.6. As a metter of fact, for some y, we compute all w k and substitute them in equation (4.9). Thus, the update of y is performed by taking the dominant eigenvector of

K k=1 w * k T k + w k T H
k , since γ must be maximized. Some remarks about the CPS algorithm:

input : T ∈ C n+1×n+1× K . output: y, w k , with k ∈ {1, K}. Initialize y; repeat V ← 0; for k = 1 to K do w k ← y H T k y; V ← V + w * k T k + w k T H k ; end
Set y as the eigenvector of V whose eigenvalue is maximum; until some stopping criterion is satisfied;

Algorithm 6: CPS algorithm

• As described in Remark 3, if the definition of p(y) were p(y) = 2m j=1 (y H Q j y) 2 , (P9) would be a partially symmetric problem. However, the Lagrangian would be

L = µ y 2 -1 + K k=1 T k -w k yy T 2 ,
in which case the stationary equation (4.9) would not be an eigevalue problem on y. The problem could still be adapted to an eigenvalue problem but with a more complex matrix;

• The complexity per iteration of CPS algorithm is given by O{2( K + ρ)(n + 1) 2 }, where ρ is the number of steps3 used in Lanczo's algorithm to compute the dominant eigenvector. For large n, K is also large, so that ρ can be dropped of the complexity, since it is generally a very small value;

• The CPS algorithm can also be used to tackle the partially symmetric rank-1 approximation problem (P4). In this case, the equations on w k and γ becomes w k = y T T k y and K k=1 w k (T k + T T k )y = γy.

Examples on real and complex fields

In this section, we present some examples of solving multivariate quadratic systems. In the real case, we use the standard ALS algorithm whereas for complex systems the algorithm used is the proposed CPS. For both algorithms the stopping criterion is satisfied when either the maximum number of iterations (itr max = 10000) or the threshold

| T -T [itr] -T -T [itr-1] | ≤ 10 -6
for an estimate T at iteration itr, is attained.

Since CPS is also an initialization-dependent algorithm, we test four initialization algorithms: SeROAP, THOSVD and ST-HOSVD (described in Chapter 2), and a conjugated random initialization method. The latter consists of the selection of the tensor that minimizes the approximation error from a sample of rank-1 tensors with conjugated factors, generated according to a uniform probability distribution with entries in [-1, 1]. Let L be the size of this sample and y (l) , w (l) , l ∈ {1, 2, . . . , L}, be the factors of the l th -generated tensor. The factor y is initialized in CPS algorithm as follows y = arg min y (l) Ty (l) ⊗ y (l) * ⊗ w (l) .

The initializations can also be applied to the ALS algorithm, so that the conjugated random initialization becomes a symmetric random initialization. The experiments have performed with Matlab c in an MAC OSX 10.8.5 with processor 3.2 GHz Intel core i5, and memory 8G 1600MHz DDR3.

Example 1: Let the overdetermined real system with 5 unknowns and 7 equations

                           x 2 1 + x 2 x 3 + x 5 = 3 x 2 2 + x 2 x 5 -2x 3 x 4 + x 2 5 = 5 x 1 x 3 + x 1 x 5 -x 2 4 + x 2 -x 4 = 0.25 x 1 x 2 + x 2 3 + x 3 x 5 -3x 2 = -1 x 1 x 4 -3x 1 x 5 + 2x 2 2 + x 1 -x 4 = -5 x 2 1 + x 2 3 + x 2 4 -x 2 5 + x 4 x 5 = -0.75 x 2 2 + x 1 + x 3 = 0 , (4.10)
where x = [1 0 -1 0.5 2] T is the real solution. The tensor T ′ constructed from the system has dimensions 6 × 6 × 21. The slices k ′ such that T k ′ ≤ 10 -6 , for k ′ ≤ 36 were ignored.

Table below shows the error of the estimation xx , and the time of execution to compute the solution of (4.10), using ALS for four different initialization methods. The time to initialize the algorithm is not considered in the table. Table 4.1: Error estimation for different initializations using ALS algorithm for a (5, 7) system.

For all initializations, ALS delivers the approximative solution x = [1 0 -1.0003 0.4998 2] T .

Example 2: For the same overdetermined real system, consider the following diagram representation:

x y f + + ỹ ε where f = [f 1 f 2 . . . f 5 ]
T , and y = [y 1 y 2 . . . y 5 ] T , such that f i (x) = y i is the i-th equation in (4.10). We estimate x for a given noisy output ỹ = y + ε, where ε is a perturbation vector distributed according to a gaussian measure with variance σ 2 . We use SeROAP to initialize ALS for all simulations.

Table below shows the average estimation xx for a sample of 1000 realizations with three different values of σ. We notice that the average estimation is approximately equal to σ, which evinces that our rank-1 tensor approach to solve the system is robust under perturbations. Table 4.2: Average error estimation using ALS algorithm for a (5, 7) noisy system.

Example 3: Let the under-determined real system with 6 unknowns and 3 equations

       3x 2 1 + x 1 x 2 -3x 2 2 -x 2 x 6 = 1 -x 2 3 + 2x 4 x 5 + x 1 + x 2 = 3 x 2 1 + 3x 2 4 + 2x 2 5 + x 2 6 -2x 2 = 4 , (4.11)
Table 4.3 shows that one solution is extracted for ALS for all initializations. Indeed, max j |x T Q j x| ≤ 1.027•10 -4 , 1 ≤ j ≤ 2m, for all initializations. The solution x = [1.1843 1.8156 0 0 0 -2.4957] T was obtained with the finite algorithms SeROAP, ST-HOSVD, and THOSVD. While with the random initialization we have extracted the solution x = [1.6030 1.5943 -1.1678 0.9735 0.5991 1.0286] T . Initialization max j |x T Q j x| Time ALS (seconds) Random (L = 50000) 1.8113e-5 0.2187 THOSVD 8.8285e-5 0.8268 ST-HOSVD 9.1748e-5 0.8159 SeROAP 1.0270e-4 0.8632 Table 4.3: Error estimation for different initializations using ALS algorithm for a (6, 3) system.

Example 4: Let the overdetermined complex system with 3 unknowns and 4 equations below

           2|x 1 | 2 + x * 1 x 2 + ix * 3 = 1 + i -2x 1 x * 3 + x * 1 x 3 + |x 2 | 2 -x 2 = 1 -2i -i|x 3 | 2 + x * 2 x 3 -(1 + i)x 1 + x * 2 = 1 -|x 1 | 2 + 3|x 3 | 2 + 3ix * 2 = -1 , (4.12)
where x = [ii 1] T is the complex solution. As mentioned in Section 4.4, one way to solve this system is to substitute each unknown by its complex form a + bj and solve a real system with twice as many equations. However, this transformation introduces additional calculations, thereby we go along to solve the problem directly using the rank-1 tensor approach.

The tensor T constructed from the system has dimensions 4 × 4 × 16. As in the real case, some slices of T were not taken into account because their norms are negligible. The estimation and time performance using the CPS algorithm are presented in table below. For all initializations, the algorithm extracts the solution [ii 1] T .

Initialization

xx Time CPS (seconds)

Random (L = 50000) 9.5133e-6 1.3620 THOSVD 9.5666e-6 1.8400 ST-HOSVD 9.9688e-6 1.8066 SeROAP 1.2052e-5 1.6464 Table 4.4: Error estimation for different initializations using CPS algorithm for complex systems.

Performance evaluation of rank-1 tensor approach

In this section, we present some numerical experiments at which a solution of generic quadratic systems are extracted using the standard ALS and CPS algorithms for real and complex rank-1 tensor approximations, respectively. For the real case, a total of 500 systems for each scenario (n, m) was generated with the entries of A j , b j , c j , 1 ≤ j ≤ m, distributed according to a uniform measure in [-1, 1]. For complex systems, the entries of the real and imaginary parts of A j , b j , c j , d j , 1 ≤ j ≤ 2m, are also uniformly distributed in [-1, 1]. We have employed the symmetric (or conjugated symmetric) random initialization.

Table 4.5 summarizes the percentage p of real systems in which a solution was successfully obtained, and the average computational time for initialization, construction of T ′ , and the rank-1 approximation itself. We assume that a solution satisfies a system if

max(|y T Q 1 y|, |y T Q 2 y|, . . . , |y T Q m y|) ≤ 0.001.
All time measurements are denoted in seconds.

The results confirms that a solution of an MQ system can be extracted in a few seconds with ALS algorithm. Notice that the computational time of the ALS algorithm is larger for systems having variables as many as equations. For instance, compare scenarios (4,4) and (6, 6) with scenarios (5, 4) and (7, 9), respectively. Table 4.6 shows the performance of CPS for some of the aforementioned scenarios but in the complex field. We use the same metrics to evaluate our algorithm, except the condition in which a solution satisfies a system that now is given by max(|y The performance of CPS algorithm to compute a conjugated partially symmetric rank-1 approximation and thus extract one solution of System (P6) is not good as in the real case. We eliminate the scenarios where p < 70% so that only four scenarios are presented. Moreover, the initialization and CPS times are large compared to the real scenarios. As before, the square system also presents larger computational time among all scenarios.

H Q 1 y|, |y H Q 2 y|, . . . , |y H Q 2m y|) ≤ 0.

Chapter Summary and Directions

We conclude this chapter by presenting the main points observed in the results and also we given some directions for future work.

• Reduction of multivariate quadratic systems to best rank-1 three-way tensor approximation problems. We showed that a general real MQ system can be reduced to a best rank-1 approximation of a three-way real tensor. The proof was delineated along four propositions in Section 4.2. This approach arises as a new mathematical tool to tackle MQ systems when one solution is suitable. Some discussions about the rank-1 tensor approach were outlined in Section 4.3. Therein, we saw that the rank-1 approximation itself does not depend on the number of equations, which means that the dimensions of a tensor T , constructed from a system of m equations, is independent of m. We also saw that there is no constraints of using our approach, contrary to some other methods, such as the Newton-based one in [Dennis Jr 1996], where the system must be square and nonsingular, or [Verschelde 1999], where the system must be squared before applying the homotopy continuation. On the other hand, the main limitation of the rank-1 tensor approach is the nonequivalence between problems (P1) and (P5). Although ALS extracts one solution of an MQ system in general, we do not have yet a complete answer for the following questions: (i) which cases y n+1 = 0 in (P2)? (ii) What is the set of tensors at which S P 4 = S P 5 ? For (i), if we map all non trivial solutions of (P1) when c j = 0, 1 ≤ j ≤ m, then all solutions of (P2) such that y n+1 = 0 are also mapped. We intend to study and perform this mapping in a future work. For (ii), we claim that the set of partially symmetric tensors whose rank-1 approximations are also partiality symmetric has volume measure equal to 1, but no proof was presented here. Nevertheless, the iterative ALS algorithm applied to a partially symmetric tensor delivers a partially symmetric rank-1 approximation in general. Thus, from an algorithmic point of view, the nonequivalence between problems (P1) and (P5) does not impose any issue.

• Extension to complex quadratic systems. we extended our rank-1 tensor approach to complex MQ systems with no monomials of the form x i x j or x * i x * j , for 1 ≤ i, j ≤ 2m. Indeed, we saw that such complex systems can be reduced to conjugated partially symmetric rank-1 three-way approximation problems. As discussed in Section 4.4, any complex MQ system is equivalent to a real MQ system with twice the number of variables and equations. However, with our rank-1 tensor approach, we kept the same number of variables as in the real case.

• Proposition of the CPS algorithm. In order to deal with the conjugated partially symmetric rank-1 approximation problem, we proposed an alternating algorithm called CPS, based on the stationary equations of (P9), which can also be applied to real tensors.

• Performance of the rank-1 tensor approach. Numerical experiments showed that ALS algorithm extract a solution even for real systems of 20 variables/unknowns in a reasonable time, which is not possible with standard algebraic geometry methods. In the complex case, we saw that CPS algorithm is also capable of obtaining conjugated partially symmetric rank-1 approximations but the performance is poorer than that obtained for real systems. We presented some examples of real with/without perturbations, an example of complex system, and also a set of Monte Carlo experiments for different sizes of systems to evince the efficiency of the rank-1 tensor approximation approach.

CHAPTER 5

Alternating projections on orthogonal CP tensor decomposition

This chapter concerns the exact decomposition of tensors with one semi-unitary factor matrix. Firstly, we present a brief discussion on the alternating projection method and the orthogonal tensor decomposition. We emphasize some applications and theoretical results in the literature.

Next, the description of our algorithm, called CAPD, combining the alternating projection method and the deflation procedure is presented in details. We develop the algorithm based on the set (manifold) of unit rank-1 matrices and a linear subspace. In the sequel, we perform a convergence study of the CAPD algorithm, whose main features are the extraction of the right factors, and the convergence under the transversality and non-tangential concepts. In this part, we resort to the main theorems of [Lewis 2008, Andersson 2013] pointing out the convergence to the intersection of the manifolds. These results are contextualized to the column-wise tensor decomposition, so that the convergence to the exact decomposition is ensured.

Finally, we draw some numerical experiments evincing the efficiency of our proposed algorithm. We finish this chapter presenting some conclusions and putting out future works.

Main contributions of this chapter:

• CAPD algorithm. We propose an algorithm, called CAPD, that combines the alternating projection method and the deflation of rank-1 tensors to solve the exact decomposition of a three-way tensor with one semi-unitary factor matrix. The alternating projection part consists of iterative projections onto the manifold of unitary rank-1 matrices and a linear subspace, aiming the computation of one column of each factor matrix. Thus, we can construct a rank-1 component of the tensor and deflated this component from the original tensor, so that the rank is reduced by 1. The process is repeated until all components can be extracted.

• Extraction of the right rank-1 component. We proof that a right rank-1 component is extracted by our proposed algorithm if the alternating projection part converges to a point in the intersection of the manifolds. This is shown in Theorem 5.4.3.

• Convergence under the transversality concept. Under some conditions, we prove that if points in the intersection of the manifolds involved in the column-wise orthogonal tensor decomposition are transversal, and the starting point of CAPD algorithm is close to that intersection, then the linear convergence of our algorithm is ensured. This results is presented in Proposition 5.4.7. However, the transversality concept applied to our problem imposes some limitations on the dimensions of the manifolds. Actually, apart from a special case, only a few components of the tensor can be extracted.

• Convergence under the non-tangential concept. We prove a stronger result based on the nontangential concept. For the tensor decomposition problem with an orthogonal factor matrix, we show that CAPD algorithm converges to a point in the intersection of the manifolds, when the starting point is close enough to this intersection, and the manifolds are non-tangential at any point in the intersection. This approach mitigates the limitation of the transversality concept since now there is no restriction on the dimensions of the manifolds. This is depicted in Proposition 5.4.9.

• Performance of CAPD algorithm. Finally, we show by simulations that the distance of iterates of the CAPD algorithm converges log-linearly. We also show the drawbacks of estimating the factor matrices all-at-once using the alternating projection method. Yet, the results evince that CAPD is suitable for solving the aforementioned tensor problem.

An introduction on alternating projections

The method of alternating projections consists of projecting iteratively a point onto manifolds, one after the other, in order to find a point in their intersection. To describe the method, let M 1 , M 2 , . . . , M n be n manifolds in some space E. The alternating projection is given as follows

x (k+1) = π Mn • • • π M 2 π M 1 (x (k) ), ∀k = 0, 1, . . . (5.1)
where π M i , 1 ≤ i ≤ n is the orthogonal projection operator onto M i , and k is the iteration.

This method was originally proposed by von Neumann in 1930's for two manifolds, and he also proved that the sequence of alternating projections of two linear manifolds converges to a point in the intersection [START_REF] Neumann | [END_REF].

Method 5.1 has also been much emphasized in the literature when the problem has two convex manifolds. We highlight, for instance, the articles [START_REF] Combettes | Method of successive projections for finding a common point of sets in metric spaces[END_REF], Higham 2002, Gubin 1967, Breg 1965] for some discussion and application of the alternating convex projection method.

The extension of alternating projections to nonconvex manifolds plays an important role in several applications. We cite, for instance, those in control design [START_REF] Grigoriadis | [END_REF]] and in image processing [Bauschke 2002]. For n = 2, some important theoretical results on convergence of the alternating nonconvex projection method are highlighted in the following.

In [Zangwill 1967], it is shown that if the sequence {x (k) } ∞ k=0 is bounded and the distance to M 1 ∩ M 2 is strictly decreasing, then there exists a convergent subsequence to a point in M 1 ∩ M 2 . This results is known as the Zangwill's global conergence theorem. However, it is a weak result since it does not provide any information about the convergence of the complete sequence.

The authors in [Lewis 2008] introduce the concept of transversal manifolds and show that a linear convergence of the alternating projection method can be ensured if the manifolds are transversal at some point in M 1 ∩ M 2 and the starting point x (0) is close enough to M 1 ∩ M 2 . However, this approach presents a restriction on the dimensions of the manifolds. In order to drop the limitation of the traversal concept, the authors in [Andersson 2013] proposed a new concept called non-tangencial intersection point, which generalizes that of transversality. Both the transversality and non-tangential concepts will be studied in the context of orthogonal tensor decomposition. Some results of these two last publications are presented here and contextualized to our specific problem.

Column-wise orthogonal tensor decomposition

A tensor is orthogonal if its rank-1 decomposable components are orthogonal to one another, i.e., X i , X j = 0, ∀i = j, for the rank-R tensor T = X 1 + X 2 + • • • + X R , where X i are rank-1 tensors. The notion of orthogonality between a couple of decomposable tensors depends on the level of coupling between their factors. To be clearer, consider two three-way decomposable tensors

X 1 = a 1 ⊗ b 1 ⊗ c 1 and X 2 = a 2 ⊗ b 2 ⊗ c 2 . The tensors X 1 and X 2 are completely orthogonal if a 1 ⊥ a 2 , b 1 ⊥ b 2 and c 1 ⊥ c 2 .
Actually, all factors do not need to be orthogonal to one another to ensure the orthogonality of the tensor. For instance, if only a 1 ⊥ a 2 , then X 1 , X 2 = 0. Indeed,

X 1 , X 2 = a 1 , a 2 b 1 , b 2 c 1 , c 2 .
Herein, for a rank-R tensor, we assume that the factor matrix A of a tensor T is semi-unitary1 , which implies that a i , a j = δ(ij), ensuring the orthogonality of T . We call column-wise orthogonal tensor decomposition the decomposition of a tensor whose at least one of the factors matrices is semi-unitary. Notice that if all factor matrices are semi-unitary then the tensor is of unit norm and completely orthogonal. For other notions of orthogonality, see the article [Kolda 2001].

The column-wise orthogonal tensor decomposition plays an important role in several applications, such as blind source separation [De Lathauwer 2007], array processing [Miron 2008],and wireless systems [START_REF] Sørensen | [END_REF]]. Some important results are drawn in the literature. In [Wang 2015], the authors show that the best orthogonal low rank approximation of a tensor, with one semi-unitary factor matrix, always exists. Additionally, for almost all tensors, they also show that the ALS algorithm converges globally. In [START_REF] Sørensen | [END_REF]], the authors discuss the uniqueness and low-rank approximation properties of a CP decomposition with a semi-unitary matrix factor. In the following, we apply the results of alternating projection methods to solve this kind of decomposition.

Combined deflation and alternating projection algorithm

Let A ∈ K I×R , B ∈ K J×R and C ∈ K K×R be the factor matrices of a tensor T with entries in some field K. We assume that rank(T ) = R is known, and that one of the factor matrices is column-wise orthogonal, say A. This implies R ≤ I. Actually, without loss of generality, A can be viewed as a semi-unitary matrix due to scalar elementary indeterminacies (unitary if R = I).

Let T = A (C ⊙ B) T be the mode-1 unfolding of tensor T . Since A is semi-unitary, A H A = I. Therefore, after elementary operations

T H A = C * ⊙ B * (5.2)
For a specific column of A, namely a r , we derive from equation (5.2) that

T H a r = c * r ⊠ b * r , (5.3)
where b r and c r are the r-th column of the matrices B and C, respectively.

We propose a combined alternating projection and deflation algorithm, called CAPD, that recovers all factor matrices of a given orthogonal tensor T with a semi-unitary factor. The idea is to solve first the equation 5.3 with an alternating projection algorithm in order to estimate the factors a r , b r and c r , and to perform the deflation of the rank-1 component a r ⊗ b r ⊗ c r from T , obtaining thus a rank-(R -1) tensor Ta r ⊗ b r ⊗ c r . The process is repeated2 R -2 times until all rank-1 components are estimated.

Let T [r] be the mode-1 unfolding of a rank-(Rr + 1) tensor in K I×J×K , with 1 ≤ r ≤ R. We define the following linear subspace (5.4) and the manifold (5.5) which are the column space of T H [r] T H+ [r] and the subset of unit Kronecker vectors in K JK , respectively. By applying the Unvec operator on x, N can be actually viewed as the set of J × K unit-norm rank-1 matrices. In the proposed CAPD algorithm, the factors a r , b r and c r are computed by using the alternating projection onto these two manifolds. We set initially T [1] = T . CAPD is described in Algorithm ( 7

M r = {x ∈ K JK : T H [r] T H+ [r] -I x = 0},
N = {x ∈ K JK : x = c ⊠ b, x = 1, ∀b ∈ K J \ {0} and ∀c ∈ K K \ {0} },
(0) r ⊠ b (0) r ); k ← 0; repeat k ← k + 1; x (k) r ← π Mr π N (x (k-1) r
); until some stopping criteria are satisfied;

a r ← T H+ [r] x (k)
r ; b r ← conjugate of the left singular vector of Unvec(x

(k) r ); c r ← right singular vector of Unvec(x (k) r ); A(:, r) ← a r ; B(:, r) ← b r ; C(:, r) ← c r ; T [r+1] ← T [r] -a r ( c r ⊠ b r ) T ; end T [R]
is the mode-1 unfolding of the last rank-1 component of T .

Algorithm 7: CAPD algorithm

The algorithm works as follows. For every r ∈ {1, 2, . . . , R -1}, the vectors b

T H a r = c * r ⊠ b * r =⇒ R i=1 α i T H a i = c * r ⊠ b * r =⇒ c r ⊠ b r = R i=1 α * i (c i ⊠ b i ).
Now, let α j = 0 for some j ∈ {1, 2, . . . , R}, and write

c j ⊠ b j = 1 α * j    c r ⊠ b r - R i=1 i =j α * i (c i ⊠ b i )    .
Yet, a decomposition of T along the mode-1 is given by

T = a 1 (c 1 ⊠ b 1 ) T + • • • + a R (c R ⊠ b R ) T .
Thus, plugging the expression of c j ⊠ b j into that of T and reorganizing the factors we obtain

T = a 1 - α * 1 α * j a j (c 1 ⊠ b 1 ) T + • • • + a j-1 - α * j-1 α * j a j (c j-1 ⊠ b j-1 ) T + 1 α * j a j ( c r ⊠ b r ) T + + a j+1 - α * j+1 α * j a j (c j+1 ⊠ b j+1 ) T + • • • + a R - α * R α * j a j (c R ⊠ b R ) T .
It is important to mention that c r ⊠ b r is not in the column space of the vectors {c 1 ⊠ b 1 , . . . ,

c j-1 ⊠ b j-1 , c j+1 ⊠ b j+1 , . . . , c R ⊠ b R }, otherwise rank(T ) < R.
Since the first factor matrix is column-wise orthogonal and the decomposition is essentially unique, the vectors

{a 1 - α * 1 α * j a j , . . . , a j-1 - α * j-1 α * j a j , 1 α * j a j , a j+1 - α * j+1 α * j a j , . . . , a R - α * R α * j a j }
are orthogonal to one another. In particular, for i ∈ {1, 2, . . . , R} -{j}

a i - α * i α * j a j H 1 α * j a j = 0 =⇒ α i = 0,
which implies that a r = α j a j and c r ⊠ b r = α * j (c j ⊠ b j ). As c r ⊠ b r is a unit vector, it follows that

|α j | = 1/ c j ⊠ b j .
To eliminate scalar indeterminances, assume c j ⊠ b j = 1, which leads to |α j | = 1. Thus,

a r ⊗ b r ⊗ c r = |α j |a j ⊗ b j ⊗ c j = a j ⊗ b j ⊗ c j ,
and the proof is complete.

Corollary 5.4.4 x ⋆ r ∈ M r ∩ N =⇒ x ⋆ r is a column of C * ⊙ B * .
Proof: Up to scalar elementary indeterminances, it follows from Theorem 5.4.3 that

x ⋆ r = c * r ⊠ b * r = α j (c * j ⊠ b * j )
, for some j ∈ {1, 2, . . . , R}. Theorem 5.4.3 ensures that one rank-1 component of tensor T can be extracted if CAPD converges to a limit point x ⋆ r ∈ M r ∩ N. Actually, the convergence is not globally ensured since N is a non-convex manifold, which means that CAPD can get stuck if the starting point is badly chosen. Therefore, we focus our convergence study of CAPD method on starting points close to M r ∩ N, in which case we can draw some important results in the real field. Before, we introduce two basic definitions on manifolds.

Definition 1 [Lewis 2008] In some space E, let A and B be two C k -manifolds around a point x ∈ A ∩ B. The manifolds A and B are transverse at x if

T A (x) + T B (x) = E,
where T A (x) and T B (x) are the tangent spaces to A and B at x, respectively.

Definition 2 [Andersson 2013] In some space E, let A and B be two

C k -manifolds around a point x ∈ A ∩ B.
x is a non-tangential point if and only if

T A (x) ∩ T B (x) = T A∩B (x),
where T A (x), T B (x) and T A∩B (x) are the tangent spaces to A, B and A ∩ B at x, respectively According to [Lewis 2008], if two smooth manifolds are transverse, and the initialization is close enough to their intersection, then the alternating projection method converges linearly to a point at the intersection of the manifolds. This is stated in the following theorem Theorem 5.4.5 [Lewis 2008] In some space E, let A and B be two transverse manifolds around a point x ∈ A ∩ B. If the initial point x 0 ∈ E is close to x, then the method of alternating projections

x (k+1) = π A π B (x (k) ), (k = 0, 1, 2, . . .)
is well-defined, and the distance d A∩B (x k ) from the iterate x k to the intersection A ∩ B decreases Q-linearly3 to zero.

In our case, M r are linear subspaces and N is locally smooth in the Euclidean space (see Example 2 in [Lewis 2008]). Thus, in particular conditions, the manifolds M r and N are transverse for every r ∈ {1, 2, . . . , R -1}. This is shown in the following.

Lemma 5.4.6 Let R = {x ∈ R JK : x = c ⊠ b, ∀b ∈ R J \ {0} and ∀c ∈ R K \ {0} }. Then ∀x ∈ R it follows that T R (x) = T N (x/ x ).
Indeed, to ensure the convergence of CAPD algorithm when r = R -1, we must have R -1 ≤ R -(J -1)(K -1) =⇒ (J -1)(K -1) ≤ 1, which is true only when J ≤ 2, K ≤ 2. In other cases, Proposition 5.4.7 can only ensure, at best, the convergence for a few values of r. For some scenarios, Table 5.1 shows how many rank-1 components can be ensured by CAPD using the transversality concept. Notice that both the dimension I (omitted in the table) and the rank R must be much larger than J and K in general, in order to estimate some rank-1 components with convergence guarantees. Even for rank-100 tensors with dimensions I × 20 × 20, with I ≥ 100 for ensuring the column-wise orthogonality of factor A, the convergence to any component cannot be ensured by our algorithm. Thus, the table evinces the limitation of Proposition 5.4.7.

(J, K) R # of rank-1 components (2, 2) 2 2 (3, 3) 2 0 (3, 3) 3
Instead of considering transversality between manifolds, we can evaluate the convergence of CAPD under the non-tangential concept presented in Definition 2. According to [Andersson 2013], non-tangential points at the intersection of two manifolds is less restrictive than the transversality concept. Actually, the latter requires that the individual tangent spaces of the manifolds in a space must generate the whole space, which is not necessary for non-tangential points. The authors show that the convergence to a point at the intersection of two smooth manifolds can be ensured if the initialization of the alternating projection algorithm is close to a non-tangential point in the intersection. This result is shown in theorem below.

Theorem 5.4.8 [Andersson 2013] Let A, B and A ∩ B be C 2 -manifolds, and let x ∈ A ∩ B be a non tangential intersection point of A and B. Given ε > 0 and cos(α(x)) < c < 1, where α(x) is the minimal angle between A and B at x, there exists a ρ > 0 such that for any x (0) ∈ B ρ (x), where B ρ (x) is the open ball of radius ρ around x, the sequence of alternating projections

x (k+1) = π A π B (x (k) ), (k = 0, 1, 2, . . .) (i) converges to a point x ⋆ ∈ A ∩ B, (ii) x ⋆ -π A∩B (x (0) ) ≤ ε x (0) -π A∩B (x (0) ) , (iii) x ⋆ -x (k) ≤ cte •c k x (0) -π A∩B (x (0) ) .
We are mainly interested in the result (iii). Indeed, x (k) → x ⋆ as k → ∞. In order to apply the result of [Andersson 2013] to our case, we must first prove under some conditions that there exists at least a non-tangential point x r ∈ M r ∩ N, for every r ∈ {1, 2, . . . , R -1}. This is shown in the following proposition.

Proposition 5.4.9 Let T ∈ R I×J×K be a rank-R tensor with factor A column-wise orthogonal. If T has an essentially unique decomposition then any point x r ∈ M r ∩ N, r ∈ {1, 2, . . . , R -1}, is non-tangential.

Proof: Since the decomposition is essentially unique, the columns of B and C are unique up to scalar factors and not collinear, otherwise the Kruskal's rank of some factor matrix would equal to 1, implying non-uniqueness. Moreover, from Collorary 5.4.4, all points in M r ∩ N are columns of C ⊙ B. Thus, it turns out that ∃ρ > 0 :

B ρ (x r ) ∩ M r ∩ N = {βx r }, ∀β ∈ R such that |β| < ρ.
Since the straight along x r is a linear space, we can conclude that T Mr∩N (x r ) = {γx r }, ∀γ ∈ R.

On the other hand, the uniqueness of x r ∈ M r ∩ N within B ρ (x r ) ensures that T N (x r ) = {γx r }, ∀γ ∈ R, which leads to T N (x r )∩T Mr (x r ) = {γx r }, ∀γ ∈ R, since T Mr (x r ) = M r . Hence, T Mr (x r ) ∩ T N (x r ) = T Mr∩N (x r ), at any point x r ∈ M r ∩ N, for every r ∈ {1, 2, . . . , R -1}.

Proposition 5.4.9 brings up a strong result that ensures the convergence of CAPD algorithm. The only imposed constraint is the uniqueness of the decomposition, which is desired in practical applications.

Results on the orthogonal CP decomposition

In this section, we evaluate the performance of CAPD algorithm by numerical experiments. First, we show how the distance between iterates of the alternating projection part of our algorithm decreases as the number of iterations increases. Second, we will see that the convergence of a random tensor with a column-wise orthogonal factor matrix is log-linear. Finally, we present the advantage of CAPD with respect to the alternating projection method that computes the factor matrices all-at-once.

Iterates of CAPD algorithm for real tensors

We consider a sample of 1000 real tensors with the factor matrices generated as follows

• The entries of B and C are distributed according to a uniform measure in [-1, 1];

• For A, we first generate a random matrix with entries obeying a uniform distributed in [-1, 1],

and then set A as the left singular matrix of that random matrix, ensuring the orthogonality of the factor.

Let ∆ k = x (k+1) r -x (k) r
be the norm of the difference between two successive iterates in CAPD algorithm for some 1 ≤ r ≤ R -1 at iteration k ≥ 0 . Since we have 1000 tensors, we replace E{∆ k } by the average value at iteration k. We assume the tensors have rank 4, so that we evaluate the algorithm by estimating 3 rank-1 components (the last component is obtained directly by performing the deflation, as shown in Algorithm 7) for four scenarios.

Figure 5.1 shows for the scenarios where K = 3 that the average distance between iterates decreases as the rank of the tensor is reduced by deflation. For the scenario 20 × 20 × 20, the same phenomenon is noted for 10 -2 < E{∆ k } < 10 -6 , which covers most of the range of the other scenarios. We also note from the last scenario that for tensors with all dimensions much larger than the rank, the CAPD algorithm converges in a few iterations. Indeed, the iterates converge to zero4 in 60, 55, and 40 iterations approximately,for r equal to 1, 2 and 3, respectively. We choose the scenario 5×4×3. In Figure 5.2, we see that the distance of iterates of CPAD algorithm converges log-linearly to zero. Although we do not show other examples here, this behavior can be noted for most of the tensors of our sample. Moreover, the log-linear convergence also depends on the starting point x , for some k, does not decrease when the rank decreases. Indeed, the convergence is slower when the rank of the tensor is only 2 (r = 3).

Performance for complex tensors

Herein, we evaluate the performance of CAPD in terms of computational time and the percentage of successful decompositions denoted by %p. We assume that CAPD delivers a correct decomposition for a tensor T if T -T ≤ 10 -6 , where T is the estimated tensor. We generate a sample of 500 complex tensors. The real and imaginary parts are uniformly distributed in [-1, 1], and the columnwise factor A is the left singular matrix of a random complex matrix with uniform entries generated as before. We simulate the same four scenarios as in the previous subsection. Notice that the average time of simulations for 20 × 20 × 20 tensors corroborates with the convergence performance of iterates depicted in Figure 5.1, where we have shown that CAPD algorithm converges very fast when the dimensions are much larger than the rank. Here, it results in a small computational time compared to the other scenarios.

Simultaneous estimation of factors

As mentioned in Remark 5, the estimation of the factor matrices all-at-once might be performed by using an alternating projection algorithm based on equation (5.2). Instead of initializing the vector

c (0) r ⊠ b (0)
r , we could initialize directly a Khatri-Rao matrix C (0) ⊙ B (0) , and thus we would have the following alternating projection problem

X (k+1) = π M π N R (X (k) ), (k = 0, 1, 2, . . .), where M = X ∈ K JK×R : T H T H+ -I X = 0 , and N R = {X ∈ K JK×R : X = C ⊙ B, X = 1, ∀B ∈ K J×R \ {0} and ∀C ∈ K K×R \ {0} }.
If the algorithm converged to a point in M ∩ N R so that the right factors B and C could be properly estimated, then the factor A would be estimated by solving the matrix linear equation (5.2), whose the least square solution is A = T H+ (C ⊙ B).

This approach presents two crucial drawbacks:

• The complexity of this all-at-once alternating projection algorithm is larger than that of the CAPD algorithm. Indeed, the employment of an alternating projection algorithm to solve equation (5.2) is equivalent to employ R parallel alternating projection algorithms to solve R equations (5.3). Notice, however, that CAPD solves only R -1 times the equation ( 5.

3) .

• There is no guarantee that the estimated columns of the factors A, B, and C be the same of A, B and C. As the initialization

C (0) ⊙ B (0) is actually R rank-1 initializations c (0) r ⊠ b (0)
r (of R parallel problems), some of those initializations can converge to the same point, which means that the estimated matrices might have collinear columns. Indeed, this happens more often that one can imagine.

In order to show the limitation of computing the factors all-at-once using the alternating projection method, we count how many column factors are properly estimated for the matrices A, B and C. We evaluate the following real and complex scenarios for a sample of 100 tensors: 15×10×5, 20×10×5, 20 × 20 × 5 and 20 × 20 × 20 tensors. For each tensor, we perform a random initialization.

Figure 5.3 shows that the estimation of all factors all-at-once using the alternating projection method happens in a very few cases. In fact, the best scenario is the one with 20 × 20 × 5 real tensors, in which only for 7 tensors it was possible to extract completely the factor matrices. In all scenarios, the algorithm recovers 3 or 4 column factors of A, B and C in general. To increase the chances of estimating the right columns factors, we could test more initializations. However we would run again the algorithm, increasing though the number of operations to be computed. Therefore, CAPD is more efficient than the all-at-once alternating projection computation.

Chapter Summary and Directions

This chapter finishes with the following conclusions and further directions: • CAPD algorithm. The alternating projection method combined with the deflation of rank-1 decomposable tensors arose as an alternative method to compute the exact CP decomposition of an orthogonal tensor with a semi-unitary factor matrix.

• Extraction of factors with CAPD algorithm. We proved that CAPD algorithm extracts one column of each factor matrix when the iterates of our proposed algorithm converges to a point in the intersection of two manifolds.

• Transversality of points in CAPD algorithm. Under some conditions on the factor matrices B and C, we showed that the convergence of the alternating projection part of the CAPD algorithm to a point in M r ∩ N is linear when the initialization of the algorithm is close to M r ∩ N and the transversality at points in the intersection of the manifolds M r and N is satisfied. In that case, we could evaluate how many columns factors can be extracted for different sizes and rank of tensors.

• Non-tangential concept in CAPD algorithm. In order to drop the restriction on the dimensions of tangent spaces of the manifolds under the transversality point of view, we used the concept of non-tangential points. In this case, we proved that CAPD algorithm converges to the exact decomposition when the initialization is close enough to M r ∩ N, for every 1 ≤ r ≤ R -1.

• Performance of CAPD algorithm. From numerical simulations, we saw that the distance between iterates in the alternating projection part of CAPD decreases log-linearly to zero. We also showed that the computation of factors using CAPD performs better than the computation of factors all-at-once using an other alternating projection method applied directly to the whole factor matrices. Indeed, it was shown that the complexity of CAPD is smaller and the all-at-once method does not ensure the exact estimation of the factor matrices.

• Low approximation with CAPD. For future work, we intend to evaluate CAPD algorithm under noisy scenarios. We already note that the deflation procedure compromises the estimation of the factor matrices since one of the estimated rank-1 factors is not an exact component of the noiseless tensor. We do not know exactly in which scenarios we can take advantage of CPDA when compared with other standard algorithms, such as ALS and those described in [START_REF] Sørensen | [END_REF]].

Conclusions Main conclusions

We summarize the main conclusions of this thesis.

Chapter 2

• The computation of the best rank-1 approximation using the moment approach due to Lasserre was compromised by the exponential growth of variables due to the problem relaxations, in spite of the method has provided a certificate of global optimality. This variable drawback limited simulations for real tensors of order 3 and dimension 2 × 2 × 2. A test with tensors with dimension 3×3×3 showed that the computational time to compute a best rank-1 approximation is approximately 30 minutes, which show the inefficency of the method.

• The proposed SeROAP algorithm is always at least as good as the ST-HOSVD algorithm (and consequently as THOSVD) for three-way tensors. The mathematical proof was confirmed by numerical experiments. Moreover, the computational complexity of SeROAP is smaller than that of THOSVD and it is of the same order of the complexity of the ST-HOSVD algorithm.

• The proposed CE is an iterative algorithm to compute rank-1 three-way tensor approximations. It presents small computational complexity than ALS algorithm when one dimension is sufficiently larger than the other two dimensions. We saw that the approximate error and the number of iterations needed to converge is smaller for CE in this scenario. We also presented some theoretical results ensuring that the objective function of the rank-1 approximation problem converges to a stationary value when CE is employed.

Chapter 3

• We proposed an iterative deflation algorithm to compute low rank tensor approximations whose mainly idea is to update rank-1 components using successive rank-1 approximations. This algorithm, called DCPD, is different to the method propose in [Cichocki 2009c] that updating the factors of the tensors using ALS. With our approach, we saw by simulations that residuals are quickly reduced, so that DCPD converges in a few iterations.

• We performed a study on the convergence of DCPD considering the best rank-1 approximation of the components within. Three main results were presented: (i) the monotonic decrease of a sequence of residuals; (ii) The recovering of an exact decomposition for DCPD when the residuals obeys some geometric conditions; (iii) a conjecture about the existence of continuous probability measures ensuring the convergence of DCPD with high probability.

expect to improve the computational time to compute the best rank-1 approximation taking into account their sparsity.

• SeROAP arose as a simple and efficient non-iterative algorithm to compute rank-1 approximations. In a future work, we intend to develop new strategies to tackle higher order tensors outperforming SeROAP. We also have some ideas to improve SeROAP, such as finding better approximations by looking for other rank-1 tensors from that one delivered by SeROAP.

• For the CE algorithm, we do not prove that it converges to a stationary point, but only that the objective of the best rank-1 approximation problem converges to some value. The proof is still a challenge.

Chapter 3

• The proof of Conjecture 3.2.8 is an open problem, even for some constrained cases (e.g. orthogonal tensors). We think that some tools in theory measure and algebraic geometry can gives us some directions to overcome this problem.

Chapter 4

• There are some unanswered questions about the reduction of polynomial quadratic systems into rank-1 approximation problems: (i) What are the conditions ensuring the equivalence between problems (P1) and (P2)? (ii) What is the set of tensors at which the equality S P 4 = S P 5 holds?

We intend to tackle these two question soon.

Chapter 5

• We do not evaluate the CAPD algorithm under noisy scenarios. Actually, the theoretical results only work for the exact CP decomposition. However, we intend to evaluate the impact of noise on the estimation of factor matrices. An idea for mitigating the errors due to noise is to approximate the unfolding matrix T [r] to a rank Rr + 1 matrix by computing the SVD. We do not know yet whether this will be able to improve the estimations.

input : T ∈ K I 1 ×I 2 ו••×I N : input data, R: rank parameter. output: A (n) ∈ K In×R , for n = 1, . . . , N : factor matrices Initialize A (1) , . . . , A (N ) ; a

(n) r ← a (n) r / a (n) r , ∀r, n = 1, 2, . . . , N -1; E = T -R r=1 a (1) r ⊗ a (2) r ⊗ • • • ⊗ a (N ) r ; repeat for r = 1 to R do T [r] ← E + a (1) r ⊗ a (2) r ⊗ • • • ⊗ a (N ) r ; for n = 1 to N do a (n) r ← T (n) [r] a (N ) * r ⊗ a (N -1) * r ⊗ a (n+1) * r ⊗ • • • ⊗ a (n-1) * r ⊗ • • • ⊗ a (2) * r ⊗ a (1) * r ; if n = N then a (n) r ← a (n) r / a (n) r ; end E ← T [r] -a (1) r ⊗ a (2) r ⊗ • • • ⊗ a (N )
r ; end end until some stopping criterion is satisfied;

Algorithm 9: HALS algorithm

A.2 Moment relaxation: a brief overview

Let the following optimization problem

p min = inf x∈K p(x), (A.1) where K = {x ∈ R n : g 1 (x) ≥ 0, g 2 (x) ≥ 0, . . . , g m (x) ≥ 0}, with p(x), g i (x) ∈ R[x], for i ∈ {1, 2, . . . , m}.
According to Lasserre's approach [Lasserre 2001], Problem (A.1) is equivalent to

p min = inf µ K p(x)µ(dx), (A.2)
where µ is a probability measure on R n supported by the set K. Note that p(x) = α p α x α , α ∈ N n , so that it turns out that p(x)µ(dx) = p T y, where p is the vector of coefficients p α abiding by some monomial order indexed by α, and y the sequence of moments of µ such that y α = x α µ(dx). For more details about monomial orders refer to [START_REF] Cox | [END_REF]], and on theory of moments see [Berg 1987[START_REF] Curto | [END_REF][START_REF] Putinar | [END_REF]. To compute V C (ker(M s (y))), we can use the eigenvalue method in the multivariate case applied to the ideal constructed from a flat extension of M s (y). More details is presented in [Henrion 2005].

Problem (A.2) can be formulated as

The software GloptiPoly c [START_REF] Henrion | [END_REF][START_REF] Henrion | [END_REF] implements the moment problems (A.5) with the lexicographic basis, tests optimality and extracts global optimizers.

Introduction

L'approximation tensorielle de rang faible joue ces dernières années un rôle important dans plusieurs applications, telles que la séparation aveugle de source, les télécommunications, le traitement d'antennes, les neurosciences, la chimiométrie, et la fouille de données. La décomposition tensorielle Polyadique Canonique (CP) est très attractive comparativement à des outils matriciels classiques, notamment pour l'identification de systèmes. En fait, l'approche tensorielle présente des propriétés intéressantes d'unicité sous des faibles contraintes, tandis qu'avec la SVD (un exemple d'outil matriciel largement utilisé) la contrainte d'orthogonalité doit être imposée pour assurer l'unicité de la décomposition. L'objectif de la décomposition CP est de trouver toutes les composantes λ r et a (r) i , 1 ≤ i ≤ N, 1 ≤ r ≤ R, à partir du tenseur T en utilisant un algorithme comme la méthode des moindres carrées (ALS), par exemple.Cependant, dans la pratique, les tenseurs sont corrompus par du bruit, si bien que nous devons résoudre un problème d'approximation de rang faible, une fois que le rang du tenseur devient inconnue. Ainsi, le problème d'approximation peut être posé de la façon suivante:

arg min λr, a (r) i =1 T - R r=1 λ r (a (r) 1 ⊗ • • • ⊗ a (r) N ) . (B.2)
Cette thèse présente différents objectifs et est divisée en quatre autres chapitres hors cette introduction. Dans le chapitre suivant, notre objective est calculer l'approximation tensorielle de rang-1 en utilisant des méthodes itératifs et non-itératifs. Dans le troisième chapitre, nous utilisons des algorithmes non-itératifs de rang-1 dans le contexte de la déflation itérative, avec le but de résoudre le problème d'approximation (B.2). Le quatrième chapitre établit un lien entre la meilleure approximation de rang-1 et les systèmes quadratiques multivariés. Finalement, nous traitons le cas particulier de la décomposition CP exacte pour des tenseurs avec l'une des matrices facteurs semi-unitaire.

Chapitre: Approximation de Rang-1

Nous étudions des alternatives pour résoudre le problème (B.2) lorsque R = 1. En fait, ce problème est NP-difficile, de sorte que de nouveaux algorithmes pour le traiter ont besoin d'être proposés.

Ainsi, nous présentons les principales contributions de ce chapitre:

• Adaptation du problème de meilleure approximation de rang-1 pour assurer l'estimation de la solution global en utilisant la méthode de Lasserre [Lasserre 2001]. Dans le problème

p min = min λ,a i ∈K I i T -λ • a 1 ⊗ • • • ⊗ a N 2 s.t. a i = 1. , (B.3)
nous établissons une valeur maximale pour λ afin d'assurer l'existence du minimum global des problèmes hiérarchiques de Lasserre (voir l'appendice A.2 du manuscrit en anglais).

• Algorithme non-itératif pour calculer une approximation de rang-1. Nous proposons un algorithme qui s'appelle SeROAP, dont l'idée est basée sur des successives SVDs et projections. Comparé à la méthode de HOSVD tronquée (T-HOSVD) [De Lathauwer 2000], l'algorithme SeROAP présente une complexité numérique plus petite. Par rapport à l'algorithme ST-HOSVD [START_REF] Vannieuwenhoven | [END_REF]], SeROAP est un algorithme compétitif en termes de complexité numérique pour des tenseurs de grandes dimensions et petit ordre.

• Erreur d'approximation de SeROAP. Nous démontrons algébriquement que l'erreur d'approximation obtenue avec SeROAP est au moins aussi bon que celle calculée avec l'algorithme ST-HOSVD (et par conséquent que HOSVD tronquée), pour des tenseurs d'ordre 3.

• Limitation de SeROAP. Les résultats numériques ont montré que l'algorithme SeROAP présente une performance faible de l'erreur d'approximation pour des tenseurs d'ordre supérieur à quatre, par rapport aux autres méthodes non-itératifs de rang-1 déjà mentionnées. En fait, même pour l'ordre 4, la performance de ST-HOSVD sous ce critère surpasse celle de notre algorithme. En revanche, SeROAP est statistiquement mieux que la HOSVD tronquée avec une probabilité élevée pour cette même ordre.

• Algorithme d'approximation de rang-1 itératif. 1 , Cela nous permet de mettre à jour la deuxième composante appelée X

(1) 2 . Nous continuons avec la même procédure jusqu'à la composante R, de manière à obtenir un nouveau résidu E (1) tel que E (1) ≤ E (0) . La déflation itérative continue jusqu'à la convergence (ou un autre critère d'arrêt), représentée dans la figure par l'itération k.

Les contributions de ce chapitre sont les suivantes:

• Déflation itérative + algorithme itératif d'approximation de rang-1. Avec le but de réduire le plus rapidement possible la norme des résidus E (k) , on introduit les algorithmes SeROAP, ST-HOSVD et T-HOSVD pour mettre à jour les composantes X (j) r , ce qui entraîne plus d'intelligence dans la procédure de mise à jour des composantes en comparaison à la méthode de déflation itérative HALS [Cichocki 2009c].

• Étude théorique sur la convergence de l'algorithme DCPD. Nous montrons que la norme des résidus E (k) , k ≥ 1, est une séquence décroissante. Nous montrons aussi que l'algorithme DCPD estime la décomposition CP exacte d'un tenseur donné lorsque les résidus n'appartient pas à un cône de volume arbitrairement petit. Dans une seconde étape, sous les contraintes de probabilité absolument continue de la distribution des tenseurs et de la continuité d'une certaine fonction des résidus, on montre que DCPD + meilleure approximation de rang-1 peut réduire le résidu initial d'un facteur inférieur à une constante connue (qui dépend de l'angle du cône) après L itérations avec une probabilité élevée. Finalement, on conjecture (Conjecture 3.2.8) l'existence des mesures de probabilités qui assurent la convergence de DCPD + meilleure approximation de rang-1 vers la décomposition exacte avec une haute probabilité.

• Performance de DCPD avec les algorithmes d'approximation de rang-1 non-itératifs. Notre conjecture est renforcée par des résultats numériques. On génère des tenseurs dont ces composantes suivent une distribution uniforme. Nous montrons également la réduction rapide du résidu en utilisant DCPD + SeROAP par rapport à la méthode de déflation itérative HALS. Nous présentons aussi un ensemble de résultats numériques pour des scénarios bruités et non-bruités, avec la variation du rang et de l'ordre tensoriels.

Chapitre: Systèmes Quadratiques Multivariés et la Meilleure Approximation de Rang-1

Dans ce chapitre, l'objectif est de montrer que le système quadratique P1:

       x T A 1 x + b T 1 x + c 1 = 0 . . . x T A m x + b T m
x + c m = 0, où x ∈ R n , A j ∈ R n×n sont des matrices symétriques, b j ∈ R n et c j ∈ R, 1 ≤ j ≤ m, peut être réduit à un problème de meilleure approximation de rang-1 d'un tenseur d'ordre 3. La généralisation pour la cas complexe est faite à partir du système P2 ci-après

• Extraction des bonnes composantes de rang-1. Dans le cas de convergence de l'étape de projection alternée de l'algorithme CAPD, nous montrons qu'une correcte composante de rang-1 du tenseur est toujours obtenue.

• Convergence vs. transversalité. Sous des certaines conditions, nous montrons que s'il existe des points transverses à l'intersection des variétés qui définissent la décomposition CP orthogonal de facteur semi-unitaire. Ainsi, pour une initialisation de l'algorithme CAPD proche des variétés impliquées, la transversalité assure une convergence linéaire [Lewis 2008] vers un point à l'intersection des variétés. Cependant, le concept de transversalité est limité dans notre cas, une fois que l'estimation correcte de seulement quelques composantes de rang-1 du tenseur peut être assurée.

• Convergence vs. non-tangentialité. Nous démontrons un résultat plus fort basé sur le concept de la non-tangentialité entre variétés. Pour l'algorithme CAPD, nous prouvons que l'estimation de toutes les composantes de rang-1 peuvent être assurée lorsque des points à l'intersection des variétés impliquées sont non-tangentielles. Effectivement, le concept de tranversalité est plus restreint que celui de la non-tangetialité [Andersson 2013].

• Performance de l'algorithme CAPD. Nos présentons des résultats numériques montrant que la distance entre les estimations à chaque itération de l'algorithme CAPD converge log-linéairement. Nous montrons aussi la problématique existante lors de l'estimation de toutes les matrices facteurs à la fois, avec la méthode de projection alternée. En fait, sans la déflation de tenseurs de rang-1, l'estimation de colonnes proportionnelles dans le matrices facteurs est possible, de façons que n'existe aucune assurance de trouver les facteurs correctement. Les résultats numériques ont montré que cela est un problème courant.

Perspectives

Nous présentons maintenant les principales perspectives de mon travail de thèse. Les détails des activités envisagées sont décrits par chapitre.

Chapitre: Approximation de Rang-1

• Développer d'autres méthodes dans la Géométrie Algébrique pour résoudre le problème de meilleure approximation de rang-1. Malgré les techniques décrites en [Bucero 2014, Nie 2014] présentent une meilleure performance que la méthode de Lasserre en termes de complexité numérique, elles subissent à une augmentation exponentielle du nombres de variables lors de la relaxation des problèmes d'optimisation. Ainsi, une fois que le nombre de monômes de l'objectif du problème (B.3) est parcimonieux, nous espérons réduire la complexité en tenant compte de ce fait.

• Nous envisageons soit une amélioration de SeROAP, soit la proposition d'un nouvel algorithme non-itératif de rang-1 pour traiter des tenseurs d'ordre supérieur, afin de réduire l'erreur d'approximation. Dans le cas de SeROAP, une idée d'amélioration c'est de choisir l'ordre de dépliage par rapport le rang des versions dépliées du tenseur. Plus précisément, on choisit l'ordre de dépliage dont les matrices aient le rang le plus petit. Ainsi, on garde plus d'information du tenseur déplié.

• Pour algorithme CE, nous n'avons pas encore prouvé sa convergence vers un point stationnaire, mais seulement que l'objectif du problème de meilleure approximation de rang-1 converge vers une valeur limite. La démonstration de convergence vers la solution optimale est encore un défi.

Chapitre: Déflation Itérative

• La preuve de la Conjecture 3.2.8 est un problème ouvert, même pour des cas particulières (e.g. tenseurs orthogonaux). Nous envisageons d'utiliser les outils mathématiques de la théorie de la mesure et de la géométrie algébrique pour traiter ce problème.

Chapitre: Systèmes Quadratiques Multivariés et la Meilleure Approximation de Rang-1

• Il y a des questions autour de la réduction des systèmes quadratiques à des problèmes d'approximation de rang-1 pour lesquelles nous n'avons pas encore de réponses. Par exemple: quelles sont les conditions assurant l'équivalence entre les problèmes P1 et la meilleure approximation de rang-1? Nous envisageons de traiter cette question bientôt.

Chapitre: Projection Alternée Appliquée aux Tenseurs Orthogonaux

• Nous n'avons pas encore évalué la performance de l'algorithme CAPD sous des scénarios bruités. Les résultats théoriques marchent seulement pour le cas de la décomposition CP exacte. Cependant, nous avons l'intention d'évaluer l'impact du bruit sur l'estimation des matrices facteurs. Une idée pour mitiger des erreurs d'estimation est d'approximer les formes dépliées des tenseurs à des matrices de rang désiré, avant la déflation. D'ailleurs, avec cette approche, nous ne savons pas si l'estimation des facteurs sera satisfaisante comparée à d'autres méthodes telles que celles décrites en [START_REF] Sørensen | [END_REF]].

Abstract: Low rank tensor decomposition has been playing for the last years an important role in many applications such as blind source separation, telecommunications, sensor array processing, neuroscience, chemometrics, and data mining. The Canonical Polyadic tensor decomposition is very attractive when compared to standard matrix-based tools, manly on system identification. In this thesis, we propose: (i) several algorithms to compute specific low rank-approximations: finite/iterative rank-1 approximations, iterative deflation approximations, and orthogonal tensor decompositions. (ii) A new strategy to solve multivariate quadratic systems, where this problem is reduced to a best rank-1 tensor approximation problem. (iii) Theoretical results to study and proof the performance or the convergence of some algorithms. All performances are supported by numerical experiments.

Keywords: CP decomposition, rank-1 approximation, deflation, quadratic system, orthogonal tensor.

Résumé : L'approximation tensorielle de rang faible joue ces dernières années un rôle important dans plusieurs applications, telles que la séparation aveugle de source, les télécommunications, le traitement d'antennes, les neurosciences, la chimiométrie, et l'exploration de données. La décomposition tensorielle Canonique Polyadique est très attractive comparativement à des outils matriciels classiques, notamment pour l'identification de systèmes. Dans cette thèse, nous proposons (i) plusieurs algorithmes pour calculer quelques approximations de rang faible spécifique: approximation de rang-1 itérative et en un nombre fini d'opérations, l'approximation par déflation itérative, et la décomposition tensorielle orthogonale; (ii) une nouvelle stratégie pour résoudre des systèmes quadratiques multivariés, où ce problème peut être réduit à la meilleure approximation de rang-1 d'un tenseur; (iii) des résultats théoriques pour étudier les performances ou prouver la convergence de quelques algorithmes. Toutes les performances sont illustrées par des simulations informatiques.

Mots clés : décomposition CP, approximation de rang-1, déflation, système quadratique, tenseur orthogonal.

Resumo : A aproximação tensorial de baixo posto desempenha nestes últimos anos um papel importante em várias aplicações, tais como separação cega de fontes, telecomunicações, processamento de antenas, neurociênca, quimiometria e exploração de dados. A decomposição tensorial canônica é bastante atrativa se comparada às técnicas matriciais clássicas, principalmente na identificação de sistemas. Nesta tese, propõe-se (i) vários algoritmos para calcular alguns tipos de aproximação de posto: aproximação de posto-1 iterativa e em um número finito de operações, a aproximação por deflação iterativa, e a decomposição tensorial ortogonal; (ii) uma nova estratégia para resolver sistemas quadráticos em várias variáveis, em que tal problema pode ser reduzido à melhor aproximação de posto-1 de um tensor; (iii) resultados teóricos visando estudar o desempenho ou demonstrar a convergência de alguns algoritmos. Todas os desempenhos são ilustrados através de simulações computacionais.

Palavras-chave : decomposição CP, aproximação de posto-1, deflação, sistema quadrático, tensor ortogonal.

  introduction to canonical polyadic tensor decomposition . . . . . . . . . . . . . 15 1.2 Why study the canonical polyadic tensor decomposition? . . . . . . . . . . . . . . . 16 1.3 Different objectives and the scope of the thesis . . . . . . . . . . . . . . . . . . . . . 19 1.4 Structure of the thesis and outline . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 1.5 Publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Rank-1 tensor approximation 21 2.1 Rank-1 tensor approximation: description and algorithms . . . . . . . . . . . . . . . 24 2.2 An algebraic geometric approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.3 Finite rank-1 approximation algorithms . . . . . . . . . . . . . . . . . . . . . . . . 28 2.4 Iterative rank-1 approximation methods . . . . . . . . . . . . . . . . . . . . . . . . 40 2.5 Chapter Summary and Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 Iterative deflation 47 3.1 Deflation on CP decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.2 Study on iterative deflation and best rank-1 tensor approximations . . . . . . . . . . 51 3.3 Estimation of the CP decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.4 Chapter Summary and Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 Multivariate quadratic systems and the best rank-1 tensor approximation 65 4.1 An introduction to multivariate quadratic systems . . . . . . . . . . . . . . . . . . . 68 4.2 From quadratic systems to best rank-1 approximations . . . . . . . . . . . . . . . . 69 4.3 Discussion on the real tensor approach . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.4 Extension to the complex field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.5 Conjugated partially symmetric rank-1 approximation . . . . . . . . . . . . . . . . . 77 4.6 Examples on real and complex fields . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.7 Performance evaluation of rank-1 tensor approach . . . . . . . . . . . . . . . . . . . 81 5 Nomenclature K the set of numbers in a generic field.

  spanthe space spanned by a set of vectors,matrices or tensors.card the cardinality of a set.ker the kernel of a matrix.trace the trace of a matrix.OLandau's notation for computational complexity. Only multiplications are counted. introduction to canonical polyadic tensor decomposition

Figure 1

 1 Figure 1.1: CP decomposition of a three-way tensor.

Figures 2 Figure 2 . 3 :

 223 Figures 2.1, 2.2 and 2.3 depict the better performance of SeROAP over THSOVD for p 1 , p 2 , and p 3 , respectively.

  Figure 2.4: SeROAP vs ST-HOSVDp 1 = [1 2 3].

Figure 2 . 7 :

 27 Figure 2.7: SeROAP vs ST-HOSVD and THOSVD -4-th and 5-th order tensors.

  eigenvalues whose eigenvectors are x t and y t , respectively.The eigenpair (λ (y) t+1 , y t+1 ) obtained by solving equation (2.15) with x = x t , is solution of the maximization problem λ

Corollary 3 .

 3 2.2 In DCPD algorithm, the inequality E[r, l] ≤ sin(γ[r, l]) E [r -1, l] holds for any 1 < r ≤ R.Proof: By replacing X , E and γ in Lemma 3.2.1 with X [r, l -1], E [r -1, l] and γ[r, l] respectively, the result follows directly. Corollary 3.2.3 For any l > 1 and c

Figure 3 .

 3 Figure 3.2 presents the performance of the algorithms in terms of the average of the normalized norm of residuals per iteration, namely E{ E[R, l] / T }, for different values of the signal-to-noise (SNR) ratio for 5 × 5 × 5-rank-3 tensors. Again, a sample of 500 complex tensors with real and imaginary parts of the entries uniformly distributed in [-1, 1] is considered. Both ALS and HALS algorithms are randomly initialized. The same stopping criteria of Section 3.3.2 is established here. Additive Gaussian noise is considered in our simulations.

Figure 3 .Figure 3 . 2 :

 332 Figure 3.3 shows that ALS performs better than the other algorithm in all scenarios. In the scenario with SNR of 10 dB, while approximately 60% of simulations finished up to 100 iterations for ALS,

Figure 3 . 3 :

 33 Figure 3.3: Percentage of tensors at which the simulations has finished up to l iterations.

Figure 3 .

 3 Figure 3.5 shows the average normalized residual µ under tensor order variation. We evaluate the performance of the algorithms under a noiseless scenario, and under scenarios with SNR = 20 dB and SNR = 0 dB. A sample of 500 complex tensors with dimensions I 1 = I 2 = • • • = I N = 3 and rank R = 3 is generated according to a uniform distribution as before.

Figure 3 . 4 :

 34 Figure 3.4: Mean of the normalized residual µ for different ranks and 5 × 5 × 5 complex tensors.

Figure 3 . 5 :

 35 Figure 3.5: Mean of the normalized residual for different orders and I 1 = I 2 = • • • I N = 3.

  . y = 1.

Lemma

  ). The upper indices denote the k-th iteration of the algorithm, and π Mr and π N are the projection operators onto M r and N, respectively. Notice that π Mr = T H [r] T H+ [r] . input : T ∈ K I×J×K : input data output: A, B and C: factor matrices T ← mode-1 unfolding of T ; T [1] ← T ; Define N; for r = 1 to R -1 do Define the linear subspace M r from T [r]

r

  r ∈ ker(T [r] ), otherwise x (0) r = 0 and the algorithm would not extract the r-th rank-1 component. The vector x (0) r is introduced into a repeat loop that consists of the alternating projection phase. If x (k) r converges to a point in M r ∩ N , then the estimated factor a r is computed by solving in the least square sense the equation (5.3) for a given x ∈ ker(T [r] ), which implies that a r = 0, b r = 0, and c r = 0. From Lemma 5.4.2, a r can be written as a linear combination of the columns of the factor A.That is, a r = R i=1 α i a i . Since x (k) r → x ⋆ r ∈ M r ∩ N,it follows that a r , b r and c r satisfy equation (5.3), so that x ⋆ r = c * r ⊠ b * r . Substituting a r in (5.3)

  Number of estimated components ensured by CAPD under the transversality concept.

  Figure 5.1: Convergence of the iterates of CAPD algorithm for rank-4 tensors.

  Figure 5.2: Convergence of the iterates of CAPD algorithm for a single rank-4 tensor.

Figure 5 . 3 :

 53 Figure 5.3: Tensors in which n column factors are correctly extracted with the all-at-once algorithm.

  Theorem A.2.2[Laurent 2009] Let t ≥ max(d p , d K ) and K p = {x ∈ K|p(x) = p min }. Let also y ∈ N n 2t be a global minimum of the problem (A.5). Assume card(K p ) < ∞ and ∃s s.t. max(d p , d K ) ≤ s ≤ t and rank(M s (y)) = rank(M s-d K (y)). Then, p mom t = p min and V C (ker(M s (y))) ⊆ K p .

  Soit T ∈ K I 1 ×I 2 ו••×I N un tenseur d'ordre N avec ses composantes définies dans un corps K. La décomposition CP est représentée par où R est le nombre minimal de composantes de rang 1, aussi connue comme le rang du tenseur[Hitchcock 1927]. La Figure B.1 décrit la décomposition CP pour des tenseurs d'ordre 3.

Figure

  Figure B.1: Décomposition CP d'un tenseur d'ordre 3.

Table 2 .

 2 3: Performance of CE method compared to ALS initialized with SeROAP.

	3 × 4 × 5 3 × 4 × 20 3 × 20 × 20 20 × 20 × 20	91.00 % 87.00 % 86.67 % 62.67 %	27 44 71 102	13 14 27 52
	Table 2.2: Performance of CE method compared to ALS with a random initialization.
	Scenario 3 × 4 × 5 3 × 4 × 20 3 × 20 × 20 20 × 20 × 20	% ∆φ [n] ≥ 0 mean iteration (ALS) mean iteration (CE) 97.33 % 23 10 96.00 % 36 11 91.33 % 57 19 78.00 % 86 42

Table 3 .

 3 1: Number of operations (multiplications) per iteration of tensor algorithms.

	Scenario	ALS HALS DCPD-	DCPD-	DCPD-
			SeROAP STHOSVD THOSVD
	81 375 1029 1715 135 625 2187 3645 243 486 729 1701 3 × 3 × 3 × 3 × 3 × 3 2187 5832 3 × 3 × 3 5 × 5 × 5 7 × 7 × 7 9 × 9 × 9 3 × 3 × 3 × 3 3 × 3 × 3 × 3 × 3	360 1500 3920 8100 1170 3600 10890	342 1450 3822 7938 1098 3366 10170	702 3250 8918 18954 2754 10206 36450

Table 3 .

 3

2: Complexity per iteration /R of tensor algorithms.

Table 3 .

 3 3: Percentage of tensors in which the exact CP decomposition is succeeded.

	Algorithm ALS	3 × 3 × 3 5 × 5 × 5 7 × 7 × 7 9 × 9 × 9 77.8 % 99.6 % 99.0 % 99.6 %
	HALS	48.8 %	92.0 %	98.4 %	96.6 %
	DCPD-SeROAP	63.8 %	97.8 %	99.4 %	99.6 %
	DCPD-ST-HOSVD	67.0 %	98.0 %	99.6 %	99.8 %
	DCPD-THOSVD	70.0 %	99.0 %	99.6 %	99.8 %

  table below.

		complexity	percentage of	convergence residual residual residual
	Algorithm	per	successful	versus	versus	versus	versus
		iteration	decompositions	iteration	iteration	rank	order
	ALS	1	1	1	3	1	1
	HALS	2	5	5	5	5	3
	DCPD-SeROAP	4	4	2	1	2	5
	DCPD-STHOSVD	3	3	3	2	3	2
	DCPD-THOSVD	5	2	4	4	4	4

Table 4 .

 4 5: Performance of ALS algorithm to extract one solution of generic quadratic systems.

	Scenario (n, m) Initialization Time Time for constructing T ′ Time ALS Total Time p (%)
	(2, 3)	1.2707	8.5767e-04	0.0852	1.3567	100
	(4, 4)	1.5190	0.0014	2.0137	3.5340	100
	(5, 4)	1.5288	0.0024	0.4613	1.9925	100
	(6, 6)	1.5800	0.0035	3.5453	5.1288	100
	(7, 9)	1.6444	0.0051	0.9939	2.6433	98.6
	(10, 8)	2.7514	0.0133	2.0391	4.8037	100
	(10, 10)	2.7705	0.0133	8.8520	11.6358	100
	(20, 20)	11.8718	0.1993	51.4871	63.5582	100

Table 4 .

 4 002. 6: Performance of CPS algorithm to extract one solution of generic quadratic systems.

	Scenario (n, m) Initialization Time Time for constructing T ′ Time CPS Total Time p (%)
	(2, 3)	7.2546	9.9674e-04	1.0424	8.2980	96
	(4, 4)	8.1409	0.0020	30.9283	39.0712	70.2
	(5, 4)	8.7274	0.0030	9.0716	17.8020	97.2
	(10, 8)	27.2700	0.0217	36.4495	63.7413	79

Table below

 below

			summarizes
	the results.
	Scenario	Average time (seconds) %p
	5 × 4 × 3 20 × 4 × 3 20 × 20 × 3 20 × 20 × 20	0.7231 0.7131 0.7942 0.2723	99.6 99 100 100

Table 5 .

 5 2: Time and percentage of correct decompositions of CAPD for rank-4 tensors.

  Nous proposons un algorithme appelé CE dont l'idée est construite à partir d'un problème de valeur propre couplée. L'objectif de cet algorithme est de résoudre le problème d'approximation de rang-1 pour le cas d'ordre 3. Nous décrivons également des résultats théoriques assurant l'amélioration et la convergence de l'objectif du problème (B.3). • Performance de l'algorithme CE. La méthode CE présente d'avantage par rapport à ALS. D'abord, en termes de complexité numériques, CE est plus attractive si l'une des dimensions du tenseur est largement supérieur aux autres dimensions. En plus, nous montrons par des simulations que, pour deux types d'initialisations, l'erreur d'approximation obtenue par CE est plus petite que celle calculée par ALS. De même pour le nombre moyen d'itération pour la convergence. du tenseur T est mise à jour. Ensuite, on combine le résidu E

	pour construire le tenseur Y	(1) 1 avec la deuxième composante X 2 , lequel est décomposé de la même manière que Y (1) (1)	(0) 2

The CP decomposition is also called PARAFAC or CANDECOMP in the literature.

We also use equivalently the term essentially unique to indicating that the decomposition is unique with the exception of the elementary indeterminacies of scaling and permutation of factors.

This thesis is composed of a brief introduction, four chapters with results (chapter 2 to 5), and a summary of conclusions and perspectives. We also have included two appendices. The main points of each part of the manuscript is presented in the following.

The computational time of ALS and SeROAP algorithms is of around few milliseconds for

× 2 × 2 tensors.

See the explanation in next bulet to understand why the scenario 3 × 3 × 3 cannot be used to infer that our conjecture fails.

Tensor-based and rank-1 tensor algorithms are not the same techniques.

In chapters 1 and 2, we have used the variable k instead. Here k is already used as the index of the third mode of T .

The upcoming reasoning is identical if B or C were semi-unitary.

After the last deflation, the algorithm yields a rank-1 tensor, whose factors are easily obtained by computing an economic SVD.

The distance decreases Q-linearly to zero is lim k→∞ x kx ⋆ / x k+1x ⋆ = 0, for some x ⋆ ∈ A ∩ B.

It is the virtual zero or the floating-point relative accuracy of Matlab, whose value is 2 -52 .

Acknowledgements

I would like to thank the Conselho Nacional de Desenvolvimento Científico e Tecnológico in Brazil for funding this thesis, the director of the GIPSA-lab Jean-Marc Thiriet who gladly received me into the laboratory, and the specialist responsible of the doctorate school EEATS Gang Feng for helping me with important guidance during my PhD progress. I would like to express my gratitude to my thesis advisor, Pierre Comon, to whom I owe my progress as a researcher, for pointing carefully my mistakes and suggesting ideas to develop all this work. I thank my co-advisor André Lima Férrer de Almeida for all fruitful discussions and remarks during my thesis. I could not forget Professor João Cesar who guided me in the beginning of my research still in Brazil.

A special thanks to DECODA team, Souleymen, Yang, Rodrigo, Jeremy, Kostya, Cheng Feng, Miguel, Francesca, for all contributions that they made in the meetings. The discussions about some points of my work were really helpful.

I am grateful to my friends and family to encourage me in hard times, specially my wife, Débora, who was patient and was always by my side supporting me in all circumstances.

Main contributions of this chapter:

• Reduction of multivariate quadratic systems to best rank-1 tensor approximations. The best rank-1 tensor approximation is an optimization problem than can be applied to solve multivariate polynomial quadratic systems. We proof that a general quadratic system in the real field can be reduced to a best rank-1 approximation problem of a three-way real tensor. This tensor is constructed from the parameters of the system. The extension to specific complex systems is also done using the same approach. In this case, we show that the system is reduced to a conjugated partially symmetric rank-1 three-way tensor approximation problem, which can be tackled by an alternating algorithm.

• Rank-1 approximation algorithm. We propose an algorithm to compute the (conjugated) partially symmetric rank-1 three-way tensor approximation problem. This algorithm, called CPS, is developed from the stationary equations of the least squares problem min y,w Ty ⊗ y * ⊗ w 2 .

• Performance of the rank-1 tensor approach. ALS and CPS algorithms present satisfactory performance to extract a solution of a quadratic system in a reasonable time (e.g. a few seconds for a real square system of 10 variables/equations). We draw some examples and perform Monte Carlo experiments to evince the efficiency of our approach for real and complex systems.

Clearly, the factors b r and c r are directly obtained from the singular vectors of Unvec(x

r ) as shown in Algorithm (7).

After the estimation of a r , b r , and c r , we deflate the original tensor by the rank-1 component a r ( c r ⊠ b r ) T , and perform again R -2 times the alternating projection procedure in order to find the other components of the factor matrices A, B and C.

Remark 4 Neither the convergence of the CAPD algorithm to a point x ⋆ r ∈ M r ∩N nor the estimation of the right factors a r , b r and c r can be ensured. However, theoretical results in the next section guarantee local convergence and the extraction of the suitable factors.

Remark 5 Instead of estimating one-by-one the rank-1 components, we could try to estimate all components all-at-once from equation (5.2). However, the orthogonality of the estimated factor A would not be ensured anymore unless a constraint is imposed. Actually, all columns of the matrix factors could be collinear, as we will show in Section 5.5. Moreover, the estimation of all factors all-at-once requires more computations.

Convergence study on CAPD algorithm

This section deals with important results on the convergence of CAPD algorithm. We start with the proof of some essential lemmas, which will be important to demonstrate that the extraction of one rank-1 component of the tensor is ensured if CAPD method converges to some point in M r ∩ N for every r ∈ {1, 2, . . . , R -1}. Finally, based on the ideas within [Lewis 2008], we show that CAPD algorithm has locally linear convergence properties. ) be the unit rank-1 approximation of vector x

). The same reasoning can be applied to the next iterations, and the proof is complete.

r ) be the minimal norm solution at k-th iteration in CAPD algorithm. Then for all k ≥ 0, a (k) r is a linear combination of the columns of A.

Proof: It follows directly from the mode-1 unfolding T = A(C ⊙ B) T that the columns of A span Ran(T ). Yet, a

Theorem 5.4.3 Let T be a tensor whose factor A is a column-wise orthogonal matrix. Assume that T has an essentially unique decomposition, and c

Proof: Note that R can be viewed as the set of J × K rank-1 matrices. According to [Lewis 2008], the tangent space of R is given by

where {u 1 , u 2 , . . . , u J } and {v 1 , v 2 , . . . , v K } are the sets of left and right singular vectors of the matrix Unvec

Since the tangent space does not depend on the singular value σ 1 , and

Proposition 5.4.7 Let T ∈ R I×J×K be a rank-R tensor with factor A column-wise orthogonal. Assume neither C nor B have collinear columns. In CAPD algorithm, M r and N are transverse manifolds at any point

According to Lemma 5.4.6, the tangent space of N at x is given by

for left and right singular vectors

Note that all vectors in K(x) are orthogonal to one another so that they span a

is the only vector of K(x) in T Mr (x) so that any linear combination of at least two vectors in K(x) does not lie in T Mr (x). In other words, Span{K(x

Thus, in order to ensure transversality between the manifolds, we should have

Although Proposition 5.4.7 ensures the transversality between the manifolds M r and N and, thus, the linear convergence of our algorithm, the condition r ≤ R -(J -1)(K -1) is very restrictive.

Chapter 4

• In our thesis, we dedicated a chapter for tackling multivariate quadratic systems. We showed that this problem can be reduced to the best rank-1 approximation problem. Thus, some tensor tools, such as iterative algorithms, can be employed to extract solutions of a particular system. The reduction was proved for the real field. We also gave a list of advantages and drawbacks of our tensor approach.

• We generalized the tensor approach of the real case for a specific but more general complex system. We showed that our complex system can be reduced to a conjugated partially symmetric rank-1 approximation problem.

• To compute a conjugated partially symmetric rank-1 approximation of a tensor, we proposed an algorithm, called CPS, based on the stationary equations of (P9). We saw by numerical experiments, that this algorithm can extract one complex solution of the system.

Chapter 5

• we proposed an algorithm, called CAPD, combining the alternating projection and deflation methods to compute the exact decomposition of tensors with a column-wise orthogonal matrix.

• We showed that CAPD extracts all rank-1 components of the tensor if the iterates of the algorithm converge to points at the intersection of two known manifolds.

• If some transversality conditions are satisfied, we showed that the iterates of CAPD converge to a point at the intersection of manifolds. However, we saw that the transversality concept applied to our problem does not ensure the extraction of all rank-1 components of a tensor.

• To overcome the restrictive condition imposed by transversality, we used the non-tangential concept to show that the convergence of the CAPD algorithm can be ensured when starting points (initialization of the algorithm) are close to the intersection of predefined manifolds.

Thus, the extraction of all rank-1 components is ensured.

• We also showed that the extraction of rank-1 components one-by-one is less computationally complex than the estimation of the factor matrices all-at-once.

Perspectives

In the following, we draw the main directions for future works:

Chapter 2

• The employment of other algebraic geometry tools such as [Bucero 2014] to compute best rank-1 approximations can outperform the Lasserre's relaxation method, but it still concerns tensors with small dimensions. The monomials of the objective of Problem 2.2 are sparse, so that we A

Appendices

A.1 Iterative CP decomposition algorithms

This appendix presents the description of two iterative tensor decomposition algorithms used in our simulations: ALS [Smilde 2005, Comon 2009b] and HALS [Cichocki 2009c]. These algorithms are used to compute low rank CP tensor approximations.

In the following, the input parameter R denotes the rank of the output tensor. Assuming R 0 is the rank of the input tensor T , if R 0 ≤ R, then the algorithms perform an exact decomposition. On the other hand, if R 0 > R, a lower rank-R approximation is computed.

Alternating least squares

The most commonly used algorithm for solving the CP decomposition is ALS [Smilde 2005]. The goal is to update alternately each factor matrix in each iteration by solving a least squares problem conditioned on previous updates of the other factor matrices. The implementation is quite simple and it is detailed in Algorithm 8.

)V + end until some stopping criterion is satisfied;

Algorithm 8: ALS algorithm

Hierarchical ALS

The Hierarchical ALS algorithm was originally proposed in [Cichocki 2009c] in the context of nonnegative matrix factorizations. Herein, we present its higher order version for complex tensors. 107 p min = inf y p T y s.t. y 0 = 1, y is a sequence of moments with a representing measure on K.

For a sequence of moments, an infinite-dimensional moment matrix can be set up with entries indexed by N n , videlicet each element (α, β) of the matrix is given by y α+β , for α, β ∈ N n [Lasserre 2001].

Let M (y) be the moment matrix for y and define g (y) := M (y)g, with g the sequence of coefficients of a polynomial g(x) whose entries are sorted according to the same monomial order of y. The α-entry of g (y) is given by (g (y) ) α = β g β y α+β . From Lemma 4.2 within [Laurent 2009], the last constraint in Problem (A.3) can be relaxed by the weak conditions M (y) 0 and M (g

0, ∀i ∈ {1, 2, . . . , m}, which leads to the following relaxed optimization problem whose global solution is a lower bound for p min .

Since it is not clear how to proceed with infinity matrices, the previous problem must be relaxed by truncating the involved moment matrices. For this, define M t (y) the moment matrix for y ∈ N n 2t (the sequence of moments up to degree 2t). Let d p = ⌈deg(p(x))/2⌉, d g i = ⌈deg(g i (x))/2⌉ and d K = max(d g 1 , d g 2 , . . . , d gm ) (set d K = 1 if m = 0). For any integer t ≥ max(d p , d K ), it follows the sequence of SDP problems

One clearly notes that p mom t ≤ p mom ≤ p min . Problems (A.5) can be viewed as moment hierarchical problems: as t gets larger, p mom t is closer to p mom . The number of variables increases to n+2t

2t , as long as the lexicographic monomial order is considered.

Next, we present two known theoretical results on moment relaxations whose main ideas boil down to ensure an attainable infimum of Problems (A.5), and a certificate of optimality allowing to extract global minimizers of Problem (A.1).

Proposition A.2.1 [Laurent 2009] If the ball constraint R 2 -n i=1 x 2 i ≥ 0 is present in the description of the semialgebraic set K, then the feasible region of the relaxed problems (A.5) is bounded and its infimum is attained.

Chapitre: Déflation Itérative

Ce chapitre présente un étude détaillé sur l'utilisation d'algorithmes d'approximation de rang-1 dans le contexte de la déflation itérative (appelé DCPD dans notre manuscrit). L'objective est de calculer une approximation de rang faible pour un tenseur donné.

La déflation itérative combinée avec la meilleure approximation de rang-1 est illustrée dans la

1 est la meilleure approximation de rang-1 de T . Une fois que cette procédure de déflation ne réduit pas le rang en général [Stegeman 2007], il existe un résidu E (0) tel que T = E (0) + R r=1 X

(0) r . À la itération 1, on combine le résidu E (0) avec la première composante X (0) 1 de façon à obtenir le tenseur Y

(1) 1

1 comme la somme de sa meilleure approximation de rang-1, i.e. X

P2:

Dans ce cas, nous montrons que P2 peut être réduit à un problème de meilleure approximation de rang-1 conjugué et partialement symétrique, lequel est représenté par le problème suivant P3: Les contributions du chapitre sont les suivantes:

• Réduction du système quadratique au problème de meilleure approximation de rang-1. Nous montrons qu'une solution d'un système polynomial quadratique dans le corps des réels peut être obtenue à partir de la solution du problème (B.3) pour le cas d'ordre 3. De la même façon, on démontre que le cas le plus général P2 est réduit au problème P3. Ainsi, nous pouvons appliquer un nouvel approche pour traiter des systèmes quadratiques multivariés.

• Algorithme pour l'approximation de rang-1. Nous proposons un algorithme pour résoudre le problème de meilleure approximation de rang-1 conjugué et partiellement symétrique. Cet algorithme s'appelle CPS (Conjugated Partially Symmetric) et son idée est basée sur les équations des points stationnaires du problème P3.

• Performance de l'approche de rang-1. Les algorithmes itératifs ALS et CPS présentent des performances satisfaisantes lors de l'extraction d'une solution des systèmes quadratiques et dans un temps de simulation raisonnable. Nous présentons quelques exemples des systèmes et aussi un ensemble des simulations pour prouver l'efficacité de l'approche dans les cas réel et complexe.

Chapitre: Projection Alternée Appliquée aux Tenseurs Orthogonaux

L'objectif de ce chapitre est de résoudre la décomposition CP exacte dans le cas où une des matrices facteurs est semi-unitaire. Pour cela, on combine la méthode de déflation avec celle de projection alternée. Les concepts des transversalité et non-tangentialité sont également étudiés dans le contexte tensoriel pour évaluer la convergence d'un algorithme. Les principaux résultats sont les suivants:

• Algorithme CAPD. Nous proposons un algorithme appelé CAPD que combine la méthode des projections alternées avec la déflation de tenseurs de rang-1. Cela nous permet de résoudre la décomposition CP exacte d'un tenseur d'ordre 3 avec l'une des matrices facteurs semi-unitaire.