Solving Multilinear Systems via Tensor Inversion

. Higher order tensor inversion is possible for even order. This is due to the fact that a tensor group endowed with the contracted product is isomorphic 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 decomposition and multilinear SVD provided that some constraints are satisﬁed. Moreover, within this group 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 ﬁt the framework in the least-squares sense. These are cases when there is an odd number of modes or when each mode has distinct dimension. Numerically we solve multilinear systems using iterative techniques, namely, biconjugate gradient and Jacobi methods.

PDEs.Since the inversion of tensor impinges upon a tensor-tensor multiplication definition, the contracted product for tensor multiplication was chosen since it provides a natural setting for multilinear systems and high-dimensional eigenvalue problems considered here.It is also an intrinsic extension of the matrix product rule.Still other choices of multiplication rules should be considered for particular application in hand.For example, in the matrix case, there is the alternative multiplication of Strassen [50] which improves the computational complexity by using block structure format and the optimized matrix multiplication based on blocking for improving cache performance by Demmel [20].In a recent work of Ragnarsson and Van Loan [44,45], the idea of blocking is extended to tensors.In the work of Braman [5], an alternative tensor multiplication is used in image processing applications.Our choice of the standard canonical tensor-tensor multiplication provides a useful setting for algorithms for decompositions, inversions, and multilinear iterative solvers.
Associated with tensors are multilinear systems.Multilinear systems model many phenomena in engineering and science.For example, in continuum physics and engineering, isotropic and anisotropic elasticity are modelled [41] as multilinear systems.Approximating solutions to PDEs in high dimensions amounts to solving multilinear systems.Current tensor-based methods for solving PDEs require a reduction of the spatial dimensions and some applications of tensor decomposition techniques; here we focus on tensor iterative methods for solving high-dimensional 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 Anderson 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 conjectures 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 structure and evolution.In addition, numerical (multi)linear algebra techniques are becoming useful tools in understanding complicated models and difficult problems in quantum statistical mechanics.For example, Bai et al. [2] have developed numerical linear algebra methods for the manyelectrons Hubbard model and quantum Monte Carlo simulations.
We develop a tensor-based visualization tool for localization 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 Anderson models are even order tensors which satisfy the symmetry requirement in the tensor SVD that we describe in section 3.Moreover, the numerical results provide 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.First, we define the tensor group which provides the framework for formulating multilinear systems and tensor inversion.Second, we discuss tensor decompositions derived 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,52,16] provided that some conditions are satisfied.These decompositions appear in many signal processing applications, e.g., see [10,18] and the references therein.Moreover, we describe multilinear systems in PDEs and quantum models while providing numerical methods for solving multilinear systems.Multilinear systems which do not fit in the framework are addressed by pseudoinversion methods.

Preliminaries.
We denote the scalars in R with 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 = a ijk , (A) ij = a ij , (a) i = a i unless noted otherwise.The superscripts indicate the length of the vector or the size of the matrices.For example, b K is a vector with length 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 any N , the Einstein product is defined by the operation * N via For example, if T , S ∈ R I×J×I×J , the operation * 2 is defined by the following: 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 theory of relativity [22].Notice that the Einstein product * 1 is the standard matrix multiplication since Definition 2.3 (Tucker mode-n product).Given a tensor T ∈ R I×J×K and some matrices A ∈ R Î×I , B ∈ R Ĵ×J , and C ∈ R K×K , the Tucker mode-n products are the following: Notice that the Tucker product • n is the Einstein 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].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 ∈ R I×J is partitioned into blocks of kl matrices of size k×l.If the partition is (I, 1), then matrix A ∈ R I×J is partitioned into blocks of 1 × J row vectors of size J, denoted by a Definition 2.6 (tensor blocks of third order tensor).Given a tensor A of size I × J × K with partition (l, m, n) on I, J, and K, a third order tensor 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 matrix block (top-bottom matrix slice) is denoted by Moreover, the standard flattenings of third order tensor are concatenations of the matrix slices forming the mode-one, mode-two, and mode-three matricizations [16]: Definition 2.7 (tensor blocks of N th order tensor).Given a tensor A ∈ R I1×I2×•••×IN with N ≥ 4 with partition (J 1 , J 2 , . . ., J N ) on the index sets Here we discuss some relevant tensor blocks in this paper.Given a fourth order tensor A ∈ R I1×I2×I3×I4 with partition (1, 1, I 3 , I 4 ), A ∈ R I1×I2×I3×I4 is partitioned into matrix blocks of size I 1 × I 2 .Each block is denoted by A (3,4) i3,i4 = A(:, :, i 3 , i 4 ) ∈ R I1×I2 with i 3 = 1, . . ., I 3 and i 4 = 1, . . ., I 4 .Similarly, given a fourth order tensor A ∈ R I1×I2×I3×I4 with partition (I

Tensor group set-up.
Recall that the group of invertible N × N matrices is called the general linear group of degree n, denoted GL(n, R), where the binary operation is the matrix multiplication.Let's denote this group M N,N (R) consisting of all invertible N × N matrices with the matrix multiplication.
Next in section 3, we will show that T N,M,N,M with the Einstein product (2.2) is a group and the transformation f is an isomorphism between T N,M,N,M and M N •M,N •M .So even order tensor inverses exist through the mapping to the general linear group.Consequently, we will show that tensors with different mode dimensions have no inverses under the contraction product.These are tensor analogues to rectangular matrices.In addition, we also discuss the notion of pseudoinverses in multilinear systems for odd and even orders with distinct mode lengths.

Standard tensor decompositions.
In 1927, Hitchcock [29,30] introduced the idea that a tensor is decomposable into a sum of a finite number of rank-one tensors.Today, we refer to this decomposition as CP tensor decomposition (also known as CANDECOMP [8] or PARAFAC [28]).CP is a linear combination of rankone tensors, i.e., where The column vectors a r , b r , c r , and d r 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 [39] proved that for third order tensor, is the sufficient condition for uniqueness of T = R r=1 a r • b r • c r up to permutation and scalings where Kruskal's rank k A 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]: Another decomposition called higher order SVD (also known as Tucker and multilinear SVD) was introduced by Tucker [51,52] in which a tensor is decomposable into a core tensor multiplied (Definition 2.3) by a matrix along each mode, i.e., where T , S ∈ R I×J×K×L are fourth order tensors with four orthogonal factors A ∈ R I×I , B ∈ R J×J , C ∈ R K×K , and D ∈ R L×L .The Tucker decomposition is not unique, i.e., if there exist invertible matrices R ∈ R I×I , S ∈ R J×J , T ∈ R K×K , and U ∈ R L×L , then where its core tensor S CP ∈ R I×J×K×L 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 decompositions built from the isomorphic map.

Tensor group, multilinear systems and decompositions.
For the sake of clarity, the main discussion in this section is limited to fourth order tensors, although the definitions and theorems presented are easily extended to the general case of even ordered tensors.Here we define a group structure on a set of fourth order tensor through a push-forward map on the general 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 , where • refers to the usual matrix multiplication.
(1) According to the definition of f , we can define a map h Clearly, the map h is a bijection since f is a bijection.
(2) Since f is a bijection, for some 1 ≤ i, j ≤ I 1 I 2 , there exist some unique For every 1 ≤ r ≤ I 1 I 2 , there exist some unique u, v such that (u − 1) It follows from the properties of f that the Einstein product (2.2) can be defined through the transformation: Consequently, the inverse map f −1 satisfies

any bijection. Then we can define a group structure on T by defining
Downloaded 11/17/14 to 129.72.129.221.Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.phpfor 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 isomorphism.Proof.The generalization of the transformation f (2.5) on the set T I1,...,IN ,J1,...,JN (R) with the binary operation * N easily provides the extension for this case.
Theorem 3.4.The ordered pair Theorem 3.4 implies that odd order tensors have no inverses with respect to the operation * N , although such binary operation may exist in which the set of odd order tensors exhibits a group structure.Lemma 3.1 and Theorem 3.2 show that the transformation f (2.4) is an isomorphism between groups T and M. From Corollary 3.3, it follows that these structural properties are preserved for any ordered pair (T I1,...,IN ,J1,...,JN (R), * N ) for any N .

Multilinear systems. A linear system is conveniently expressed as
The bilinear map has the linearity properties for some scalars c, c, d, d and vectors x, x 1 , x 2 ∈ R N , y, y 1 , y 2 ∈ R M .In general, we define multilinear transformations for the following systems: See [25] for more discussions on multilinear transformation.
Multilinear systems model many phenomena in engineering and sciences.For example, in continuum physics and engineering, isotropic and anisotropic elastic models are multilinear systems [41] like where T and E are second order tensors modeling stress and strain, respectively, and the fourth order tensor C refers to the elasticity tensor.Multilinear systems are also prevalent in solving PDEs numerically.In section 4, we derive and solve multilinear systems in discretized Poisson problems and eigenvalue problems in quantum models by using tensor-based iterative methods.There are modern applications in emerging fields like system biology, e.g., sparse multilinear systems model interactions among genes and proteins in living cells [43].

Solving multilinear systems.
Our first approach for solving multilinear systems numerically is to use the Gauss-Newton algorithm for approximating A −1 through the function 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 biconjugate gradient (BiCG) and the higher order Jacobi methods 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 method requires symmetric and positive definite matrix so that the multilinear system is premultiplied by its transpose A T which is later defined in Definition 3.5.The BiCG method solves the multilinear system by searching along It follows that φ( X ) attains a minimum iteratively and precisely at an optimizer X , where A * 2 X = B.The higher order Jacobi method is also implemented for comparison.The Jacobi method for tensors is based on split-Downloaded 11/17/14 to 129.72.129.221.Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.phpting the tensor A into its diagonal entries from the lower and upper diagonal entries of the tensor.
We approximate the solution to the multilinear systems considered in section 4 using these two multilinear iterative methods; see Figure 4.2 for the pseudocodes of the algorithms.We also discuss the advantages of these methods.
3.1.2.Overdetermined multilinear systems.According to Theorem 3.4, odd order tensors have no inverses with respect to the operation * N .Nonrectangular tensors also do not have inverses.Multilinear systems with tensors of odd order or that are nonrectangular (distinct mode sizes) are referred to as overdetermined multilinear systems, for example, 1.
, and B ∈ R I×J×K×L .A −1 does not exist for either case.We extend the concepts of pseudoinversion for odd order tensors and nonrectangular tensors.The optimization formulations, are considered to find multilinear least-squares solutions of the systems.Note that the Frobenius norm, • F , is defined as IN .The linear least-squares (LLS) method is a well-known method for overdetermined linear systems.Often the number of observations b exceed the number of unknown parameters x in LLS, forming an overdetermined system, e.g., Ax = b, (3.4) where where r = b − Ax, the overdetermined system (3.4) is solved.This is called the LLS method.The minimizer, the vector x * , is called the 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.Definition 3.5 (transpose).A transpose of S ∈ R I×J×I×J is a tensor T which has entries t ijkl = s klij .We denote the transpose of S as More generally, tensor transposition can be described through the permutation of indices.Let permutation σ be defined as σ(i In Definition 3.5, S T σ = S klij , where σ(ijkl) = klij.A transpose of a third order tensor A ∈ R I×J×K is A T σ = A kij with σ(ijk) = kij.Recall that (A T ) T = A in the matrix case.In the following lemmas, we describe some similar properties for higher order tensors.
Lemma 3.6.Let A ∈ R I×J×K and ρ = (ρ 1 , ρ 2 , ρ 3 ) be a cyclic permutation of length three on the index set {ijk}.Then ((A T ρ1 ) T ρ2 ) T ρ3 = A. Downloaded 11/17/14 to 129.72.129.221.Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.phpProof.For the index set {ijk}, there are two cyclic permutations: 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 permutations on the index set {ijst} and six cyclic permutations of length four.
For N th order tensors, the number of tensor transposes is N !with (N − 1)! cyclic permutations on the index set {i 1 i 2 . . .i N }.
Henceforth, we drop the subscript σ in the transposes.

Consider the multilinear system
where A ∈ R I×J×K , x ∈ R K , and B ∈ R I×J and define Lemma 3.9.Any minimizer x ∈ R K of φ 1 satisfies the following system: where A T ∈ R K×I×J with entries A T kij = (A) ijk .Proof.We expand the objective function, and calculate ∂φ1 ∂x (x).It follows that Clearly, the minimizer x of φ 1 satisfies For the problem where A ∈ R I×J×S×T , X ∈ R S×T ×K×L , and B ∈ R I×J×K×L and the objective function we have the following lemma.
Lemma 3.10.Any minimizer X ∈ R S×T ×K×L of φ 2 satisfies the system Remark 3.11.We omit the proof for Lemma 3.10 since it is similar to that of Lemma 3.9.The critical points, x = (A T * 3.2.Decompositions via isomorphic group structures.Theorem 3.2 implies that (T, * 2 ) is structurally similar to (M, •).Thus we endow (T, * 2 ) with the group structure such that (T, * 2 ) and (M, •) are isomorphic as groups.Here we discuss some of the definitions, theorems, and decompositions preserved by the transformation.
Definition 3.12 (diagonal tensor).A tensor D ∈ R I×J×I×J is called diagonal if d ijkl = 0 when i = k and j = l.Definition 3.13 (identity tensor).The identity tensor I is a diagonal tensor with entries where It generalizes to for the 2N th order identity tensor.

4). The SVD for tensor A has the form
where U ∈ R I×J×I×J and V ∈ R I×J×I×J are orthogonal tensors and D ∈ R I×J×I×J is a diagonal tensor with entries σ ijij called the singular values.Moreover, (3.11) can be written as a sum of fourth order tensors where the left and right singular matrices, U ( From the isomorphic property (3.2) and Theorem 3.2, we have where U and V are orthogonal matrices and D is a diagonal matrix.In addition, 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 Cardoso [18] when I = J.The fourth order tensor in [18] is a quadricovariance tensor in the blind identification of underdetermined mixtures problems in signal processing.Also, see Figure 3.2 for the comparison of the core tensors from CP decomposition (2.6) and tensor SVD (3.18).
Lemma 3.20.Let T ∈ R I×J×I×J and R = rank(f (T )), where f is the transformation in (2.4).The tensor SVD (3.11) in Theorem 3.17 is equivalent to CP (2.6) Proof.Define r = k + (l − 1)I.Then where σ klkl = σrr .Since u ijkl = (U Then, Moreover, the factor matrices A ∈ R I×IJ , B ∈ R J×IJ , C ∈ R I×IJ , and D ∈ R J×IJ are built from concatenating the vectors a r , b r , c r , and d r , respectively. Remark 3.21.To satisfy existence and uniqueness 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, f (T ) is sufficiently low rank, that is, R = rank(f (T )) < IJ for some I and J, then (2.7) holds.Downloaded 11/17/14 to 129.72.129.221.Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.phpLemma 3.22.Let T ∈ R I×J×I×J and R = rank(f (T )), where f is the transformation in (2.4).The tensor SVD (3.11) in Theorem 3.17 is equivalent to multilinear SVD (2.8) Proof.From (3.12), we have Remark 3.23.Typically, the core tensor of a multilinear 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 I 2 J 2 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 Corollary 3.24.Let T ∈ R I×J×I×J be symmetric and R = rank(f (T )), where f is the transformation in (2.4).The tensor EVD (3.14)   to be rank-one matrices.In Corollary 3.25, the partial symmetry t jikl = t ijkl implies that P (3,4) kl is symmetric (as well as rank-one).Thus, (P σrr a ir a jr a îr a ĵr .This decomposition is known as symmetric CP decomposition [11]. where Ω = {(x, y) : 0 < x, y < 1} with boundary Γ, f is a given function, and We compute an approximation of the unknown function v(x, y) in (4.1).Several problems in physics and mechanics are modeled by (4.1), where 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-direction.From the standard central Downloaded 11/17/14 to 129.72.129.221.Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.phpdifference approximations, the difference formula, (4.2) is obtained.If we assume Δx = Δy, then the difference equation (4.2) is equivalent to where , and (4.5) The entries of V and F are the values on the mesh on the unit square where, ( Here the Dirichlet boundary conditions are imposed so the values of V are zero at the boundary of the unit square, i.e., v i0 = v iN +1 = v 0,j = v N +1j = 0 for 0 < i, j < N + 1.Typically, V and F are vectorized which gives the following linear system: Poisson's equation in two dimensions is expressed as a sum of Kronecker products [19], i.e., Moreover, the discretized problem in three dimensions is where A N ∈ R N ×N ×N ×N and matrices V ∈ R N ×N and F ∈ R N ×N are the discretized functions v and f on a unit square mesh defined in (4.5).The nonzeros entries of the matrix block A N (2,4) k=α,l=β ∈ R N ×N are the following: where ĀN ∈ R N ×N ×N ×N ×N ×N and V, F ∈ R N ×N ×N .Both V and F are discretized on the unit cube.The entries on the tensor block ( ĀN ) The tensor representation of Poisson problems consists 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 Poisson's equation is near one [19].Also, the approximation in Figure 4.2 is first order accurate.
There are some advantages in computing in a tensor structured domain.First, without vectorization, the solution and the computational grid have a one-to-one correspondence so that sophisticated boundary conditions are easily integrated in the multilinear system.Second, the Laplacian tensor preserves the low bandwidth since the main nodal points sit on the tensor diagonal entries and the rest of the stencil points lie on the off-diagonal positions.Although 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 bandwidth than the Laplacian matrix.In fact, reducing the bandwidth of these sparse matrices directly and substantially improves the number of operations and storage locations; see the reordering methods of Cuthill and McKee [13] and George and Liu [24].

Multilinear systems in eigenvalue problems: The Anderson model and localization properties.
The Anderson model is the most studied model for understanding spectral and transport properties of an electron in a disordered medium.In 1958, Anderson [1] described the behavior of electrons in a crystal with impurities, that is, when electrons can deviate from their sites by hopping from atom to atom and are constrained to an external random potential modeling the random environment.He argued heuristically that electrons in such 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 Z d .More specifically, the Anderson model is a random Hamiltonian H ω on 2 (Z d ), d ≥ 1, defined by  The random Schrödinger operator models disordered 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 away from the lattice positions.This phenomena can be attributed to imperfections in the crystallization, glassy materials, 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) = e itHω ψ 0 .
Thus the spectral properties of H ω are studied to extract valuable information.In this case, the localization properties of the Anderson model are of interest.For instance, the localization properties are characterized by the spectral properties of the Hamiltonian H ω ; see the references [31,38,49]  Suppose f (x) is the Dirac delta function, i.e., Then V f(x) = v(x 0 )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 localization in an interval I if for all q > 0 and all φ ∈ 2 (Z d ) with compact support, where χ I is an indicator function and X is a multiplicative operator from 2 (Z d ) → Dynamical localization in this form implies that all moments of the position operator 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 localization for d = 1 for all energies and arbitrary disorder λ; see Kunz and Souillard [40] and Carmona, Klein, and Martinelli [7].For any d, it is known that the localized states 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 amount 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 beams and X-ray micro diffraction (e.g., see [32,33]) have been used to study elemental composition to understand how the atoms are arranged and to identify some defects.Powerful computer simulations and mathematical modeling are being applied increasingly as visualization tools for understanding crystal structure and evolution.Here we develop a tensor-based visualization tool for localization and for verifying some conjectures in high dimension (d ≥ 3) for all disorder λ and at various distributions.

Approximation of eigenvectors.
To approximate the eigenvectors of the multidimensional Anderson model, the eigenvalue decomposition 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 that the center nodes are around zero and have random entries, i.e., (H ω where σ and τ are random numbers with uniform distribution on [−1, 1] accounting for the random diagonal potential V ω .Higher order 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 representation of the Hamiltonian is necessary before the appropriate Einstein product rules and mappings are applied.
In       Schrödinger model in one dimension which are consistent with the results in [31] for the Anderson model in one dimension.Observe that for large amount of disorder (e.g., λ = 1), the localized states are apparent.However this is not true for a smaller amount of disorder (e.g., λ = .1).The localization is not so apparent for λ = .1 for Downloaded 11/17/14 to 129.72.129.221.Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.phpN = 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 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 than 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 in 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 localization is detected in the plots in Figure 4.5 when λ = 1.The plots in Figure 4.6 are generated by applying the higher order orthogonal iteration algorithm [17] to the Hamiltonian tensors (4.14) and (4.15).
The numerical results provide some validations that these localizations exist for large disorder for dimensions d > 1 for a sufficient amount of atoms.

Conclusion.
We have shown that even order tensors with a rectangular form 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 invertible, we have extended the notion of pseudoinversion for these cases.Alternative multilinear SVD and EVD decompositions arise from these isomorphic properties.Notably, these decompositions give a natural framework for the eigenvalue problems of quantum models like the discrete Schrödinger and Anderson models as well as solving Poisson problems in high dimensions.Moreover, CP and multilinear SVD can be computed through these alternative decompositions provided that some symmetry constraints hold.These alternative decompositions also provide a simple way to factorize quadricovariance 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 requirement low and (2) a oneto-one correspondence between the solution and the computational grid facilitate the integration of complicated boundary conditions.Moreover, the tensor representation of eigenvalue problems in quantum statistical mechanics yields visualization tools for studying localization properties of the Anderson model in high dimension.
To improve the applicability of tensor-based methods in solving high-dimensional eigenvalue problems and multilinear systems, we plan to further develop these proposed methods in several research lines.We plan to make the iterative solvers more efficient.Starting with the high-dimensional Poisson problems and eigenvalue problems, a novel method which can operate on each sparse 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 applications in analyzing genes and proteins interaction.Further developments of tensor SVD and tensor EVD for partial symmetric fourth order tensors are needed to provide direct 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 operation on a set G is a rule that assigns to each ordered pair (A, B) of elements of G some element of G. Downloaded 11/17/14 to 129.72.129.221.Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.phpThus E * 2 A = A * 2 E = A for ∀A ∈ T. Therefore E i1i2j1j2 = δ i1j1 δ i2j2 is the identity element for * 2 on T.
Finally, we know that f −1 (I I1I2×I1I2 ) = E. (A3) Show that for each A ∈ T, there exists an inverse A such that A * 2 A = A * 2 A = E.
It follows that for each A ∈ T, there exists an inverse A such that A * 2 A = A * 2 A = E.
Therefore, the ordered pair (T, * 2 ) is a group where 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.
×I 3 on the (i 2 , i 4 ) entry of the main block matrix (b) Matrix block A (3,4) i 3 i 4 ∈ R I 1 ×I 2 on the (i 3 , i 4 ) entry of the main block matrix
and b ∈ R m with m > n.Through minimization of the residual, min x r 2 , are unique minimizers for (3.6) and (3.8), respectively, since φ 1 and φ 2 are quadratic functions.Equations (3.7) and (3.9) are called the high order normal equations.See Tables 3.1-3.2for specific tensor transposes in multilinear systems.
kl ) îĵ = (V (3,4) r ) îĵ = c îr d ĵr , it follows that a ir b jr a îr b ĵr with identical factors, A = C and B = D.As in Remark 3.17, the existence of the factors A and B requires the matrix blocks P (3,4) kl

(4. 8 )Fig. 4 . 1 .
Fig. 4.1.Stencils for higher order tensors.The main node of the five-point and seven-point stencils sits on the diagonal entries of fourth order and sixth order Laplacian tensors, respectively.

( 4 . 12 )
Fig. 4.2.A solution to the Poisson equation in two dimensions with Dirichlet boundary conditions.

Remark 4 . 1 .
The random potential V ω is a multiplication operator on 2 (Z d ) with matrix elements V ω (x) = v x (ω), where (v x (ω)) x∈Z d is a collection of i.i.d.random variables with distribution ρ indexed by Z d .
-4.7, the eigenvectors of the Anderson model are compared against the eigenvectors of the discrete Schrödinger operator model for dimension d = 1, 2, 3. Recall that the discrete Schrödinger operator models perfect crystals without disorder, that is, where λ = 0.In Figure 4.5, the 2D eigenvectors of the Anderson and the discrete Schrödinger are calculated through the tensor SVD (3.18) based on the isomorphic map and the standard multilinear SVD (2.8).

1 Fig. 4 . 7 .
Fig. 4.7.Two views (first column and second column) of the three-dimensional eigenvectors of the Anderson model (left) and the discrete Schrödinger operator (right) for varying disorder.

Table 3 . 1
Third order tensor transposes in multilinear systems.
in Theorem 3.18 is equivalent to CP (2.6) if there exist A ∈ R I×I×J , B ∈ R J×I×J such that a ikl b jkl = p ijkl .Corollary 3.25.Let T ∈ R I×J×I×J with symmetries t ijkl = t klij and t jikl = t ijkl with R = rank(f (T )).The tensor EVD (3.14) in Theorem 3.18 is equivalent to CP (2.6) if there exist A ∈ R I×I×J , B ∈ R J×I×J such that a ikl b jkl = p ijkl .Remark 3.26.From Corollary 3.24, we have T ij îĵ = . The Hamiltonian H ω exhibits spectral localization if H ω has almost surely pure point spectrum with exponentially decaying eigenfunctions.= σ p (H) ∪ σ ac (H) ∪ σ sc (H) corresponding to the invariant subspaces H p of point spectrum, to H ac of absolutely continuous spectrum, and to H sc of singular continuous spectrum.The localization properties of the Anderson model can be described by spectral or dynamical properties.Let I ⊂ R. Definition 4.3.We say that H ω exhibits spectral localization in I if H ω almost surely has pure point spectrum in I (with probability one), that is, σ(H ω ) ∩ I ⊂ σ p (H ω ) with probability one.Moreover, the random Schrödinger operator H ω has exponential spectral localization in I, and the eigenfunctions corresponding to eigenvalues in I decay exponentially.Thus if for almost all ω, the random Hamiltonian H ω has a complete set of eigenvectors (ψ ω,n ) n∈N in the energy interval I satisfying Downloaded 11/17/14 to 129.72.129.221.Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Downloaded 11/17/14 to 129.72.129.221.Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php