GLOBALLY CONVERGENT JACOBI-TYPE ALGORITHMS
FOR SIMULTANEOUS ORTHOGONAL
SYMMETRIC TENSOR DIAGONALIZATION∗
JIANZE LI† , KONSTANTIN USEVICH‡ , AND PIERRE COMON‡
Abstract. In this paper, we consider a family of Jacobi-type algorithms for simultaneous orthogonal diagonalization problem of symmetric tensors. For the Jacobi-based algorithm of [SIAM J.
Matrix Anal. Appl., 2(34):651–672, 2013], we prove its global convergence for simultaneous orthogonal diagonalization of symmetric matrices and 3rd-order tensors. We also propose a new Jacobi-based
algorithm in the general setting and prove its global convergence for sufficiently smooth functions.

arXiv:1702.03750v2 [math.NA] 27 Jul 2017

Key words. orthogonal tensor diagonalization, Jacobi rotation, global convergence, Lojasiewicz
gradient inequality, proximal algorithm
AMS subject classifications. 15A69, 49M30, 65F99, 90C30

1. Introduction. Higher-order tensor decompositions have attracted a lot of
attention in the last two decades because of the applications in various disciplines,
including signal processing, numerical linear algebra and data analysis [5, 9, 20].
The most popular decompositions include the Canonical Polyadic and Tucker decompositions, where additional constraints are often imposed, such as symmetry [10],
nonnegativity [22] or orthogonality [19].
An important class of tensor approximation problems is approximation with orthogonality constraints. In particular, the orthogonal symmetric tensor diagonalization problem for 3rd and 4th-order cumulant tensors is in the core of Independent
Component Analysis [6, 7, 8], and is a popular way to solve blind source separation
problems in signal processing [11]. In the same context, simultaneous orthogonal matrix diagonalization [4] is widely used; simultaneous orthogonal tensor diagonalization
for slices of 4th-order cumulants is also used to solve source separation problems [16].
def

Main notations. Let Rn1 ×···×nd = Rn1 ⊗ · · · ⊗ Rnd denote the space of dth-order
tensors. In the paper, the tensors are typeset with a bold calligraphic font (e.g, T ,
W), and matrices are in bold (e.g, Q, U ); the elements of tensors or matrices are
typeset as Tijk (or sometimes Ti,j,k ) or Qij respectively.
For a tensor T ∈ Rn1 ×···×nd and a matrix M ∈ Rm×nk , their k-mode product
[20, subsection 2.5] is the tensor (T •k M ) ∈ Rn1 ×···×nk−1 ×m×nk+1 ×···×nd defined as
def

(T •k M )i1 ,...,id =

nk
X

Ti1 ,...,ik−1 ,j,ik+1 ...,id Mi1 ,j

j=1

A tensor T ∈ Rn×···×n is called symmetric if its entries do not change under any
permutation of indices; its diagonal is, by definition, the vector
def 
diag{T } = T1...1
∗ Submitted

···

Tn...n

T

.

to the editors DATE.
Funding: This work was supported by the ERC project “DECODA” no.320594, in the frame
of the European program FP7/2007-2013. The first author was partially supported by the National
Natural Science Foundation of China (No.11601371).
† School of Mathematics, Tianjin University, Tianjin 300072, China (lijianze@tju.edu.cn).
‡ GIPSA-Lab, CNRS and Univ. Grenoble Alpes, France (firstname.lastname@gipsa-lab.fr).
1

2

J. LI, K. USEVICH, AND P. COMON

We will always denote by k · k the Frobenius norm of a tensor or a matrix, or the
Euclidean norm of a vector. Finally, let On ⊂ Rn×n denote the orthogonal group,
that is, the set of orthogonal matrices. Let SOn ⊂ On denote the special orthogonal
group, the set of orthogonal matrices with determinant 1.
Problem statement. In this paper, we consider the following simultaneous orthogonal symmetric tensor diagonalization problem. Let {A(`) : 1 ≤ ` ≤ m} ⊂ Rn×···×n
be a set of symmetric tensors. We wish to maximize
(1)

Q∗ = arg max

Q∈SOn

m
X

k diag{W (`) }k2 ,

`=1

where W (`) = A(`) •1 QT · · · •d QT for 1 ≤ ` ≤ m. Problem (1) has the following
well-known problems as special cases:
• orthogonal tensor diagonalization problem if m = 1 and d > 2;
• simultaneous orthogonal matrix diagonalization problem if m > 1 and d = 2.
Several algorithms have been proposed in the literature to solve special cases of
problem (1). The first were the Jacobi CoM (Contrast Maximization) algorithm for
orthogonal diagonalization of 3rd and 4th-order symmetric tensors [6, 7, 8] and the
JADE (Joint Approximate Diagonalization of Eigenmatrices) algorithm for simultaneous orthogonal matrix diagonalization [4]. An algorithm for simultaneous orthogonal
3rd-order tensor diagonalization was proposed in [14]. These Jacobi-type algorithms
have been very widely used in applications [11], and have the advantage that Jacobi
rotations can be computed by rooting low-order polynomials. Nevertheless, up to our
knowledge, the convergence of these methods was not proved, although it was often
observed in practice [8, 13].
Contribution. In this paper, we consider several Jacobi-type algorithms to solve
problem (1), and study their global convergence properties. By global convergence, we
mean that, for any starting point, the whole sequence of iterations produced by the
algorithm always converges to a limit point1 . Note that the global convergence does
not guarantee convergence to a global maximum, since the cost function is multimodal.
First, we consider the Jacobi-based algorithm proposed in [17] for best low multilinear rank approximation of 3rd-order symmetric tensors. The algorithm uses a
gradient-based order of Jacobi rotations, hence we call this algorithm Jacobi-G in our
paper. For the Jacobi-G algorithm the convergence of subsequences of iterations to
stationary points was established in [17]. We prove that, for problem (1), Jacobi-G
algorithm converges globally and the limit point is always a stationary point in the
cases d = 2, 3. The proof is based on a variant Lojasiewicz theorem developed in [25]
for analytic submanifolds of Euclidean space.
Second, we propose a new Jacobi-based algorithm inspired by proximal methods
[24], which is called the Jacobi-PC algorithm in this paper. We show that JacobiPC algorithm always converges globally to a stationary point in the general setting.
In particular, for cases d = 3, 4 of problem (1), Jacobi-PC algorithm allows for a
simple algebraic solution to find the optimal Jacobi rotation, as in the Jacobi CoM
algorithm [6, 7, 8]. In addition, this algorithm does not need the order of Jacobi
rotations introduced in Jacobi-G algorithm. Finally, the global convergence of JacobiPC algorithm is proved for sufficiently smooth functions. Therefore, these results may
be applied to other tensor approximation problems (e.g., nonsymmetric orthogonal
diagonalization [23],or Tucker approximation of symmetric [17] or antisymmetric [3]
1 This

was also called single-point convergence in [26].

GLOBALLY CONVERGENT JACOBI-TYPE ALGORITHMS

3

tensors).
Organization. The paper is organized as follows. In sections 2 to 4, we present
some more or less known results or easy extensions. In sections 5 to 7, we show
our main results of this paper. In section 2, we introduce the abstract optimization
problem on SOn and the general Jacobi algorithm to solve this abstract problem. In
section 3, we recall the Jacobi-G algorithm proposed in [17] and its local convergence
properties in the general setting. In section 4, we list some properties that are specific to orthogonal tensor diagonalization problem. In section 5, we prove the global
convergence of Jacobi-G algorithms for simultaneous orthogonal diagonalization of
matrices and 3rd-rder tensors. In section 6, we propose the Jacobi-PC algorithm
and prove its global convergence. We derive formulas for optimal Jacobi rotations
in the case of 3rd and 4th-order tensors. In section 7, we provide some numerical
experiments.
2. Problem statement and Jacobi rotations.
2.1. Abstract optimization problem. Let f : SOn → R be a function. The
abstract optimization problem is to find Q∗ ∈ SOn such that
Q∗ = arg max f(Q).

(2)

Q∈SOn

Example 1. Problem (1) in section 1 is a special case of problem (2). Because
kA(`) k = kW (`) k for any 1 ≤ ` ≤ m, problem (1) is equivalent to find
Q∗ = arg min

Q∈SOn

where offdiag{W

(`)

m
X

k offdiag{W (`) }k2 ,

`=1

} is the vector of elements in W (`) except diagonal elements.

Example 2. In [17], the best low multilinear rank approximation problem for 3rd
order symmetric tensors was formulated as a special case of problem (2). In fact,
based on [15, Theorem 4.1], the cost function has the form
f(Q) = kA •1 M QT •2 M QT •3 M QT k2


Ir 0
where r < n and M =
∈ Rn×n .
0 0

(3)

Example 3. Under some assumptions, it was shown that the orthogonal tensor
diagonalization problem could be solved in the sense of maximization of the trace of a
tensor [12], and is also a special case of problem (2).
2.2. Givens rotations and the general Jacobi algorithm. Let θ ∈ R be
an angle and (i, j) be a pair of indices with 1 ≤ i < j ≤ n. We denote the Givens
rotation matrix by


1


..


.
0




cos θ
− sin θ




(i,j,θ)
..
G
=
,
.




sin θ
cos θ




.
.


.
0
1

4

J. LI, K. USEVICH, AND P. COMON

i.e., the matrix defined by

(G(i,j,θ) )k,l



1,





cos θ,
= sin θ,



− sin θ,



0,

k = l, k 6∈ {i, j},
k = l, k ∈ {i, j},
(k, l) = (j, i),
(k, l) = (i, j),
otherwise

for 1 ≤ k, l ≤ n. We summarize the general Jacobi algorithm in Algorithm 1.
Algorithm 1 General Jacobi algorithm
Input: A function f : SOn → R, a starting point Q0 .
Output: Sequence of iterations {Qk }k≥1 .
• For k = 1, 2, . . . [until a stopping criterion is satisfied] do
•
Choose the pair (ik , jk ) according to a certain pair selection rule.
•
Compute the angle θk∗ that maximizes the function hk (θ) defined as
def

hk (θ) = f (Qk−1 G(ik ,jk ,θ) ).

(4)
def

∗

•
Set U k = G(ik ,jk ,θk ) , and update Qk = Qk−1 U k .
• End for
Note that in Algorithm 1 (and all the following algorithms, unless mentioned
explicitly) we do not specify a stopping criterion, since our goal is to analyse the whole
sequence of iterations {Qk }k≥1 produced by the algorithm. The stopping criterion
applies only to software implementation of algorithms.
Algorithm 1 performs the sequence of iterations Qk by multiplicative updates:
(5)

Qk = U 1 · · · U k−1 U k ,

where each U k is an elementary rotation. The advantage of this Jacobi-type algorithm
is that each update is a one-dimensional optimization problem, which can be solved
efficiently.
Remark 2.1. Note that the choice of pair (ik , jk ) in every iteration is not specified
in Algorithm 1. One of the most natural rules is in cyclic fashion as follows.
(1, 2) → (1, 3) → · · · → (1, n) →
(6)

(2, 3) → · · · → (2, n) →
· · · → (n − 1, n) →
(1, 2) → (1, 3) → · · · .

We call the Jacobi algorithm with this cyclic-by-row rule the Jacobi-C algorithm.
In Jacobi-C algorithm, the choice of pairs is periodic with the period n(n − 1)/2.
Each set of n(n − 1)/2 iterations is called a sweep. The pair selection rule (6) was
used in the Jacobi CoM algorithm [6, 7, 8] and JADE algorithm [4].
Remark 2.2. If several equivalent maximizers are present in (4), also in all the
following algorithms, we choose the one with the angle of smaller magnitude.

GLOBALLY CONVERGENT JACOBI-TYPE ALGORITHMS

5

3. Jacobi-G algorithm. In this section, we recall the Jacobi-based algorithm in
[17] and its properties. For simplicity, this algorithm is called the Jacobi-G (gradientbased Jacobi-type) algorithm in our paper.
3.1. Technical definitions. Define the matrix

0

..

.


0

d (i,j,θ)

..
=
δi,j =
G
.

dθ
θ=0

1



0
and introduce the notation


0
−1
0
..

.













0

def

di,j (Q) = Qδi,j
for Q ∈ On . Next, for a differentiable function f : On → R, we define the projected
gradient [17, Lemmma 5.1] as
def

Proj∇ f (Q) = QΛ(Q),

(7)
where

QT ∇f(Q) − (∇f(Q))T Q
2
and ∇f(Q) is the Euclidean gradient of f as a function of the matrix argument. The
projected gradient Proj∇ f (Q) is exactly the Riemannian gradient if On is viewed as
an embedded submanifold of Rn×n [2].
def

(8)

Λ(Q) =

3.2. Jacobi-G algorithm. In [17], a modification of Algorithm 1 was proposed,
that choose a pair (i, j) at each iteration that satisfies
|hProj∇ f (Q), di,j (Q)i| ≥ εk Proj∇ f (Q)k,

(9)

where ε is a small positive constant.
Algorithm 2 Jacobi-G algorithm
Input: A function f : SOn → R, a small positive 0 < ε ≤ n2 , a starting point Q0 .
Output: Sequence of iterations {Qk }k≥1 .
• For k = 1, 2, . . . [until a stopping criterion is satisfied] do
•
Choose the pair (i, j) = (ik , jk ) satisfying inequality (9) at Q = Qk−1 .
•
Compute the angle θk∗ that maximizes the function hk (θ) defined in (4).
def

∗

•
Set U k = G(ik ,jk ,θk ) , and update Qk = Qk−1 U k .
• End for
Note that, as shown in [17, Proof of Lemma 5.3] the left hand side of (9) can be
also written in terms of the function hk :
(10)

|h0k (0)| = |hProj∇ f (Qk−1 ), dik ,jk (Qk−1 )i|.

The following lemma (an easy generalisation of [17, Lemma 5.2]) shows that it is
always possible to choose such a pair (ik , jk ).

6

J. LI, K. USEVICH, AND P. COMON

Lemma 3.1. For any differentiable function f : SOn → R, Q ∈ SOn and 0 < ε ≤
2/n, it is always possible to find (i, j), with i < j, such that (9) holds.
Proof. First of all, thanks to the representation (7), rotation invariance of the
Euclidean norm and Λ(Q) being skew-symmetric, we have that the condition (9) is
equivalent to
2|Λi,j (Q)| ≥ εkΛ(Q)k.

(11)

Now we choose the pair (i, j) that maximizes the left-hand side of (11)
(i, j) = arg max |Λk,l (Q)|,
1≤k,l≤n

Since Λ(Q) is skew-symmetric, we can choose such a pair with i < j. Finally, by the
well known inequality between matrix norms,
2|Λi,j (Q)| = 2 max |Λk,l (Q)| ≥
1≤k,l≤n

2
kΛ(Q)k ≥ εkΛ(Q)k.
n

Remark 3.2. Lemma 3.1 was proved in [17, Lemma 5.2] only for the cost function (3), although it is valid for any differentiable function. Hence, the convergence
results of [17] are valid for a wide class of functions (see the following theorem).
Theorem 3.3 ([17, Theorem 5.4] combined with Lemma 3.1). Let f be a C ∞
function. Then every accumulation point2 Q∗ of the sequence {Qk }k≥1 produced by
Algorithm 2 is a stationary point of the function f ( i.e, Proj∇ f (Q∗ ) = 0).
3.3. Variants of the Jacobi-G algorithm. Note that the description of Algorithm 2 does not say precisely how to select the pairs (ik , jk ) that satisfy (9). The
first option is suggested by the proof in Lemma 3.1: set ε = n2 and select the maximal
element in (8). We summarize this option in Algorithm 3.
Algorithm 3 Jacobi-G-max algorithm
Input: A function f : SOn → R, a starting point Q0 .
Output: Sequence of iterations {Qk }k≥1 .
• For k = 1, 2, . . . until a stopping criterion is satisfied do
•
Choose the pair (ik , jk ) = (i, j) that maximizes |Λi,j (Qk−1 )| in (8).
•
Compute the angle θk∗ that maximizes the function hk (θ) defined in (4).
def

∗

•
Set U k = G(ik ,jk ,θk ) , and update Qk = Qk−1 U k .
• End for
However, choosing the maximal element in (8) requires a search over all the
elements, which may take additional time. In [17], it was suggested to take ε  1. If
we choose (ik , jk ) in the cyclic order and ε is very small, then it is natural to expect
that the inequality (9) will be often satisfied, and thus the behavior of the algorithm
will be very close to the behavior of the Jacobi-C algorithm.
To make this idea more rigorous, we construct a modification of Jacobi-C algorithm that skips the rotations if the magnitude of the directional derivative is below
a given threshold. The modification is described in Algorithm 4.
2 i.e.,

the limit of every convergent subsequence.

GLOBALLY CONVERGENT JACOBI-TYPE ALGORITHMS

7

Algorithm 4 Jacobi-C-threshold algorithm
Input: A function f : SOn → R, a threshold δ > 0, a starting point Q0 .
b such that k Proj∇ f (Q)k
b ≤ δ.
Output: A point Q,
• Set k = 1.
• Do (cyclic sweep)
•
For i = 1, . . . , n, j = i + 1, . . . , n
•
If |Λi,j (Qk−1 )| > nδ
•
Set (ik , jk ) = (i, j).
•
Compute the angle θk∗ that maximizes hk (θ) defined in (4).
∗

def

•
Set U k = G(ik ,jk ,θk ) , and update Qk = Qk−1 U k .
•
Set k ← k + 1.
•
End If
•
End For
• While there was progress in the last sweep.
b = Qk−1 .
• Return Q

Note that the output of Algorithm 4 is different from other algorithms in terms
of the output, because the rule of skipping the rotations also yields a well-defined
stopping criterion. Also, note that, compared to Jacobi-G, the inequality (9) is not
needed. But, Algorithm 4 can be viewed as a special case of Algorithm 2, as shown
by the following remark.
Remark 3.4. There exists ε = ε(δ) such that the iterates Qk produced by Algorithm 4 are the first elements of the sequence produced by a Jacobi-G algorithm
(Algorithm 2). Indeed, consider c = maxQ∈SOn k Proj∇ f (Q)k, which is finite if f is
δ
is the required value of ε.
smooth (since SOn is compact). Hence ε = nc
Also note that, Algorithm 4 always terminates if the Jacobi-G algorithm converges
to a stationary point.

4. Derivatives in the orthogonal tensor diagonalization problems .

4.1. Projected gradient: the case of a single tensor (m = 1) . In this
subsection, we derive a concrete form of projected gradient (7), which will be used
in section 5. Let A ∈ Rn×···×n be a dth-order symmetric tensor and Q ∈ On be an
orthogonal matrix. Let W = A •1 QT · · · •d QT and

(12)

f(Q) =

n
X
j=1

2
Wjj...j
=

n
X

(

X

j=1 p1 ,p2 ,...,pd

Ap1 ,p2 ,...,pd Qp1 ,j Qp2 ,j . . . Qpd ,j )2 .

8

J. LI, K. USEVICH, AND P. COMON

First, we calculate the Euclidean gradient of f at Q. Let us fix i and j. Then


h
X
d
∂f
d−2
= 2Wjj...j dQd−1
A
+
(d
−
1)
Qi,j
(
Qk1 ,j Ai,...,i,k1 ) + . . .
i,i,...,i
i,j
∂Qi,j
d−1
k1 6=i
 
X
d
+q
Qq−1
Qk1 ,j · · · Qkd−q ,j Ai,...,i,k1 ,...,kq ) + . . .
i,j (
q
k1 ,...,kd−q 6=i
i
X
+d
Qk1 ,j · · · Qkd−1 ,j Ai,k1 ,...,kd−1 )
k1 ,...,kd−1 6=i

= 2dWjj...j

n
X

Qk1 ,j · · · Qkd−1 ,j Ai,k1 ,...,kd−1 = 2dWjj...j Vij...j ,

k1 ,...,kd−1 =1

where V = A •2 QT · · · •d QT . Noting that V = W •1 Q, we get that



V11...1 V12...2 · · · V1n...n W1...1 0 · · ·
0
 V21...1 V22...2 · · · V2n...n  
.. 
..
..


.
.
0
. 

 ..
.


..  
∇f(Q) = 2d  .
·
·
·
·
·
·

.
.
.
 .

.
.
.
.
.
0 


0
···
0 Wn···n
Vn1...1 Vn2...2 · · · Vnn...n



W1...1 0 · · ·
0
W11...1 W12...2 . . . W1n...n
 W21...1 W22...2 . . . W2n...n  
.. 
..
..

 0
.
.
. 
.



...
...
...  .
= 2dQ  . . .

.
.
..
..

  ..
0 
Wn1...1 Wn2...2 . . . Wnn...n
0
···
0 Wn···n
After projecting ∇f(Q) onto the tangent space at Q to the manifold On , we get
Proj∇ f (Q) = QΛ(Q),
where Λ(Q) is the matrix with


0,
(13)
Λk,l (Q) = d(Wkl...l Wl...l -Wk...k Wk...kl ),


−Λk,l (W),

k = l;
k < l;
k > l.

for any 1 ≤ k, l ≤ n. This is a special case of (8) for function (12).
Remark 4.1. Let A be a 3rd order symmetric tensor. Then Proj∇ f (Q) =


0
W122 W222 −W111 W112 W133 W333 −W111 W113
W
W
−W
W
0
W
W
−W
W
3Q
.
111
112
122
222
233
333
222
223
W111 W113 −W133 W333 W222 W223 −W233 W333

0

4.2. The cost function at each iteration. Now consider a single iteration in
Algorithm 1 for the cost function (12) . For simplicity3 , we consider only the case
def

(ik , jk ) = (1, 2). We also denote G(θ) = G(1,2,θ) without loss of generality. Then
Qk = Qk−1 G(θk∗ ).
3 other

cases follow by substitution of indices

9

GLOBALLY CONVERGENT JACOBI-TYPE ALGORITHMS

In this subsection, we take Q = Qk−1 , so that
W = A •1 QT •2 QT · · · •d QT ,
is the rotated tensor before the kth iteration. We also use notation T for the candidate
tensors after k-th iteration, i.e.
T = T (θ) = W •1 GT (θ) •2 GT (θ) · · · •d GT (θ),
so that hk (θ) = k diag{T (θ)}k2 . Note that we omit k in the notation for W and T ,
but this will not lead to confusion within one iteration.
Although hk (θ) is defined on the whole real line, it is periodic. Apart from the
obvious period 2π, it has a smaller period π/2. Indeed,




0 −1 0 · · · 0
− sin(θ) − cos(θ) 0 · · · 0
1 0 0 · · · 0
 cos(θ) − sin(θ) 0 · · · 0






 0
0
1 · · · 0 = G(θ) 
G (θ + π/2) = 
0 0 1 · · · 0 .
 ..

..
.. . .
.. 
..
..
.. . .
.. 
.

. .
. .
.
.
.
.
.
0

0 ···

0

1

0

0

0 ···

1

Hence, the tensor T (θ + π/2) differs from T (θ) by permutations of the indices 1
and 2 and change of signs, which does not change the cost function. Hence, due to
Remark 2.2, we are, in fact, maximizing hk (θ) on the interval [−π/4, π/4] . In fact,
in all the algorithms we choose θk∗ ∈ [−π/4, π/4] with hk (θk∗ ) = max hk (θ).
θ∈R

It is often convenient to rewrite the function hk (θ) in (17) in a polynomial form.
Consider the change of variables θ = arctan(x). Then minimization of hk (θ) on
[−π/4, π/4] is equivalent to minimization of
def

τk (x) = hk (arctan(x))
on [−1, 1]. After the change of variables, since x = tan θ, we have




1
P (θ)
1 −x
(14)
G(θ) =
, where P (θ) = √
.
I n−2
x2 + 1 x 1
Hence, the cost function is a rational function in x
(15)

τk (x) =

ρ(x)
,
(1 + x2 )d

where ρ(x) is a polynomial of degree 2d (note that we dropped the index k for simplicity).
In [8, eqn. (22)-(23)] it was shown that finding critical points of τk (x) can be
reduced to finding roots of a quadratic polynomial if d = 3 or a quartic polynomial if
d = 4. We do not provide these expressions in this section, but their generalizations
can be found in subsection 6.2.
4.3. Derivatives for matrices and third-order tensors. We recall the expressions for derivatives of the cost function (12) in the case d ∈ {2, 3} from [8], but
give proofs for completeness.
Lemma 4.2 ([8, eqn. (22)-(23)]). In the notations of subsection 4.2, the derivatives of hk have the following form.

10

J. LI, K. USEVICH, AND P. COMON

• For d = 3:
0

hk (θ) = 6(T111 T112 − T122 T222 ),
00

2
2
2
2
hk (θ) = −6(T111
+ T222
− 3T112
− 3T122
− 2T111 T122 − 2T112 T222 ).

• For d = 2:
0

hk (θ) = 4(T11 T12 − T12 T22 ),
00

2
2
2
hk (θ) = −4(T11
+ T22
− 2T11 T22 − 4T12
).
0

00

The expressions for hk (0) and hk (0) can be found by substituting T to W in the above
expressions.
Proof. From (13) we have
0

0

0

hk (θ) = hProj∇ f(Qk−1 G(θ)), Qk−1 G (θ)i = hQk−1 G(θ)Λ(Qk−1 G(θ)), Qk−1 G (θ)i
T

0

= hΛ(Qk−1 G(θ)), G(θ) G (θ)i = −2Λ1,2 (Qk−1 G(θ))
(
6(T111 T112 − T122 T222 ), d = 3,
=
4(T11 T12 − T12 T22 ),
d = 2,
where Λ is defined in (13).
After straightforward differentiation, we have that for d = 3
dT111
dT112
dT122
dT222
= 3T112 ,
= 2T122 − T111 ,
= T222 − 2T111 ,
= −3T122 ,
dθ
dθ
dθ
dθ
and for d = 2
dT11
= 2T12 ,
dθ

dT12
= (T22 − T11 ),
dθ

dT22
= −2T12 .
dθ

00

The equation for hk (θ) follows by substitution.
4.4. The general case (m > 1). In the general case, the cost function in
problem (1) is
(16)

f (Q) =

m
X

f (`) (Q),

`=1

where f (`) (Q) = k diag{W (`) }k2 for any 1 ≤ ` ≤ m.
By linearity, the cost function hk (θ) in (4) in every iteration can be conveniently
written as
(17)

hk (θ) =

m
X

(`)

hk (θ),

`=1

where
(`)

def

hk (θ) = f (`) (Qk−1 G(ik ,jk ,θ) ) = k diag{T (`) }(θ)k2 ,
and T (`) (θ) is defined in the same way as T in subsection 4.2. Therefore, the deriva(`)
tives of hk can be obtained in the same way as in subsection 4.3

GLOBALLY CONVERGENT JACOBI-TYPE ALGORITHMS

11

Finally, as in subsection 4.2, we can use a change of variables
(18)

τk (x) = hk (arctan(x)),

which leads to
τk (x) =

m
X

(`)

τk (x)

`=1

where, for 1 ≤ ` ≤ m,
(`)

(`)

τk (x) = hk (arctan(x)) =

ρ(`) (x)
.
(1 + x2 )d

5. Global convergence of Jacobi-G algorithm for symmetric low-order
tensors.
5.1. Lojasiewicz gradient inequality. In this subsection, we recall some important results about the convergence of iterative algorithms. The discrete-time analogue of classical Lojasiewicz’s theorem was proposed in [1], and make it possible to
prove the convergence of many algorithms [26, 27]. In [25], to prove the convergence
of projected line-search methods on the real-algebraic variety of real m × n matrices
of rank at most k, the optimization problem
min f(x)

x∈M

on a closed subset M ⊆ Rn was considered. Suppose that the tangent cone Tx M at
x ∈ M is a linear space. Let Proj∇ f (x) be the projection of the Euclidean gradient
∇f(x) on the tangent space at x. We first introduce the definition of an analytic
submanifold in Rn .
Definition 5.1 ([21, Def. 2.7.1]). A set M ⊆ Rn is called an m-dimensional real
analytic submanifold if, for each p ∈ M, there exists an open subset U ⊆ Rm and a
real analytic function f : U → Rn which maps open subsets of U onto relatively open
subsets of M and which is such that
p ∈ f (U) and rank J f (u) = m, ∀u ∈ U,
where J f (u) is the Jacobian matrix of f at u.
The following results were proved in [25].
Lemma 5.2. Let M ⊆ Rn be an analytic submanifold. Then any point x ∈ M
satisfies a Lojasiewicz inequality for Proj∇ f (x), that is, there exist δ > 0, σ > 0 and
ζ ∈ (0, 1/2] such that for all y ∈ M with ky − xk < δ, it holds that
|f(x) − f(y)|1−ζ ≤ σk Proj∇ f (x)k.
Theorem 5.3 ([25, Theorem 2.3]). Let M ⊆ Rn be an analytic submanifold and
{xk : k ∈ N} ⊂ M be a sequence. Suppose that f is real analytic and, for large enough
k,
(i) there exists σ > 0 such that
|f(xk+1 ) − f(xk )| ≥ σk Proj∇ f (xk )kkxk+1 − xk k;
(ii) Proj∇ f (xk ) = 0 implies that xk+1 = xk .
Then any accumulation point of {xk : k ∈ N} ⊆ M is the only limit point.

12

J. LI, K. USEVICH, AND P. COMON

Now we apply Theorem 5.3 to the compact orthogonal group On ⊂ Rn×n and
get Corollary 5.4, which will allow us to prove the global convergence of Jacobi-G
algorithm in subsection 5.2.
Corollary 5.4. Let f be real analytic in Algorithm 1. Suppose that, for large
enough k,
(i) there exists σ > 0 such that
|f(Qk ) − f(Qk−1 )| ≥ σk Proj∇ f (Qk−1 )kkQk − Qk−1 k;

(19)

(ii) Proj∇ f (Qk−1 ) = 0 implies that Qk = Qk−1 .
Then the iterations {Qk : k ∈ N} converge to a point Q∗ ∈ On .
Remark 5.5. Under the same assumptions as in Corollary 5.4, Theorem 3.3 tells
us that Algorithm 2 converges to a stationary point.
5.2. Global convergence of Jacobi-G algorithm for matrices and 3rdorder tensors. In this section, we consider the case d ∈ {2, 3}, that is, one of the
following options:
• for a set of 3rd-order symmetric tensors {A(`) : 1 ≤ ` ≤ m} ⊆ Rn×n×n , the
cost function is
(20)

f (Q) =

m
X

k diag{A(`) •1 QT •2 QT •3 QT }k2 ;

`=1

• for a set {A(`) : 1 ≤ ` ≤ m} ⊆ Rn×n of symmetric matrices, the cost function
is
(21)

f (Q) =

m
X

k diag{QT A(`) Q}k2 .

`=1

Theorem 5.6. For the cost function (20) or (21), Algorithm 2 converges to a
stationary point of f in On , for any starting point Q0 .
Remark 5.7. In the case m = 1 and d = 3, this is the Jacobi-G algorithm for
orthogonal diagonalization of 3rd order symmetric tensors. Theorem 5.6 shows the
global convergence of this algorithm.
Before giving the proof of Theorem 5.6, we formulate several lemmas.
Lemma 5.8. In the case d ∈ {2, 3}, for the cost function τk (x) defined in (18),
the following identities hold true
(22)

(23)

τk (x) − τk (0) =

0

τk (x) =

0
1 00
1
(h (0)(x − x3 ) + hk (0)x2 ),
(1 + x2 )2 k
2

0
00
1
(hk (0)(1 − 6x2 + x4 ) + hk (0)(x − x3 )).
2
3
(1 + x )

Proof. First, by linearity of the expressions (22) and (23), and from subsection 4.4
we can prove the identities only for the case of a single tensor or matrix (i.e., the
cost function (12)). Second, the equality (23) follows from (22) by straightforward
differentiation and the fact that
0

τk (x) = (τk (x) − τk (0))0 .

GLOBALLY CONVERGENT JACOBI-TYPE ALGORITHMS

13

Hence, we are left to prove (22) for the cost function (12).
Recall the notation of subsection 4.2, and consider the case d = 3. By substitution
(14), and due to the fact that the rotation affects only first two elements on the
diagonal, we get that
2
2
2
2
τk (x) − τk (0) = hk (θ) − hk (0) = T111
+ T222
− W111
− W222
1
[6(W111 W112 − W122 W222 )(x − x5 )
=
(1 + x2 )3
2
2
2
2
− 3(W111
+ W222
− 3W112
− 3W122
− 2W111 W122 − 2W112 W222 )(x2 + x4 )]
0
1
1 00
=
[hk (0)(x − x3 ) + hk (0)x2 ],
2
2
(1 + x )
2

where the last equality follows from Lemma 4.2.
The case d = 2 is analogous: we have
2
2
2
2
τk (x) − τk (0) = T11
+ T22
− W11
− W22
2
[(2W11 W12 − 2W12 W22 )(x − x3 )
=
(1 + x2 )2
2
2
2
+ W22
− 2W11 W22 − 4W12
)x2 ]
− (W11
0
1
1 00
=
[hk (0)(x − x3 ) + hk (0)x2 ],
2
2
(1 + x )
2

where the last equality follows again from Lemma 4.2.
Lemma 5.9. In each iteration of Algorithm 2 for the cost function (20) or (21),
the following inequality holds true
√
2ε
k Proj∇ f (Qk−1 )k · kQk − Qk−1 k
|f (Qk ) − f (Qk−1 )| ≥
4
for any k ∈ N.
Proof. At each iteration, from (9) and (10) we have
(24)

0

|hk (0)| ≥ εk Proj∇ f(Qk−1 )k.

Next, for an optimal angle θ∗ = θk∗ , we have
(25)

√
kQk − Qk−1 k = kG(ik ,jk ,θ∗ ) − I n k = 2 2 sin



θ∗
2


≤
0

√

2|θ∗ |.

Note that the tangent x∗ = tan(θ∗ ) should satisfy equation τk (x∗ ) = 0.
0
If x∗ = 0 or ±1, then hk (0) = 0 from (23), hence k Proj∇ f(Qk−1 )k = 0 from (24)
and the result is obvious. Consider the case 0 < |x∗ | < 1. Then from (23) we get
0

(26)

00

hk (0) =

hk (0)(x4∗ − 6x2∗ + 1)
,
x∗ (x2∗ − 1)

and thus
0

(27)

τk (x∗ ) − τk (0) =

x∗ hk (0)
2(1 − x2∗ )

14

J. LI, K. USEVICH, AND P. COMON

by substituting (26) into (22). Finally, by combining (24)–(27), we get
0

|f (Qk ) − f (Qk−1 )| = |hk (θ∗ ) − hk (0)| =
0
1
≥ | tan(θ∗ )hk (0)| ≥
2

√

x∗ hk (0)
2(1 − x2∗ )

2ε
k Proj∇ f(Qk−1 )k · kQk − Qk−1 k
4

Proof of Theorem 5.6. Lemma 5.9 guarantees that condition (19) in Corollary 5.4
holds true. Since the cost function is analytic, by Corollary 5.4, the sequence Qk
converges to a stationary point Q∗ . Finally, by Theorem 3.3, Q∗ is a stationary point
of f in (20).
6. Jacobi-PC algorithm and its global convergence. The Jacobi-G algorithm has several disadvantages: the convergence for 4th-order tensors is currently
unknown, and the parameter ε needs to be chosen in a proper way. In this section,
we propose an new Jacobi-based algorithm, which is inspired by proximal algorithms
in convex [24] and nonconvex [?] optimization .
6.1. Jacobi-PC algorithm and its global convergence. Suppose that we
are given a twice continuously differentiable function f : SOn → R, such that
f (QG(i,j,θ) ) = f (QG(i,j,θ+π/2) )

(28)

for any Q ∈ SOn and 1 ≤ i < j ≤ n (i.e., it is π/2-periodic along any geodesic). Then
we propose the Jacobi-PC algorithm (Jacobi-C algorithm with a proximal term) in
Algorithm 5.
Algorithm 5 Jacobi-PC algorithm
Input: A smooth function f : SOn → R, and a positive δ0 , a starting value Q0 .
Output: Sequence of iterations {Qk }k≥1 .
• For k = 1, 2, . . . until a stopping criterion is satisfied do
•
Choose the pair (ik , jk ) according to the pair selection rule (6).
•
Compute the angle θk∗ that maximizes the function
h̃k (θ) = f(Qk−1 G(ik ,jk ,θ) ) − δ0 γ(θ),
where
γ(θ) = 2 sin2 (θ) cos2 (θ).
def

∗

•
Set U k = G(ik ,jk ,θk ) , and update Qk = Qk−1 U k .
• End for
Remark 6.1. The periodicity condition (28) is not necessary for the global convergence of the algorithm in Theorem 6.2, but we add it due to its presence in the
orthogonal tensor diagonalization problem. If the condition (28) does not hold, another proximal term γ(θ) may be needed. Finally, other pair selection rules than (6)
can be used.
Theorem 6.2. The sequence produced by Algorithm 5 converges to a stationary
point Q∗ ∈ On for any starting point Q0 .

GLOBALLY CONVERGENT JACOBI-TYPE ALGORITHMS

15

Proof. We first prove the convergence. Since
f (Qk ) − f (Qk−1 ) − δ0 γ(θk∗ ) = h̃k (θk∗ ) − h̃k (0) ≥ 0,
we get that
f (Qk ) − f (Qk−1 ) ≥ δ0 γ(θk∗ ) ≥ 0.

(29)

Note that f (Qk ) is bounded since On is compact. Then f (Qk ) → c < +∞ and thus
f (Qk ) − f (Qk−1 ) → 0.
By (29), we have that γ(θk∗ ) → 0. Note that γ(θ) ≥ 8|θ|2 /π 2 for θ ∈ [−π/4, π/4].
Then θk∗ → 0 and thus there exists Q∗ ∈ On such that Qk → Q∗ .
0
Next we prove that h̃k (0) → 0, that is Λik ,jk (Qk−1 ) → 0. Define
h̄(θ, Q) = f(QG(i,j,θ) ) − δ0 γ(θ)
for θ ∈ R, Q ∈ On and 1 ≤ i < j ≤ n. Let
def

M1 =

max

Q∈On ,θ∈R,
1≤i<j≤n

∂ 2 h̄
(θ, Q) .
(∂θ)2

Then M1 < +∞ since f is C2 smooth, h̄ is periodic with respect to θ and On is
compact. Therefore, we have that
0

0

0

|h̃k (0)| = |h̃k (θk∗ ) − h̃k (0)| ≤ |θk∗ |M1
0

for any Qk−1 ∈ On , θk∗ ∈ R and 1 ≤ ik < jk ≤ n, and thus h̃k (0) → 0.
Finally we prove that Q∗ ∈ On is a stationary point of f, that is Λ(Qk−1 ) → 0.
We have proved that Λik ,jk (Qk−1 ) → 0 in the above part. Now we prove other
entries of Λ(Qk−1 ) also converge to 0. For simplicity, take (ik , jk ) = (1, 2) and
(ik+1 , jk+1 ) = (1, 3) for instance. It is enough to prove that Λ1,2 (Qk ) → 0. In fact, if
we define
φ : R → R, θ 7→ Λ1,2 (Qk−1 G(θ)),
then Λ1,2 (Qk ) = φ(θk∗ ) and Λ1,2 (Qk−1 ) = φ(0). Define
φ̄(θ) = Λ1,2 (QG(i,j,θ) )
for θ ∈ R, Q ∈ On and 1 ≤ i < j ≤ n. Let
def

M2 =

max

Q∈On ,θ∈R,
1≤i<j≤n

0

|φ̄ (θ)|.

Then M2 < +∞ since φ̄ is smooth and periodic. Therefore
|Λ1,2 (Qk ) − Λ1,2 (Qk−1 )| = |φ(θk∗ ) − φ(0)| ≤ |θk∗ |M2 ,
and thus Λ1,2 (Qk ) → 0.

16

J. LI, K. USEVICH, AND P. COMON

6.2. Elementary rotations for orthogonal tensor diagonalization. The
cost function (16) in simultaneous orthogonal tensor diagonalizatiom has the property
(28), hence the Jacobi-PC algorithm is guaranteed to converge. Moreover, it allows
for finding the update using an algebraic algorithm in the cases d = 3, 4.
Let us show how to find θk∗ in every iteration of Algorithm 5. Let τek (x) =
e
hk (arctan x) be as in subsection 4.2. Then we obtain that
τek (x) =

x2
ρ(x)
−
2δ
,
0
(1 + x2 )d
(1 + x2 )2
0

where ρ(x) is the polynomial defined in (15). Then τek (x) = 0 is equivalent to ω(x) = 0,
where
0

0

ω(x) = (1 + x2 )d+1 τek (x) = ρ (x)(1 + x2 ) − 2dxρ(x) − 4δ0 x(1 − x2 )(1 + x2 )d−2
is a polynomial of degree 2d.
Note that from subsection 4.2 and π2 -periodicity of γ, we have that e
hk (θ) = e
hk (θ+
π/2) for any θ ∈ R, hence τek (x) = τek (−1/x). Now we represent the algebraic solutions
of ω(x) by this property, that is, ω(x) = 0 has the same solutions as ω(−1/x) = 0
except the possible roots at the origin. Let ξ = x − 1/x. Then
ω(x) = x2(d−d1 −d2 ) (1 + x2 )d2

d1
Y

(x − xj )(x + 1/xj )

j=1

= x2(d−d1 −d2 ) (1 + x2 )d2

d1
Y

(x2 − ξj x − 1) = x2d−d1 −2d2 (1 + x2 )d2

d1
Y

(ξ − ξj )

j=1

j=1

for some 0 ≤ d1 , d2 ≤ d. Now we have that ω(x) = 0 if and only if
Ω(ξ) =

d1
Y

(ξ − ξj ) = 0,

j=1

except the possible roots at the origin. If the algebraic roots ξj can be calculated then
the roots (xj , −1/xj ) could be deduced by rooting the polynomials x2 − ξj x − 1 = 0.
Now we restrict ourselves to the case of a single tensor (i.e. the cost function
(12)). In fact, if A is of 3rd or 4th-order, it can be shown that Ω(ξ) has algebraic
solutions and thus (xj , −1/xj ) can be determined. The following Lemma 6.3 provides
the specific form of Ω(ξ) in these cases, and is a direct generalisation of the results
the ordinary Jacobi algorithm in [8, Appendix].
Lemma 6.3. (i) Let A ∈ Rn×n×n be a 3rd order symmetric tensor and
a = 6(W111 W112 − W122 W222 );
2
2
2
2
b = 6(W111
+ W222
− 3W112
− 3W122
− 2W111 W122 − 2W112 W222 ) + 4δ0 .

Then
ω(x) = a(1 − 5x2 − 5x4 + x6 ) + b(x5 − x) = x2 (1 + x2 )[aξ 2 + bξ − 4a];
Ω(ξ) = aξ 2 + bξ − 4a.

GLOBALLY CONVERGENT JACOBI-TYPE ALGORITHMS

17

(ii) Let A ∈ Rn×n×n×n be a 4th order symmetric tensor and
a = 8(W1111 W1112 − W1222 W2222 );
2
2
2
2
b = 8(W1111
− 3W1122 W1111 − 4W1112
− 4W1222
+ W2222
− 3W1122 W2222 )+4δ0 ;

c = 8(18W1112 W1122 − 7W1111 W1112 + 3W1111 W1222
− 18W1122 W1222 − 3W1112 W2222 + 7W1222 W2222 );
d = 8(9W1111 W1122 − 32W1112 W1222 − 2W1111 W2222
2
2
2
+ 9W1122 W2222 + 12W1112
− 36W1122
+ 12W1222
) + 4δ0 ;

e = 80(6W1122 W1222 − W1111 W1222 − 6W1112 W1122 + W1112 W2222 ).
Then
ω(x) = a(x8 + 1) + b(x7 − x) + c(x6 + x2 ) + d(x5 − x3 ) + ex4
1
1
1
1
= x4 [a(x4 + 4 ) + b(x3 − 3 ) + c(x2 + 2 ) + d(x − ) + e]
x
x
x
x
= x4 (aξ 4 + bξ 3 + (4a + c)ξ 2 + (3b + d)ξ + 2a + 2c + e);
Ω(ξ) = aξ 4 + bξ 3 + (4a + c)ξ 2 + (3b + d)ξ + 2a + 2c + e.
Note that if we set δ0 = 0, we obtain exactly the expressions from [8, Appendix].
Remark 6.4. The expressions for Ω(ξ) in the case of simultaneous orthogonal
diagonalization problem can be also easily found in the same way as in Lemma 6.3,
by exploiting the additivity of the corresponding expressions in subsection 4.4.
7. Numerical results. In this section, we present numerical experiments in
order to compare the presented algorithms in the case of orthogonal diagonalization
problems for symmetric tensors. The algorithms were implemented in MATLAB
and the codes are available on request.
The setup of all the experiments is as follows:
• A diagonal tensor D is chosen. (For convenience, we choose the tensors such
that kDk = 1.)
• A random rotation matrix Q is applied to obtain
A0 = D •1 QT · · · •d QT .
• The test tensor is constructed as A = A0 + E, where E is the symmetrization
of a tensor containing realization of i.i.d. Gaussian noise with variance σ 2 .
To each test example we apply the following algorithms:
• Jacobi-C : Algorithm 1 with the order of pairs (6).
• Jacobi-G-max : Algorithm 3.
• Jacobi-G: Algorithm 2 for various values of ε.
• Jacobi-PC : Algorithm 5 for various values of δ0 (shortened to “Jacobi-P” in
the plots).
The stopping criterion is chosen to be the maximum number of iteration.
In each of the plots, we plot kAk2 − f (Qk ), which is exactly the squared norm
of the off-diagonal elements. In all the plots, the markers correspond to the places
where the new sweep starts.

18

J. LI, K. USEVICH, AND P. COMON

7.1. Test 1: equal values on the diagonal. In this subsection, we consider
10 × 10 × 10 and 10 × 10 × 10 × 10 tensors where the diagonal values are given by
1
Di...i = √ .
10
We plot the results in Figures 1 and 2.
Te st 1, d = 3, n = 10, σ = 10 − 4

0

10
J ac ob i - C
J ac ob i - G- m ax

−1

10

J ac ob i - G, ε =
J ac ob i - G, ε =

2
100n

10

−2

J ac ob i - PC , δ 0 = 0. 0001
J ac ob i - PC , δ 0 = 0. 001
J ac ob i - PC , δ 0 = 0. 01

−3

10

−4

10

10

J ac ob i - G, ε =

2
10n

J ac ob i - G, ε =

2
100n

J ac ob i - PC , δ 0 = 0. 0001
J ac ob i - PC , δ 0 = 0. 001
J ac ob i - PC , δ 0 = 0. 01

−3

10

−4

10

−5

−5

10

10

−6

0

J ac ob i - C
J ac ob i - G- m ax

−1

2
10n

kAk 2 − f (Q k )

kAk 2 − f (Q k )

−2

10

10

Te st 1, d = 4, n = 10, σ = 10 − 4

0

10

−6

50

100

150
Iteration (k)

200

250

10

300

0

50

100

(a) 3-rd order

150

200
250
Iteration (k)

300

350

400

(b) 4-th order

Fig. 1. Equal values on the diagonal, small noise.

Te st 1, d = 3, n = 10, σ = 10 − 2

0

Te st 1, d = 4, n = 10, σ = 10 − 2
0

10

10

J ac ob i - C
J ac ob i - G- m ax

J ac ob i - G, ε =

2
10n

J ac ob i - G, ε =

2
10n

J ac ob i - G, ε =

2
100n

J ac ob i - G, ε =

2
100n

J ac ob i - PC , δ 0 = 0. 0001

kAk 2 − f (Q k )

kAk 2 − f (Q k )

J ac ob i - C
J ac ob i - G- m ax

J ac ob i - PC , δ 0 = 0. 001
J ac ob i - PC , δ 0 = 0. 01

−1

10

−2

10

0

J ac ob i - PC , δ 0 = 0. 0001
J ac ob i - PC , δ 0 = 0. 001
J ac ob i - PC , δ 0 = 0. 01

−1

10

−2

50

100
Iteration (k)

(a) 3-rd order

150

200

10

0

50

100

150
Iteration (k)

200

250

300

(b) 4-th order

Fig. 2. Equal values on the diagonal, higher noise.

As we see in Figures 1 and 2, in all the examples all the methods converge to the
same cost function value. We observe that the behavior of the Jacobi-PC algorithm
it not too different from the behavior of the Jacobi-C algorithm.
The convergence of Jacobi-G-max is the fastest, but the difference is marginal.
Also, the Jacobi-G-max is typically slower in the beginning, but accelerates when
the algorithm is closer to the local maximum. Finally, if ε is small, the behavior of
Jacobi-G is almost indistinguishable from Jacobi-C, as pointed out in Remark 3.4.
7.2. Test 2: different values on the diagonal. In this subsection, we consider
10 × 10 × 10 and 10 × 10 × 10 × 10 tensors where the diagonal values are given by
i
Di...i = √
.
385

19

GLOBALLY CONVERGENT JACOBI-TYPE ALGORITHMS

We plot the results in Figures 3 and 4.
Te st 2, d = 3, n = 10, σ = 10 − 4

0

10
J ac ob i - C
J ac ob i - G- m ax

−1

10

2
10n

J ac ob i - G, ε =

2
100n

10

−2

J ac ob i - PC , δ 0 = 0. 0001
J ac ob i - PC , δ 0 = 0. 001
J ac ob i - PC , δ 0 = 0. 01

−3

10

−4

10

10

J ac ob i - G, ε =

2
10n

J ac ob i - G, ε =

2
100n

J ac ob i - PC , δ 0 = 0. 0001
J ac ob i - PC , δ 0 = 0. 001
J ac ob i - PC , δ 0 = 0. 01

−3

10

−4

10

−5

−5

10

10

−6

0

J ac ob i - C
J ac ob i - G- m ax

−1

J ac ob i - G, ε =

kAk 2 − f (Q k )

kAk 2 − f (Q k )

−2

10

10

Te st 2, d = 4, n = 10, σ = 10 − 4

0

10

−6

50

100

150
Iteration (k)

200

250

10

300

0

50

100

(a) 3-rd order

150

200
250
Iteration (k)

300

350

400

(b) 4-th order

Fig. 3. Different values on the diagonal, small noise.

Te st 2, d = 3, n = 10, σ = 10 − 2

0

Te st 2, d = 4, n = 10, σ = 10 − 2
0

10

10

J ac ob i - C
J ac ob i - G- m ax

J ac ob i - G, ε =

2
10n

J ac ob i - G, ε =

2
10n

J ac ob i - G, ε =

2
100n

J ac ob i - G, ε =

2
100n

J ac ob i - PC , δ 0 = 0. 0001

kAk 2 − f (Q k )

kAk 2 − f (Q k )

J ac ob i - C
J ac ob i - G- m ax

J ac ob i - PC , δ 0 = 0. 001
J ac ob i - PC , δ 0 = 0. 01

−1

10

−2

10

0

J ac ob i - PC , δ 0 = 0. 0001
J ac ob i - PC , δ 0 = 0. 001
J ac ob i - PC , δ 0 = 0. 01

−1

10

−2

50

100
Iteration (k)

(a) 3-rd order

150

200

10

0

50

100

150
Iteration (k)

200

250

300

(b) 4-th order

Fig. 4. Different values on the diagonal, higher noise.

In Figures 3 and 4 we see that this scenario is less favorable for Jacobi-PC: if
the value of δ0 is too high, then it slows down the convergence of the algorithm. We
also see that typically the Jacobi-G algorithms are the fastest, but the difference with
Jacobi-C is not significant again. Also for small values of ε, the behavior of Jacobi-G
resembles the behavior of Jacobi-C.
7.3. High noise and local minima. In this subsection, we consider the case
of high noise. We repeat only the 4th-order experiments (for a single tensor) from
subsection 7.2 except with σ = 10−1 . We take two different realizations of E and plot
the results in Figure 5.
In Figure 5, we see that the behavior of the algorithms is more erratic, and they
may converge to different cost function values. This is explained by the non-convexity
of the problem and presence of different local minima, which is typical for the tensor
approximation problems [18]. Next, the Jacobi-G-max algorithm here has the worst
performance. This is also explained well by the non-convexity of the problem, because
the compatibility of the Jacobi rotation with the gradient (eqn. (9)) may not be

20

J. LI, K. USEVICH, AND P. COMON

Te st 2, d = 4, n = 10, σ = 10 − 1
J ac ob i - C
J ac ob i - G- m ax

0.93

10

0.91

10

2
10n

J ac ob i - G, ε =

2
100n

0.86

10

J ac ob i - PC , δ 0 = 0. 0001
J ac ob i - PC , δ 0 = 0. 001
J ac ob i - PC , δ 0 = 0. 01

0.89

10

J ac ob i - C
J ac ob i - G- m ax

0.87

10

J ac ob i - G, ε =

kAk 2 − f (Q k )

kAk 2 − f (Q k )

Te st 2, d = 4, n = 10, σ = 10 − 1

0.87

J ac ob i - G, ε =

2
10n

J ac ob i - G, ε =

2
100n

J ac ob i - PC , δ 0 = 0. 0001

0.85

10

J ac ob i - PC , δ 0 = 0. 001
J ac ob i - PC , δ 0 = 0. 01

0.84

10

0.83

10

10

0.82

10
0.85

0.81

10

10
0

100

200

300
400
500
Iteration (k)

600

700

800

0

100

200

(a) 1-st realization

300
400
500
Iteration (k)

600

700

800

(b) 2-nd realization

Fig. 5. Different values on the diagonal, high noise, different noise realizations.

optimal.
7.4. Simultaneous diagonalization. We conclude the section by a small example of simultaneous diagonalization. We take 4th order 10 × 10 × 10 × 10 tensor A
generated as in subsections 7.1 to 7.2 (for the noise level σ = 10−2 ), and consider its
10 slices B(1) , . . . , B(10) ∈ Rn×n×n along the last dimension, i.e.
(i)

Bk,l,s = Ak,l,s,i .
Then, we perform the joint diagonalization of tensors B(1) , . . . , B(m) (for m = 10) and
run the same algorithms as in the previous experiments, but for the cost function in
the case of simultaneous diagonalization. The results are plotted in Figure 6.
Te st 1, d = 3, n = 10, m = 10, σ = 10 − 2

0

Te st 2, d = 3, n = 10, m = 10, σ = 10 − 2

0

10

10

J ac ob i - C
J ac ob i - G- m ax

J ac ob i - G, ε =

2
10n

J ac ob i - G, ε =

2
10n

J ac ob i - G, ε =

2
100n

J ac ob i - G, ε =

2
100n

J ac ob i - PC , δ 0 = 0. 0001

kAk 2 − f (Q k )

kAk 2 − f (Q k )

J ac ob i - C
J ac ob i - G- m ax

J ac ob i - PC , δ 0 = 0. 001
J ac ob i - PC , δ 0 = 0. 01

−1

10

−2

10

0

J ac ob i - PC , δ 0 = 0. 0001
J ac ob i - PC , δ 0 = 0. 001
J ac ob i - PC , δ 0 = 0. 01

−1

10

−2

50

100

150
Iteration (k)

(a) Test 1

200

250

300

10

0

50

100

150
Iteration (k)

200

250

300

(b) Test 2

Fig. 6. Simultaneous diagonlization of tensor slices.

The results in Figure 6 exhibit a similar behavior to the results in subsections 7.1
to 7.2. When comparing the results of single tensor diagonalization for the same
tensors, (see Figure 2 and Figure 4, subfigures (b)), we can see that the results are
comparable, and even the simultaneous diagonalization may yield a slightly higher
cost function value. But the cost function in this case is different because the tensor
is not rotated along the last mode.

GLOBALLY CONVERGENT JACOBI-TYPE ALGORITHMS

21

8. Conclusions. We showed that by modifying the well-known Jacobi CoM
algorithm [6, 8] for orthogonal symmetric tensor diagonalization problem, it is possible
to prove its global convergence. The global convergence of Jacobi-G algorithm [17]
is proved for the case of simultaneous orthogonal symmetric matrix (or 3rd-order
tensor) diagonalization. The global convergence for 4th-order case is still unknown.
Our new proximal-type algorithm Jacobi-PC is globally convergent for a wide range of
optimization problems, and shows a good performance in the numerical experiments.
REFERENCES
[1] P. A. Absil, R. Mahony, and B. Andrews, Convergence of the iterates of descent methods
for analytic cost functions, SIAM Journal on Optimization, 16 (2005), pp. 531–547.
[2] P.-A. Absil, R. Mahony, and R. Sepulchre, Optimization Algorithms on Matrix Manifolds,
Princeton University Press, Princeton, NJ, 2008.
[3] E. Begovic and D. Kressner, Structure-preserving low multilinear rank approximation of
antisymmetric tensors, ArXiv e-prints, (2016), https://arxiv.org/abs/1603.05010.
[4] J. Cardoso and A. Souloumiac, Blind beamforming for non-gaussian signals, IEE Proceedings F (Radar and Signal Processing), 6 (1993), pp. 362–370.
[5] A. Cichocki, D. Mandic, L. D. Lathauwer, G. Zhou, Q. Zhao, C. Caiafa, and H. A.
PHAN, Tensor decompositions for signal processing applications: From two-way to multiway component analysis, IEEE Signal Processing Magazine, 32 (2015), pp. 145–163.
[6] P. Comon, Independent Component Analysis, in Higher Order Statistics, J.-L. Lacoume, ed.,
Elsevier, Amsterdam, London, 1992, pp. 29–38.
[7] P. Comon, Independent component analysis, a new concept ?, Signal Processing, 36 (1994),
pp. 287–314.
[8] P. Comon, Tensor Diagonalization, A useful Tool in Signal Processing, in 10th IFAC Symposium on System Identification (IFAC-SYSID), M. Blanke and T. Soderstrom, eds., vol. 1,
Copenhagen, Denmark, July 1994, IEEE, pp. 77–82.
[9] P. Comon, Tensors : A brief introduction, IEEE Signal Processing Magazine, 31 (2014),
pp. 44–53.
[10] P. Comon, G. Golub, L.-H. Lim, and B. Mourrain, Symmetric tensors and symmetric tensor
rank, SIAM Journal on Matrix Analysis and Applications, 30 (2008), pp. 1254–1279.
[11] P. Comon and C. Jutten, eds., Handbook of Blind Source Separation, Academic Press, Oxford,
2010.
[12] P. Comon and M. Sorensen, Tensor diagonalization by orthogonal transforms, Report ISRN
I3S-RR-2007-06-FR, (2007).
[13] L. De Lathauwer, Algebraic methods after prewhitening, in Handbook of Blind Source Separation, Independent Component Analysis and Applications, P. Comon and C. Jutten, eds.,
Academic Press (Elsevier), 2010, pp. 155–177.
[14] L. De Lathauwer, B. De Moor, and J. Vandewalle, Blind source separation by simultaneous third-order tensor diagonalization, in 1996 8th European Signal Processing Conference
(EUSIPCO 1996), 1996, pp. 1–4.
[15] L. De Lathauwer, B. De Moor, and J. Vandewalle, On the best rank-1 and rank-(r1 ,r2
,. . .,rn) approximation of higher-order tensors, SIAM Journal on Matrix Analysis and
Applications, 21 (2000), pp. 1324–1342.
[16] L. De Lathauwer, B. De Moor, and J. Vandewalle, Independent component analysis and
(simultaneous) third-order tensor diagonalization, IEEE Transactions on Signal Processing, 49 (2001), pp. 2262–2271.
[17] M. Ishteva, P.-A. Absil, and P. Van Dooren, Jacobi algorithm for the best low multilinear
rank approximation of symmetric tensors, SIAM J. Matrix Anal. Appl., 2 (2013), pp. 651–
672.
[18] M. Ishteva, P.-A. Absil, S. Van Huffel, and L. De Lathauwer, Tucker compression and
local optima, Chemometrics and Intelligent Laboratory Systems, 106 (2011), pp. 57–64,
https://doi.org/10.1016/j.chemolab.2010.06.006.
[19] T. G. Kolda, Orthogonal tensor decompositions, SIAM Journal on Matrix Analysis and Applications, 23 (2001), pp. 243–255.
[20] T. G. Kolda and B. W. Bader, Tensor decompositions and applications, SIAM Review, 51
(2009), pp. 455–500.
[21] S. Krantz and H. Parks, A Primer of Real Analytic Functions, A Primer of Real Analytic
Functions, Birkhäuser Boston, 2002.

22

J. LI, K. USEVICH, AND P. COMON

[22] L.-H. Lim and P. Comon, Nonnegative approximations of nonnegative tensors, Journal of
Chemometrics, 23 (2009), pp. 432–441.
[23] C. D. M. Martin and C. F. V. Loan, A jacobi-type method for computing orthogonal tensor
decompositions, SIAM Journal on Matrix Analysis and Applications, 30 (2008), pp. 1219–
1232.
[24] N. Parikh and S. Boyd, Proximal algorithms, Found. Trends Optim., 1 (2014), pp. 127–239.
[25] R. Schneider and A. Uschmajew, Convergence results for projected line-search methods on
varieties of low-rank matrices via lojasiewicz inequality, SIAM Journal on Optimization,
25 (2015), pp. 622–646.
[26] A. Uschmajew, A new convergence proof for the higher-order power method and generalizations, Pac. J. Optim., 11 (2015), pp. 309–321.
[27] L. Wang and M. T. Chu, On the global convergence of the alternating least squares method
for rank-one approximation to generic tensors, SIAM Journal on Matrix Analysis and
Applications, 35 (2014), pp. 1058–1072.

