METHODS FOR COMPUTING THE CANONICAL POLYADIC DECOMPOSITION
MIKAEL SØRENSEN ∗§ ‡, LIEVEN DE LATHAUWER ∗§ ‡, PIERRE COMON ¶, SYLVIE ICART †, AND LUC DENEIRE †
Abstract. The aim of this paper is to compute a Canonical Polyadic Decomposition (CPD) of a tensor using a Simultaneous Generalized Schur Decomposition (SGSD). First we propose new methods to compute a CPD with triangular factor matrices. This leads to more relaxed and robust SGSD methods. Next, we propose new ways to deal with partially (Hermitian) symmetric CPDs via the SGSD approach. For the case when the CPD contains a partial Hermitian symmetry we also provide a new uniqueness result. Finally, numerical experiments are reported.
Key words. tensor, polyadic decomposition, canonical decomposition (CANDECOMP), paral-lel factor (PARAFAC), simultaneous generalized Schur decomposition, partially symmetric tensor, triangular matrices.
1. Introduction. A Simultaneous Generalized Schur Decomposition (SGSD) approach for computing a third-order Canonical Polyadic Decomposition (CPD) with a partial symmetry was proposed in [20]. It was later adapted to the general unsym-metric third-order CPD case in [5]. The SGSD method has been applied in blind separation of communication signals [20, 7], in multidimensional harmonic retrieval [11, 1] and in blind identification of underdetermined mixtures [8]. The SGSD method requires that two of the factor matrices of the CPD have full column rank. However, a CPD problem with one full column rank factor matrix can often be converted to a CPD problem with all factor matrices having full column rank [6].
The original SGSD method [20], here called SGSD1-CP, for computing a third-order CPD requires that all three of the involved factor matrices have full column rank, as will be explained in section 4. Moreover, it does not fully take the struc-ture among the involved triangular matrices into account, as will also be explained in section 4. We present new methods for computing a CPD with triangular fac-tor matrices. Based on these results we first propose a new SGSD method, called SGSD2-CP, for the computation of a third-order CPD with three full column rank factor matrices. It attempts to take the structure among the involved triangular ma-trices into account and can be understood as an improved and generalized version of the procedure for computing the parameters of a structured tensor [4]. Second, we present a more relaxed SGSD method, called SGSD3-CP, which only requires that
∗KU Leuven - E.E. Dept. (ESAT) - SCD-SISTA, Kasteelpark Arenberg 10, B-3001 Leuven-Heverlee, Belgium.
§Group Science, Engineering and Technology, KU Leuven Kulak, E. Sabbelaan 53, 8500 Kortrijk, Belgium, {Mikael.Sorensen, Lieven.DeLathauwer}@kuleuven-kulak.be.
¶GIPSA-Lab, CNRS UMR5216, Grenoble Campus, BP.46 F-38402 St Martin d’Heres Cedex, France, pierre.comon@gipsa-lab.grenoble-inp.fr
†Laboratoire I3S - UMR7271 - UNS CNRS, Les Algorithmes - Euclide-B, 06903 Sophia Antipolis, France, {icart, deneire}@i3s.unice.fr.
‡Research supported by: (1) Research Council KU Leuven: GOA-Ambiorics, GOA-MaNet, CoE EF/05/006 Optimization in Engineering (OPTEC), CIF1, STRT1/08/23, (2) F.W.O.: (a) project G.0427.10N, (b) Research Communities ICCoS, ANMMM and MLDM, (3) the Belgian Federal Sci-ence Policy Office: IUAP P7 (DYSCO II, Dynamical systems, control and optimization, 2012-2017), (4) EU: ERNSI.
The work of Mikael Sørensen was supported by the EU by a Marie-Curie Fellowship under contract No MEST- CT-2005-021175.
two of the matrix factors of the given CPD have full column rank. The proposed SGSD2-CP and SGSD3-CP methods mainly differ from the existing methods [20, 5] in the way the involved triangular matrices are computed. A part of this work has been presented in the conference paper [18]. CPDs containing partial (Hermitian) symmetries are common in many practical problems such as blind source separation based on simultaneous matrix diagonalization [3, 15, 8]. Hence, the third contribu-tion of this paper are new SGSD methods for computing a CPD containing partial (Hermitian) symmetries. For the case when the CPD contains a partial Hermitian symmetry we also provide a new uniqueness result.
The paper is organized as follows. The rest of the introduction will present our notation followed by a quick review of the CPD and SGSD. In section 2 we discuss the relation between the SGSD and CPD. Next, in section 3 we present a generalized version of the extended QZ method for computing the SGSD. Section 4 will review the original SGSD1-CP method. Next, in sections 5 and 6 we propose the two new SGSD2-CP and SGSD3-CP methods for computing a CPD. Section 7 will discuss new SGSD approaches for the computation of a CPD containing partial (Hermitian) symmetries while section 8 presents a new uniqueness result for a CPD with a partial Hermitian symmetry. Section 9 will explain how to extend the SGSD methods to Nth-order tensors. Finally, in section 10 we conclude the paper.
1.1. Notation. Vectors, matrices and tensors are denoted by lower case bold-face, upper case boldface and upper case calligraphic letters, respectively. The symbols ⊗ and ⊙ denote the Kronecker and Khatri-Rao products, defined as
A ⊗ B , a11B a12B . . . a21B a22B . . . .. . ... . .. , A ⊙ B , a1⊗ b1 a2⊗ b2 . . . ,
in which (A)mn = amn and ar and br denote the rth column vector of A and B,
respectively. Throughout the paper ar denotes the rth column vector of the matrix
A. Furthermore, the outer product of N vectors a(n)∈ CIn is denoted by a(1)◦ a(2)◦
· · ·◦a(N )∈ CI1×I2×···×IN and it satisfies the relation a(1)◦ a(2)◦ · · · ◦ a(N )
i1i2...iN =
a(1)i1 a(2)i2 · · · a(N )iN .
Further, (·)T, (·)∗, (·)H, (·)†, k · kF, |·|, Col (·), Null (·), (·)⊥, Re {·} and Im {·}
denote the transpose, conjugate, conjugate-transpose, Moore-Penrose pseudo-inverse, Frobenius norm, determinant, column space, null space, complementary subspace, real part and imaginary part of a matrix, respectively.
Let δij be the Dirac delta function which is zero except when the two indices are
equal, in which case δii = 1. The identity matrix and all-ones vector are denoted
by IR ∈ CR×R and 1R = [1, . . . , 1]T ∈ CR, respectively. Let 0M,N ∈ CM×N and
0M ∈ CM denote an all-zero matrix and an all-zero vector, respectively. The vector
e(R)p ∈ CR is equal to one at entry p and zero elsewhere. Matlab index notation will
be used to denote submatrices of a given matrix. For example, A(1 : k, :) denotes the submatrix of A consisting of the rows from 1 to k.
Let Blockdiag (X1, . . . ,XR) denote a block-diagonal matrix holding the matrices
{Xr} on its block-diagonal. The notation triu (·) and diag (·) is used to denote the
operator that sets the strictly lower triangular elements and off-diagonal elements of a matrix equal to zero, respectively. Let A ∈ CJ×K with J, K ≥ R and x ∈ R, then
triuR(A) ∈ CR×R and sign (x) are equal to
triuR(A) = triu (A (1 : R, 1 : R)) , sign (x) =
1 , x > 0 0 , x = 0 −1 , x < 0 .
Moreover, let A ∈ CI×J, then Vec (A) ∈ CIJ denotes the column vector defined
by (Vec (A))i+(j−1)I = (A)ij. Let a ∈ CIJ, then the reverse operation is Unvec (a) =
A ∈ CI×J such that (a)
i+(j−1)I = (A)ij. Let A ∈ CI×I, then Vecd (A) ∈ CI
denotes the column vector defined by (Vecd (A))i = (A)ii. Let A ∈ CI×J, then
Dk(A) ∈ CJ×J denotes the diagonal matrix holding row k of A on its diagonal.
The k-rank of a matrix A is denoted by k (A). It is equal to the largest integer k(A) such that every subset of k (A) columns of A is linearly independent.
1.2. Canonical Polyadic Decomposition (CPD). A N th-order rank-1 ten-sor X ∈ CI1×···×IN is defined as the outer product of non-zero vectors a(n) ∈ CIn,
n∈ {1, . . . , N}, such that Xi1...iN =
QN n=1a
(n)
in . The rank of a tensor X is equal to
the minimal number of rank-1 tensors that yield X in a linear combination. Assume that the rank of X is R, then it can be written as
X =
R
X
r=1
a(1)r ◦ · · · ◦ a(N )r , (1.1)
where a(n)r ∈ CIn. In this paper decomposition (1.1) will be referred to as the
Canon-ical Polyadic Decomposition (CPD) of X . Let us stack the vectors {a(n)r } into the
matrices
A(n)=h a(n)1 , · · · ,a(n)R i
∈ CIn×R, n∈ {1, . . . , N} . (1.2)
The matrices A(n)in (1.2) will be referred to as the CPD factor matrices. In this paper we say that the CPD of the tensor X in (1.1) is partially symmetric if A(m) = A(n) or A(m)= A(n)∗for some m 6= n.
The CPD of X in (1.1) is said to be unique if all the N-tupletsA(1), . . . ,A(N ) satisfying (1.1) are related via A(n) = A(n)Π∆A(n), ∀n ∈ {1, . . . , N} , where Π is a
permutation matrix and {∆A(n)} are diagonal matrices satisfying
QN
n=1∆A(n)= IR.
The following three matrix representations of a CPD of a third-order tensor X ∈ CI1×I2×I3 will be used throughout the paper. Let X(i1··)∈ CI2×I3 denote the matrix
such that X(i1··) i2i3 = Xi1i2i3, then X (i1··)= A(2)D i1 A(1)A(3)T and CI1I2×I3∋ X (1) , X(1··) .. . X(I1··) = A(2)D1 A(1) .. . A(2)DI1 A(1) A (3)T =A(1) ⊙ A(2)A(3)T.
Let the matrices X(·i2·) ∈ CI3×I1 be constructed such that X(·i2·)
i3i1 = Xi1i2i3, then X(·i2·)= A(3)D i2 A(2)A(1)T and CI2I3×I1∋ X (2) , X(·1·) .. . X(·I2·) = A(3)D1 A(2) .. . A(3)DI2 A(2) A (1)T =A(2) ⊙ A(3)A(1)T.
Let X(··i3)∈ CI1×I2 satisfy X(··i3) i1i2 = Xi1i2i3, then X (··i3)= A(1)D i3 A(3)A(2)T and CI1I3×I2∋ X (3) , X(··1) .. . X(··I3) = A(1)D1 A(3) .. . A(1)DI3 A(3) A (2)T =A(3) ⊙ A(1)A(2)T.
In practice, we are usually only interested in a low-rank approximation of X . Typically the least-squares cost function
fA(1),A(2),A(3)= I3 X i3=1 X(··i3) − A(1)Di3 A(3)A(2)T 2 F (1.3)
is used. Unfortunately a minimizer of (1.3) may not exist [10, 14]. That is, instead of a minimizer the cost function has an infimum and the minimization may yield diverging components [14].
1.3. Simultaneous Generalized Schur Decomposition (SGSD). Given a tensor X ∈ CI1×I2×I3 of rank R and I
1, I2≥ R, a SGSD of X is defined by
X(··i3)= QR(i3)ZT,
∀i3∈ {1, . . . , I3} , (1.4)
where Q ∈ CI1×I1 and Z ∈ CI2×I2 are unitary matrices and
R(i3)= " R(i3) 0R,I2−R 0I1−R,R 0I1−R,I2−R # ∈ CI1×I2 (1.5) in which R(i3)
∈ CR×R are upper triangular matrices. In practice, the rank of X may
be larger than R and an approximate SGSD is of interest. The unknown matrices Q, Z and {R(i3)
} will be found by minimizing the least-squares cost function
gQ, Z, {R(i3)}= I3 X i3=1 X(··i3)− QR(i3)ZT 2 F , (1.6)
where Q ∈ CI1×I1 and Z ∈ CI2×I2 are unitary matrices and R(i3)∈ CI1×I2 is of the
form (1.5).
It can be shown in a similar way as in [1], that the problem of minimizing the cost function (1.6) is equivalent to the maximization of the cost function
h(Q, Z) = I3 X i3=1 triuR QHX(··i3)Z∗ 2 F . (1.7)
Let us stack the entries of the upper triangular part in r(i3)=hR(i3)(1, 1), R(i3)(1 : 2, 2), . . . , R(i3)(1 : R, R)i
T
∈ CR(R+1)2 .
We have r(i3)= WHVec
X(··i3)with W given by
Q =q1, . . . ,qI1 ∈ CI1×I1 Z = [z1, . . . ,zI2] ∈ C I2×I2 W = [z1⊗ q1,z2⊗ q1,z2⊗ q2, . . . ,zR⊗ q1, . . . ,zR⊗ qR] ∈ CI1I2× R(R+1) 2 .
Due to the reformulation of the SGSD problem from (1.6) into (1.7), as suggested in [1], it is clear that an optimal solution for the SGSD of X exists. Indeed, due to the continuity of the function (1.7) and the compactness of its domain one can conclude that an optimal solution to (1.7) always exists which implies that an optimal solution to (1.6) also always exists.
2. Link between SGSD and CPD. The computation of CPD via SGSD was proposed in [5, 20]. In this section we will formulate more explicitly than in [5, 20] how the two decompositions are connected. The derivation will make use of the following lemma.
Lemma 2.1. Consider A ∈ CI×Rwith k (A) ≥ 2. Then there exists a nonsingular matrix N ∈ CI×I for which A = NA is such that D
j1 A
−1
Dj2 A
exists and has distinct non-zero diagonal elements for some j1, j2∈ {1, . . . , I}.
Proof. A proof can for instance be found in [19].
In the following proposition a generalized eigenvalue is a scalar λ ∈ C which reduces the rank of a pencil, say X(··1)− λX(··2)∈ CI1×I2.
Proposition 2.2. Given a rank-R tensor X ∈ CI1×I2×I3 of which the CPD
involves the factor matrices A(n) ∈ CIn×R, n ∈ {1, 2, 3}, such that it has matrix
representations
X(··i3)= A(1)D i3
A(3)A(2)T, i3∈ {1, . . . , I3} . (2.1)
Assume that A(1) and A(2) have full column rank and kA(3) ≥ 2. If there exist unitary matrices Q ∈ CI1×I1 and Z ∈ CI2×I2 such that
R(··i3)= QHX(··i3)Z∗ = " R(··i3) 0R,I2−R 0I1−R,R 0I1−R,I2−R # ∈ CI1×I2, i 3∈ {1, . . . , I3},(2.2)
where R(··i3)∈ CR×R are upper triangular matrices, then there exists a permutation
matrix P ∈ RR×R and a diagonal matrix D ∈ CR×R such that
A(1)P= QR, A(2)P= ZL, A(3)PD= VecdR(··1) T VecdR(··2)T .. . VecdR(··I3) T , (2.3) in which R= R 0I1−R,R ∈ CI1×R, R∈ CR×R is upper triangular, (2.4) L= L 0I2−R,R ∈ CI2×R, L∈ CR×R is lower triangular. (2.5)
The other way around, if A(1) admits the QR-factorization A(1)= QR and A(2) the QL-factorization A(2) = ZL, then (2.2) holds with {R(··i3)
} upper triangular. Proof. We first show the second part of the proposition. Substitution of A(1) = QR and A(2)= ZL in (2.1) immediately yields (2.2) with R(··i3)= RDi3
A(3)LT upper triangular.
We now show the first part of the proposition. From Lemma 2.1 we know that kA(3) ≥ 2 implies that there exists a nonsingular matrix N ∈ CI3×I3 for which
A(3)= NA(3) is such that Dj1 A(3) −1 Dj2
A(3)exists and has distinct non-zero diagonal entries for some j1, j2∈ {1, . . . , I3}. We denote X(1)= X(1)NT. We have
X(··j1)= A(1)Dj1 A(3)A(2)T (2.6) X(··j2)= A(1)Dj2 A(3)A(2)T. (2.7)
The distinct diagonal entries of Dj1
A(3)
−1
Dj2
A(3) are the generalized eigen-values of the pencil X(··j1)
− λX(··j2)
. The matrices of right and left eigenvectors are
A(2)T† and A(1)†, respectively.
Unitary matrices Q ∈ CI1×I1 and Z ∈ CI2×I2 satisfying (2.2) yield a Generalized
Schur Decomposition of the pair (X(··j1),X(··j2)): R(··j1)= QHX(··j1)Z∗, R(··j1)= triu R R(··j1) (2.8) R(··j2)= QHX(··j2)Z∗, R(··j2)= triu R R(··j2), (2.9) where R(··j1) and R(··j2)
are upper triangular. The diagonal entries of the diagonal matrix diagR(··j1)
diagR(··j2)−1
also correspond to the generalized eigenvalues of the pencil X(··j1)− λX(··j2). There exists a permutation matrix P ∈ RR×R such
that PTDj1 A(3) −1 Dj2 A(3)P = diagR(··j1)diagR(··j2) −1 . The pencil of upper triangular matrices R(··j1)
− λR(··j2)
has upper triangular matrices of right and left eigenvectors, which we denote by ER ∈ CR×R and EL ∈ CR×R, respectively.
Consider ER= ER 0I1−R,R ∈ CI1×R, E L= EL 0I2−R,R ∈ CI2×R,
then ERand ELare matrices of right and left eigenvectors, respectively, of the pencil
R(··j1)− λR(··j2). Combining (2.6)–(2.7) with (2.8)–(2.9) we obtain as in (2.3) that
PTA(2)T†= Z∗E R⇔ A(2)P = Z E†RT = ZL , (2.10) PTA(1)†= ELQH ⇔ A(1)P = Q (EL)†= QR , (2.11) where R = R 0I1−R,R = (EL)†, L = L 0I2−R,R =E†R T
in which R and LT are upper triangular matrices. Since ER and EL have full
col-umn rank, the matrices L and R must also have full colcol-umn rank. By inserting (2.10) and (2.11) into (2.2) we get R(··i3) = RPD
i3
A(3)PTLT. Let R(··i3)
R(··i3)(1 : R, 1 : R) and D = diag Rdiag L and recall that PTD i3 A(3)P = Di3 A(3)P, then
diagR(··i3)= diagRPTDi3
A(3)PLT = PTDi3 A(3)Pdiag Rdiag L = Di3 A(3)PD, ∀i3∈ {1, . . . , I3} .
We finally obtain DPTA(3)T =hVecdR(··1), . . . ,VecdR(··I3)i.
Recall that the factor matrices A(1), A(2) and A(3) of the (approximate) CPD of the tensor X ∈ CI1×I2×I3 are typically found by minimizing the cost function
fA(1),A(2),A(3)= I3 X i3=1 X (··i3) − A(1)Di3 A(3)A(2)T 2 F. (2.12)
Let A(1) = QR be the QR factorization of A(1) where Q ∈ CI1×I1 is a unitary
matrix and R ∈ CI1×R is a matrix of the form (2.4). Similarly, let A(2) = ZL
be the QL factorization of A(2) where Z ∈ CI2×I2 is a unitary matrix and L ∈
CI2×R is a matrix of the form (2.5). Due to Proposition 2.2 we can replace the
matrices {RDi3
A(3)LT} in (2.12) by the upper triangular matrices {R(··i3)} and
still preserve the uniqueness of the solution. This means that the SGSD method will attempt to minimize the cost function
gQ, Z, {R(··i3)}= I3 X i3=1 X(··i3)− QR(··i3)ZT 2 F. (2.13)
We know from subsection 1.3 that minimizing the cost function (2.13) is equivalent to maximizing the cost function
h(Q, Z) = I3 X i3=1 triuR QHX(··i3)Z∗ 2 F. (2.14)
Thus, the SGSD-CP method first computes the unitary matrices Q and Z that make the matrices {X(··i3)} as upper triangular as possible by maximizing (2.14). This means that we have gone from the original non-unitary optimization problem (2.12), which may not have an optimal solution, to the unitary optimization problem (2.14) in which an optimal solution always exists.
3. Generalized version of the extended QZ method. Various methods have been proposed to numerically maximize the cost function (2.14). For instance, a two-sided Jacobi method for the case when the CPD is real valued has been proposed in [5]. Moreover, a one-sided Jacobi method has been proposed in [21] which is also valid for complex data. In [20] the extended QZ method was proposed for the case when I1= I2= R. In this section we generalize the extended QZ method to the case
I1, I2≥ R. The unitary matrices Q and Z are generated as
Q = S Y s=1 IY1−1 r=0 Q[r, s] , Z = S Y s=1 IY2−1 r=0 Z[r, s],
where S denotes the number of executed sweeps, Q[r, s]H = Ir−1 0r−1,I1−r+1 0I1−r+1,r−1 Us ∈ CI1×I1 Z[r, s] = 0R−r,rAs 0Ir,R−rR−r 0R−r,ICs2−R Bs 0I2−R,R−r Ds ∈ CI2×I2,
where Us ∈ C(I1−r+1)×(I1−r+1) is unitary and As ∈ Cr×r, Bs ∈ C(I2−R)×r, Cs ∈
Cr×(I2−R) and D
s∈ C(I2−R)×(I2−R) are such that Z[r, s] is unitary. In the extended
QZ method the matrices Q and Z are updated in alternating manner. The updating rules are X(··i3) ← Q[r, s]HX(··i3), Q ← Q · Q[r, s] and X(··i3) ← X(··i3)Z[r, s]∗,
Z ← Z · Z[r, s].
The matrix Q[r, s]H is designed to minimize the below-diagonal part of the rth
column of the matrices {X(··i3)}. Let eX(r)∈ C(I1−r+1)×I3 be defined by
e
X(r)=hX(··1)(r : I1, r), · · · , X(··I3)(r : I1, r)
i
, r∈ {1, . . . , min (R, I1− 1)} ,
and let eX(r) = UΣVH denote the compact SVD of eX(r). A unitary matrix Us that
minimizes the Frobenius norm of
UsXe (r)
(2 : I1− r + 1, :) is Us= UH.
The matrix Z[r, s] is designed to minimize the part of the matrices {X(··i3)}
that corresponds to the vectors {X(··i3)(r, 1 : r − 1)} and {X(··i3)(r, R + 1 : I
2)}. Let b X(r)∈ CI3×(r+I2−R)be defined by b X(r) = X(··1)(r, 1 : r) , X(··1)(r, R + 1 : I2) .. . X(··I3)(r, 1 : r) , X(··I3)(r, R + 1 : I 2) , r∈ {2 − sign (I2− R) , . . . , R} ,
and let bX(r)= UΣVH denote the compact SVD of bX(r). A unitary matrix Z∗s=
Vs Xs
Ws Ys
that minimizes the Frobenius norm of b X(r) Vs Xs Ws Ys (:, [1 : r − 1, r + 1 : r+ I2− R]) is Z∗s= [V(:, 2 : r), V(:, 1), V(:, r + 1 : r + I2− R)].
An outline of the generalized version of the extended QZ procedure is presented as Algorithm 1.
The following sections 4, 5 and 6 will discuss methods to find A(1), A(2)and A(3) once Q and Z have been found.
4. SGSD1-CP. The original SGSD1-CP method proposed in [20] will be re-viewed in this section. The working assumption is that the factor matrices A(n) ∈ CIn×R, n ∈ {1, 2, 3} have full column rank. Let X ∈ CI1×I2×I3 be a tensor of rank
Rconstructed from the full column rank factor matrices A(n)∈ CIn×R, n ∈ {1, 2, 3},
such that
X(··i3)= A(1)D i3
Algorithm 1 Outline of the extended QZ method when I1, I2≥ R.
Initialize: Q and Z Repeat until convergence for r = 1 to min (R, I1− 1) do e X(r)=hX(··1)(r : I1, r), · · · , X(··I3)(r : I1, r) i UΣVH= svd e X(r) X(··i3)(r : I 1,:) ← UHX(··i3)(r : I1,:) Q(:, r : I1) ← Q(:, r : I1)U end for for r = R to max (1, R − I2+ 2) do b X(r)= X(··1)(r, 1 : r) , X(··1)(r, R + 1 : I2) .. . X(··I3)(r, 1 : r) , X(··I3)(r, R + 1 : I 2) UΣVH= svd b X(r) V ← [V(:, 2 : r), V(:, 1), V(:, r + 1 : r + I2− R)] X(··i3)(:, [1 : r, R + 1 : I 2]) ← X(··i3)(:, [1 : r, R + 1 : I2])V Z(:, [1 : r, R + 1 : I2]) ← Z(:, [1 : r, R + 1 : I2])V∗ end for
Assume that Q and Z have been found and let R(··i3)
= triuR
QHX(··i3)Z∗. We
know from Proposition 2.2 that A(3) = hVecdR(··1), . . . ,VecdR(··I3)iT
up to the inherent column scaling and permutation ambiguities. Assume that A(3) has been found and compute F = X(1)
A(3)T†. In the exact case, we know that F = A(1)⊙ A(2) and therefore A(1) and A(2) follow from the R decoupled best rank-1 approximation problems min a(1)r ,a(2)r fr− a(1)r ⊗ a(1)r 2 F = mina(1)r ,a(2)r Unvec (fr) − a(2)r a(1)Tr 2 F , r∈ {1, . . . , R} .
In the presence of noise a refinement of the factor matrices {A(n)} is often necessary. The subsequent refinement step should be implemented by a method which fully takes the CPD structure of the problem into account. For optimization methods, see [17] and references therein.
5. SGSD2-CP. In this section the SGSD2-CP method for computing a third-order CPD will be presented. This method attempts to take the relations among the triangular matrices {R(··i3)
} into account when computing A(3) whereas the SGSD1-CP method only uses the diagonals of the triangular matrices {R(··i3)} when
comput-ing A(3). As in the SGSD1-CP method, the working assumption of the SGSD2-CP method is that all three factor matrices A(n)∈ CIn×R, n ∈ {1, 2, 3} have full column
rank.
column rank factor matrices A(n)∈ CIn×R, n ∈ {1, 2, 3}, such that X(··i3)= A(1)D i3 A(3)A(2)T ∈ CI1×I2, ∀i 3∈ {1, . . . , I3} .
Assume that Q and Z have been found. Construct the tensor R ∈ CR×R×I3 such
that Ri1i2i3 = R (··i3) i1i2 = triuR QHX(··i3)Z∗ i1i2 , then R(··i3) = RD i3 A(3)LT. Moreover, let R(i1··) ∈ CR×I3 satisfy the relation R(i1··)
i2i3 = Ri1i2i3, then R (i1··) =
LDi1 R
A(3)T. Stack the matrices {R(i1··)
} into the matrix
CR2×I3 ∋ R(1)= R(1··) .. . R(R··) = LD1 R .. . LDR R A(3)T = R ⊙ LA(3)T.
We will now present an algorithm to compute the factor matrices from Col R(1)
. Let R(1)= UΣWHdenote the compact SVD of R(1), where U ∈ CR
2×R
and W ∈ CI3×R
are columnwise orthonormal matrices and Σ ∈ CR×R is a positive diagonal matrix.
The Khatri-Rao product of two full column rank matrices has full column rank and consequently the matrix R⊙ L has full column rank. Since the matrix A(3) has full column rank by assumption we can now conclude that there exists a nonsingular matrix M ∈ CR×R such that
W∗M−T = A(3) (5.1)
UΣM = R⊙ L ⇔ UΣmr= rr⊗ lr, r∈ {1, . . . , R} . (5.2)
Let ar= rr(1 : r) and br= lr(r : R), then (5.2) is equivalent to
UΣmr= X 1≤p≤r X r≤q≤R e(R)p ⊗ e(R)q (ar⊗ br) = E(r)(ar⊗ br) , r∈ {1, . . . , R} , (5.3) where E(r) = [e1⊗ er, . . . ,e1⊗ eR, . . . ,er⊗ er, . . . ,er⊗ eR] ∈ CR 2×r(R−r+1) . , (5.4) Equation (5.3) can be written as
UΣmr= E(r)(ar⊗ br) ⇔ h UΣ, −E(r)i a mr r⊗ br = 0R2, r∈ {1, . . . , R} . (5.5) Proposition 5.2 will show that (5.5) has a unique solution. Moreover, it will provide us with a method to solve the system (5.5). In order to prove it we need the following lemma.
Lemma 5.1. Let R ∈ CR×R be a nonsingular upper triangular matrix and let L∈ CR×R be a nonsingular lower triangular matrix. Then the matrix S = R ⊙ L ∈
CR2×R has the following property
Proof. The proof is straightforward, noting that S (t + (t − 1)R, u) = (R)tu(L)tu (5.6) and (R)tu∈ C , u < t (R)tu6= 0 , u= t (R)tu= 0 , u > t (5.7) (L)tu= 0 , u < t (L)tu6= 0 , u= t (L)tu∈ C , u > t (5.8)
By combining (5.6), (5.7) and (5.8) we can conclude that
S (r + (r − 1)R, :) = αre(R)Tr , αr∈ C, ∀r ∈ {1, . . . , R} .
Proposition 5.2. Let R = [r1, . . . , rR] ∈ CR×R be a nonsingular upper
triangu-lar matrix and let L = [l1, . . . , lR] ∈ CR×R be a nonsingular lower triangular matrix.
Let E(r) be of the form (5.4). If the column vectors of U ∈ CR2×R
constitute a basis for Col (R ⊙ L), then the matrix
h
U,−E(r)i∈ CR2×(R+r(R−r+1)) (5.9)
has rank R + r(R − r + 1) − 1.
Proof. Consider the basis {rm⊗ lm} for Col (R ⊙ L). On one hand, rr⊗ lr is a
linear combination of the columns of E(r). On the other hand, by Lemma 5.1, vector rn⊗ ln, with n 6= r has a non-zero at position n + (n − 1)R, while vectors rp⊗ lp,
n6= p 6= r, and the columns of E(r) are zero at that position. Hence, matrix h
R ⊙ L, −E(r)i
has rank r(R −r +1)−1. In conclusion, the matrix (5.9) has rank R+r(R−r +1)−1. Proposition 5.2 tells us that we may ignore the Kronecker product structure ar⊗
br when solving (5.5) since the matrix
h
UΣ, −E(r)i has a one-dimensional kernel. We now reduce the number of equations by projecting on the orthogonal complement of ColE(r).
Corollary 5.3. Let R = [r1, . . . , rR] ∈ CR×R and L = [l1, . . . , lR] ∈ CR×R
be a nonsingular upper and lower triangular matrix, respectively. Define the diagonal projection matrix
CR2×R2 ∋ D(r) = IR2− E(r)E(r)T, (5.10)
where E(r)∈ CR2×r(R−r+1) is of the form (5.4). If the column vectors of U ∈ CR2×R
constitute a basis for Col (R ⊙ L), then the matrix
has rank R − 1, ∀r ∈ {1, . . . , R}.
Proof. Consider the columns {rs⊗ ls} of R ⊙ L. On one hand, rr⊗ lr is a linear
combination of the columns of E(r) and is therefore set equal to zero by D(r). On the other hand, by Lemma 5.1, vector rn⊗ ln, with n 6= r has a non-zero at position
n+ (n − 1)R, while vectors rp⊗ lp, n 6= p 6= r, are zero at that position. The entries at
position n + (n − 1)R, 1 ≤ n 6= r ≤ R, are not set equal to zeros by D(r). We conclude that D(r)(R ⊙ L) has rank R − 1. Hence, any matrix U of which the columns form a basis for Col (R ⊙ L) yields a matrix D(r)U of rank R − 1.
By projecting (5.5) on the orthogonal complement of E(r) we obtain
U(r)mr= 0R2, r∈ {1, . . . , R} (5.12)
with U(r)= D(r)UΣ and D(r)defined in (5.10). According to Corollary 5.3, U(r)still has a one-dimensional kernel. In practice, we first remove the all-zero row-vectors of U(r) to obtain eU(r)∈ C(R2−r(R−r+1))×R, and we solve
e
U(r)mr= 0R2−r(R−r+1), r∈ {1, . . . , R} . (5.13)
Note that the set of R2 equations in R + r(R − r + 1) unknowns (5.5) has been
reduced to a set of R2− r(R − r + 1) homogeneous equations in R unknowns. In
the case of an approximate tensor decomposition, we set mr equal to the singular
vector associated with the smallest singular value of eU(r). Once M = [m1, . . . ,mR]
has been obtained from (5.12), then from (5.1) we get A(3) = W∗M−T. Compute F = X(1)
A(3)T† = X(1)MTWT, then as in the SGSD1-CP method, the matrices
A(1) and A(2)now follow from the R decoupled best rank-1 approximation problems min a(1)r ,a(2)r Unvec (fr) − a(2)r a(1)Tr 2 F, r∈ {1, . . . , R} .
Finally, we obtain A(3) via relation A(3)T =A(1)⊙ A(2)†X(1).
6. SGSD3-CP. In this section we present the SGSD3-CP method for computing a third-order CPD where only two of the factor matrices are required to have full column rank. This is in contrast to the SGSD1-CP and SGSD2-CP methods which both require that all three of the involved factor matrices have full column rank. In subsection 6.1 we first consider the case when kA(3)≥ 2 while in subsection 6.2 we relax this constraint to kA(3)≥ 1.
6.1. SGSD3-CP with kA(3)≥ 2. Let the third-order tensor X ∈ CI1×I2×I3
of rank R be constructed from the two full column rank factor matrices A(1)∈ CI1×R,
A(2)∈ CI2×Rand the factor matrix A(3)∈ CI3×R with k
A(3)≥ 2 such that X(··i3)= A(1)D i3 A(3)A(2)T, i3∈ {1, . . . , I3} .
Assume that Q and Z have been found and denote R(··i3) = triu R
QHX(··i3)Z∗.
As in the SGSD1-CP method, we find A(3) from the relation A(3)T =hVecdR(··1), . . . ,Vecd (R)(··I3)i
Let the tensor R ∈ CR×R×I3 be constructed such that R i1i2i3 = R (··i3) i1i2 , then R (··i3)= RDi3
A(3)LT. Moreover, let R(·i2·)∈ CI3×R satisfy the relation R(·i2·)
i3i1 = Ri1i2i3,
then R(·i2·)= A(3)D i2 L
RT. Stack the matrices {R(·i2·)} into the matrix
CRI3×R∋ R (2) = R(·1·) .. . R(·R·) = A(3)D1 L .. . A(3)DR L RT =L ⊙ A(3)RT.
Let R(2) = UΣVH denote the compact SVD of R(2). Since L has full column rank,
the matrix L ⊙A(3)has full column rank. This implies that there exists a nonsingular matrix M ∈ CR×R such that
UM = L⊙ A(3)⇔ Umr= lr⊗ a(3)r , r∈ {1, . . . , R} . (6.1) Let br= lr(r : R) ∈ CR−r+1 and E(r)=he(R)r ,e (R) r+1, . . . ,e (R) R i ∈ CR×(R−r+1), (6.2) then lr= E(r)br and lr⊗ a(3)r = E(r)⊗ a(3)r br. (6.3)
By inserting (6.3) into (6.1) we obtain Umr= E(r)⊗ a(3)r br⇔ h U, −E(r)⊗ a(3)r i mr br = 0R2, r∈ {1, . . . , R} .(6.4)
Due to the following Proposition 6.1 the linear system (6.4) has a unique solution. The least-squares solution to (6.4) is given by the right singular vector of the matrix h
U, −E(r)⊗ a(3)r
i
associated with its smallest singular value. Proposition 6.1. Given A(3) ∈ CI3×R with k
A(3)≥ 2 and let L ∈ CR×R be
a nonsingular lower triangular matrix. Let E(r) be of the form (6.2). If the column vectors of U ∈ CRI3×R constitute a basis for Col
L⊙ A(3), then the matrix h
U,−E(r)⊗ a(3)r i∈ CRI3×(2R−r+1) (6.5)
has rank 2R − r.
Proof. Consider the columns {lr⊗ a(3)r } of Col
L ⊙ A(3), then U = [u1, . . . ,uR] = h l1⊗ a(3)1 , . . . ,lR⊗ a(3)R i = l11a(3)1 0I3 · · · 0I3 l21a(3)1 l22a(3)2 . .. ... .. . ... . .. 0I3 lR1a(3)1 lR2a(3)2 · · · lRRa(3)R . We also have E(r)⊗ a(3)r = h e(R)r ⊗ a(3)r , . . . ,e (R) R ⊗ a (3) r i .
Since ur= lr⊗ a(3)r = R X q=r lqre(R)q ⊗ a(3)r
the problem of determining the rank of h
U, −E(r)⊗ a(3)r
i
∈ CRI3×(2R−r+1)
reduces to finding the rank of h u1, . . . ,ur−1,ur+1, . . . ,uR,e(R)r ⊗ a(3)r ,e (R) r+1⊗ a(3)r , . . . ,e (R) R ⊗ a (3) r i ∈ CRI3×(2R−r). (6.6) By reorganizing the column vectors of (6.6) we obtain
h u1 . . . ur−1 e(R)r ⊗ a(3)r ur+1 e(R)r+1⊗ a (3) r . . . uR e(R)R ⊗ a (3) r i = l1,1a(3)1 0I3 0I3 0I3 0I3 0I3 0I3 .. . . .. ... ... ... ... ... ... lr−2,1a(3)1 0I3 0I3 0I3 0I3 0I3 0I3 lr−1,1a(3)1 lr−1,r−1a(3)r−1 0I3 0I3 0I3 0I3 0I3 lr,1a(3)1 lr,r−1a(3)r−1 a (3) r 0I3 0I3 0I3 0I3 lr+1,1a(3)1 lr+1,r−1a(3)r−1 0I3 lr+1,r+1a (3) r+1 a (3) r 0I3 0I3 lr+2,1a(3)1 lr+2,r−1a(3)r−1 0I3 lr+2,r+1a (3) r+1 0I3 0I3 0I3 .. . ... ... ... . .. ... ... lR,1a(3)1 lR,r−1a(3)r−1 0I3 lR,r+1a (3) r+1 0I3 lR,Ra (3) R a (3) r . (6.7) The matrix (6.7) is lower block-triangular. Since kA(3)≥ 2 and L is nonsingular the matrices on the block-diagonal of (6.7) have full column rank which means that the matrix (6.7) has full column rank. Hence, the matrix (6.5) has rank 2R − r. Consequently any matrix U of which the columns form a basis for ColL ⊙ A(3) yields a matrix of rank 2R − r.
After L ∈ CR×R has been computed, then via b
r in (6.4) we obtain L =
[LT,0T I2−R,R]
T ∈ CI2×R and A(2) = Z∗L. Finally, we compute A(1) via relation
A(1)T =A(2)⊙ A(3)†X(2).
6.2. SGSD3-CP with kA(3)= 1. Assume now that the third-order tensor X ∈ CI1×I2×I3 of rank R is constructed from the two full column rank factor matrices
A(1) ∈ CI1×R and A(2) ∈ CI2×R but the factor matrix A(3) ∈ CI3×R now only has
kA(3)= 1. In this case A(3) can be written as
A(3) =h1TL1⊗ c1, . . . ,1 T LQ⊗ cQ
i
up to a permutation and scaling of the columns of A(3) and R =PQq=1Lq. We can
write the decomposition of X as follows
X = Q X q=1 AqBTq ◦ cq, (6.9) where A(1)= [A1, . . . ,AQ] , Aq ∈ CI1×Lq A(2)= [B1, . . . ,BQ] , Bq ∈ CI2×Lq.
The decomposition (6.9) is known as the (Lq, Lq,1)-Block Term Decomposition (BTD)
of X [9]. Since AqBTq = AqFqF−1q B T
q for any nonsingular matrix Fq ∈ CLq×Lq we
can from X in (6.9) only find A(1) and A(2)up to a column permuted block-diagonal matrix. We now explain that the SGSD3-CP method is still valid even if kA(3)= 1. It will based on generalizations of Propositions 2.2 and 6.1.
Proposition 6.2. Given a rank-R tensor X ∈ CI1×I2×I3 of which the (L
q, Lq,
1)-BTD involves the factor matrices A(n)∈ CIn×R, n ∈ {1, 2, 3}, such that it has matrix
representations
X(··i3)= A(1)D i3
A(3)A(2)T, i3∈ {1, . . . , I3} . (6.10)
Assume that A(1) and A(2) have full column rank and A(3) is of the form (6.8) with k(C) ≥ 2. If there exist unitary matrices Q ∈ CI1×I1 and Z ∈ CI2×I2 such that
R(··i3)= QHX(··i3)Z∗= " R(··i3) 0R,I2−R 0I1−R,R 0I1−R,I2−R # ∈ CI1×I2, i 3∈ {1, . . . , I3}, (6.11)
where R(··i3)∈ CR×R are upper triangular matrices, then there exists a permutation
matrix P ∈ RR×R, a diagonal matrix D ∈ CR×R and nonsingular matrices F q ∈ CLq×Lq such that A(1)FP= QR , A(2)F−1P= ZL , A(3)PD= VecdR(··1) T VecdR(··2) T .. . VecdR(··I3) T ,(6.12)
in which F = Blockdiag (F1, . . . , FQ), R ∈ CI1×R and L ∈ CI2×R are of the forms
(2.4) and (2.5), respectively. Conversely, if A(1) admits the QR-factorization A(1)= QRand A(2) the QL-factorization A(2)= ZL, then (6.11) holds with {R(··i3)} upper triangular.
Proof. The second part of the proposition is obvious. Indeed, by substitut-ing A(1) = QR and A(2) = ZL in (2.1) immediately yields (6.11) with R(··i3) = RDi3
A(3)LT upper triangular.
We now show the first part of the proposition. Since k (C) ≥ 2 we know from Lemma 2.1 that there exists a nonsingular matrix N ∈ CI3×I3 for which C = NC
is such that Dj1 C
−1
Dj2 C
exists and has distinct non-zero diagonal entries for some j1, j2∈ {1, . . . , Q}. We denote X(1)= X(1)NT. We have
X(··j1) = A(1)Dj1 A(3)A(2)T (6.13) X(··j2) = A(1)Dj2 A(3)A(2)T. (6.14)
The diagonal entry of Dj1
A(3)
−1
Dj2
A(3)denoted by λqis the generalized
eigen-values of the pencil X(··j1)
− λqX (··j2)
with algebraic multiplicity Lq. The matrices
of right and left eigenvectors areA(2)T† and A(1)†, respectively. Since Lq ≥ 1 the
right and left eigenvectors can only be unique up to a column permuted block-diagonal matrix of the form of F = Blockdiag (F1, . . . ,FQ) in which Fq ∈ CLq×Lq.
Unitary matrices Q ∈ CI1×I1 and Z ∈ CI2×I2 satisfying (6.11) yield a Generalized
Schur Decomposition of the pair (X(··j1)
,X(··j2) ): R(··j1)= QHX(··j1)Z∗, R(··j1)= triu R R(··j1) (6.15) R(··j2)= QHX(··j2)Z∗, R(··j2)= triu R R(··j2) (6.16) where R(··j1) and R(··j2)
are upper triangular. The diagonal entries of the diagonal matrices diagR(··j1)
diagR(··j2)−1
also correspond to the generalized eigenvalues of the pencil X(··j1)− λX(··j2). Combining (6.13)–(6.14) with (6.15)–(6.16) we obtain the first part of the proposition in the same way as explained in Proposition 2.2.
Proposition 6.3. Given A(3)=h1TL1⊗ c1, . . . , 1 T LQ⊗ cQ i ∈ CI3×R, C= [c 1, . . . , cQ] ∈ CI3×Q (6.17)
with k (C) ≥ 2 and R =PQq=1Lq. Let L ∈ CR×R be a nonsingular lower triangular
matrix and let E(r) be of the form (6.2). Let the column vectors of U ∈ CRI3×R
constitute a basis for ColL⊙ A(3). Denote L(q)=Pqp=1Lp, then the matrix
G(r) =hU,−E(r)⊗ a(3) r
i
∈ CRI3×(2R−r+1) (6.18)
has rank 2R − r + Lq when r ∈ {L(q−1)+ 1, . . . , L(q)}. Furthermore, when r ∈
{L(q−1)+ 1, . . . , L(q)}, then the column vectors of V(r) ∈ C(2R−r+1)×(L(q)−r+1) given
by V(r)= mL(q−1)+r bL(q−1)+r , mL(q−1)+r+1 0 bL(q−1)+r+1 , · · · , mL(q) 0L(q)−r bL(q) (6.19)
constitute a basis for NullG(r), where the vectors mp∈ CRensures that G(r)V(r)=
0RI3,(L(q−1)−r+1) and bp= rp(1 : p).
Proof. The proof is analogous to Proposition 6.1. Consider the columns {lr⊗a(3)r }
of L ⊙ A(3), then U = [u1, . . . ,uR] = h l1⊗ a(3)1 , . . . ,lR⊗ a(3)R i .
Since a(3)L
(q−1)+1 = a (3)
L(q−1)+2 = · · · = a (3)
L(q) we know that when r ∈ {L(q−1) +
1, . . . , L(q)} then the column vectors of
ur,ur+1, . . . ,uL(q)
can be written as lin-ear combinations of the column vectors of E(r)⊗ a(3)r . Hence, when r ∈ {L(q−1)+
1, . . . , L(q)} the problem of determining the rank of
h
U, −E(r)⊗ a(3)r
i
∈ CRI3×(2R−r+1) (6.20)
reduces to finding the rank of H(r)∈ CRI3×(2R−r+Lq)given by
H(r) =hu1, . . . ,ur−1,uL(q)+1, . . . ,uR,E (r)
⊗ a(3)r
i
. (6.21)
As observed in the proof of Proposition 6.1, the matrix (6.21) is a column permuted lower block-triangular with nonsingular matrices on its block-diagonal. Hence, the matrix (6.20) has rank 2R − r + Lq when r ∈ {L(q−1)+ 1, . . . , L(q)}. Consequently
any matrix U of which the columns form a basis for ColL ⊙ A(3)yields a matrix G(r)given by (6.18) of rank 2R − r + Lq when r ∈ {L(q−1)+ 1, . . . , L(q)}.
It is also clear from (6.21) that the column vectors of V(r) in (6.19) form a basis for NullG(r).
Due to Propositions 6.2 and 6.3 we can still compute the CPD of a tensor via the SGSD3-CP procedure up to its inherent ambiguities even if kA(3)= 1.
7. CPD with a partial symmetry. In many applications the CPD of a tensor has a partial symmetry. This is for instance the case when dealing with a tensor in which covariance matrices are stacked [3, 15, 8]. As a further motivation, let us consider the unsymmetric CPD in which one of the factor matrices, say A(3), has full column rank. For instance, let X ∈ CI1×I1×I3 of rank R be constructed from the
factor matrices A(1) ∈ CI1×R, A(2) ∈ CI2×R and full column rank A(3) ∈ CI3×R. It
has the following matrix representation CI1I2×I3 ∋ X
(1) =
A(1)⊙ A(2)A(3)T. (7.1) A mild uniqueness condition that applies in this case is stated next.
Theorem 7.1. Consider the tensor X ∈ CI1×I2×I3 with rank R and matrix
representation X(1)= A(1)⊙ A(2)A(3)T. Define C ∈ CI12I 2 2×R(R−1)/2 by c(i
1−1)(I1I22)+(j1−1)I22+(i2−1)I2+j2,(r2−2)(r2 −1)2 +r1 =
a(1)i1,r1 a(1)i1,r2 a(1)j1,r1 a (1) j1,r2 · a(2)i2,r1 a(2)i2,r2 a(2)j2,r1 a (2) j2,r2 , where 1 ≤ i1, j1≤ I1, 1 ≤ i2, j2 ≤ I2 and 1 ≤ r1< r2 ≤ R. If A(3) and C have full
column rank, then the CPD of X is unique [6, 13]. Generically1, this is true if [6] 2R(R − 1) ≤ I1(I1− 1)I2(I2− 1) and R ≤ I3. (7.2)
Condition (7.2) is also sufficient for generic uniqueness when A(1) = A(2)∗, with A(1) generic complex. The condition R ≤ I3 is meant to ensure that A(3) has full column
rank.
1A CPD property is called generic if it holds with probability one when the entries of the factor matrices are drawn from absolutely continuous probability density functions.
Assume that the conditions stated in Theorem 7.1 are satisfied. Let X(1) =
UΣVH denote the compact SVD of X(1), then there exists a nonsingular matrix
M ∈ CR×R such that A(1)
⊙ A(2) = UΣM and A(3) = M−1V∗. It was explained in [6] that the CPD problem (7.1) can be transformed into the partially symmetric CPD problem
CR2×R∋ Y(1)= (M ⊙ M) DT, (7.3) where D ∈ CR×R, see [6] for the construction of Y
(1). Hence, the problem of
comput-ing a CPD with a partial symmetry must be considered as a core problem in numerical multilinear algebra.
In subsection 7.1 we review how a partial symmetry was dealt with in the SGSD1-CP method in [20]. Next, in subsections 7.2 and 7.3 we explain how to incorporate a partial symmetry into the SGSD2-CP and SGSD3-CP methods. We limit the discus-sion to the case where A(1)= A(2), but the results are also valid when A(1)= A(2)∗. 7.1. SGSD1-CP with a partial symmetry. Consider a tensor X ∈ CI1×I1×I3
with rank R and matrix representation CI12×I3 ∋ X
(1)=
A(1)⊙ A(2)A(3)T =A(1)⊙ A(1)A(3)T. (7.4) Assume that the factor matrices A(1) and A(3) in (7.4) have full column rank. Then in [20] the computation of a CPD with the partial symmetry A(1)= A(2) was done via the SGSD1-CP method. This was done by ignoring the relation between the factor matrices A(1) and A(2) and by just treating the problem as an unsymmetric CPD problem.
7.2. SGSD2-CP with a partial symmetry. Assume that the factor matrices A(1) and A(3) in (7.4) have full column rank. Instead of just ignoring the partial symmetry, we present a partially symmetric version of the SGSD2-CP method.
Let A(1) = QR be the QR factorization of A(1) where Q ∈ CI1×I1 is a unitary
matrix and R ∈ CI1×R is a matrix of the form (2.4). Similarly, let A(2)= ZL be the
QL factorization of A(2) where Z ∈ CI1×I1 is a unitary matrix and L ∈ CI1×R is a
matrix of the form (2.5). Since A(1) = QR = ZL = A(2) we can express L in terms of R and vice-versa:
L = WR , R = WHL , (7.5)
in which W = V(1 : R, 1 : R) and V = QHZ. We first determine Q and Z by maximizing
h(Q, Z) = I3 X i3=1 triuR QHX(··i3)Z∗ 2 F.
Assume Q and Z have been found. Let the tensor R ∈ CR×R×I3 be constructed such
that Ri1i2i3 = R (··i3) i1i2 = triuR QHX(··i3)Z∗ i1i2 , then R(··i3) = RD i3 A(3)LT. Moreover, let R(i1··) ∈ CR×I3 satisfy the relation R(i1··)
i2i3 = Ri1i2i3, then R (i1··) =
LDi1 R
A(3)T. Stack the matrices {R(i1··)
} into the matrix
CR2×I3 ∋ R (1)= R(1··) .. . R(R··) = LD1 R .. . LDR R A(3)T = R ⊙ LA(3)T. (7.6)
Let us first use (7.5) to express (7.6) in terms of R only. We left-multiply (7.6) by (IR⊗ W) and work with (IR⊗ W) R(1)instead of R(1). This in turn means that (7.6)
can be expressed in terms of R as follows S(1) = (IR⊗ W) R(1) = R ⊙ R
A(3)T. We will now present an algorithm to compute the factor matrices from Col S(1)
. Let S(1)= UΣVH denote the compact SVD of S(1), then there exists a nonsingular
matrix M ∈ CR×R such that
UM = R ⊙ R ⇔ Umr= rr⊗ rr, r∈ {1, . . . , R} . (7.7)
Let ar= rr(1 : r), then (7.7) is equivalent to
Umr= X 1≤p≤r X 1≤q≤r e(R)p ⊗ e(R)q (ar⊗ ar) = E(r)(ar⊗ ar) , r ∈ {1, . . . , R}, (7.8) where E(r)∈ CR2×r2 is given by E(r)=he(R)1 ⊗ e (R) 1 , . . . ,e (R) 1 ⊗ e(R)r , . . . ,e(R)r ⊗ e (R) 1 , . . . ,e(R)r ⊗ e(R)r i . (7.9) Equation (7.8) can be written as
h
UΣ, −E(r)i amr
r⊗ ar
= 0R+r2, r∈ {1, . . . , R}. (7.10)
Proposition 7.3 tells us first that hUΣ, −E(r)i in (7.10) has an r-dimensional kernel. Second, it states that if we impose in (7.10) that an, 1 ≤ n ≤ R, are the
columns of a nonsingular upper triangular matrix, then we can determine ar in a
unique manner. In order to prove it we need the following lemma.
Lemma 7.2. Let R ∈ CR×R be a nonsingular upper triangular matrix, then the matrix S = R ⊙ R ∈ CR2×R
has the following property
S(r + (r − 1)R, 1 : r) = αre(r)Tr , αr∈ C, r ∈ {1, . . . , R} . (7.11)
Proof. The proof is straightforward, noting that
S (t + (t − 1)R, u) = (R)tu(R)tu (7.12) and (R)tu∈ C , u < t (R)tu6= 0 , u= t (R)tu= 0 , u > t (7.13)
By combining (7.12) and (7.13) we get (7.11).
Proposition 7.3. Let R = [r1, . . . , rR] ∈ CR×R be a nonsingular upper
trian-gular matrix and let E(r) ∈ CR2×r2
be of the form (7.9). If the column vectors of U∈ CR2×R
constitute a basis for Col (R ⊙ R), then the matrix G(r)=hU,−E(r)i∈ CR2×(R+r2)
has rank R + r(r − 1). Furthermore, let the column vectors of V(r) ∈ C(R+r2)×r constitute a basis for NullG(r), then W(r) = V(r) R+ (r − 1)r + 1 : R + r2,:∈
Cr×r has rank 1 and ColW(r)= Col (a
r), in which ar= rr(1 : r) ∈ Cr.
Proof. Let us first prove that G(r)in (7.14) has an r-dimensional kernel. Consider the columns {rn⊗ rn} of Col (R ⊙ R). On one hand, the vectors {rn⊗ rn}n≤r are
linear combinations of the columns of E(r). On the other hand, by Lemma 7.2, vector rn⊗rn, with n > r has a non-zero entry at position n+(n−1)R, while vectors rp⊗rp,
p≤ r, and the columns of E(r) are zero at that position. Furthermore, due to the property (7.11) we also know that the vectors {rp⊗ rp}p>r are linearly independent.
Consequently, the matrixhR ⊙ R, −E(r)ihas rank R+r(r −1). In conclusion, matrix (7.14) has rank R + r(r − 1).
Let us now prove that W(r) has rank 1 and ColW(r) = Col (ar), 1 ≤ r ≤
R. Note that E(r)(bp⊗ bp) = rp⊗ rp when bp =
ap 0r−p , and therefore xp = mp bp⊗ bp
∈ NullG(r), 1 ≤ p ≤ r, where the vectors mp ∈ CR are chosen
such that G(r)xp = 0R2, 1 ≤ p ≤ r. Due to Lemma 7.2 we know that the vectors
{bp⊗ bp}rp=1and hence the vectors {xp}rp=1 are linearly independent. Consequently,
{xp}rp=1 is a basis for Null
G(r). Let the column vectors of V(r) ∈ C(R+r2)×r constitute an alternative basis for NullG(r), then there exists a nonsingular matrix N ∈ Cr×r such that V(r)= [x1, . . . ,xr] N = m1 ,· · · , mr b1⊗ b1 br⊗ br N . (7.15) From (7.15) we observe that the column vectors of
W(r) = V(r)((r − 1) R + 1 : (r − 1) R + r, :) ∈ Cr×r
are all proportional to ar= br, i.e., the sth column vector of W(r) satisfies w (r)
s =
αsarfor some αs∈ C. Thus, we can conclude that any matrix U of which the columns
form a basis for Col (R ⊙ R) yields a rank-1 matrix W(r) of which the column space ColW(r)= Col (ar), r ∈ {1, . . . , R}.
The following Corollary 7.4 tells us that by projecting on the orthogonal comple-ment of ColE(r)we may reduce the number of equations in (7.10).
Corollary 7.4. Let R = [r1, . . . , rR] ∈ CR×Rbe a nonsingular upper triangular
matrix and define the diagonal projection matrix
CR2×R2 ∋ D(r) = IR2− E(r)E(r)T, (7.16)
where E(r) ∈ CR2×r2
is of the form (7.9). If the column vectors of U ∈ CR2×R
constitute a basis for Col (R ⊙ R), then the matrix
CR2×R∋ U(r)= D(r)U (7.17) has rank R − r, r ∈ {1, . . . , R}. Moreover, let the column vectors of V(r) ∈ CR×r
constitute a basis for NullU(r)and CR2×r
W(r)((r − 1) R + 1 : (r − 1) R + r, :) ∈ Cr×r has rank 1 and ColZ(r) = Col (r r),
r∈ {1, . . . , R}.
Proof. Let us first prove that U(r) = D(r)U has rank R − r, r ∈ {1, . . . , R}. Consider the columns {rn⊗ rn} of R ⊙ R. On one hand, the vectors in the set
{rn⊗ rn}n≤r can be written as a linear combination of the column vectors E(r). On
the other hand, by Lemma 7.2, vector rn⊗ rn, with n > r has a non-zero entry at
position n + (n − 1)R, that is not set to zero by D(r). Furthermore, due to property (7.11) we know that the vectors {rn⊗ rn}n>r are linearly independent. Hence, we
can conclude that any matrix U of which the columns form a basis for Col (R ⊙ R) yields a matrix U(r)= D(r)U of rank R − r, r ∈ {1, . . . , R}.
Let us now prove that W(r) has rank 1 and ColW(r)= Col (ar), 1 ≤ r ≤ R.
First note that U(r) = D(r)U = D(r)(R ⊙ R) F for some nonsingular matrix F ∈ CR×R. Since D(r) sets the non-zero entries of {rp⊗ rp}p≤r equal to zero while the vectors {rp⊗ rp}p>r are still linearly independent we conclude that {e(R)p }p≤r is a
basis for NullU(r). This means that
W(r)= (R ⊙ R) FE(r)N
= ([r1, . . . ,rr] ⊙ [r1, . . . ,rr]) M . (7.18)
for some nonsingular matrices M, N ∈ Cr×r, where E(r)=he(R)1 , . . . ,e (R) r
i
∈ CR×r. From (7.18) we observe that the column vectors of
Z(r)= W(r)((r − 1) R + 1 : (r − 1) R + r, :)
are all proportional to rr, that is z(r)s = αsrrfor some αs∈ C. Thus, we can conclude
that any matrix U of which the columns form a basis for Col (R ⊙ R) yields a matrix Z(r) with rank 1 and ColZ(r)= Col (rr), r ∈ {1, . . . , R}.
By projecting (7.10) on the orthogonal complement of E(r) we obtain
U(r)mr= 0R2, r∈ [1, . . . , R} (7.19)
with U(r)= D(r)UΣ and D(r)defined in (7.16). Due to Corollary 7.4, U(r)still has an r-dimensional kernel. Corollary 7.4 also tells us that we can find rrfrom Null
U(r). We first remove the all-zero row-vectors of U(r), obtaining eU(r) ∈ C(R2−r2)×R, so that (7.19) reduces to
e
U(r)mr= 0R2−r2, r∈ [1, . . . , R} . (7.20)
Note that set (7.10) of R2 equations in R + r2 unknowns has been reduced to a set
of R2− r2 equations in R unknowns. Next, we set the column vectors of V(r) equal
to the right singular vectors associated with the r smallest singular values of eU(r). Construct W(r) = U(r)V(r) and set ar equal to the singular vector associated with
the largest singular value of Z(r)= W(r)((r − 1) R + 1 : (r − 1) R + r, :). The vector rr is now given by rr =
aT
r,0TR−r
T
∈ CR. Once the vectors {r
r} have been found,
So far, we have been working with S(1) = R ⊙ R
A(3)T. As an alterna-tive, we may also use (7.5) to express (7.6) in terms of L only. We form T(1) =
WH⊗ IR
R(1) = L⊙ L
A(3)T and work as follows. Let T(1) = UΣVH denote
the compact SVD of T(1), then there exists a nonsingular matrix N ∈ CR×R such
that Unr= lr⊗ lr= X r≤p≤R X r≤q≤R e(R)p ⊗ e(R)q (br⊗ br) = E (r) (br⊗ br) , (7.21) where br= lr(r : R), and E (r) ∈ CR2×(R−r+1)2 is given by E(r)=he(R)r ⊗ e(R)r , . . . ,e(R)r ⊗ e (R) r+1, . . . ,e (R) R ⊗ e(R)r , . . . ,e (R) R ⊗ e (R) R i . (7.22) As in the upper triangular case, we first determine nrby solving a homogeneous linear
system. Let CR2×R
∋ U(r) = D(r)U in which CR2×R2
∋ D(r) = IR2 − E (r)
E(r)T. Let the column vectors of V(r) ∈ CR×(R−r+1) constitute a basis for NullU(r)
, W(r)= UV(r) and
Z(r)= W(r)(1 : (R − r + 1), :) ∈ C(R−r+1)×(R−r+1). (7.23) As in the upper triangular case, the vector br now follows from the best rank-1
approximation of Z(r) in (7.23).
We overall find a column of A(1) from the corresponding column of R (via S(1))
or L (via T(1)), depending on which has the least number of nonzero entries. This
reduces the computational complexity. Consequently, if a(1)r denotes the rth column
vector of A(1), then when r ≤ ⌊R+1
2 ⌋ we construct a (1)
r from R. On the other hand,
when r ≥ ⌊R+1
2 ⌋ + 1 we construct a (1)
r from L. Thus, we obtain
A(1)= Q · R :, 1 : R+ 1 2 ,Z · L :, R+ 1 2 + 1 : R .
The matrix A (3) now follows via relation A (3)T =A(1)⊙ A(1)†X(1).
7.3. SGSD3-CP with a partial symmetry. In this subsection we adapt the SGSD3-CP method to the case where a CPD has a partial symmetry. The method assumes that A(1)has full column rank and kA(3)≥ 1. This is in contrast with the previously discussed SGSD1-CP and SGSD2-CP methods for computing a CPD with a partial symmetry which require that both A(1) and A(3) have full column rank.
Consider the QR factorization A(1)= QR, where Q ∈ CI1×I1 is a unitary matrix
and R ∈ CI1×R is a matrix of the form (2.4). Consider also the QL factorization
A(1) = ZL, where Z ∈ CI1×I1 is a unitary matrix and L ∈ CI1×R is a matrix
of the form (2.5). Due to Proposition 2.2 or 6.2, the first step is to compute the unitary matrices Q and Z that make the matrices {X(··i3)} as upper triangular as
possible by maximizing f (Q, Z) = PI3 i3=1 triuR QHX(··i3)Z∗ 2 F. Assume that
Q and Z have been found and construct the tensor R ∈ CR×R×I3 according to
Ri1i2i3 = R (··i3) i1i2 = triuR QHX(··i3)Z∗ i1i2
partial symmetry we find A(3) from A(3)T = hVecdR(··1), . . . ,VecdR(··I3)i.
Next, we determine A(1). Let R(·i2·) ∈ CI3×R satisfy the relation R(·i2·)
i3i1 = Ri1i2i3,
then R(·i2·)= A(3)D i2 L
RT. Stack the matrices {R(·i2·)} into the matrix
CRI3×R∋ R (2) = R(·1·) .. . R(·R·) = A(3)D1 L .. . A(3)DR L RT =L ⊙ A(3)RT.
Let R(2)= UΣVH denote the compact SVD of R(2)and let M ∈ CR×R be a
nonsin-gular matrix such that
UM = L ⊙ A(3)⇔ Umr= lr⊗ a(3)r , r∈ {1, . . . , R} . (7.24)
Furthermore, let br = lr(r : R) ∈ CR−r+1, then due to Proposition 6.1 or 6.3 we
know that we can obtain L by solving the homogeneous linear systems h U, −E(r)⊗ a(3)r i mr br = 0I3R, r∈ {1, . . . , R} , (7.25)
where E(r)∈ CR×(R−r+1)is of the form (6.2).
Alternatively, let R(··i3) ∈ CR×R satisfy the relation R(··i3)
i1i2 = Ri1i2i3, then
R(··i3)= RD i3
A(3)LT. Stack the matrices {R(··i3)
} into the matrix
CRI3×R∋ R (3) = R(··1) .. . R(··I3) = RD1 A(3) .. . RDI3 A(3) L T =A(3)⊙ RLT.
Compute the compact SVD R(3)= SΣVHand let N ∈ CR×Rbe a nonsingular matrix
such that
SN = A(3)⊙ R ⇔ Snr= a(3)r ⊗ rr, r∈ {1, . . . , R} . (7.26)
Furthermore, let ar = rr(1 : r) ∈ Cr, then again due to Proposition 6.1 or 6.3 we
know that we can obtain R by solving the homogeneous linear systems h S, −a(3) r ⊗ F(r) i nr ar = 0I3R, r∈ {1, . . . , R} , (7.27) where F(r)=he(R)1 ,e (R) 2 , . . . ,e(R)r i ∈ CR×r. (7.28)
The influence of perturbation noise can be reduced by to simultaneously make use of relations (7.25) and (7.27). Recall from (7.5) that L = WR and R = WHL, in which W = V(1 : R, 1 : R) and V = QHZ. From (7.26) and relation L = WR we obtain (II3⊗ W) Snr= a (3) r ⊗ lr= a(3)r ⊗ E(r) br, (7.29)
where br = lr(r : R) ∈ CR−r+1 and E(r) ∈ CR×(R−r+1) is given by (6.2). By
com-bining (7.25) and (7.29) we get " U 0I3R,R −E (r) ⊗ a(3)r 0I3R,R (II3⊗ W) S −a (3) r ⊗ E(r) # mnrr br = 02I3R, r∈ {1, . . . , R} .(7.30)
Similarly, from (7.24) and relation R = WHL we obtain WH⊗ II3 Umr= rr⊗ a(3)r = F(r)⊗ a(3)r ar, (7.31)
where ar= rr(1 : r) ∈ Cr and F(r)∈ CR×r is given by (7.28). By combining (7.27)
and (7.31) we get " WH⊗ II3 U 0I3R,R −F (r) ⊗ a(3)r 0I3R,R S −a (3) r ⊗ F(r) # mnrr ar = 02I3R, r∈ {1, . . . , R}.(7.32)
We overall find a column of A(1)from the corresponding column of R (via (7.30)) or L (via (7.32)), depending on which has the least number of nonzero entries. Con-sequently, if a(1)r denotes the rth column vector of A(1), then when r ≤ ⌊R+12 ⌋ we
construct a(1)r from R. On the other hand, when r ≥ ⌊R+12 ⌋ + 1 we construct a(1)r
from L. Thus, we get A(1)= Q · R :, 1 : R+ 1 2 ,Z · L :, R+ 1 2 + 1 : R .
8. CPD with a partial Hermitian symmetry. In this section we consider the partial Hermitian symmetry A(2) = A(1)∗. This is for instance relevant for ap-plications involving covariance matrices of complex data [3, 8, 2]. Let us define the tensor Y ∈ CI1×I1×I3 by y
i1i2i3 = x ∗
i2i1i3 for all indices. Let us now stack X and Y
in Z ∈ CI1×I1×2I3 in such a way that Z (1) =
X(1),Y(1)
∈ CI2
1×2I3. In the case of
partial Hermitian symmetry, the first factor matrices of X and Y are both equal to A(1), and the second factor matrices of X and Y are both equal to A(1)∗. We have
Z(1)=
A(1)⊙ A(1)∗ hA(3),A(3)∗iT. (8.1) We may now study the uniqueness of the CPD of Z in order to obtain uniqueness results for X . We formulate as Theorem 8.3 a relaxation of Theorem 7.1 for the case with partial Hermitian symmetry. The proof makes use of Lemma 8.2 which in turn makes use of the following result.
Lemma 8.1. Given an analytic function f : Cn → C. If there exists a vector
x∈ Cn such that f (x) 6= 0, then the set { x | f (x) = 0 } is of Lebesgue measure zero.
Proof. For a proof, see for instance [12]. Lemma 8.2. Let P ∈ CI×R, then hPT, PHi
T
∈ C2I×R generically has rank
min (2I, R).
Proof. The proof is straightforward, simply notice that P P∗ = II √ −1 · II II −√−1 · II Re {P} Im {P} = T Re {P} Im {P} . (8.2)
If 2I ≥ R, then we let hRe {P}T,Im {P}TiT = hITR,0T2I−R,R
iT
and if 2I < R, then we let hRe {P}T,Im {P}TiT =hIT2I,0T2I,R−2I
iT
. Since T ∈ C2I×2I in (8.2) is
nonsingular we know that for this particular case the rank ofhPT,PHiT is equal to min (2I, R). By invoking Lemma 8.1 we can now conclude thathPT,PHiT generically has rank min (2I, R).
Theorem 8.3. Consider the tensor X ∈ CI1×I1×I3 with rank R and matrix
representation CI12×I3∋ X (1)=
A(1)⊙ A(1)∗A(3)T.Define C ∈ CI14×R(R−1)/2 by
c(i
1−1)I13+(j1−1)I12+(i2−1)I1+j2,(r2−2)(r2 −1)2 +r1 =
a(1)i1,r1 a(1)i1,r2 a(1)j1,r1 a(1)j1,r2 · a(1)∗i2,r1 a(1)∗i2,r2 a(1)∗j2,r1 a(1)∗j2,r2 , where 1 ≤ i1, j1, i2, j2≤ I1 and 1 ≤ r1< r2≤ R. If h A(3)T, A(3)Hi T
and C have full column rank, then the CPD of X is unique. Generically, this is true if
2R(R − 1) ≤ I12(I1− 1)2 and R≤ 2I3. (8.3)
Proof. By exploiting the partial symmetry X(··i3) = A(1)
Di3
A(3)A(1)H we can build the I3 additional matrix representations
Y(··i3)= X(··i3)H = A(1)D i3 A(3)∗A(1)H, i3∈ {1, . . . , I3} and obtain Z(1)∈ CI 2 1×2I3 given by Z(1) = h X(1),Vec Y(··1), . . . ,VecY(··I3)i=A(1) ⊙ A(1)∗ hA(3)T,A(3)Hi. The theorem now follows directly from Theorem 7.1 and Lemma 8.2.
In the case of complex Hermitian symmetry one may also resort to Theorem 7.1. Generic bound (7.2) also applies, with I1= I2. Comparing condition (7.2) with (8.3)
it is clear that Theorem 8.3 provides a more relaxed bound. In applications such as [3, 8, 2] this means that fewer covariance matrices are needed.
Under the conditions in Theorem 8.3 we may compute the CPD of X through the CPD of Z as in [6] and briefly described in section 7. Let Z(1) = UΣVH denote the
compact SVD of Z(1)in (8.1), then there exists a nonsingular matrix M ∈ CR×R such
that A(1)⊙ A(1)∗ = UΣM andhA(3)T,A(3)HiT = M−1V∗. The CPD of X follows
from the partially symmetric CPD (7.3) also in case where A(3) is rank deficient. Once A(1) and A(1)∗ have been obtained from M, we obtain A(3) via relation
A(3)T = A(1)⊙ A(1)∗ A(1)∗⊙ A(1) † X(1) Y∗(1) .
Moreover, if the CPD contains the partial Hermitian symmetry A(1)= A(2)∗, then due to (8.1) we know that as long as A(1) andhA(3)T,A(3)Hi
T
have full column rank the SGSD1-CP and SGSD2-CP methods can still be used, i.e., the rank constraint on A(3) has been relaxed to a rank constraint onhA(3)T,A(3)HiT.
9. Extension to N th-order tensors. In this section we briefly explain how to extend the SGSD-CP methods to N th-order tensors. Given a tensor X ∈ CI1×···×IN
with rank R constructed from the factor matrices A(n) ∈ CIn×R, n ∈ {1, . . . , N},
such that X = R X r=1 a(1)r ◦ · · · ◦ a(N )r . (9.1)
By appropriate reshaping the tensor (9.1) has the following matrix representation
X = (A ⊙ B) CT, (9.2) where A = A(1)⊙ · · · ⊙ A(P )∈ CI×R, I= P Y n=1 In B = A(P +1)⊙ · · · ⊙ A(Q)∈ CJ×R, J = Q Y n=P +1 In C = A(Q+1)⊙ · · · ⊙ A(N )∈ CK×R, K= N Y n=Q+1 In
in which 1 ≤ P < Q ≤ N. By ignoring the structural relations among the matrices in A, B and C, then (9.2) reduces to a third-order CPD problem. Assume that A and B have full column rank and k (C) ≥ 1, then in order to compute A and B from (9.2) the discussed SGSD-CP methods can be applied. Let A = QR be the QR factorization of A, where Q ∈ CI×I is a unitary matrix and R ∈ CI×R is a matrix
of the form (2.4) and A = ZL be the QL factorization of A where Z ∈ CJ×J is a
unitary matrix and L ∈ CJ×R is a matrix of the form (2.5), then cost function (2.14)
is equal to g(Q, Z) = K X k=1 triuR QHX(··k)Z∗ 2 F , (9.3) where X(··k)= ADk(C) BT, k∈ {1, . . . , K} .
Assume that A, B and C have been found via one of the proposed SGSD methods. The factor matrices now follow from best rank-1 tensor approximation problems. For instance, from A we obtain
min a(1)r ,..., a(P )r ar− a(1)r ⊗ · · · ⊗ a(P )r 2 F, r∈ {1, . . . , R} . (9.4)
Notice that when P = 2, then (9.4) reduces to best rank-1 matrix approximation problems which can be solved by standard numerical linear algebra methods. By sim-ilar procedure we obtain {A(n)}Qn=P +1 and {A
(n)
}N
n=Q+1 from rank-1 approximation
If A = B∗, then we work with the tensor Z = (A ⊙ A∗) C C∗ T ∈ CI2×2K (9.5)
obtained via the procedure described in section 8.
To summarize, we can compute the CPD of an N th-order tensor by treating it as a third-order tensor. Note that a third-order tensor is usually highly overdetermined and therefore we do not expect to lose any significant accuracy by combing the modes of an N th-order tensor.
10. Conclusion. Based on a link between the SGSD and CPD it is possible to address some CPD problems as SGSD problems. We reviewed the original SGSD-CP method which was called SGSD1-CP in this paper. It was originally designed to only jointly diagonalize a set of square matrices. Hence, we first generalized the extended QZ method to the non-square case.
The SGSD1-CP method does not fully take the relations among the involved triangular matrices into account. This motivated us to propose the SGSD2-CP which attempts to take the relations among the involved triangular matrices into account. In the derivation of SGSD2-CP we showed that it is possible to deduce the involved triangular matrices from the column space of a matrix representation of the tensor. Moreover, we provided a numerical method to compute the triangular factor matrices from the knowledge of the column space of the matrix representation of the tensor of interest.
Another drawback of the original SGSD1-CP method is that it requires that all three of the involved factor matrices have full column rank. This motivated us to develop the SGSD3-CP method which only requires that two of the involved factor matrix to have full column rank. This is in fact the necessary condition for the SGSD approach to the CPD problem. Again, we provided a numerical method to compute the involved triangular factor matrices.
In many practical problems the CPD has a partial (Hermitian) symmetry. We explained how to incorporate such a constraint into the SGSD2-CP and SGSD3-CP methods. For the case when the CPD having a partial Hermitian symmetry we also proved a new uniqueness result. Finally, we explained how to extend the SGSD-CP methods for computing a N th-order SGSD-CPD with possibly also partial (Hermitian) symmetries.
REFERENCES
[1] K. Abed-Meraim and Y. Hua, A Least-Squares Approach to Joint Schur Decomposition, Proc. ICASSP, May 12-15, Seattle, USA, 1998.
[2] K. Abed-Meraim and Y. Xiang and J. H. Manton and Y. Hua, Blind Source Separation Using Second-Order Cyclostationarity Statistics, IEEE Trans. Signal Processing, 49(2001), pp. 694–701.
[3] A. Belouchrani and K. Abed-Meraim and J.-F. Cardoso and E. Moulines, A Blind Source Separation Technique Using Second-Order Statistics, IEEE Trans. Signal Process-ing, 45(1997), pp. 434–444.
[4] P. Comon and M. Sørensen and E. Tsigaridas, Decomposing Tensors with Structured Matrix Factors Reduces to Rank-1 Approximations, Proc. ICASSP, March 14-19, Dallas, USA, 2010
[5] L. De Lathauwer and B. De Moor and J. Vandewalle, Computation of the Canonical Decomposition by means of a Simultaneous Generalized Schur Decomposition, SIAM J. Matrix Anal. Appl., 26(2004), pp. 295–327.
