{'fraction_non_alphanumeric': 0.027526690998396833, 'fraction_numerical': 0.010778729065047583, 'mean_word_length': 4.961366409109394, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 31, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solving Tensor Structured Problems with Computational Tensor Algebra

Since its introduction by Gauss, Matrix Algebra has facilitated understanding of scientific problems, hiding distracting details and finding more elegant and efficient ways of computational solving. Today's largest problems, which often originate from multidimensional data, might profit from even higher levels of abstraction. We developed a framework for solving tensor structured problems with tensor algebra that unifies concepts from tensor analysis, multilinear algebra and multidimensional signal processing. In contrast to the conventional matrix approach, it allows the formulation of multidimensional problems, in a multidimensional way, preserving structure and data coherence; and the implementation of automated optimizations of solving algorithms, based on the commutativity of all tensor operations. Its ability to handle large scientific tasks is showcased by a real-world, 4D medical imaging problem, with more than 30 million unknown parameters solved on a current, inexpensive hardware. This significantly surpassed the best published matrix-based approach.

. The standard way of dealing with such problems is not intrinsically multidimensional, and assumes reduction of problem formulation to the classical matrix formalism, which is built on two-dimensional (matrix) and one-dimensional (vector) data structures. This is achieved by the reordering of multidimensional terms from the original problem formulation, into matrices and vectors. The derived matrix formulation is then treated using well-known techniques and algorithms of matrix algebra. Finally the solution is reconstituted into its natural multidimensional form.
While this approach reduces a multidimensional problem to the well-known standard form of matrix algebra, it introduces certain limitations. The resulting formulation is no longer explicitly multidimensional, -a fact that may create difficulty in identifying and understanding important properties inherent to the particular problem [20]. Vectorization of multidimensional objects may lead to loss of spatial data coherence [6,13], which can adversely affect the performance of solving algorithms. In many cases, the derived problem formulation does not have a straightforward, intuitive connection with the process of generating efficient solving algorithms.
Currently, a growing interest in objects with more then two dimensionstensors, has become evident in the scientific community. Introduced in Tensor Analysis and Multilinear Algebra, tensors gained the attention of practitioners of diverse fields (see the excellent review by Kolda [14]). However, translating tensor mathematics to a convenient computational framework raises many issues [19]; including, convenient notation, ease of algorithm implementation, performance.
Here we propose a comprehensive framework for solving tensor structured problems by tensor algebra that allows natural and elegant formulation of multidimensional problems using multidimensional data structure directly. The proposed framework is based on a generalization of the concepts of matrix algebra to multiple dimensions; it incorporates and unifies existing approaches from multidimensional signal processing [5,10,2] , recent developments in the field of tensor analysis [12], and multilinear algebra [27,17,16,11,14,1,19]. It is based on ordered vector spaces and a generalized definition of tensor multiplication, based on the concept of tensor equations and the tensor inverse, and induces the design of fast computational solving algorithms with built-in expression analysis for automated generation of efficient implementations on serial and parallel computers. In our companion paper we show how to computationally solve large real-world problems using our framework, and claim that this novel approach offers a natural, higher level abstraction to solving a broad range of very large computational problems which arise in multidimensional signal processing and other scientific fields.

Tensor nomenclature, basic tensor operations
In our framework, multidimensional data and multidimensional transformations are described based on a formalism originating in Tensor Analysis [12] and Physics [7]. Each object in the formulation represents a multidimensional entity -a tensor. Tensors describing multidimensional data belong to a space that is the outer product of a number of vector spaces (see A.3). A tensor can be represented by a multidimensional array of its components 1 . For example, a scalar field defined on a discrete uniform three-dimensional grid with extents [N X , N Y , N Z ] can be expressed as a third order tensor T xyz , where tensor order is determined by the number of indices. Note that in our tensor notation, direct correspondence between designation of tensor indices (x, y, z) and dimensions of the physical space (X, Y, Z) is typically used. Tensors indices can be either covariant (subscripts) or contravariant (superscripts), depending on the way the tensor components transform in respect to a change of basis (see A.1). This property of a tensor index is called variance. For the common case of Euclidean spaces with orthonormal bases, there is no difference between covariant and contravariant indices; thus the choice of index variance is influenced by the context and consistency of tensor expressions.
Consider a transformation applied to tensor T xyz which results in a tensor in the same tensor product space. The transformation can be expressed by A x1y1z1 xyz . Its application to T implies the use of classical tensor operations, such as the outer product, which expands tensor order and the contraction, which reduces tensor order For designation of the contraction, we use the Einstein convention [7] which assumes implicit summation over a pair of indices with equal designation but different variance. From the practical point of view, it is convenient to combine the outer product and the contraction into a single operation, in analogy with the matrix multiplication. Thus by allowing simultaneous contraction over multiple corresponding index pairs we get Our framework relies on the fact that a tensor is an object which can implicitly contain all its properties as indices and component values. Due to this object nature in a practical implementation (2.3) can be reduced to D := A*T, where "*" denotes tensor multiplication.

Concept of ordered spaces, commutativity of tensor multiplication
In some cases, multidimensional transformations can be decomposed to a product of one-dimensional transformations. Such transformations are called sepa-1 Further in the document the term tensor is identified with tensor components rable. They are reduced to sequential application of one-dimensional transformations from the decomposition. Using tensor formalism in three dimensions this can be expressed as In conventional tensor formalism (and likewise in matrix formalism), the result of expression (2.4) depends on the order of multiplicands (lack of commutativity). This is due to concatenation of non contracted indices, which is used for assurance of a unique result of the tensor product 2 . In contrast, separable data handling operations, like filtering or resampling in signal processing, lead to the same result independently of the order in which separable transformations are applied. To avoid this formal, but practically important mismatch, and to add more flexibility to handling of tensor expressions, we introduce the following constraint: in analogy to physical space, which is ordered (right hand rule), we require vector spaces used for construction of the tensor product to have a unique predefined order. This constraint 3 leads to a fundamental new property of the framework, compared to conventional approaches for matrix and tensors. The tensor multiplication becomes commutative, in addition to being associative. Formally this means that indices of the resulted tensor have unique positions determined by a predefined order of vector spaces. Thus complex tensor products can be evaluated in arbitrary order but lead to the same result. For example, for ordered spaces X, Y, Z we can write The commutativity of the tensor multiplication, achieved hereby, considerably increases the ease of handling tensor expressions. One of the important practical values of this property consists in simplification of semantic analysis of tensor expressions performed for the purpose of reduction of computational and storage requirements. For example, expression A z1t z ·B y1 y ·C x1 x ·T xyz t , (N X < N Y < N Z < N T , N X N Y > N T ), which according to conventional tensor and matrix formalisms can be evaluated only in a single chain order, in our framework can be computed by reordering C x1 x · B y1 y · (A z1t z · T xyz t ) . This gives the best computational performance with minimal memory requirements for storing intermediate results. At this point, it is possible to use the tensor multiplication, in combination with tensor addition and subtraction, defined from multilinearity of tensor product space (see A.4) for formulating multidimensional problems in a natural and elegant multidimensional representation. Note that formulations expressed in tensor terms are also convenient for differentiation in respect to multidimensional terms (see A.8), an important prerequisite for multidimensional optimization problems.
2.2 Extension of the tensor notation According to Einstein notation it is not permitted to have more than one index with same designation and variance. However, in our settings we do allow this, with agreement on some consistency rules. As an example, consider the following expression This operation is a tensor analogue of the Khatri-Rao product [26] which is a matching elementwise (over i) Kronecker product of two sets of vectors. Note that two formally equivalent indices are merged into a single one, without dependency on outer product or contraction applied to other indices. But the constraint is, that indices can be contracted only once after all possible merges resulting in a single pair of superscript and subscript. Here is an example which represents the most general form of the tensor product, combining outer product, contraction and elementwise product Notice how neatly the proposed notation unifies different types of products which exist separately in Matrix Algebra (Matrix product, Kronecker product, Khatri-Rao product).

Tensor equations and tensor inverse
Expression A x1y1z1 xyz · U xyz = B x1y1z1 represents a relation between a multidimensional input U xyz and a multidimensional output B x1y1z1 . This relation is determined by a transformation A x1y1z1 xyz . When the input is unknown, the expression describes an inverse problem. In cases where a unique solution exists, the existence of a tensor inverse is implied -a transformation which maps A to the identity tensor x · δ y2 y · δ z2 z is the identity transformation in three dimensions expressed in terms of Kronecker delta (see A.5). Note that equation 3.1 can be represented in an equivalent form as a matrix mapping from a vector space to a vector space, with proper reshaping of tensor terms. In cases when system tensor A has some structure, treatment of the problem using the tensor representation can be more advantageous than using the conventional matrix formalism, that is explained in the next sections.

Direct tensor solvers, iterative tensor solvers, Krylov subspace tensor solvers
Here we introduce tensor solvers which can be explicitly applied to solve such tensor equations computationally. The most straightforward way to find the solution tensor and, if necessary, a tensor inverse -is based on use of a tensor extension of Gauss elimination or a tensor analogue of LU decompositionboth belonging to the class of direct tensor solvers. These methods can benefit from the explicit tensor structure of equation coefficients, which is typically manifested by handling multidimensional sub-blocks of nonzero coefficients. In many practical problems, however, tensor equations are so large that solving them using direct methods becomes prohibitive. In this instance, we propose the class of iterative tensor solvers. For example, one can use a tensor Jacobi solver which is based on the iterative computation of tensor multiplication and belongs to the class of stationary iterative solvers. The pseudocode of the algorithm and an actual software implementation are shown on Fig. 4.1.
Note that the presented object based implementation of the Jacobi solver does not use explicit indexing of tensor components. All information about indices is contained in tensor objects which are defined once at the construction stage. This is sufficient for performing further computations with automatic semantic analysis for optimizing the performance, and without complicating the solving process by explicit use of indices.
From the point of view of computational and storage requirements, the presented iterative solver is particularly attractive for application to tensor structured problems. As a basic example, consider the Poisson equation in three dimensions ∇ 2 u(x, y, z) = f (x, y, z) which after discretization on a uniform grid can be expressed by the following tensor equation Here U is an unknown tensor, F is the observation tensor, A, B, C are one-dimensional finite difference transformations which are equivalent to filters with impulse response dependent on the order of the Laplacian approximation. This decomposition allows very efficient computation of the tensor multiplication by performing sequential separable transformations along individual dimensions of the tensor U. The separability renders the computation of the whole set of equation coefficients unnecessary, thus dramatically reducing storage requirements of the presented iterative algorithm (obviously, for such simple example, exploiting the tensor characteristics of the problem is also possible using the standard matrix approach with careful index bookkeeping).
Another type of iterative algorithms which are frequently used for solving large inverse problems, is the family of Krylov subspace solvers, which includes GMRES, CG, BICGSTAB and other solvers [23]. In our extension of this solver class to tensors, most of them can be expressed as an iterative application of the tensor multiplication (see Fig. 4.2) and thus profit from the · · · END (c) Fig. 4.1 (a) Pseudocode of tensor Jacobi iterator for three dimensions (E is the inverse of diagonal tensor constructed from main diagonal of A ). (b) Actual implementation of solver for any number of dimensions using our tensor library designed in Tensor Oberon programming language, which offers algebraic operators and optimized implementations for the basic tensor operations. Here "*" and "+*" stand for tensor multiplication and inner product respectively. (c) A simple code example of using the solver. same benefits discussed above for the tensor Jacobi iterator. In contrast to Jacobi solver which only converges well for a limited class of problems, Krylov solvers perform significantly better for the large majority of problems. Fig. 4.2 presents a pseudocode and Oberon implementation of the tensor Conjugate Gradients (CG) algorithm applied to solve a 3D problem A x1y1z1

Tensor Multigrid algorithms unify algebraic multigrid and geometric multiresolution
In the matrix domain, multigrid methods for solving linear system of equations are known for their fast convergence and computational efficiency. These methods perform iterative improvement of a solution approximation, based on smoothing of the solution error on multiple scales of the problem. Basically, there are two types of multigrid/multiresolution strategies: geometric multiresolution preserves spatial coherence in multidimensions by working on datasets at different resolutions but its formulation is often problem specific. Algebraic matrix multigrid algorithms (AMG) offer a 'black box' approach for solving matrix equations, but they rely on values of matrix coefficients only and risk neglecting data coherence, at a potential cost of degraded performance.
We have used tensor equations to express both approaches in a unified fashion. The proposed tensor multigrid algorithm (TMG) has an important advantage over AMG: in contrast to the matrix algorithm, which neglects multidimensionality and the structure of the problem, TMG unifies both concepts of geometric and algebraic multiresolution, by preserving and exploiting multidimensionality and spatial data coherence (see Fig. 4.3).

Fig. 4.3 shows pseudocode for a Tensor Algebraic
Multigrid V-cycle algorithm applied to solve a 3D tensor problem A x1y1z1 xyz · C xyz = B x1y1z1 . Problems in higher dimensions are handled in a similar way. Fig. 4.4 compares visually how AMG and TMG are applied to a problem of reconstruction of a twodimensional image from a few arbitrarily taken samples, using non-uniform spline interpolation [3]. This example shows that performance of the AMG algorithm is significantly degraded due to data distortion, which is introduced at coarse matrix scales, whilst the TMG converges to a good approximation of the solution within a few iterations, due to its preservation of spatial coherence.
% get preliminary constructed system scale and % corresponding restriction/prolongation operators % separable restriction of the residual

Automatic tensor expression analysis for high-performance implementations
Because the high level tensor formulation contains the complete tensorial structure information of a complex mathematical problem, automatic expression analysis can be performed, resulting in automated software optimization and code generation. We have experimentally verified this aspect in our tensor library and have found that such automatic analysis is capable of optimizing code, with regard to computational and memory requirements, in both a problem-specific and hardware-specific manner. This result profits strongly from the commutativity of all tensor operations the significant differentiating feature over standard matrix formulations which are non-commutative with respect to matrix multiplication. Some aspects of such automated code generation have also been discussed in [4]. This particular topic is a promising area for future research.

Application of the framework to large real-world problems
Our framework was successfully applied to solve a large real-world computational problem -a spline-based global variational reconstruction [3] of a multidimensional signal from incomplete and spatially scattered measurements in four dimensions. The problem was originally formulated in terms of matrix algebra. In four dimensions, the matrices involved in the computations become extremely large, essentially prohibiting computation of the problem on current standard hardware. For example, for a moderate size of the reconstruction grid 128 × 128 × 128 × 16 we have to deal with a matrix which is represented by 7'514'144'528 nonzero elements (in compressed block-band diagonal storage, 30 GBytes in single precision). Using our framework, we reformulated the problem in terms of tensors. The resulting tensor formulation helped to analyze the mathematical structure of the problem and to derive its decomposition in the form of a 1-rank tensor decomposition known as Canonical Decomposition (CANDECOMP) [11,14,1,19]. More details about proposed tensor formulation and solving algorithm can be found in our companion paper. The identified tensorial structure allowed efficient computation of the tensor multiplication which therefore can be efficiently used within a tensor Krylov solver. Another important implication of the identified property of the tensor formulation, is a dramatic reduction of storage requirements. Our approach does not require explicit storage either of the system tensor or its decomposition, and thus the problem of size, discussed above (33'554'432 unknowns parameters), for about 9'000'000 measurements, can be computed on current inexpensive multi core computer equipped with only 2 GBytes of physical memory in less than 60 minutes, surpassing the capability of a published, matrix-based solving algorithm [3] for the same problem by far. The decomposability of the tensor multiplication allows convenient and complete parallelization of computations using all currently available parallelism paradigms, including Single Instruction Multiple Data (SIMD), Multi Core, Clusters of PCs and General Purpose GPU (GPGPU) computing technology.
Note that parallelizability of tensor operations, in general, has proven to be very beneficial: on Multi Core CPUs we have observed a speed increase almost proportional to the number of cores; and on a recent graphics card (ATI Radeon HD 4870 X2) an additional increase has been measured for specific examples, of more than 10, compared to a Quad Core CPU.

Discussion
In this work we introduced a framework for computational solving of large tensor structured problems. The proposed approach leads to a natural, dimensionalitypreserving formulation of tensor structured problems, and by maintaining the multidimensional structure of the data, and the commutativity properties of the framework, directly induces computationally efficient solving strategies that can profit from automatic expression analysis, separability properties of the tensor formulation, and the preservation of spatial coherence of the data, to speed up convergence in tensor solvers as compared to their matrix counterparts.
The chosen formalism based on Einstein notation together with introduced commutativity of tensor multiplication, makes the framework well suited for implementing a tensor programming language, offering self-optimized computations of tensor expressions by semantic analysis of their terms while avoiding explicit use of tensor indexing in the actual implementations of solving programs. This was verified by implementing our own tensor library built on the principles of the proposed framework.
The properties presented above distinguish our framework from one often used in the field of tensor decompositions [27,17,16,11,14,1,19], but do not limit its use for solving problems studied in that field.
The ongoing research shows that the framework is well suited for solving problems of signal processing in higher dimensions, which have inherent tensor structure. We are currently investigating its usefulness in problems from image reconstruction, computer vision and fluid dynamics.
In summary, this high-level approach fits well into the landscape of increasing interest in tensors as a workhorse for elegant and efficient solving of those scientific problems which arise from multidimensional data.   28] Computer Vision -Data compression [13] Data Mining -Independent Component Analysis [6,15] Chemistry -Signal filtering [20] Neuroscience -Modeling of fluorescence data [1] -Social Network Analysis [1] Solving inverse multidimensional problems: Signal Processing -Least squares problems [3,21] Computer Vision -Differential equations (Poisson, Navier-Stokes, Finite Machine Learning Element Methods) [18,9] Fluid Dynamics -Integral equations [8] Chemistry Electromagnetism Modeling of properties of complex systems consisting of many Computational Physics interacting elements expressed by tensor expressions: Quantum Chemistry -Modeling of electronic and optical properties of Computational Chemistry molecules and their interactions [25,22] Material Science

A.1 Dual vectors spaces, contravariant and covariant mechanisms of transformation
Let V be a real vector space with finite number of dimensions N V . We designate by V * unique vector space of same dimensionality that is dual in respect to V .
For a given base ev of V transformations to a new baseêv 1 and vice versa are defined byê where P v v 1 and S v 1 v 2 are direct and indirect base transformations respectively. A vector from V given by Components of such vectors which transform indirectly in respect to transformation of the base are called contravariant.
For the base ev of V and dual base e v 1 of V * the following equation is satisfied where δ v 1 v is Kronecker delta. Thus the dual base e v 1 transforms according tô In contrast to contravariant vector components, components of vectors from dual space V * transform in the same way as base of V Components of such vectors are called covariant.

A.2 Generalized definition of a tensor
Definition A.1 Let ordered set of real vector spaces U 1 , . . . , U P , V 1 , . . . , V Q with finite dimensions I 1 , . . . , I P , J 1 , . . . , J Q and their dual vector space U * 1 , . . . , U * P , V * 1 , . . . , V * Q be given. By identifying (U * i ) * with U i we define every vector u i ∈ U i to be a linear map U * i → ℜ given by for arbitrary u * i ∈ U * i , where , denotes the inner product. With (P + Q) vectors where a * i , b j are arbitrary vectors in U * i and V j respectively. The element T is termed a (P +Q)-th order decomposed tensor, P -times contravariant and Q-times covariant. The space generated by all linear combinations of decomposed tensors is termed the tensor product space. Any arbitrary tensor can be represented as a weighted sum of decomposed tensors [11,1,19,15]. This definition is compatible with, and extends [17,21].

A.3 Tensor notation
To each dimension of the physical space we assign a vector space with a given dimensionality. For example for 3-dimensional physical space we define vector spaces X, Y, Z with sizes N X , N Y , N Z . In such vector spaces one can define multiple coordinate systems related to each other by transformations discussed above. Components of a tensor in some tensor product space are designated as an array with multiple indices. Each index in the designation has direct correspondence to a dimension of the physical space and thus to one of the defined vector spaces. For example for 3-dimensional case a tensor which is an element of X × Y × Z can be designated by T x i y j z k . Subindices i, j, k are used for distinguishing between different coordinate systems.

A.4 Elementwise tensor operations
Linearity of the defined tensor product space implies definition of the following elementwise tensor operations such as tensor addition, subtraction and multiplication by a scalar: With introduced constraint of ordered vector spaces all properties of these operations such as associativity and commutativity remain unchanged.

A.5 Kronecker delta
Kronecker delta presented in expression (A.4) is a second order tensor which is equivalent to the identity transformation. Mixed version of the tensor applied to components of a contravariant or covariant first order tensor (vector) does not introduce changes Covariant or contravariant versions of the Kronecker delta change variance of tensor indices, but do not change values of tensor components δx 1 x · T x = Tx 1 ≡ Tx (A.14) Note that here indices x and x 1 correspond to identical coordinate systems. In the same way Kronecker delta can be applied along some dimension of a higher order tensor δ x 1 x · T xyz ≡ δ y 1 y · T xyz · δ z 1 z · T xyz ≡ T xyz (A.16) The tensor multiplication of Kronecker deltas for different vector spaces forms multidimensional identity transformation δ x 1 y 1 z 1 xyz = δ x 1 x · δ y 1 y · δ z 1 z (A.17) δ x 1 y 1 z 1 xyz · T xyz = T x 1 y 1 z 1 ≡ T xyz (A. 18) A.6 Inner product the inner product of two third order Euclidean tensors A xyz and B xyz is defined by A, B = (A xyz · δxx 1 yy 1 zz 1 ) · B xyz · δ x 1 y 1 z 1 xyz =Â x 1 y 1 z 1 ·Bx 1 y 1 z 1 (A. 19) At the same time using the language of elementwise products the inner product can be expressed as A, B = (A xyz · B xyz ) · δxx 1 yy 1 zz 1 · δ x 1 y 1 z 1 xyz = C x 1 y 1 z 1 · Ix 1 y 1 z 1 (A.20) where Ix 1 y 1 z 1 is a tensor with all components equal to one, which can be easily checked based on the properties of Kronecker delta.
Note that in general case of non-Euclidean spaces the inner product is defined with the use of the metric tensor which for Euclidean case reduces to the Kronecker delta.

A.7 Tensor transposition
In classical settings transposition of a tensor is defined by changing positions of tensor indices. But in our framework the order of indices is unique according to a predefined order of vector spaces. Thus patterns like A T A from matrix normal equations correspond to the tensor multiplication with change of variance by the Kronecker delta in Euclidean spaces or the metric tensor in general case.

A.8 Differentiation in respect to a tensor
When dealing with an optimization problem, differentiation in respect to the optimization parameter may be required. Consider an example of differentiating a tensor valued function of the form f (U ) = A x 1 y 1 z 1 xyz · U xyz ∂ ∂U f (U ) = ∂ ∂U x 2 y 2 z 2 A x 1 y 1 z 1 xyz · U xyz = A x 1 y 1 z 1 xyz · ∂U xyz ∂U x 2 y 2 z 2 (A. 21) where indices (x, y, z) and (x 2 , y 2 , z 2 ) correspond to identical coordinate systems. Using tensor transformational properties it is straightforward to show that ∂U xyz ∂U x 2 y 2 z 2 = δ xyz In a similar way one can compute derivative of more complex scalar and tensor valued functions.