SIAM J. MATRIX ANAL. APPL. ©c 2013 Societ y for Industrial and Applied Mathematics
Vol. 34, No. 2, pp. 542–570
SOLVING MULTILINEAR SYSTEMS VIA TENSOR INVERSION∗
M. BRAZELL† , N. LI‡ , C. NAVASCA§ , A ND C. TAMON
¶
Abstract. Higher order tensor inversion is possible for eve n order. This is due to the fact that
a tensor group endowed with the contracted product is isomo rphic to the general linear group of
degree n. With these isomorphic group structures, we derive a tensor SVD which we have shown
to be equivalent to well-known canonical polyadic decomposi tion and multilinear SVD provided
that some constraints are satisfied. Moreover, within this gro up structure framework, multilinear
systems are derived and solved for problems of high-dimensional PDEs and large discrete quantum
models. We also address multilinear systems which do not fit the framework in the least-squares
sense. These are cases when there is an odd number of mod es or when each mode has distinct
dimension. Numerically we solve multilinear systems using itera tive techniques, namely, biconjugate
gradient and Jacobi methods.
Key words. tensor and matrix inversions, multilinear s ystem, tensor decomposition, least-
squares method polynomial
AMS subject classifications. 15A69, 15A18, 65F15, 65F 99
DOI. 10.1137/100804577
1. Introduction. Tensor decompositions have b een successfully applied across
many fields which include, among others, chemometr ics [54, 48, 6], signal process-
ing [10, 15], and computer vison [53]. More recent applications are in large-scale
PDEs through a reduced rank representation of opera tors with applications to quan-
tum chemistry [34] and aerospace engineering [21]. S tate-of-the-art tensor methods
have been applied to problems in quantum chemistry: Khoromskij, Khoromskaia, and
Flad [35] solved the fundamental Hatree–Fock equation and Beylkin and Mohlenkamp
[3, 4] worked on multidimensional operators in quan tum models. Hackbusch and
Khoromskij [26] and Hackbusch, Khoromskij, and Tyr tyshnikov [27] have solved mul-
tidimensional boundary and eigenvalue problems usi ng a reduced low-dimensional
tensor-product space through separated representati on and hierarchical Kronecker
tensor from the underlying high spatial dimensions. Se e the survey papers [15, 34, 36]
and the references therein for more applications and te nsor-based methods. Extensive
studies (e.g., [11, 14, 16, 37]) have exposed many asp ects of the differences between
tensors and matrices despite the fact that tensors are m ultidimensional generalizations
of matrices.
In this paper, we continue to investigate the rela tionship between matrices and
tensors. Here we address the following question: w hen is it possible to matricize
(tensorize) and apply matrix-(tensor) based methods to high-dimensional problems
and data with inherent tensor (matrix) structure. S pecifically, we address tensor
inversion through group theoretic structures and by providing numerical methods
for specific multilinear systems in quantum mechanica l models and high-dimensional
∗ Received by the editors August 5, 2010; accepted for publi cation (in revised form) by B. Hen-
drickson February 20, 2013; published electronically May 14, 2 013. The second and third authors
were both supported in part by National Science Foundation gr ant DMS-0915100.
http://www.siam.org/journals/simax/34-2/80457.html
† Department of Mechanical Engineering, University of Wyoming, Laramie, WY 82071
(michaeljbrazell@gmail.com).
‡ Department of Mathematics, Clarkson University, Potsdam , NY 13699 (nali@clarkson.edu).§Corresponding author. Department of Mathematics, Univ ersity of Alabama at Birmingham,
Birmingham, AL 35294 (cnavasca@uab.edu).
¶ Department of Computer Science, Clarkson University, Potsdam, NY 13699 (tino@clarkson.edu).
542
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
SOLVING MULTILINEAR SYSTEMS VIA TE NSOR INVERSION 543
PDEs. Since the inversion of tensor impinges upon a te nsor-tensor multiplication def-
inition, the contracted product for tensor multiplicatio n was chosen since it provides
a natural setting for multilinear systems and high-dim ensional eigenvalue problems
considered here. It is also an intrinsic extension of the matrix product rule. Still
other choices of multiplication rules should be consid ered for particular application
in hand. For example, in the matrix case, there is th e alternative multiplication of
Strassen [50] which improves the computational comp lexity by using block structure
format and the optimized matrix multiplication based o n blocking for improving cache
performance by Demmel [20]. In a recent work of Rag narsson and Van Loan [44, 45],
the idea of blocking is extended to tensors. In the wor k of Braman [5], an alternative
tensor multiplication is used in image processing applic ations. Our choice of the stan-
dard canonical tensor-tensor multiplication provides a useful setting for algorithms
for decompositions, inversions, and multilinear iterativ e solvers.
Associated with tensors are multilinear systems. M ultilinear systems model many
phenomena in engineering and science. For example, i n continuum physics and engi-
neering, isotropic and anisotropic elasticity are model led [41] as multilinear systems.
Approximating solutions to PDEs in high dimensions amounts to solving multilinear
systems. Current tensor-based methods for solving P DEs require a reduction of the
spatial dimensions and some applications of tensor dec omposition techniques; here we
focus on tensor iterative methods for solving high-dim ensional Poisson problems in
the multilinear system framework.
Tensor representations are also common in large discrete quantum models like
the discrete Schrödinger and Anderson models. The A nderson model1 [1] is the most
studied model for spectral and transport properties of an electron in a disordered
medium. The study of spectral theory of the Anderson model is a very active research
topic, but there are still many open problems and co njectures for high-dimensional
d ≥ 3 cases; see [31, 38, 49] and the references therein .
Powerful computer simulations and mathematical modeling are increasingly used
as visualization tools for understanding crystal structu re and evolution. In addition,
numerical (multi)linear algebra techniques are becomin g useful tools in understanding
complicated models and difficult problems in quantu m statistical mechanics. For
example, Bai et al. [2] have developed numerical linear algebra methods for the many-
electrons Hubbard model and quantum Monte Carlo s imulations.
We develop a tensor-based visualization tool for loc alization and for verifying some
conjectures in high dimension (d ≥ 3) for all disorder λ and at various distributions.
The Hamiltonians of the discrete Schrödinger and An derson models are even order
tensors which satisfy the symmetry requirement in th e tensor SVD that we describe
in section 3. Moreover, the numerical results provi de some validation that these
localizations exist for large disorder for dimension d > 1 for a sufficient amount of
atoms; see the conjectures in [31, 38, 49].
The contributions of this paper are three-fold. Fi rst, we define the tensor group
which provides the framework for formulating multili near systems and tensor inver-
sion. Second, we discuss tensor decompositions deriv ed from the isomorphic group
structure and show that they are special cases of the well-known canonical polyadic
(CP) decomposition [8, 28] and multilinear SVD [51, 5 2, 16] provided that some con-
ditions are satisfied. These decompositions appear in m any signal processing applica-
tions, e.g., see [10, 18] and the references therein. Mo reover, we describe multilinear
1P.W. Anderson received the Nobel Prize in Physics in 19 77 for his work on the spectral and
transport properties of an electron in a disordered medium.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
544 M. BRAZELL, N. LI, C. NAVASCA, AN D C. TAMON
systems in PDEs and quantum models while providin g numerical methods for solv-
ing multilinear systems. Multilinear systems which d o not fit in the framework are
addressed by pseudoinversion methods.
2. Preliminaries. We denote the scalars in R w ith lowercase letters (a, b, . . . )
and the vectors with bold lowercase letters (a,b, . . . ) . The matrices are written as
bold uppercase letters (A,B, . . . ) and the symbol for tensors are calligraphic letters
(A,B, . . . ). The subscripts represent the scalars (A) ijk = aijk, (A)ij = aij , (a)i = ai
unless noted otherwise. The superscripts indicate the length of the vector or the size
of the matrices. For example, bK is a vector with len gth K and B
N×K is a N ×K
matrix.
The order of a tensor refers to the cardinality of the index set. A matrix is a
second order tensor and a vector is a first order tensor .
Definition 2.1 (even and odd tensors). Given an N th tensor T ∈ R
I1×I2×···×IN .
If N is even (odd), then T is an even (odd) N th order tensor.
Definition 2.2 (Einstein product [22]). For a ny N , the Einstein product is
defined by the operation ∗N via ∑
(2.1) (A ∗N B)i1...iNkN+1...kM = ai1i2...iNk1...k N bk1...kNkN+1kN+2...kM ,
k1...k N
where A ∈ RI1×···×IN×K1×···×KN and B ∈ RK1×···×
KN ×KN+1×···×KM .
For example, if T ,S ∈ RI×J×I×J , the∑op∑eration ∗
2 is defined by the following:
I J
(2.2) (T ∗2 S)ijîĵ = tijuvs uvîĵ .
u=1 v=1
The Einstein product is a contracted product that it is widely used in the area
of continuum mechanics [41] and in the study of the t heory of relativity [22]. Notice
that the Einstein product ∗1 is the sta∑ndard matrix m ultiplication since K
(2.3) (M ∗1 N)ij = miknkj = (M N)ij
k=1
for M ∈ RI×K ,N ∈ RK×J .
Definition 2.3 (Tucker mode-n product). Give n a tensor T ∈ R
I×J×K and
some matrices A ∈ RÎ×I , B ∈ RĴ×J , and C ∈ RK̂×K , the Tucker mode-n products
are the following: ∑ I
(T •1 A)̂i,j,k = tijkaîi ∀î, j, k (mo de-1 product),
∑i=1 J
(T •2 B) ĵ,i,k = tijkbĵj ∀ĵ, i, k (mode-2 product),
∑j=1 K
(T •3 C)k̂,i,j = tijkck̂k ∀k̂, i, j (mo de-3 product).
k=1
Notice that the Tucker product •n is the Einste in product ∗1 when the mode
summation is specified.
Below we define matrix and tensor blocks. These blockings have been described
in the papers of Ragnarsson and Van Loan [44, 45].
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
SOLVING MULTILINEAR SYSTEMS VIA TE NSOR INVERSION 545
Definition 2.4 (partition of indices). Let I be the index set with n parti∑tions and
|I| cardinality. Then there exist Ik subindex sets for k = 1, . . . , n such that k |Ik| =
|I|. Collectively, (n1, n2, . . . , nm) denotes the number of partitions on multiple index
sets I1, I2, . . . , Im.
Definition 2.5 (matrix blocks). Given a matrix A of size I × J with partition
(k, l) on I and J , respectively, a matrix A ∈ RI×J is partitioned into blocks of kl
matrices of size k×l. If the partition is (I, 1), then matr ix A ∈ RI×J is partitioned into
(1)
blocks of 1× J row vectors of size J , denoted by ak = A(k, :) ∈ JR for k = 1, . . . , I.
Similarly, if the partition is (1, J), then matrix A ∈ RI×J is partitioned into blocks
of I × (2)1 column vectors of size I, denoted by al = A I (:, l) ∈ R for l = 1, . . . , J .
Definition 2.6 (tensor blocks of third order ten sor). Given a tensor A of size
I × J ×K with partition (l,m, n) on I, J , and K, a t hird order tensor A ∈ RI×J×K
is partitioned into l ·m · n number of tensor blocks of size l ×m× n.
Here are some special cases. If the partition is (I, 1, 1), then A ∈ R
I×J×K is
partitioned into matrix blocks of size J ×K. Each ma trix block (top-bottom matrix
(1)
slice) is denoted by A = A(i, :, :) ∈ RJ×Ki for i = 1, . . . , I. If the partition is
(1, J, 1), then A ∈ RI×J×K is partitioned into matri x blocks of size I × K. Each
(2)
matrix block (left-right matrix slice) is denoted by Aj = A(:, j, :) ∈ RI×K for
j = 1, . . . , J . If the partition is (1, 1,K), then A ∈ RI× J×K is partitioned into matrix
× (3)blocks of size I K. Each matrix block (front-back ma trix slice) is denoted by Ak =S(:, :, k) ∈ RI×J for k = 1, . . . ,K. Moreover, the stan dard flattenings of third order
tensor are concatenations of the matrix slices forming the mode-one, mode-two, and
(1) (1) (1)
mode-three matricizations [16]: A(1) = [(A )T (A )T . . . (A )T1 2 I ] ∈ RK×I·J ,
(2) (2) (2) (2)A = [A A . . . A ] ∈ RI×J·K (3) (3) (3)1 2 J , andA(3) = [(A T 1 ) (A2 )T . . . (A TK ) ] ∈
RJ×K·I .
Definition 2.7 (tensor blocks of Nth order te nsor). Given a tensor A ∈
RI1×I2×···×IN with N ≥ 4 with partition (J1, J2, . . . , J N ) on the index sets I1, I2, . . . , IN ,
a tensor A ∈ RI1×I2×···×IN is partitioned into j1j2 . . . jN number of tensor blocks of
size J1 × J2 × · · · × JN .
Here we discuss some relevant tensor blocks in t his paper. Given a fourth or-
der tensor A ∈ RI1×I2×I3×I4 with partition (1, 1, I , I ), A ∈ RI1×I2×I3×I43 4 is par-
(3,4)
titioned into matrix blocks of size I1 × I 2. Each block is denoted by A =
×
i3,i4
A(:, :, i , i ) ∈ RI1 I23 4 with i3 = 1, . . . , I3 and i4 = 1, . . . , I4. Similarly, given a
fourth order tensor A ∈ RI1×I2×I3×I4 with partition (I1, I2, 1, 1), A ∈ RI1×I2×I3×I4
(1,2)
is partitioned into matrix blocks Ai ,i of size I
3 × I4. See Figure 2.1. Now given1 2
a sixth order tensor A ∈ RI1×I2×I3×I4×I5×I6 with pa rtition (1, I2, 1, I4, 1, I6), A ∈
I ×I (2,4,6)R 1 2×I3×I4×I5×I6 is partitioned into tensor blocks A i i ,i of size I1 3 5 1 × I3 × I5. See
Figure 2.2.
2.1. Tensor group set-up. Recall that the grou p of invertible N ×N matrices
is called the general linear group of degree n, denote d GL(n,R), where the binary
operation is the matrix multiplication. Let’s denote t his group MN,N (R) consisting
of all invertible N ×N matrices with the matrix mult iplication.
Definition 2.8 (transformation). Define the tra nsformation f : TN,M,N,M(R)
−→ MN ·M,N ·M (R) with f(A) = A defined component- wise as
f
(2.4) (A)nmnm −→ (A)[n+(m−1)N ][n+( m−1)N ],
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
546 M. BRAZELL, N. LI, C. NAVASCA, AN D C. TAMON
(2,4) (3,4)
(a) Matrix block A ∈ RI1×I3i i on the (b) Matrix block A I1×I22 4 i i ∈ R on the3 4
(i2, i4) entry of the main block matrix (i3, i4) ent
ry of the main block matrix
Fig. 2.1. Fourth order tensor A ∈ RI1×I2×I3×I4 with I 1 = I2 = I3 = I4 = 4 with index
partitioning: (1, I2, 1, I4) (left) and (1, 1, I3, I4) (right).
Fig. 2.2. Sixth order tensor A ∈ RI1×I2×I3×I4×I5×I6 with index partitioning (1, I2, 1, I4, 1, I6)
(2,4,6) (2,4,6)
into tensor blocks Ai i ,i of size I1 × I3 × I5. The tensor blocks A1 3 5 i i i are arranged in a third2 4 6
order tensor structure on the (i2, i4, i6) entry of the main block tensor.
where TN,M,N,M(R) = {A ∈ N×M×N×MR : det(f(A)) = 0}. In general, the transfor-
mation is defined as
f
(2.5) (A) i1i2...inj1j2...j ∑ ∏n −→ (A)[i + n (i k−1 ∑n ∏k−11 k=2 k−1) l=1 Il ][j1+ k=2(jk−1) l=1 Jl]
when A ∈ TI1,...,In,J1,...,Jn(R) and A ∈ MI1·I2...In−1·In, J1·J2...Jn−1·Jn(R).
Remark 2.9. These transformations are known as column (row) major format in
several computer languages which are typically used to enhance efficiency in accessing
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
SOLVING MULTILINEAR SYSTEMS VIA TE NSOR INVERSION 547
arrays. This mapping is also the matrix unfolding of fourth order tensors in signal
processing applications, e.g., see [18].
Next in section 3, we will show that TN,M,N,M with the Einstein product (2.2) is a
group and the transformation f is an isomorphism betw een TN,M,N,M andMN ·M,N ·M .
So even order tensor inverses exist through the mappi ng to the general linear group.
Consequently, we will show that tensors with differ ent mode dimensions have no
inverses under the contraction product. These are te nsor analogues to rectangular
matrices. In addition, we also discuss the notion of pseudoinverses in multilinear
systems for odd and even orders with distinct mode le ngths.
2.2. Standard tensor decompositions. In 192 7, Hitchcock [29, 30] introduced
the idea that a tensor is decomposable into a sum o f a finite number of rank-one
tensors. Today, we refer to this decomposition as C P tensor decomposition (also
known as CANDECOMP [8] or PARAFAC [28]). CP i s a linear combination of rank-
one tensors, i.e.,
∑ R
(2.6) T = ar ◦ br ◦ cr ◦ dr ,
r=1
where T ∈ RI×J×K×L, a ∈ RI , b ∈ RJ , c ∈ RKr r r , and dr ∈ RL. The column vectors
ar, br, cr, and dr form the so-called factor matrices A , B, C, and D, respectively.
The tensorial rank [30] is the minimum R ∈ N such that T can be expressed as a sum
of R rank-one tensors. Moreover, in 1977 Kruskal [3 9] proved that for third order
tensor,
2R+ 2 ≤ kA(A) + kB(B∑) + k C(C)
is the sufficient condition for uniqueness of T = R r=1 ar ◦ br ◦ cr up to permutation
and scalings where Kruskal’s rank kA is the maximum
number r such that any set
of r columns of A is linearly independent. Kruskal’s uniqueness condition was then
generalized for n ≥ 3 by Sidiropoulous and Bro [47]:
∑ n
(2.7) 2R+ (n− 1) ≤ kA(j) (A(j) )
∑ j=1
T R (1)for = r=1 ar ◦ · · · ◦
(n)
a . r
Another decomposition called higher order SVD ( also known as Tucker and mul-
tilinear SVD) was introduced by Tucker [51, 52] in w hich a tensor is decomposable
into a core tensor multiplied (Definition 2.3) by a mat rix along each mode, i.e.,
(2.8) T = S •1 A •2 B •3 C •4 D ,
where T ,S ∈ RI×J×K×L are fourth order tensors wit h four orthogonal factors A ∈
I×I , B ∈ J×JR R , C ∈ K×KR , and D ∈ L×LR . Th e Tucker decomposition is not
unique, i.e., if there exist invertible matrices R ∈ RI×I , S ∈ RJ×J , T ∈ RK×K , and
U ∈ RL×L, then
T = S •1 A •2 B •3 C •4 D = S̄ •1 Ā • 2 B̄ •3 C̄ •4 D̄,
where S̄ = • R−1 • S−1 • T−1 • U−11 2 3 4 , Ā = RA, B̄ = SB, C̄ = TC, and D̄ = UD.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
548 M. BRAZELL, N. LI, C. NAVASCA, AN D C. TAMON
The CP decomposition is a special case of the Tu cker decomposition framework,
i.e., ∑ R
(2.9) T = ar ◦ br ◦
cr ◦ dr = SCP •1 A • 2 B •3 C •4 D,
r=1
where its core tensor S I×J×K×L CP ∈ R is diagonal, that is, the nonzero entries are
located at (S)iiii for i = 1, . . . , R with R = min{I, J ,K, L}. Next, we discuss the
connections of CP and Tucker to the tensor decomposi tions built from the isomorphic
map.
3. Tensor group, multilinear systems and d ecompositions. For the sake
of clarity, the main discussion in this section is limit ed to fourth order tensors, al-
though the definitions and theorems presented are e asily extended to the general
case of even ordered tensors. Here we define a grou p structure on a set of fourth
order tensor through a push-forward map on the gen eral linear group. Also several
consequential results from the group structure will be discussed.
Lemma 3.1. Let f be the map defined in (2.4). Then the following properties
hold:
1. The map f is a bijection. Moreover, there exists a bijective inverse map
f−1 : MI1I2,I1I2(R) → TI1,I2,I1,I2(R).
2. The map satisfies f(A∗2 B) = f(A) · f(B), wh ere · refers to the usual matrix
multiplication.
Proof. Here we denote [n] = {1, 2, . . . , n} and its c ardinality as |n|.
(1) According to the definition of f , we can defin e a map h : [I1]× [I2] → [I1I2]
by h(i1, i2) = i1 + (i2 − 1)I1. Clearly, the map h is a b ijection since f is a bijection.
(2) Since f is a bijection, for some 1 ≤ i, j ≤ I1I2, there exist some unique
indices i1, i2, j1, j2 for 1 ≤ i1, j1 ≤ I1, 1 ≤ i2, j2 ≤ I2 su ch that (i2 − 1)I1 + i1 = i and
(j2 − 1)I1 + j1 = j. So, ∑
[f(A ∗2 B)]ij = (A ∗2 B)i1i2j1j 2 = ai1i2uvbuvj1j2 ,
∑ u,v
| I1I2|
[f(A) · f(B)]ij = [f(A)]ir[f(B)]rj.
r=1
For every 1 ≤ r ≤ I1I2, there exist some unique u, v su ch that (u − 1)I1 + v = r. So,∑ |I∑ 1I2|
ai1i2uvbuvj1j2 = [f(A)]ir [f( B)]rj.
u,v r=1
It follows from the properties of f that the Einste in product (2.2) can be defined
through the transformation:
(3.1) A ∗2 B = f−1[f(A ∗ −1 2 B)] = f [f(A) · f(B)].
Consequently, the inverse map f−1 satisfies
(3.2) f−1(A ·B) = f−1(A) ∗2 f−1 (B).
Theorem 3.2. Suppose (M, ·) is a group. Let f : T → M be any bijection. Then
we can define a group structure on T by defining
A ∗2 B = f−1[f(A) · f(B )]
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
SOLVING MULTILINEAR SYSTEMS VIA TE NSOR INVERSION 549
for all A B ∈ , T. In other words, the binary operation ∗2 satisfies the group axioms.
Moreover, the mapping f is an isomorphism.
The proof is straightforward; see the details in the appendix. As a matter of fact,
the group structure can be formulated as a ring isomo rphism.
Corollary 3.3. Let Im = Jm for m = 1, .
. . , N . Then the ordered pair
∗ (TI1,...,IN ,J1,...,JN (R), N ) is a group where the operation ∗N is the generalized Einstein
product in (2.1).
Proof. The generalization of the transformation f ( 2.5) on the set TI1,...,IN ,J1,...,JN
(R) with the binary operation ∗N easily provides the e xtension for this case.
Theorem 3.4. The ordered pair ( TI1,I2...,I2N−1(R), ∗N ) is not a group under the
operation ∗N .
Proof. TakeN = 2. Then T ∗2S ∈/ TI1,I2...,I2N−1(R) , where T ,S ∈ TI1,I2...,I2N−1(R).
It follows that TI1,I2...,I 2N−1(R) is not closed under ∗2. Thus the ordered pair
(TI1,I2...,I2N−1(R), ∗2) is not a group.
Theorem 3.4 implies that odd order tensors have no inverses with respect to the
operation ∗N , although such binary operation may ex ist in which the set of odd or-
der tensors exhibits a group structure. Lemma 3.1 an d Theorem 3.2 show that the
transformation f (2.4) is an isomorphism between g roups T and M. From Corol-
lary 3.3, it follows that these structural properties are preserved for any ordered pair
( TI1,...,IN ,J1,...,JN (R), ∗N ) for any N .
3.1. Multilinear systems. A linear system is c onveniently expressed as Ax =
b, where A ∈ RM×N ,x ∈ RN , and b ∈ RM . Recall th at A ∈ RM×N defines a linear
transformation L : RN → RM such that L(x) = Ax . A bilinear system is defined
through B : RM × RN → R with B(x,y) = yTBx = B •1 x •2 y, where B ∈ RM×N .
The bilinear map has the linearity properties
B(cx1 + dx2,y) = cB(x1,y) + d B(x2,y)
and
B(x, c̄y1 + d̄y2) = c̄B(x,y1) + d̄ B(x,y2)
for some scalars c, c̄, d, d̄ and vectors x,x1,x2 ∈ RN ,y ,y1,y2 ∈ RM . In general, we
define multilinear transformations for the following sy stems:
1. B •1 x •2 y = b, where B ∈ RI×J×K , x ∈ I , R y ∈ RJ , and b ∈ R
2. M∗ X ∗ Y = b, where M ∈ I×J×K×L, X ∈ K×L I×J2 2 R R , Y ∈ R , and b ∈ R
3. M ∗2 X ∗3 Y = B, where M ∈ RI×J×K×L×M×N , X ∈ RM×N×O, Y ∈
RK×L×O, and B ∈ RI×J .
See [25] for more discussions on multilinear transform ation.
Multilinear systems model many phenomena in en gineering and sciences. For ex-
ample, in continuum physics and engineering, isotropic and anisotropic elastic models
are multilinear systems [41] like
C ∗2 E = T,
where T and E are second order tensors modeling stre ss and strain, respectively, and
the fourth order tensor C refers to the elasticity tenso r. Multilinear systems are also
prevalent in solving PDEs numerically. In section 4, w e derive and solve multilinear
systems in discretized Poisson problems and eigenvalu e problems in quantum models
by using tensor-based iterative methods. There are m odern applications in emerging
fields like system biology, e.g., sparse multilinear syst ems model interactions among
genes and proteins in living cells [43].
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
550 M. BRAZELL, N. LI, C. NAVASCA, AN D C. TAMON
(a) Higher Order Biconjugate Gradient
(b) Higher Order Jacobi
Fig. 3.1. Pseudocodes for iterative solvers.
3.1.1. Solving multilinear systems. Our first a pproach for solving multilinear
systems numerically is to use the Gauss–Newton algo rithm for approximating A−1
through the function
g(X ) = A ∗2 X − I = 0 ,
where I is the identity fourth order tensor defined in (3 .10) and tensor X is unknown.
This method is highly inefficient because calculating the inversion of a Jacobian is
very expensive. To save memory and operational costs , we consider iterative methods
for solving multilinear systems. The pseudocodes in Figure 3.1 describe the bicon-
jugate gradient (BiCG) and the higher order Jacobi m ethods for solving multilinear
system A ∗2 X = B. Note that the algorithms solely use tensor computations, i.e.,
no matricizations are involved. Recall that the BiCG m ethod requires symmetric and
positive definite matrix so that the multilinear system i s premultiplied by its transpose
AT which is later defined in Definition 3.5. The BiCG method solves the multilinear
system by searching along Xk = Xk−1 + αk−1Pk−1 wit h a line parameter αk−1 and a
search direction Pk−1 while minimizing the objective fu nction φ(Xk+αk−1Pk), where
φ(X ) = 12X T ∗ T2A∗2X −X ∗2B. It follows that φ(X̂ )
attains a minimum iteratively
and precisely at an optimizer X̂ , where A∗2 X̂ = B. T he higher order Jacobi method
is also implemented for comparison. The Jacobi metho d for tensors is based on split-
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
SOLVING MULTILINEAR SYSTEMS VIA TE NSOR INVERSION 551
ting the tensor A into its diagonal entries from the low er and upper diagonal entries
of the tensor.
We approximate the solution to the multilinear s ystems considered in section 4
using these two multilinear iterative methods; see Fig ure 4.2 for the pseudocodes of
the algorithms. We also discuss the advantages of the se methods.
3.1.2. Overdetermined multilinear systems . According to Theorem 3.4,
odd order tensors have no inverses with respect to t he operation ∗N . Nonrectan-
gular tensors also do not have inverses. Multilinear sys tems with tensors of odd order
or that are nonrectangular (distinct mode sizes) are referred to as overdetermined
multilinear systems, for example,
1. A • x = B, where A ∈ RI×J×K , x ∈ RK , an d B ∈ RI×J3 ,
2. A ∗ X = B, where A ∈ RI×J×S×T ,X ∈ RS×T ×K×L, and B ∈ RI×J×K×L.
A−1 does not exist for either case. We extend the conc epts of pseudoinversion for odd
order tensors and nonrectangular tensors. The optimi zation formulations,
(3.3) min ‖A •3 x−B‖F and min ‖A ∗ X − B‖F ,
x X
are considered to find multilinear least-squares∑solution s of the systems. Note that the
Frobenius norm, ‖ · ‖F , is defined as ‖A‖2 2F = i i ...i |ai1i2...,iN | for AI1×I2×...IN .1 2 N
The linear least-squares (LLS) method is a well-kno wnmethod for overdetermined
linear systems. Often the number of observations b e xceed the number of unknown
parameters x in LLS, forming an overdetermined syst em, e.g.,
(3.4) Ax = b,
where A ∈ Rm×n, x ∈ Rn, and b ∈ Rm with m > n. Through minimization of the
residual,
min ‖r‖ 2,
x
where r = b − Ax, the overdetermined system (3.4 ) is solved. This is called the
LLS method. The minimizer, the vector x∗, is called t he least-squares solution of the
linear system (3.4).
We now describe tensor transposition to facilitate the discussions on the normal
equations for multilinear systems.
S Definition 3.5 (transpose). A transpose of ∈ RI×J×I×J is a tensor T
which has entries tijkl = sklij . We denote the tra nspose of S as T = ST . If
S ∈ IR 1×I2···×IN×J1×···×JN , then (T ) = (STi1,i2,...,in,j1,j2,...,jn )j1,j2,...,jn,i1,i2,...,in is
the transpose of S.
More generally, tensor transposition can be descr ibed through the permutation
of indices. Let permutation σ be defined as σ(i1i2 . . . i n) = im(1)im(2) . . . im(n), where
m(j) ∈ {1, 2, . . . , n}. Then the transpose of A ∈ RI1×I2···×I N with respect to σ is
ATσ = Aim(1)i ...i .m(2) m(n)
In Definition 3.5, STσ = Sklij , where σ(ijkl) = klij . A transpose of a third order
tensor A ∈ I×J×K TR is Aσ = Akij with σ(ijk) = kij.
Recall that (AT )T = A in the matrix case. In the following lemmas, we describe
some similar properties for higher order tensors.
Lemma 3.6. Let A ∈ RI×J×K and ρ = (ρ1, ρ2, ρ3) be a cyclic permutation of
length three on the index set {ijk}. Then ((ATρ )Tρ )Tρ = A.1 2 3
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
552 M. BRAZELL, N. LI, C. NAVASCA, AN D C. TAMON
Proof. For the index set {ijk}, there are two c yclic permutations: ρ1(ijk) =
jki, ρ2(jki) = kij, ρ3(kij) = ijk and ρ̄1(ikj) = k ji, ρ̄2(kji) = jik, ρ̄3(jik) =
ikj. It follows that ((AT )T Tρ ρ )ρ = (A )T = A s ince AT = A and (AT )T1 2 3 kij ρ3 ijk ρ1 jki ρ1 ρ2
= (A T jki)ρ .2
Although there are six permutations on the index set {ijk}, there are only two
set of cyclic permutations. For the fourth order tensor , there are twenty-four permu-
tations on the index set {ijst} and six cyclic permuta tions of length four.
Lemma 3.7 (property of fourth order tensor trans pose). Let A ∈ RI×J×S×T and
ρ be a cyclic permutation ρ = (ρ1, ρ2, ρ3, ρ4) of length
four on the index set {ijst}.
Then (((AT )Tρ ρ )T Tρ )ρ = A. For Nth order tensors, th e number of tensor transposes1 2 3 4
is N ! with (N − 1)! cyclic permutations on the index set {i1i2 . . . iN}.
Henceforth, we drop the subscript σ in the transp oses.
Definition 3.8 (critical point). Let φ : Rn → R be a continuously differentiable
function. A critical point of φ is a point x̄ ∈ nR such that
∇ φ(x̄) = 0.
Consider the multilinear system
(3.5) A •3 x = B,
where A ∈ RI×J×K , x ∈ RK , and B ∈ RI×J and defin e
‖A • − ‖2 (3.6) φ1(x) = 3 x B F .
Lemma 3.9. Any minimizer x̄ ∈ RK of φ1 satisfie s the following system:
(3.7) AT ∗2 A •3 x = AT ∗2 B ,
where AT ∈ RK×I×J with entries AT kij = (A)ijk .
Proof. We expand the objective function,
φ1(x) = 〈A •3 x−B,A •3 x−B〉 = (A •3 x)T (A • x)− 2BT3 (A •3 x) +BTB
and calculate ∂φ1∂x (x). It follows that
∂φ 1
(x) = 2AT ∗2 A • T3 x− 2A ∗2 B
∂x
since
∂ [ ] ∂
(A •3 x)T (A •3 x) = ⎣
⎡∑(∑ )⎤
xkA kijAijlxl ⎦∂x ∂x
ij kl
∂
= ⎣
⎡∑⎝⎛∑ ⎞⎤ xkA kijAijlxl⎠⎦∂x [
∂ ∑kl ij ]
= x (ATk ∗ A )klx = 2(ATl ∗ 2 A)x∂x
kl
= 2AT ∗2 A •3 x
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
SOLVING MULTILINEAR SYSTEMS VIA TE NSOR INVERSION 553
and
[ ] ⎡∑(∑ )⎤
⎣ ⎡∂ ∂ ∂ ∑⎛∑ ⎞⎤2 (A •3 x)TB = 2 xkAkijBij ⎦ = 2 ⎣ ⎝ xkAkijBij⎠⎦
∂x ∂x [ ∂x∑ij k ] k ij ∂
= 2 xk(AT ∗2 B)k
∂x
k
= 2AT ∗2 B.
Clearly, the minimizer x̄ of φ1 satisfies
AT ∗2 A •3 x = AT ∗2 B .
Furthermore, the critical point is x̄ = (AT ∗ A)−1∗ AT 2 2 ∗2B provided that (A
T ∗2A)−1
exists.
For the problem
A ∗2 X = B,
where A ∈ I×J×S×TR ,X ∈ S×T×K×LR , and B ∈ I×J×K×L R and the objective
function
(3.8) φ 2 2(X ) = ‖A ∗2 X − B‖F ,
we have the following lemma.
Lemma 3.10. Any minimizer X̄ ∈ RS×T×K×L of φ2 satisfies the system
(3.9) AT ∗ T2 A ∗2 X = A ∗2 B ,
where AT ∈ RS×T×I×
J denotes the transpose of A ∈ RI×J×S×T . Moreover, the
critical point of φ2 is X̄ = (AT ∗ A)−12 ∗ T2A ∗2B provid ed that (AT ∗2A) is invertible.
Remark 3.11. We omit the proof for Lemma 3.1 0 since it is similar to that of
Lemma 3.9. The critical points, x̄ = (AT ∗2A)−1 ∗ AT2 ∗2B and X̄ = (A
T ∗2A)−1 ∗2
AT ∗2 B, are unique minimizers for (3.6) and (3.8), respectively, since φ1 and φ2
are quadratic functions. Equations (3.7) and (3.9) ar e called the high order normal
equations. See Tables 3.1–3.2 for specific tensor transp oses in multilinear systems.
Table 3.1
Third order tensor transposes in multilin ear systems.
A x B A T
RI×J×K RK RI×J RK×I×J
RJ×K×I RI RJ×K RI× J×K
RK×I×J RJ RK×I RJ× K×I
RI×K×J RJ RK×I RJ× I×K
RK×J×I RI RK×J RI×K×J
× × RJ I K RK RJ×I RK×J×I
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
554 M. BRAZELL, N. LI, C. NAVASCA, AN D C. TAMON
Table 3.2
Fourth order tensor transposes in multili near systems.
A X B AT
RI×J×S×T RS×T×K×L RS×T×K×L RS×T×I×J
RJ×S×T×I RT×I×K×L RJ×S×K×L RT×I×J×S
RS×T×I×J RI×J×K×L RS×T×K×L RI×J×S×T
RT×I×J×S RJ×S×K×L RT×I×K×L RJ×S×T×I
3.2. Decompositions via isomorphic group structures. Theorem 3.2 im-
plies that (T, ∗2) is structurally similar to (M, ·). Thus w e endow (T, ∗2) with the group
structure such that (T, ∗2) and (M, ·) are isomorphic a s groups. Here we discuss some
of the definitions, theorems, and decompositions prese rved by the transformation.
Definition 3.12 (diagonal tensor). A tensor D ∈ RI×J×I×J is called diagonal
if dijkl = 0 when i = k and j =
l.
Definition 3.13 (identity tensor). The identity tensor I is a diagonal tensor
with entries
Ii1i2j1j2 = δi1j 1δi2j2 ,
where {
1, l = k,
δlk =
0, l = k.
It generalizes to
∏ N
(3.10) (I )i1i2...iN j1j2...jN = δi kjk
k=1
for the 2N th order identity tensor.
Remark 3.14. The diagonal core tensor D ∈ RI1×···×I N of CP (2.6) has nonzero
entries di1,i2,...,in when i1 = · · · = iN .
Definition 3.15 (orthogonal tensor). A tensor U ∈ RI×J×I×J is orthogonal if
UT ∗2 U = I, where I is the identity tensor under the binary operation ∗2.
Definition 3.16 (symmetric tensor). A tensor S ∈ RI1×I2···×IN×J1×···×JN is
symmetric if S = ST , that is, si1,i2,...,in,j1,j2,...,jn = sj1 ,j2,...,jn,i1,i2,...,in .
Theorem 3.17 (SVD). Let A ∈ RI×J×I×J with R = rank(f(A)), where f is the
transformation in (2.4). The SVD for tensor A has th e form
(3.11) A = U ∗ T 2 D ∗2 V ,
where U ∈ RI×J×I×J and V ∈ RI×J×I×J are orthogon al tensors and D ∈ RI×J×I×J
is a diagonal tensor with entries σijij called the singula r values. Moreover, (3.11) can
be written as ∑∑
A (3,4) (3,4)(3.12) = σklkl(Ukl )ij ◦ (V kl )̂iĵ ,
kl ij,̂iĵ
(3,4)
a sum of fourth order tensors where the left and right singular matrices, U
ij
and
(3,4)
Vij , are matrix blocks of U and V, respectively.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
SOLVING MULTILINEAR SYSTEMS VIA TE NSOR INVERSION 555
The symbol ◦ denotes the outer product where A ijkl = Bij ◦Ckl = BijCkl.
Proof. Let A = f(A). From the isomorphic prope rty (3.2) and Theorem 3.2, we
have
−1
A = U ·D ·VT −f−→ A = U ∗ 2 D ∗2 V
T ,
where U and V are orthogonal matrices and D is a diagonal matrix. In addition,
f−1
U ·UT = I and V ·VT = I −−→ UT ∗2 U = I and VT ∗2 V = I.
Theorem 3.18 (eigenvalue decomposition (EVD ) for symmetric tensor). Let
Ā ∈ I×J×I×JR and R = rank(f(A)), where f is the t ransformation in (2.4). Ā is a
real symmetric tensor if and only if there is a real ort hogonal tensor P ∈ RI×J×I×J
and a real diagonal tensor D̄ ∈ RI×J×I×J such that
(3.13) Ā = P ∗ D̄ ∗ PT , 2 2
where P ∈ I×J×I×J I×J×I×JR is an orthogonal tensor and D̄ ∈ R is a diagonal tensor
with entries σ̄klkl called the e∑igen∑values. Moreover, (3 .13) can be written as
Ā (3,4) (3,4)(3.14) = σ̄klkl(P kl )ij ◦ (Pkl )̂iĵ ,
kl ijîĵ
(3,4)a sum of fourth order tensors where the eigenmatrices P ∈ RI×Jkl are the matrix
blocks of P.
Proof. Through the mapping provided through the isomorphic property (3.2) and
Theorem 3.2, there exist an orthogonal tensor P and a diagonal tensor D̄ such that
Ā P ∗ D̄ ∗ PT (3,4) (3,4)= 2 2 . Moreo∑ver, the fourth order tenso r P̂ijîĵ = (Pkl )ij ◦ (Pkl )̂iĵ∑ P̂ ∑ (3,4) ◦ (3,4) ∑ ∑is symmetric since ijîĵ = ij,̂iĵ(Pkl )ij (Pkl )̂ Tiĵ = ij,̂iĵ PijklP = P ·îĵkl s rs
PT = P ·PTr̂s s r̂s rs = ij,̂iĵ PîĵklPT
(3,4) (3,4)
ijkl = (Pkl )̂iĵ ◦ (Pkl )ij = P̂îĵij .
(3,4) (3,4)
Remark 3.19. If the eigenmatrix (Pkl ) is sy mmetric, that is, (Pkl )ij =
(3,4)
(P kl )ji, then the entries of Ā have the following symmetry: ājilk = āijkl. If
ā jilk = āijkl and āijkl = āklij , then (3.14) is exactly the tensor eigendecomposition
found in the paper of De Lathauwer, Castaing, and C ardoso [18] when I = J . The
fourth order tensor in [18] is a quadricovariance tenso r in the blind identification of
underdetermined mixtures problems in signal processi ng. Also, see Figure 3.2 for the
comparison of the core tensors from CP decomposition (2.6) and tensor SVD (3.18).
3.3. Connections to standard tensor decom positions. Here we relate ten-
sor SVD (3.11) to both CP (2.6) and Tucker (2.8) dec ompositions.
Lemma 3.20. Let T ∈ RI×J×I×J and R = rank(f (T )), where f is the transfor-
mation in (2.4). The tensor SVD (3.11) in Theorem 3.17 is equivalent to CP (2.6)
if there exist A ∈ I×I×J ,B ∈ J×I×J , C ∈ I×I×J J×I×JR R R , and D ∈ R such that
aiklbjkl = uijkl and cîkldĵkl = vîĵkl.
Proof. Define r = k + (l − 1)I. Then
∑ ∑ R
T (3,4) (3,4) (3,4) (3,4)ijîĵ = σklkl(Ukl )ij ◦ (Vkl )̂iĵ = σ̄r r(Ur )ij ◦ (Vr )̂iĵ ,
kl r=1
(3,4) (3,4)where σklkl = σ̂rr. Since uijkl = (Ukl )ij = (Ur )ij = airbjr and vîĵkl =
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
556 M. BRAZELL, N. LI, C. NAVASCA, AN D C. TAMON
(3,4)
(a) Matrix block Si i of fourth or-3 4
der core tensor from the CP dec om-
position
(3,4) (2,4)
(b) Matrix block Si i of fourth or- (c) Matrix block Si i of fourth or-3 4 2 4
der core tensor from tensor SVD der core ten sor from tensor SVD
Fig. 3.2. Matrix blocks of the core tensors A ∈ R3×3×3×3 w ith nonzero diagonal entries (purple
dots).
(3,4) (3,4)
(Vkl )̂iĵ = (Vr )̂iĵ = cîrdĵr, it follows that ∑ R
Tijîĵ = σ̂rr(U(3,4)) ◦ (V(3,4)r ij r )̂iĵ
∑r=1 ∑ R R
= σ̂rr(U
(3,4)) (3,4) r ij ◦ (Vr )̂iĵ = σ̂rrairbjrcîrd ĵr.
r=1 r=1
Then,
∑ R
(3.15) T = σ̄rrar ◦ br ◦ cr ◦ d r.
r=1
Moreover, the factor matrices A ∈ RI×IJ ,B ∈ RJ×IJ ,C ∈ RI×IJ , and D ∈ RJ×IJ
are built from concatenating the vectors ar, br, cr, an d dr, respectively.
Remark 3.21. To satisfy existence and uniquenes s of the CP decomposition, the
inequality (2.7) must hold. If R = IJ = rank(f(T ) ), then (3.15) does not satisfy
(2.7). However, if f(T ) is sufficiently low rank, that is, R = rank(f(T )) < IJ for
some I and J , then (2.7) holds.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
SOLVING MULTILINEAR SYSTEMS VIA TE NSOR INVERSION 557
Lemma 3.22. Let T ∈ RI×J×I×J and R = rank(f (T )), where f is the transfor-
mation in (2.4). The tensor SVD (3.11) in Theorem 3 .17 is equivalent to multilinear
SVD (2.8) if there exist A ∈ RI×I ,B ∈ RJ×J , C ∈ R I×I , and D ∈ RJ×J such that
a b = u ik jl ∑ijkl
and cikdjl = vijkl. ∑ (3 ,4) (3,4)
Proof. From (3.12), we have Tijîĵ = kl σklkl(Ukl )ij ◦ (Vkl )̂iĵ which implies
Tijîĵ = kl σklklaikbjlcîkdĵl.
Remark 3.23. Typically, the core tensor of a m ultilinear SVD (2.8) is dense.
However, the core tensor resulting from Lemma 3.22 is not dense (possibly sparse),
i.e., there are IJ nonzeros elements from the total I2 J2 entries in the fourth order
core tensor of size I×J×I×J . Similarly, the existence of the decomposition impinges
upon the existence of the factors A ∈ RI×I ,B ∈ RJ×J , C ∈ RI×I , and D ∈ RJ×J
such that U = A ◦B and V = C ◦D.
Corollary 3.24. Let T ∈ RI×J×I×J be symmetr ic and R = rank(f(T )), where
f is the transformation in (2.4). The tensor EVD (3.14 ) in Theorem 3.18 is equivalent
to CP (2.6) if there exist A ∈ RI×I×J ,B ∈ RJ×I×J su ch that aiklbjkl = pijkl.
Corollary 3.25. Let T ∈ I×J×I×JR with symm etries tijkl = tklij and tjikl =
tijkl with R = rank(f(T )). The tensor EVD (3.14) in Theorem 3.18 is equivalent to
CP (2.6) if there exist A ∈ RI×I×J ,B ∈ RJ×I×J such th∑at aiklbjkl = pijkl.∑ R (3,4) (3,4)Remark 3.26. From Corollary 3.24, we have Tijîĵ = r=1 σ̄rr(Pr )ij ◦(Pr )̂iĵR
= r=1 σ̄rrairbjraîrbĵr with identical factors, A = C and B = D. As in Remark
(3,4)
3.17, the existence of the factors A and B requires t he matrix blocks Pkl to be
rank-one matrices. In Corollary 3.25, the partial sym metry tjikl = tijkl implies that
∑(3,4) (3,4)Pkl is symmetric (as well as rank-one). Thus, (P kl )ij = aiklajkl =⇒ T =R
r=1 σ̄rrairajraîraĵr. This decomposition is known as symmetric CP decomposi-
tion [11].
Corollary 3.27. Let T ∈ RI×J×I×J be symmetr ic and R = rank(f(T )), where
f is the transformation in (2.4). The tensor EVD (3.14 ) in Theorem 3.18 is equivalent
to multilinear SVD (2.8) if there exist A ∈ RI×I and B ∈ R
J×J such that∑aikbjl =pijkl.
Remark 3.28. Th∑e multilinear SVD from Coro llary 3.27 is Tijîĵ = kl σklkl(3,4) (3,4)
(Pkl )ij ◦ (P )̂ kl iĵ = kl σklklaikbjlaîkbĵl following Lemma 3.22.
4. Numerical examples involving multilinea r systems.
4.1. Multilinear systems in Poisson problem s. Consider the two-dimensio-
nal (2D) Poisson problem
−∇2v = f in Ω,
(4.1)
u = 0 on Γ,
where Ω = {(x, y) : 0 < x, y < 1} with boundary Γ, f is a given function, and
∂2∇2 v ∂
2v
v = + .
∂x2 ∂y2
We compute an approximation of the unknown functio n v(x, y) in (4.1). Several prob-
lems in physics and mechanics are modeled by (4.1), w here the solution v represents,
for example, temperature, electromagnetic potential, or displacement of an elastic
membrane fixed at the boundary.
The mesh points are obtained by discretizing the unit square domain with step
sizes, Δx in the x-direction and Δy in the y-directio n. From the standard central
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
558 M. BRAZELL, N. LI, C. NAVASCA, AN D C. TAMON
difference approximations, the difference formula,
vl−1,m − 2vl,m + vl+1,m v l,m−1 − 2vl,m + v l,m+1(4.2) + = f(xl, ym),
Δx2 Δy2
is obtained. If we assume Δx = Δy, then the differen ce equation (4.2) is equivalent
to
(4.3) ANV +VAN = (Δx)
2F ,
where ⎡ ⎤
2 − 1 0
(4.4) A = ⎢⎢⎢⎢ . ⎥
⎣ −1 2
. . ⎥⎥
N . ⎥ , . .. . . −1 ⎦
⎡ 0 −1 2 ⎢ ⎤⎢ v11 v12 . . . v1N ⎢ . ⎢ v v . .. ⎥⎣ 21 22 . . ⎥(4.5) V = ⎥. , and. . . . . ⎥. . . vN− ⎦1N
⎡ vN1 . . . vNN−1 vNN ⎢ ⎤⎢ f11 f12 . . . f1N⎢ ⎢ . . ⎥⎣ f21 f
. . .. ⎥
F = 22 ⎥. .. . . . ⎥. . . . f ⎦N−1N
f . . . f − f N1 NN 1 NN
The entries of V and F are the values on the mesh on th e unit square where, (xi, yj) =
(iΔx, jΔx) ∈ [0, 1] × [0, 1]. Here the Dirichlet bound ary conditions are imposed so
the values of V are zero at the boundary of the unit sq uare, i.e., vi0 = viN+1 = v0,j =
vN+1j = 0 for 0 < i, j < N + 1.
Typically, V an⎡d F are vectorized which gives the following l⎤inea⎡r system ⎢ A + 2I ⎤
:
⎢ N N −IN 0 ⎥⎥ ⎢⎢ v⎢ 11⎣⎢ −
.
I . vN AN + 2IN . ⎥⎥⎦ ⎢⎣ 12 ⎥
⎥
(4.6) AN×N · v = . . .. ⎥. . . . ⎦ −I .N
⎡0 ⎤ −IN AN + 2I vN NN⎢ ⎥ ⎢ f11⎢⎣ f122 ⎦⎥⎥
= (Δx) . .
..
fNN
Poisson’s equation in two dimensions is expressed as a sum of Kronecker products
[19], i.e.,
(4.7) AN×N = IN ⊗AN +AN ⊗ IN.
Moreover, the discretized problem in three dimensions is
(4.8) (AN ⊗ IN ⊗ IN + IN ⊗AN ⊗ IN + IN ⊗ IN ⊗A N) · vec(V) = (Δx)3vec(F).
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
SOLVING MULTILINEAR SYSTEMS VIA TE NSOR INVERSION 559
(a) 5-point Stencil (b) 7-point Stencil
Fig. 4.1. Stencils for higher order tensors. The main no de of the five-point and seven-point
stencils sits on the diagonal entries of fourth order and sixth o rder Laplacian tensors, respectively.
The higher order tensor representation of the 2D discretized Poisson problem
(4.1) is
(4.9) AN ∗2 V = F,
where A ∈ N×N×N×NN R and matrices V ∈ N×NR and F ∈ N×NR are the discretized
functions v and f on a unit square mesh defined in (4. 5). The nonzeros entries of the
(2,4)
matrix block A N×NNk=α,l=β ∈ R⎧ are the following: ⎪ ⎪⎪⎪ (2,4)⎪⎪(A ) 4⎪ Nα,β α,β = (Δx)2 ,⎪⎪⎨ (2,4)(A ) − = −1Nα,β α 1β (Δx)2 ,
(4.10)
⎪⎪⎪⎪⎪⎪⎪
(2,4)
(ANα,β )α+1,β =
−1
(Δx)2 ,
⎪ (2,4)
⎩(A ) − =
−1
Nα,β α,β 1 (Δx)2 ,
(2,4)
(A −1Nα,β )α,β+1 = (Δx)2
for α, β = 2, . . . , N − 1. These entries form a five-poi nt stencil; see Figure 4.1. The
discretized three-dimensional Poisson equation is
(4.11) ĀN ∗3 V = F ,
where Ā ∈ RN×N×N×N×N×NN and V ,F ∈ RN×N×N . Both V and F are discretized
(2,4,6)
on the unit cube. The entries on the tensor block ( ĀN ) N×N×Nl,m,n ∈ R of Ā N
would follow a seven-point stencil, i.e.,
⎪⎧⎪⎪⎪⎪ Ā (2,4,6)
⎪((
6
N )
⎪ α,β,γ
)α,β,γ = (Δx)3 ,
⎪⎪⎪ (2,4,6)⎪((Ā −1
⎪ N )⎪ α,β,γ
)α−1,β,γ = (Δx )3 ,
⎪⎨ Ā (2,4,6) (( N )α,β,γ )α+1,β,γ = −1(Δx )3 ,
⎪⎪⎪ Ā (2,4,6)(4.12) (( ) −1 ⎪⎪ N α,β,γ )α,β−1,γ = (Δx ,⎪
)3
⎪⎪⎪⎪((Ā
(2,4,6)
) −1 N
⎪ α,β,γ
)α,β+1,γ = (Δx 3 ,
⎪
)
⎪⎪ Ā (2,4,6)⎩(( −1N ) α,β,γ )α,β,γ−1 = (Δx)3 , Ā (2,4,6)(( ) −1N α,β,γ )α,β,γ+1 = (Δx)3
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
560 M. BRAZELL, N. LI, C. NAVASCA, AN D C. TAMON
(a) Approximated Solution (b) Error vs. Iteration: Bicongugate Gradient
(blue -.-) and Jacobi (red –)
Fig. 4.2. A solution to the Poisson equation in two dimensions with Dirichlet boundary condi-
tions.
− for α, β, γ = 2, . . . , N 1 since vijk satisfies
6v 3ijk − vi−1jk − vi+1jk − vij−1k − vij+1k − vijk− 1 − vijk+1 = (Δx) fijk.
The tensor representation of Poisson problems consi sts of the multilinear systems
(4.9) and (4.11). The iterative methods used to solve (4.9) and (4.11) are described
in Figure 3.1; see Figure 4.2 for the numerical results . The convergence of Jacobi is
slow since the spectral radius with respect to the Pois son’s equation is near one [19].
Also, the approximation in Figure 4.2 is first order ac curate.
There are some advantages in computing in a ten sor structured domain. First,
without vectorization, the solution and the computa tional grid have a one-to-one
correspondence so that sophisticated boundary cond itions are easily integrated in
the multilinear system. Second, the Laplacian tensor preserves the low bandwidth
since the main nodal points sit on the tensor diagon al entries and the rest of the
stencil points lie on the off-diagonal positions. Altho ugh the Laplacian matrices in
(4.4) and (4.6) are banded, the Laplacian matrices in higher dimensions have larger
bandwidths. The Laplacian tensor has a lower bandwi dth than the Laplacian matrix.
In fact, reducing the bandwidth of these sparse matr ices directly and substantially
improves the number of operations and storage locatio ns; see the reordering methods
of Cuthill and McKee [13] and George and Liu [24].
4.2. Multilinear systems in eigenvalue probl ems: The Anderson model
and localization properties. The Anderson mod el is the most studied model
for understanding spectral and transport properties of an electron in a disordered
medium. In 1958, Anderson [1] described the behavio r of electrons in a crystal with
impurities, that is, when electrons can deviate from th eir sites by hopping from atom
to atom and are constrained to an external random pot ential modeling the random en-
vironment. He argued heuristically that electrons in suc h systems result in a loss of the
conductivity properties of the crystal, transforming it from conductors to insulators.
The Anderson model is a discrete random Schr ödinger operator defined on a
lattice dZ . More specifically, the Anderson model is a random Hamiltonian Hω on
2(Zd), d ≥ 1, defined by
(4.13) Hω = −Δ+ λVω ,
where Δ(x, y) = 1 if |x − y| = 1 and zero otherwise (the discrete Laplacian) with
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
SOLVING MULTILINEAR SYSTEMS VIA TE NSOR INVERSION 561
spectrum [−2d, 2d] and the random potential V dω = {Vω(x), x ∈ Z } consisting of
i.i.d. random variables on [−1, 1] which we assume to have bounded and compactly
supported density ρ. The disorder parameter is a non negative λ > 0. The spectrum
of Hω can be explicitly described by
− σ(Hω) = σ( Δ) + λ supp(ρ) = [−2d, 2 d] + λ supp(ρ).
Remark 4.1. The random potential Vω is a mu ltiplication operator on
d
2(Z )
with matrix elements Vω(x) = vx(ω), where (vx(ω))x∈Z d is a collection of i.i.d. random
variables with distribution ρ indexed by Zd.
The random Schrödinger operator models disorder ed solids. The atoms or nuclei
of ideal crystals are distributed in a lattice in a regular way. Since most solids are not
ideal crystals, the positions of the atoms may deviate a way from the lattice positions.
This phenomena can be attributed to imperfections in the crystallization, glassy ma-
terials, or a mixture of alloys or doped semiconductors . To model disorder, a random
potential Vω perturbs the pure laplacian Hamiltonian (−Δ) of a perfect metal. The
time evolution of a quantum particle ψ is determined by the Hamiltonian Hω, i.e.,
ψ(t) = eitH
ωψ0.
Thus the spectral properties of Hω are studied to ex tract valuable information. In
this case, the localization properties of the Anderson model are of interest. For in-
stance, the localization properties are characterized by the spectral properties of the
Hamiltonian Hω; see the references [31, 38, 49]. The Ha miltonianHω exhibits spectral
localization if Hω has almost surely pure point spectru m with exponentially decaying
eigenfunctions.
Remark 4.2. Recall that from [46] for any self-ad joint operator H , the spectral
decomposition is
σ(H) = σp(H) ∪ σac(H) ∪ σs c(H)
corresponding to the invariant subspaces Hp of point spectrum, to Hac of absolutely
continuous spectrum, and to Hsc of singular continuou s spectrum.
The localization properties of the Anderson mode l can be described by spectral
or dynamical properties. Let I ⊂ R.
Definition 4.3. We say that Hω exhibits spectra l localization in I if Hω almost
surely has pure point spectrum in I (with probability o ne), that is,
σ(Hω) ∩ I ⊂ σp(Hω) with probab ility one.
Moreover, the random Schrödinger operator Hω has e xponential spectral localization
in I, and the eigenfunctions corresponding to eigenval ues in I decay exponentially.
Thus if for almost all ω, the random Hamiltoni an Hω has a complete set of
eigenvectors (ψω,n)n∈N in the energy interval I satisfy ing
|ψω,n(x)| ≤ C −μ|x−xω ,n|ω,ne
with localization center x ω,n for μ > 0 and Cω,n < ∞, then the exponential spectral
localization hold on I.
Remark 4.4. Let V : 2(Z) → 2(Z) be a multip lication operator and suppose
v : Z → R is a function. Then, V f(x) = v(x)f(x) , and thus σ(V ) = range(v).
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
562 M. BRAZELL, N. LI, C. NAVASCA, AN D C. TAMON
(a) λ = 1, N = 50 (b) λ = .1, N = 50
Fig. 4.3. One-dimensional eigenvectors of the discrete Schrödinger operator (-x-) and the
Anderson model (black, -o-) for various modes.
Suppose f(x) is the Dirac delta function, i.e.,{
1 x = x0,
f(x) = δ(x− x 0) =
0 x = x0.
Then V f(x) = v(x0)f(x) which implies that σ(V ) = σ p(V ), i.e., V has a pure point
spectrum.
Definition 4.5. A random Schrödinger operator has strong dynamical localiza-
tion in an interval I if fo[r all q > 0 and all φ ∈ (Zd2 ) with]compact support,
‖ | |q −itHω E sup X e χI(Hω)ψ‖2 < ∞ ,t
where χI is an indicator function and X is a multipli cative operator from 2(Z
d) →
2(Z
d) defined as |X |ψ = |x|ψ(x).
Dynamical localization in this form implies that a ll moments of the position op-
erator are bounded in time.
The Anderson model is a well-studied subject. So there are numerous results in
both physics and mathematics literature; see [31] and the references therein. There
are several theoretical results on the existence of locali zation for d = 1 for all energies
and arbitrary disorder λ; see Kunz and Souillard [40] a nd Carmona, Klein, and Mar-
tinelli [7]. For any d, it is known that the localized sta tes are present for all energies
for sufficiently large disorder (λ 1). For d = 2 and with a Gaussian distribution,
it is conjectured that there is a localization for any a mount of disorder λ as in the
case for d = 1. There are many more open problems for higher dimension like the
extended state conjecture [23].
In material science, many techniques like X-ray bea ms and X-ray micro diffraction
(e.g., see [32, 33]) have been used to study elemental c omposition to understand how
the atoms are arranged and to identify some defects. Powerful computer simulations
and mathematical modeling are being applied increas ingly as visualization tools for
understanding crystal structure and evolution. Here we develop a tensor-based vi-
sualization tool for localization and for verifying some conjectures in high dimension
(d ≥ 3) for all disorder λ and at various distributions.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
SOLVING MULTILINEAR SYSTEMS VIA TE NSOR INVERSION 563
(a) λ = .1, N = 100
Fig. 4.4. One-dimensional eigenvectors of the discrete Schrödinger operator (-x-) and the
Anderson model (black, -o-) for various modes.
In Figures 4.3–4.7, the eigenvectors of the Anderso n model are compared against
the eigenvectors of the discrete Schrödinger operator m odel for dimension d = 1, 2, 3.
Recall that the discrete Schrödinger operator models perfect crystals without disor-
der, that is, where λ = 0. In Figure 4.5, the 2D eige nvectors of the Anderson and
the discrete Schrödinger are calculated through the te nsor SVD (3.18) based on the
isomorphic map and the standard multilinear SVD (2 .8).
4.2.1. Approximation of eigenvectors. To ap proximate the eigenvectors of
the multidimensional Anderson model, the eigenvalue d ecomposition in Theorem 3.18
is applied to the Hamiltonian Hω . The Hamiltonian H ω in two and three dimensions
are formed into fourth and sixth order tensors using the same stencils in Figure 4.1
with entries in (4.10) and (4.12), respectively. Note th at the center nodes are around
zero and have random entries, i.e.,
(3,4) σ
(4.14) (Hωk=α,l=β)α,β = (Δ x)
2
and
(4,5,6) τ
(4.15) (Hωl=α,m=β,n=γ)α,β,γ = (Δ
,
x)3
where σ and τ are random numbers with uniform d istribution on [−1, 1] account-
ing for the random diagonal potential Vω. Higher ord
er tensor representation easily
preserved the uniform distribution on [−1, 1] on the random potential, but this is
not necessarily true for Hamiltonians in (4.7) and (4 .8). To compute numerically
the higher-dimensional eigenvector, the tensor repres entation of the Hamiltonian is
necessary before the appropriate Einstein product rule s and mappings are applied.
In Figures 4.3, 4.4, 4.5, 4.6, and 4.7, the eigenfunc tions are approximated by the
eigenvectors from the both discrete Schrödinger and ra ndom Schrödinger (Anderson)
models. In Figures 4.3 and 4.4, the eigenvectors of t he Anderson model in one di-
mension are definitely more localized than the eigenv ectors of the discrete random
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
564 M. BRAZELL, N. LI, C. NAVASCA, AN D C. TAMON
(a) N=29
(b) N=48
Fig. 4.5. 2D eigenvectors of the discrete Schrödinger opera tor (left column) and the Anderson
model (right column) for varying disorder (λ = 10 (top), λ = 1 (middle), and λ = .1 (bottom)).
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
SOLVING MULTILINEAR SYSTEMS VIA TE NSOR INVERSION 565
Fig. 4.6. Factors of the multilinear SVD decomposition [16] of the 2D discrete Schrödinger
operator (right column) and the Anderson model (left column) for varying disorder (λ = 10 (top),
λ = 1 (middle), and λ = .1 (bottom)).
(a) λ = 10 ( b) λ = 10
(c) λ = 1 (d) λ = 1
(e) λ = 0.1 ( f) λ = 0.1
Fig. 4.7. Two views (first column and second column) of t he three-dimensional eigenvectors of
the Anderson model (left) and the discrete Schrödinger operato r (right) for varying disorder.
Schrödinger model in one dimension which are consist ent with the results in [31] for
the Anderson model in one dimension. Observe tha t for large amount of disorder
(e.g., λ = 1), the localized states are apparent. Howev er this is not true for a smaller
amount of disorder (e.g., λ = .1). The localization is not so apparent for λ = .1 for
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
566 M. BRAZELL, N. LI, C. NAVASCA, AN D C. TAMON
N = 50, but when the number of atoms is increased, that is, setting N = 100, the
localized eigenvectors are present as in the case when λ = 1; see Figures 4.3(b) and 4.4.
In the contour plots of Figures 4.5 and 4.7, the eigenvectors in two and three
dimensions of the Anderson model are more peaked th an those of the nonrandomized
Schrödinger for large disorder λ ≥ 1. As in the case for one dimension, localization
is not apparent for small disorder (λ = .1) as seen i n Figure 4.5. Moreover, as N
increases, the eigenstates of both discrete Schrödinger and Anderson models seem to
coincide for small disorder. This does not necessarily mean that the localization is
absent for this regime, but rather the localized states are harder to find for a small
amount of disorder. A larger amount of atoms has to be considered for this case. In
Figure 4.5, localization is not clearly visible for even λ = 1 in the factors calculated via
the multilinear SVD decomposition (2.8) while localiza tion is detected in the plots in
Figure 4.5 when λ = 1. The plots in Figure 4.6 are ge nerated by applying the higher
order orthogonal iteration algorithm [17] to the Hamilt onian tensors (4.14) and (4.15).
The numerical results provide some validations th at these localizations exist for
large disorder for dimensions d > 1 for a sufficient am ount of atoms.
5. Conclusion. We have shown that even order t ensors with a rectangular form
I1 × · · · × IN × J1 × · · · × JN are invertible by a transformation to the general linear
group equipped with the Einstein contracted product . While odd order tensors and
even order tensors with distinct modes are not invertib le, we have extended the notion
of pseudoinversion for these cases. Alternative multili near SVD and EVD decompo-
sitions arise from these isomorphic properties. Notabl y, these decompositions give a
natural framework for the eigenvalue problems of qua ntum models like the discrete
Schrödinger and Anderson models as well as solving P oisson problems in high dimen-
sions. Moreover, CP and multilinear SVD can be comp uted through these alternative
decompositions provided that some symmetry constrai nts hold. These alternative de-
compositions also provide a simple way to factorize qu adricovariance tensors in blind
identification of underdetermined mixtures problems.
Solving multilinear systems with the tensor-based iterative methods has several
advantages: (1) higher order tensor representation of PDEs preserve low bandwidth
thereby keeping the computational cost and memory r equirement low and (2) a one-
to-one correspondence between the solution and the co mputational grid facilitate the
integration of complicated boundary conditions. More over, the tensor representation
of eigenvalue problems in quantum statistical mechani cs yields visualization tools for
studying localization properties of the Anderson mode l in high dimension.
To improve the applicability of tensor-based metho ds in solving high-dimensional
eigenvalue problems and multilinear systems, we plan to further develop these pro-
posed methods in several research lines. We plan to m ake the iterative solvers more
efficient. Starting with the high-dimensional Poisson problems and eigenvalue prob-
lems, a novel method which can operate on each sp arse matrix (tensor) blocks as
opposed to the whole tensor structure will cut down the memory and operational
costs. Extensions of these methods in nonsymmetric and sparse multilinear systems
are important since these systems have direct applic ations in analyzing genes and
proteins interaction. Further developments of tensor S VD and tensor EVD for partial
symmetric fourth order tensors are needed to provide d irect methods in decomposing
quadricovariance tensors in prevalent signal processing applications.
Appendix: Proof of Theorem 3.2. Recall the basic definitions.
Definition 5.1 (binary operation). A binary op eration on a set G is a rule
that assigns to each ordered pair (A,B) of elements of G some element of G.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
SOLVING MULTILINEAR SYSTEMS VIA TE NSOR INVERSION 567
Definition 5.2. A group (G, ) is a set G with th e operation which satisfy the
following axioms:
(A1) The binary operation is associative, i.e., ( A B) C = A (B C) for
A,B,C ∈ G.
(A2) There is an element I ∈ G such that IX = XI for all X ∈ G. This element
I is an identity element for on G.
(A3) For each A ∈ G, there is an element Ã ∈ G with the property that Ã A =
A Ã = I.
Proof. Here we prove by checking each axiom (A1 –A3) hold in Definition 5.2.
(A1) Show that the binary operation ∗2 is associa tive.
Since we know that f is a bijective map with th e property that f(A ∗2 B) =
f(A) · f(B), we will show f−1(A ·B) = f−1(A) ∗ −12 f (B) for A,B ∈ M.
Let A,B, C ∈ T and A,B,C ∈ M, where f(A) = A, f(B) = B and f(C) = C.
Then,
(A ∗2 B) ∗ −12 C = f (A) ∗ f−12 (B) ∗ f−12 (C) = f−1( A ·B ·C) = f−1(A · (B ·C))
= f−1(A) ∗ f−1(B ·C) = A ∗ (f−1(B) ∗ f−12 2 2 (C)) = A ∗2 (B ∗2 C).
Therefore, (A ∗2 B) ∗2 C = A ∗2 (B ∗2 C).
(A2) Show that there is an identity element for ∗2 on T.
Since II1I2×I1I2 ∈ M is the identity element in the group, note that we will
suppress the superscript of I in the calculation below. Then we claim that f−1(I) is
the identity element for ∗2 on T.
For every element A ∈ T, there exists a matrix A −1 ∈ M so that f (A) = A. So,
we get
A ∗ f−1(I) = f−1(A) ∗ f−12 2 (I) = f−1(A · I) = f−1(A) = A.
Similarly,
f−1(I) ∗ A = f−1(I) ∗ f−12 2 (A) = f−1(I · A) = f−1(A) = A.
Therefore, A ∗ f−12 (I) = f−1(I) ∗2 A = A.
Define the tensor E as follows:
(E)i1i2j1j2 = δi1j1δ i2j2 ,
where {
1, l = k,
δlk =
0, l = k.
We claim that E = f−1(I). By d∑irect calculations, we have
(E ∗2 A)i i 1 2j1j2 = i1i2uvauvj1j2 = i1i2i1i2ai1i2j1j2
u,v
= δi1i1δi2i2ai1i2j1j2 = ai1 i2j1j2 = Ai1i2j1j2
∑ and
(A ∗ 2 E)i1i2j1j2 = ai1i2uvuvj1j2 = ai1 i2j1j2j1i2j1j2
u,v
= ai1i2j1j 2δj1j1δj2j2 = ai1i2j1j2 = Ai1i2j1j2 .
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
568 M. BRAZELL, N. LI, C. NAVASCA, AN D C. TAMON
Thus E ∗2 A = A ∗
2 E = A for ∀A ∈ T. Therefore Ei1i 2j1j2 = δi1j1δi2j2 is the identity
element for ∗2 on T.
Finally, we know that f−1(II1I2×I1I2) = E.
(A3) Show that for each A ∈ T, there exists an inverse Ã such that Ã ∗2 A =
A ∗2 Ã
= E.
We define Ã = f−1{[f(A)]−1} since f(A) ∈ M an d f−1 is a bijection map from
Lemma (3.1). Then,
f(Ã ∗2 A) = f(Ã) · f(A) = [f(A)]−1 · f (A) = I
I1I2×I1I2 .
From Lemma 3.1 and since f(E) = II1I2×I1I2 , we obtain Ã ∗2 A = E. Similarly,
we can get A ∗2 Ã = E.
It follows that for each A ∈ T, there exists an i nverse Ã such that Ã ∗2 A =
A ∗2 Ã = E.
Therefore, the ordered pair (T, ∗2) is a group wh ere the operation ∗2 is defined
in (2.2). In addition, the transformation f : T → M in (2.4) is a bijective mapping
between groups. Hence, f is an isomorphism.
Acknowledgments. C.N. would like to thank S hannon Starr for some fruitful
discussions on the quantum models. The authors are indebted to the referees whose
suggestions and corrections tremendously improved th is manuscript.
REFERENCES
[1] P.W. Anderson, Absence of diffusion in certain rando m lattices, Phys. Rev., 109 (1958),
pp. 1492–1505.
[2] Z. Bai, W. Chen, R. Scalettar, and I. Yamazaki, Num erical methods for quantum Monte
Carlo simulations of the Hubbard model, in Multi-Sc ale Phenomena in Complex Fluids,
T.Y. Hou, C. Liu, and J.-G. Liu, eds., Higher Edu cation Press, Beijing, China, 2009,
pp. 1–114.
[3] G. Beylkin and M.J. Mohlenkamp, Numerical operator calculus in higher dimensions, Proc.
Natl. Acad. Sci. USA, 99 (2002), pp. 10246–10251.
[4] G. Beylkin and M.J. Mohlenkamp, Algorithms for num erical analysis in high dimensions,
SIAM J. Sci. Comput., 26 (2005), pp. 2133–2159.
[5] K. Braman, Third order tensors as linear operators on a space of matrices, Linear Algebra
Appl., 433 (2010), pp. 1241–1253.
[6] D.S. Burdick, L.B. McGown, D.W. Millican, and X.M . Tu, Resolution of multicomponent
fluorescent mixtures by analysis of the excitation-emis sion-frequency array, J. Chemomet-
rics, 4 (1990), pp. 15–28.
[7] R. Carmona, A. Klein, and F. Martinelli, Anderson localization for Bernoulli and other
singular potentials, Comm. Math. Phys., 108 (1987), pp. 41–66.
[8] J.D. Carroll and J.J. Chang, Analysis of individual diff erences in multidimensional scaling
via an N-way generalization of “Eckart-Young” decom position, Psychometrika, 35 (1970),
pp. 283–319.
[9] H. Cohn, R. Kleinberg, B. Szegedy, and C. Umans, G roup-theoretic algorithms for matrix
multiplication, in Proceedings of the 46th Annual Sym posium on Foundations of Computer
Science, Pittsburgh, PA, IEEE Computer Society, 200 5, pp. 379–388.
[10] P. Comon, Tensor decompositions: State of the art and applications, in Mathematics in Signal
Processing V, J.G. McWhirter and I.K. Proudler, eds., Oxford University Press, Oxford,
2001, pp. 1–24.
[11] P. Comon, G. Golub, L.-H. Lim, and B. Mourrain, Sym metric tensors and symmetric tensor
rank, SIAM J. Matrix Anal. Appl., 30 (2008), pp. 125 4–1279.
[12] D. Coppersmith and S. Winograd, Matrix multiplicat ion via arithmetic progressions, J.
Symbolic Comput., 9 (1990), pp. 251–280.
[13] E. Cuthill and J. McKee, Reducing the bandwidth of spa rse symmetric matrices, in Proceed-
ings of the 24th ACM National Conference, New York , 1969, pp. 157–172.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
SOLVING MULTILINEAR SYSTEMS VIA TE NSOR INVERSION 569
[14] V. Da Silva and L.-H. Lim, Tensor rank and the ill-pos edness of the best low-rank approxi-
mation problem, SIAM J. Matrix Anal. Appl., 30 (200 8), pp. 1084–1127.
[15] L. De Lathauwer, A Survey of Tensor Methods, ISCAS, Taipei, 2009.
[16] L. De Lathauwer, B. De Moor, and J. Vandewalle, A multilinear singular value decom-
position, SIAM J. Matrix Anal. Appl., 21 (2000), pp. 1253–1278.
[17] L. De Lathauwer, B. De Moor, and J. Vandewalle, On the best rank-1 and rank-
(R1, R2, . . . , R
N ) approximation of higher-order tensors, SIAM J. Matrix Anal. Appl.,
21 (2000), pp. 1324–1342.
[18] L. De Lathauwer, J. Castaing, and J.-F. Cardoso, Fou rth-order cumulant-based blind iden-
tification of underdetermined mixtures, IEEE Trans. S ignal Process., 55 (2007), pp. 2965–
2973.
[19] J. Demmel, Applied Numerical Linear Algebra, SIAM, Ph iladelphia, 1997.
[20] J. Demmel, Lecture Notes on Cache Blocking, University of California Berkeley, Berkeley, CA,
http://www.cs.berkeley.edu/˜demmel/cs267 Spr99/Le ctures/Lect 02 1999b.pdf (2007).
[21] A. Doostan, G. Iaccarino, and N. Etemadi, A lea st-squares approximation of high-
dimensional uncertain systems, in Center for Turbulen ce Research, Annual Research Briefs,
Stanford University, Stanford, CA, 2007, pp. 121–132 .
[22] A. Einstein, The foundation of the general theory of re lativity, in The Collected Papers of
Albert Einstein 6, A.J. Kox, M.J. Klein, and R. Sch ulmann, eds., Princeton University
Press, Princeton, NJ, 2007, pp. 146–200.
[23] L. Erdös, M. Salmhofer, and H.-T. Yau, Towards the qu antum Brownian motion, in Mathe-
matical Physics of Quantum Mechanics, Lecture Notes in Physics 690, J. Asch and A. Joyé,
eds., Springer, Berlin, 2006, pp. 233–258.
[24] J.A. George and J.W.-H. Liu, Computer Solution of Large Sparse Positive Definite Systems,
Prentice-Hall, Englewood Cliffs, NJ, 1981.
[25] W.H. Greub, Multilinear Algebra, Springer-Verlag, Berlin , 1967.
[26] W. Hackbusch and B.N. Khoromskij, Tensor-product a pproximation to operators and func-
tions in high dimensions, J. Complexity, 23 (2007), p p. 697–714.
[27] W. Hackbusch, B.N. Khoromskij, and E.E. Tyrtyshn ikov, Hierarchical kronecker tensor-
product approximations, J. Numer. Math., 13 (2005), pp. 119–156.
[28] R.A. Harshman, Foundations of the PARAFAC proced ure: Models and conditions for an
“explanatory” multi-modal factor analysis, UCLA Wor king Papers in Phonetics, 16 (1970),
pp. 1–84.
[29] F.L. Hitchcock, The expression of a tensor or a polya dic as a sum of products, J. Math.
Phys., 6 (1927), pp. 164–189.
[30] F.L. Hitchcock, Multiple invariants and generalized ra nk of a p-way matrix or tensor, J.
Math. Phys., 7 (1927), pp. 39–79.
[31] D. Hundertmark, A short introduction to Anderson loca lization, in Analysis and Stochastics
of Growth Processes and Interface Models, P. Mörte rs, R. Penrose, H. Schwetlick, and
J. Zimmer, eds., Oxford University Press, Oxford, 200 8, pp. 194–218,.
[32] G.E. Ice, B.C. Larson, W. Yang, J.D. Budai, J.Z. Tisch ler, J.W.L. Pang, R.I. Barabash,
and W. Liu, Polychromatic X-ray microdiffraction studies of mesoscale structure and
dynamics, J. Synchrotron Rad., 12 (2005), pp. 155–162.
[33] G.E. Ice and B.C. Larson, Three-dimensional X-ray st ructural microscopy using polychro-
matic microbeams, MRS Bulletin, 29 (2004), pp. 170– 176.
[34] B.N. Khoromskij, Tensor-Structured Numerical Methods in Scientific Computing: Survey on
Recent Advances, Preprint 21/2010, MPI MIS, Leipzi g, 2010.
[35] B.N. Khoromskij, V. Khoromskaia, and H.-J. Flad, N umerical solution of the Hatree-Fock
equation in multilevel tensor-structured format, SIAM J. Sci. Comput., 33 (2011), pp. 45–
65.
[36] T. Kolda and B.W. Bader, Tensor decompositions and applications, SIAM Rev., 51 (2009),
pp. 455–500.
[37] T. Kolda, Orthogonal tensor decompositions, SIAM J. Ma trix Anal. Appl., 23 (2001), pp. 243–
255.
[38] W. Kirsch, An invitation to random Schrödinger operators , in Random Schrödinger Operators,
Panoramas et Synthèses 25, Société Mathématique de France, Paris, 2008, pp. 1–119.
[39] J.B. Kruskal, Three-way arrays: Rank and uniqueness of trilinear decompositions with appli-
cations to arithmetic complexity and statistics, Linea r Algebra Appl., 18 (1977), pp. 95–
138.
[40] H. Kunz and B. Souillard, Sur le spectre des opérateurs aux differénces finies aléatoires,
Comm. Math. Phys., 78 (1980), pp. 201–246.
[41] W.M. Lai, D. Rubin, and E. Krempl, Introduction to C ontinuum Mechanics, Butterworth-
Heinemann, Oxford, 2009.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php
570 M. BRAZELL, N. LI, C. NAVASCA, AN D C. TAMON
[42] C.D. Martin and C.F. Van Loan, A Jacobi-type meth od for computing orthogonal tensor
decompositions, SIAM J. Matrix Anal. Appl., 30 (200 8), pp. 1219–1232.
[43] B. Nazer and R.D. Nowak, Sparse interactions: Ident ifying high-dimensional multilinear
systems via compressed sensing, in Proceedings of th e 48th Annual Allerton Conference
on Communication, Control and Computation, Monticello, IL, 2010, pp. 1589–1596.
[44] S. Ragnarsson and C. Van Loan, Block tensor unfoldings, SIAM J. Matrix Anal. Appl., 33
(2012), pp. 149–169.
[45] S. Ragnarsson and C. Van Loan, Block tensor and sym metric embeddings, Linear Algebra
Appl., 438, pp. 853–874.
[46] M. Reed and B. Simon, Methods of Modern Mathemat ical Physics I: Functional Analysis,
Academic Press, New York, 1980.
[47] N.D. Sidiropoulos and R. Bro, On the uniqueness of multilinear decomposition of N-way
arrays, J. Chemometrics, 14 (2000), pp. 229–239.
[48] A. Smilde, R. Bro, and G. Giannakis, Multi-way Ana lysis. Applications in the Chemical
Sciences, John Wiley and Sons, Chichester, UK, 2004 .
[49] G. Stolz, An introduction to the mathematics of Anderso n localization, Contemp. Math., 552
(2010), pp. 71–108.
[50] V. Strassen, Gaussian elimination is not optimal, Nume r. Math., 13 (1969), pp. 354–356.
[51] L.R. Tucker, The extension of factor analysis to three-dim ensional matrices, in Contributions
to Mathematical Psychology, H. Gulliksen and N. Fr ederiksen, eds., Holt, Rinehardt, &
Winston, New York, 1963, pp. 110–127.
[52] L.R. Tucker, Some mathematical notes on three-mode factor analysis, Psychometrika, 31
(1966), pp. 279–311.
[53] M.A.O. Vasilescu and D.Terzopoulos, Multilinear subspace analysis for image ensembles, in
Proceedings of the IEEE Conference on Computer Vis ion and Pattern Recognition (CVPR
03), Madison, WI, 2003, pp. 93–99.
[54] J.P. Wold, R. Bro, A. Veberg, F. Lundby, A.N. Nilse n, and J. Moan, Active photosensi-
tizers in butter detected by fluorescence spectroscopy a nd multivariate curve resolution, J.
Argric. Food Chem., 54 (2006), pp. 10197–10204.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.
Downloaded 11/17/14 to 129.72.129.221. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php