{'fraction_non_alphanumeric': 0.045300098729340564, 'fraction_numerical': 0.013052788142221538, 'mean_word_length': 3.9982055883106895, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 32, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Post-Processing of High-Dimensional Data

Scientific computations or measurements may result in huge volumes of data. Often these can be thought of representing a real-valued function on a high-dimensional domain, and can be conceptually arranged in the format of a tensor of high degree in some truncated or lossy compressed format. We look at some common post-processing tasks which are not obvious in the compressed format, as such huge data sets can not be stored in their entirety, and the value of an element is not readily accessible through simple look-up. The tasks we consider are finding the location of maximum or minimum, or minimum and maximum of a function of the data, or finding the indices of all elements in some interval --- i.e. level sets, the number of elements with a value in such a level set, the probability of an element being in a particular level set, and the mean and variance of the total collection. The algorithms to be described are fixed point iterations of particular functions of the tensor, which will then exhibit the desired result. For this, the data is considered as an element of a high degree tensor space, although in an abstract sense, the algorithms are independent of the representation of the data as a tensor. All that we require is that the data can be considered as an element of an associative, commutative algebra with an inner product. Such an algebra is isomorphic to a commutative sub-algebra of the usual matrix algebra, allowing the use of matrix algorithms to accomplish the mentioned tasks. We allow the actual computational representation to be a lossy compression, and we allow the algebra operations to be performed in an approximate fashion, so as to maintain a high compression level. One such example which we address explicitly is the representation of data as a tensor with compression in the form of a low-rank representation.


Introduction
Many scientific and engineering computations or measurements, as well as economic or financial applications, produce large volumes of data. For simplicity, assume that we have computed or observed large quantities of one real-valued variate. Assume further that the amount of data is so large that it can not be held or stored in its entirety and has to be compressed in some way.
Frequently, the task arises for example to find the location with the maximum value of the data set, or to find the locations of all values which lie in a given interval -we call this a level set. Other, similar tasks we consider are finding the number of values in a given interval, the probability of being in a given interval, or finding the mean or the variance of the data set. These tasks are trivial for small data sets, typically performed by inspecting each datum. But if we consider truly huge amounts of data which can not be stored in full because this would exceed the capacity of any storage device, but only in a compressed manner, and where additionally it may not be possible due to time constraints to inspect each element, these tasks are not trivial any more. This is due to the fact that in the compressed representation the data values are normally not directly accessible and require additional processing. Such compression will in general be "lossy", so that not each value can be restored exactly, but only up to a certain reconstruction error.
Assuming the data as an element of some set T , the algorithms are independent of the representation of the data, as well as from the compression and reconstruction technique used, subject only to the possibility of approximately performing the operations of an Euclidean associative commutative algebra. This means that we assume that T is a vector space with an inner product, and additionally an associative commutative bilinear multiplication, making it into an algebra. The simplest example of such a structure is to envisage the data points as a vector w ∈ R N with an appropriate N ∈ N, reflecting the amount of data. The vector space operations are clear, as is the canonical Euclidean product, and for the commutative multiplication consider the point-wise or Hadamard product, i.e. the component-wise multiplication of two vectors. We shall allow that all these operations are performed only approximately, in order to maintain a high compression level.
Large volumes of data, especially when they can be thought of as samples or discrete values of some real-valued function on a high-dimensional space, can often be arranged in form of a tensor [40,41,8,31]. This offers the possibility to use approximate or compressed tensor representations. Here we will especially show one possibility, namely low-rank representations, which generalise the truncated singular value decomposition for matrices. For the sake of completeness, we shall show the possible implementation of the above mentioned algebraic operations for some of the more common low-rank formats, and the approximated, compressed, or truncated representation of their results; see also [52].
The proposed algorithms are iterative in nature, and the convergence tolerance can be adapted to the reconstruction error. The basic idea for the algorithms, which operate only on the algebraic structure, is an iterative scheme which converges to a result which solves the desired problem in some way. Such iterations typically destroy the compressed representation, so they have to be combined with re-compression or truncation.

An example
Let us give an example which motivates much of the following formulation and development. Assume that we are interested in the time evolution of some system, described by where v(t) is in some Hilbert space V and A(µ) is some parameter dependent operator; in particular A(µ) could be some parameter-dependent differential operator, for example where G ⊂ R is a domain, [0, T ] ⊂ R is the time window of interest, κ(x, µ) is a parameterised random tensor field dependent on a random parameter in some probability space µ ∈ Ω ⊂ R d with probability measure P, and one may take U = L 2 (G) [43,46]. For each µ ∈ Ω one may thus seek for solutions in L 2 ([0, T ], V) ∼ = V ⊗ L 2 ([0, T ]) = V ⊗ Z, where we have set Z := L 2 ([0, T ]). On the other hand, assume that for fixed (x, t) ∈ G × [0, T ] we are looking at the random variable v(x, t, ·) which we for simplicity assume to have finite variance, i.e. v(x, t, ·) ∈ L 2 (Ω, P) =: S. Taking into account the parametric dependence, we are hence looking for a function v(x, t, µ) which is defined on G × [0, T ] × Ω; we are thus looking for a solution in V ⊗ Z ⊗ S. This applies equally well to the abstract Eq. (1), as on the right hand side the operator depends on µ ∈ Ω, so will the solution v(t, µ) to Eq. (1), and it will lie in a similar tensor product.
Observe further that if the probability measure P on Ω in Eq. (1) or Eq. (2) is a product measure P = P 1 ⊗ · · · ⊗ P L on Ω = Ω 1 × · · · × Ω L , the space S can be split further S = L =1 S , with S := L 2 (Ω , P ). Thus in total the solution to both Eq. (1) and Eq. (2) may be viewed as an element of a high order or degree tensor space [43,46,48] If we think of samples of v(x, t, µ) at space points (x k ) K k=1 , time instances (t j ) J j=1 and parameter values (µ i ) i = (µ i 1 , . . . , µ i L ), with a multi-index i = (i 1 , . . . , i L ) ∈ N L with 1 ≤ i ≤ I , the samples of the solution v = (v(x k , t j , µ i 1 , . . . , µ i L )) = (v kji 1 form a high order or degree tensor [31]. In many chemical applications one has n = 100, and d is few hundreds or even thousands [6]. Some other high-dimensional problems from chemistry and physics like Hartree-Fock-, Schrödinger-, or Master-equations, and low-rank tensor methods of their solution are for example considered in [32,37,39,11], see also the references therein. Another example of large volumes of high-dimensional data are satellite data. Satellites collect data over a very large areas (e.g. the data collected by the National Center for Atmospheric Research (USA) [29]. Big data can also come from a computer simulator codes such as a solution of a multi-parametric equation, e.g. weather research and forecasting, and climate models [23]. Also Oil&Gas companies daily collect sensor data from multiple sources. Another source for huge volumes of data are high-energy particle accelerators like CERN [5].
For the sake of simplicity, we will treat all arguments of the function v in Eq. (4) equally, denoting them in that case by p = (x, t, µ) = (p 1 , . . . , p d ). So just consider a real-valued function p → w(p) = w(p 1 , . . . , p d ) and its samples As one may see, a tensor can be simply defined as a high-order matrix or multi-index array, where multi-indices are used instead of indices. As an example, assume that we have N = 10 16 data points. Now a vector w ∈ R N is just a tensor of first order, but one may reshape the data, say into a matrix W ∈ R 10 8 ×10 8 ∼ = R 10 8 ⊗ R 10 8 -a tensor of 2nd order. Reshaping further, one may take R 10 4 ×10 4 ×10 4 ×10 4 ∼ = 4 k=1 R 10 4 -a tensor of 4th order -as well as other combinations, such as w ∈ 16 k=1 R 10 -an array of size 10 × · · · × 10 (16 times) -which is a tensor of 16th order; and finally all the way to w ∈ 16 k=1 (R 2 ⊗ R 5 ) -a tensor of order 32. The tensors obtained in this way contain not only rows and columns, but also slices and fibres [40,41,8,31]. These slices and fibres can be analysed for linear dependencies, super symmetry, or sparsity, and may result in a strong data compression. To have a first glimpse of possible compression techniques, assume that in the above example the data has been stored in the matrix W ∈ R M ×M , where M = 10 8 = √ N , and consider its singular value decomposition (SVD) W = U ΣV T This reduced SVD is a special case of Eq. (6). For the sake of a concrete example, assume that r = 100. One may observe that the required storage has been reduced from N = M 2 = 10 16 to 2 × r × M = 2 × 10 10 for the 2 r vectors {w (1) m , w (2) m } r m=1 . The set of possible objects we want to work with will be denoted by T . We assume that this set is equipped with the operations of an associative and commutative real algebra, and carries an inner product. It is a well known result [60] that such algebras -modulo some technical details -are isomorphic via the Gel'fand representation to a function algebra -to be more precise continuous real valued functions on a compact set. Under point-wise addition and multiplication by scalars, such functions clearly form a vector space, and if one includes point-wise multiplication of functions, they form an associative and commutative real algebra. In our case this is simply the function algebra , which is obviously the same as d =1 R M . The advantage of this abstract algebraic formulation is not only that it applies to any kind of data representation on which one may define these algebraic operations, but that one may use algorithms [35] which have been developed for the algebra of real N × N matrices In our concrete examples we work directly with the data represented as tensors T := d =1 R M . Even without the identification with (M → R), it is obvious that this is naturally a vector space. As also T ∼ = R N , it is clear that T can be equipped with the canonical Euclidean inner product from R N . The associative and commutative algebra product comes from the identification with the function algebra (M → R), i.e. the pointwise product, which is also known as the Hadamard product [31]. Hence it is clear that for data w ∈ R N it is not difficult at all to define the algebraic operations, but rather how to perform them when the data is in a compressed format.
To simplify later notation on the number of operations and amount of storage needed, we will often make the assumption that M 1 = · · · = M d = n, so that the tensor in Eq. (5) represents N = n d values. As already mentioned, our example of data compression is based on low-rank approximations to elements in T . Although low-rank tensor data formats and techniques are almost unavoidable if working with large high-dimensional data sets, we would like to stress that the algorithms presented here are formulated purely in terms of the abstract algebraic structure, and are thus independent of the particular representation.
Whereas a general element w ∈ T hence has N = n d terms, a compressed representation -some low-rank versions of which will be discussed in more detail in Section 3 -will have significantly fewer terms. For example, the CP-decomposition (canonical polyadic) representation, truncated to r terms, If the rank r is reasonably small compared to N = n d independent of n and d, then we have an approximation w r ∈ T with much less storage, which also only depends on the product of n and d and is not exponential in the dimension d. But now the maximum can not be easily looked up; to get a particular element, the expression Eq. (6) has to be evaluated for that index combination. In the following we shall assume that we work with approximations such as in Eq. (6) which need much less storage. One has to make sure then that the algebraic operations of the Hadamard algebra structure on T can be performed efficiently in an approximative manner, so that they have much lower complexity than the obvious O(n d ), although the single elements of w which are usually needed for point-wise operations -and post-processing, see Subsection 1.2 -are not directly available. Thus the motivating factors for applying compression, and in particular low-rank tensor techniques, include the following: • The storage cost is reduced, depending on the tensor format, from O(n d ) to O(drn), • The low-rank tensor approximation is relatively new, but already a well-studied technique with free software libraries available. Various low-rank formats are available.
• The approximation accuracy is fully controlled by the tensor rank. The full rank gives an exact representation.
• Even more complicated operations like the Fourier transform can be performed efficiently in low-rank format. The basic fast Fourier transform on w ∈ T would have complexity O(n d log(n d )). These low-rank techniques can either be combined with the fast Fourier transform giving a complexity of O(drn), or can even be further accelerated at the price of an additional approximation error yielding a superfast Fourier transform [49,10].
On the other hand, general limitations of a compression technique are that • it could be time consuming to compute a compression, or in particular a low-rank tensor decomposition; • with sampling, it requires an axes-parallel mesh; • although many functions have a low-rank representation, in practice only some theoretical estimates exist.
But the fact still remains that there are situations where storage of all items is not feasible, and some kind of compression has to be employed.

Post-processing tasks
This work is about exploiting the compression together with the structure of a commutative algebra and the ability to perform the algebraic operations in the compressed format, at least approximately. The tensor product structure which appears in Eq. (3) is something which allows efficient calculations to be performed on a sample of the solution like Eq. (4) or Eq. (5) with the help of expressions such as Eq. (6) or other compression techniques. This tensor product structure-in this case multiple tensor product structure-is typical for such parametric problems [47]. What is often desired, is a representation which allows for the approximate evaluation of the state of Eq. (1) or Eq. (2) as function of µ ∈ M without actually solving the system again. Furthermore, one would like this representation to be inexpensive to evaluate, and for it to be convenient for certain postprocessing tasks, for example like finding the minimum or maximum value over some or all parameter values.
The tasks we consider here are • finding the location and value of maxima or minima, • finding the location and value of maxima or minima of a function of w , • finding the value in the tensor closest to some given number, • finding level sets, i.e. the indices where the value is between two levels, • finding the number of indices between two levels, • computing the probability of being between two levels, • computing the mean and the variance.
As an example, consider finding the maximum value. A naïve approach to compute the maximum would be to visit and inspect each element, but then the number of visits is O(n d ), exponential in d, which may be unacceptable for high d. Such phenomena when the total complexity/storage cost depends exponentially on the problem dimension have been called the curse of dimensionality [31].
The idea for such iterative post-processing algorithms, in the form of finding the maximum, was first presented in [14] for low-rank tensor approximations. Some further post-processing tasks for such low-rank tensor representations were investigated in [20]. To give an example from [14] of the possible savings in space and time, assume that one looks at the solution of an example of Eq. (1) in Subsection 1.1, or more specifically Eq. (2) in the stationary case with κ(ω, x) ≡ 1, a d-dimensional Poisson equation: Assume further that this is solved numerically by a standard finite-difference method with n = 100 grid-points in each direction, so that the solution vector u has N = n d data points in R N , but can also naturally be viewed as a tensor u of degree d in d k=1 R 100 . So for full storage, one needs N = n d storage locations, and if one is looking for the maximum by inspecting each element, one would have to inspect N elements.
We mention in passing that of course actually solving the discrete equation is a challenge for d > 3, but that is a different story, and for this we refer to [14]  Age of the universe: ≈ 14 × 10 9 years.
Number of hadrons (elementary particles) in the universe ≈ 10 80 . Table 1: Computing times (4th column) on 2 GHz dual-core CPU to find maximum. Now assume for the sake of simplicity that one can inspect 2 × 10 9 elements per second -on an ideal 2 GHz CPU with one inspection per cycle. Then for u ∈ (R n ) ⊗d ∼ = R n d the times needed to find the maximum for the full data-set are shown in Table 1 in the third column-assuming that it were somehow possible to store all the values indicated in the second column-whereas the actual computation with a compressed format is shown in the last and fourth column.
It is obvious that for growing d for the full representation the computational complexity and the storage requirements quickly become not only unacceptable, but totally impossible to satisfy. The second and third column behave like O(n d ) and grow exponentially with d, whereas for a low-rank representation u r of u not only can the data be stored on a modest laptop, but for the simple Algorithm 2 -to be explained later in Subsection 2.4 -the computing times in the fourth column are in terms of seconds and behave like O(n d 3 ).

State of the art
This is not a discussion of the state of the art regarding general tensor formats and their low-rank approximations, in quantum physics also known as tensor networks. For the physical motivations and numerical developments from there see [64,59,22,50,4,3]. For the mathematical and numerical view we refer to the review [41], the monographs [31,37,39], and to the literature survey on low-rank approximations [26]. In the following, we concentrate on the question of post-processing such data.
The idea of finding the largest element and the corresponding location of a tensor by solving an eigenvalue problem, where the matrix and the vectors are tensors given in the CP tensor format, was introduced in [14]. Additionally, the first numerical schemes to compute the point-wise inverse and the sign function were introduced, as well as the rank-truncation procedure for tensors given in the CP format. Later, in [20], these ideas were extended and applied to tensors which were obtained after discretisation and solution of elliptic PDEs with uncertain coefficients. Another group of authors, in [9], by combining the advances of the density matrix renormalisation group and the variational numerical renormalisation group methods, approximated several low-lying eigenpairs of large Hermitian matrices simultaneously in the block version of the TT format via the alternating minimisation of the block Rayleigh quotient sequentially for all TT cores.
In [13,12], the authors suggested methods to compute the mean, the variance, and sensitivity indices in the TT train tensor format with applications to stochastic PDEs, whereas the diagonalisation of large Toeplitz or circulant matrices via combination of the fast Fourier and the CP tensor format was shown in [49].
An investigation of approximations to eigenfunctions of a certain class of elliptic operators in R d by finite sums of products of functions with separated variables is the topic of [33]. Various tensor formats were used for a new class of rank-truncated iterative eigensolvers. The authors were able to reduce the computational cost from O(n d ) to O(n). They demonstrated the introduced algorithms by solving large-scale spectral problems from quantum chemistry: the Schrödinger, the Hartree-Fock, and the Kohn-Sham equations in electronic structure calculations.

Outline of the paper
In the following Section 2 the necessary material for abstract algebras is quickly reviewed. Then the algorithms and functions on the algebra T used to compute the post-processing tasks outlined in Subsection 1.2 are formulated in an abstract fashion, independent of the representation chosen for the data. Even as the formulation uses only the abstract algebra operations, we point out and motivate what this means in terms of the Hadamard algebra. But also for the Hadamard algebra, the presentation is independent of the tensor format to be chosen. Furthermore, for the iterative algorithm -a fixed point iteration -it is discussed how the possible truncation operation influences the convergence.
In Section 3 we present some concrete examples of compression of high-dimensional data w ∈ R N , once it is identified with a tensor w ∈ T : For such tensors of degree d we discuss the canonical polyadic (CP) representation, the Tucker representation, and the tensor-train (TT) representation, and the compression in terms of low-rank approximations based on the different tensor formats. In particular, we show how the algebra operations can be carried out in the different tensor formats, the numerical effort involved, and the effect these algebraic operations have on the compression or low-rank representation. As the compression level may deteriorate, i.e. the rank of the approximation may grow, it is important that one is able to re-compress. Pointers to such methods in the literature are included as well.
The Section 4 contains a few numerical examples and a discussion of their results used to illustrate the algorithms of Section 2 and the representations from Section 3. These examples are solutions of elliptic high-dimensional or parametric resp. stochastic partial differential equations, a domain where the data is naturally in a tensor format. But as already pointed out, for any data it is just a -often only conceptualreshape to consider it as an element of a tensor space. The conclusion is contained in Section 5.

Algorithms for post-processing
Given a data-set w in some compressed tensor format, where looking up each element is not trivial, and where additionally it may not be computationally feasible to look at each element, we still want to be able to perform tasks which essentially require looking at each datum. The tasks outlined in Subsection 1.
6. the mean (sum of all items) and variance (sum of all items squared), 7. some auxiliary functions, such as the algebraic inverse w −1 of w -in our case the Hadamard inverse -and others like the sign(w ) function which will be seen to be needed for the above tasks in Subsection 2.3.
The above tasks and auxiliary functions will be computed by finding an iteration (mapping) Φ P for each post-processing task P , such that the fixed-point v * of Φ P -i.e. Φ P (v * ) = v * -is either the sought solution for the task, or computes one of the auxiliary functions.

Iteration with Truncation
When performing computations using the operations of the algebra -like the Hadamard algebraic operations on tensors in some compressed format -after the operation the compression will typically be sub-optimal, which means that the result has to be compressed again. Thus, when performing the algebraic operations for a fixed-point iteration, the compression would either get worse and worse in each iteration, or one has to re-compress or truncate the result again to a good compression level. Therefore a re-compression after each or after a number of algebraic operations may be necessary, and we thus allow that the algebraic operations are possibly only executed approximately. In our example case the low-rank representation of tensors explained in Section 3 acts as compression, and the rank may increase through the algebraic operations, and hence we will use truncation T to low rank r with error [31,2] to avoid the problem of a deteriorating compression level.
In other words, with a general compressed representation w r of w the computation will be a truncated or perturbed iteration [34,48]. If we denote the general compression mapping by T -meaning compression with an accuracy -the iteration map is changed The general structure of the iterative algorithms for a post-processing task P or for an auxiliary function is shown in Algorithm 1. Here we want to collect some results for this kind of iteration. For such a truncated or perturbed iteration as in Algorithm 1, it is Algorithm 1 Iteration with truncation 1: Start with some initial compressed guess v 0 depending on task P .
the iterator Φ P may deteriorate the compression level 5: if compression level of z i is too bad then 6: v i+1 ← T (z i ); use truncation T to compress z i with error 7: end if 10: i ← i + 1 11: end while known that 1. if the iteration by Φ P is super-linearly convergent, the truncated iteration T • Φ P will still converge super-linearly, but finally stagnate in an -neighbourhood of the fixed point v * [34]. One could loosely say that the super-linear convergence of Φ P is stronger than the truncation by T .
2. if the iteration by Φ P is linearly convergent with contraction factor q, the truncated iteration T • Φ P will still converge linearly, but finally stagnate in an /(1 − q)neighbourhood of v * [48]. Again, one could loosely say that iteration by Φ P and truncation by T balance each other, thus resulting in a larger neighbourhood of stagnation.
We shall assume that the truncation level has thus been chosen according to the desired re-construction accuracy and taking into account the possible influence due to the convergence behaviour of the iterator Φ P .

Preliminaries and basic algebraic operations
As already pointed out, we assume in general that the set T is a vector space. Our example of the space T = d =1 R M of tensors of interest was already introduced in Subsection 1.1. This is clearly a vector space, so we are able to add two such tensors, and multiply each by some real number, in other words we may form linear combinations. The additive neutral element -the zero tensor -will be denoted as . It is important that the compressed storage format chosen allows to perform the vector space operations without going into the full representation. This will be shown for some of the familiar tensor formats in Section 3. Furthermore, as T : as vector spaces, we can carry the canonical Euclidean inner product on R N to T , which for u, v ∈ T is denoted by This makes T into a Euclidean or Hilbert space. We assume that the computation of this inner product is feasible when both u and v are in a compressed representation; again this will be demonstrated in Section 3 for the low-rank representations we shall consider as examples.
On this space T an additional structure is needed, namely a multiplication which will make it into a unital associative and commutative algebra. In our example this will be the Hadamard multiplication, which apparently goes back to Schur, which is the pointor component-wise multiplication on d It is a bi-linear operation (linear in each entry), and it is commutative. Thus this product makes T into an associative and commutative algebra. The symbol will be used both for the abstract algebra product on T , as well as for the Hadamard product on d =1 R M . As is any such algebra, for any w ∈ T it holds that w  = .
Especially, as a function of v , the product u v is a linear map L u ∈ L (T ) on T . This is the familiar canonical representation T u → L u ∈ L (T ) of an associative algebra as linear maps on itself, i.e. in this case in a commutative sub-algebra of the algebra of all linear maps L (T ) with concatenation as product.
It easy to see that this Hadamard algebra has a unit element  ∈ d =1 R M for the product, a tensor with all entries equal to unity -this makes d =1 R M into a unital algebra. Observe that by defining for all 1 ≤ k ≤ d the all-ones vectors 1 M k = [1, . . . , 1] T ∈ R M k , the Hadamard unit has rank one according to Eq. (6):  = d k=1 1 M k . Obviously one has L  = I T , the identity in L (T ).
Having defined a unit or neutral element for multiplication, we say that w ∈ T has a multiplicative inverse iff there is an element, denoted by w −1 ∈ T , such that Obviously not all elements have an inverse, e.g. the zero element  ∈ T is never invertible. In any case, not every w ∈ d =1 R M has a Hadamard inverse, for this it is necessary that all entries w m 1 ...m d = 0 are non-zero, and then Elements of the algebra which can be written as a square w = u u = u 2 are called positive, they form a convex cone T + ⊂ T ; note that obviously the zero tensor  and the multiplicative unit  are positive, and that if an invertible w is positive, so is its inverse w −1 . As usual, this defines an order relation u ≤ v ⇔ v − u ∈ T + , which implies that for any positive w ∈ T + one has  ≤ w . Further one may observe that A certain compatibility of the inner product and the algebra product is needed, as we require that i.e. that the action of the algebra product is self-adjoint, or, in other words, that the maps L w are self adjoint. This condition is satisfied for the Hadamard algebra.
The inner product with the unit element φ(w ) := w | T defines, as function of w , a positive linear functional φ, as for w ∈ T + one has It is a kind of "trace" or un-normalised state functional on the algebra, as φ() = N = dim T in the Hadamard algebra. In particular Conversely, if the algebra comes equipped with such an un-normalised positive state functional, the Eq. (11) may be taken as the definition of an inner product, which automatically satisfies Eq. (10).
We shall assume that we can compute the algebraic operations -at least approximately -in compressed representation (this will be shown for the Hadamard algebra for the low-rank formats in Section 3), and that we can compute the multiplicative inverse also in compressed representation -again at least approximately.
As in any unital algebra, if for λ ∈ C the element w − λ fails to be invertible -this means that L w − λI T is not invertible in L (T ) -we shall say that λ is in the spectrum of w ; here -as we are in finite dimensional spaces -it is an eigenvalue, i.e. there is an , there must be a multi-index m λ ∈ M such that w m λ = λ, and hence in the real Hadamard algebra T one has that always λ ∈ R. Thus the spectrum of w resp. L w is Observe that according to Eq. (10), each representing map L w from above is self-adjoint, which means that also in general all the spectra are realσ(w ) ⊂ R.
Specifically, the Euclidean Hadamard algebra T = d =1 R M which we have constructed is obviously isomorphic -as a Euclidean algebra -to R N with the canonical Euclidean inner product when equipped with the point-or element-wise Hadamard product N . Let us denote this isomorphism by V : T → R N ; it implies some ordering of the terms of each w ∈ T .
More enlightening and less obvious may be the unital algebra isomorphism with the commutative sub-algebra of diagonal matrices diag(R, N ) in the full matrix algebra with the usual matrix multiplication. Let us denote this isomorphism by M : T → diag(R, N ). If w = V (w ) ∈ R N is the vector containing all the elements of the tensor w , then M(w ) : Hence, by choosing the canonical Euclidean basis in R N and its image by V −1 as a basis in T , the canonical representation Let v ∈ T be another tensor, with associated N ), then from the definition of the isomorphism M: which contains the important equality from where one sees that which helps in understanding the following algorithms. Note that this shows that the algebra d =1 R M with the Hadamard product is already jointly diagonalised and is essentially in its Gel'fand representation [60], which is an abstract way of saying what was already stated above, namely that the spectrum σ(w ) of an element w ∈ d =1 R M are exactly the individual terms or data stored in the tensor. If in some abstract unital associative and commutative algebra the multiplication is not the point-wise multiplication -e.g. think of convolution -then it is known that modulo some technicalities [60] it is isomorphic via the Gel'fand "diagonalisation" morphism to a function algebra, and all the algorithms to follow would deal with the spectrum σ(w ) of w instead of with the "values of w ". These two notions only coincide for a function algebra.

Post-processing Tasks
Here the different post-processing tasks are explained together with how they will be computed. This may involve a number of auxiliary functions. The computation of these -again through truncated iteration Algorithm 1 -and the way in which some of the computations can be enhanced or accelerated is shown in Subsection 2.4.
Finding the maximum or minimum of w ∈ T = d =1 R M is the first task we consider. Observe, that the element of maximum modulus of the tensor w is also equal to the ∞-norm w ∞ . Finding are evidently of rank one. This shows that each datum w m of w may be found through eigenvalue computation. The eigenvalue is the data entry w m with the corresponding eigenvector e (m) indicating the location resp. the index m ∈ M.
Although quite obvious, it may be worthwhile pointing out that the linear span of any nonempty collection of eigen-resp. unit vectors T N := span{e (n) | n ∈ N } (N ⊆ M) is an invariant subspace, and the orthogonal projector on T N is L p (N ) := n∈N L e (n) with p (N ) = n∈N e (n) . With this, the spectral resolution of the identity for any and all L w ∈ L (T ) for w ∈ T is I T = m∈M L e (m) corresponding to  = m∈M e (m) , giving the spectral resolution (17) which may be seen as the trivial basis representation of w .
As the basic algorithms for eigenvalue computation like power iteration for self-adjoint linear maps converge to the eigenvalue of maximum absolute value or modulus, it is now clear that if we assume that the maximum is also the element of maximum absolute value, we only have to perform power iteration. This gives us the valueŵ = wm, and the location through its corresponding eigenvector e (m) . How to find the minimum in this case, or what to do when the element of maximum modulus is the minimum will be explained later in Subsection 2.4; it will all be accomplished with spectral shifts familiar from eigenvalue computations.
Finding the maximum or minimum of some function f of w , where it is assumed that f : σ(w ) → R, is not very difficult if one can compute f (w ), which will be simply understood as f (w ) := (f (w m )) m∈M . This definition of f (w ) is also the one which comes via the isomorphy with the algebra of diagonal matrices diag(R, N ) ⊂ gl(R, N ), where f (w ) is defined [35] through the corresponding matrix function f (M(w )) = f (W ), and is also the same as the one which comes from the abstract functional calculus on the Banach algebra T [60].
Observe that one may compute f (w ) for any real or complex function f defined on the spectrum σ(w ) -which is a finite set -at least in principle, through an interpolating polynomial p(w) = k α k w k . As usual, the number w ∈ σ(w ) is replaced by w (or L w ) in the polynomial, resulting in f (w ) = p(w ) = k α k w k . This may be highly inefficient, as in general the degree of the polynomial would be equal to the number of data, in this case the huge number N − 1. For the functions which will be needed for the post-processing tasks, the computation of the function will be done differently, namely iteratively via Algorithm 1.
Assume now that f (w ) has been computed, at least approximately. To find maxima or minima of f (w ), one now has to simply apply the considerations of the preceding paragraph to the tensor f (w ). Obviously, f (w ) has the same spectral resolution as w in Eq. (17): M . (18) Finding the index and value closest to a given number ρ ∈ R is now simply finding the eigenvector and eigenvalue of This is a special case of the preceding paragraph for the function f : t → (t − ρ) −1 . Therefore, if λ w is any eigenvalue of L w , then the corresponding eigenvalue of L y is the transformed one λ y = (λ w − ρ) −1 . Hence the value of w closest to ρ is the eigenvalue of maximum modulus of L y , and the element of largest magnitude of (w − ρ ) −1 . To find the element of smallest magnitude, or even the vanishing elements of w , one would use ρ = 0. Thus one may use the same algorithms as in the maximum search above. Observe that this operation here requires the Hadamard inverse, which will be one of the auxiliary functions.
Finding the indices in a level set requires the so-called sign function, where sign : T → T is defined component-wise as Again this is a special case of the previous one with a general function f .
The characteristic function for a subset S ⊂ R -we shall only look at intervalsχ S : T → T , is again defined component-wise: With these two auxiliary functions, it is possible to define the characteristic function of a level set, i.e. all values between ω 1 , ω 2 ∈ R. We then have with −∞ < ω 1 < ω 2 < ∞: of a characteristic function χ S as the subset of those indices where it is non-zero. Its cardinality is the number of non-zero positions in the characteristic function; requiring only the computation of one inner product with the Hadamard multiplicative unit.
Computing the probability of being in a level, as well as the mean and the variance is under the assumption that each index in w carries the same probability, as then the probability of a value of w being in a subset S is simply The mean or average and variance of w is then Given the level set function from above, one may also compute the conditional mean, conditioned on being in the set S:

Auxiliary functions and algorithmic details
Here we show how to compute the auxiliary functions used in the previous Subsection 2. 3. These are functions defined on the Euclidean Hadamard algebra T = d =1 R M , and it is worth while pointing out that the iterative algorithms will be operating on the algebra, which, as we saw in Subsection 2.2, is isomorphic to an algebra of diagonal matrices, diag(R, N ) ⊂ gl(R, N ), and these matrices are real symmetric. This means in particular that all the known algorithms for computing matrix functions [35] can be used on the Hadamard algebra as well. Of those needed here, the most basic one and the one with the simplest connection to the algebra turns out to be the inverse w −1 . Let us remark once again that the algorithms are valid in any abstract algebra, but we are mainly concerned with the example of the "function algebra" d =1 R M .
The Hadamard inverse w −1 of a Hadamard invertible w ∈ T is needed for inverse iteration and other post-processing tasks. Although the Hadamard inverse is the application of the function f : t → t −1 to w and could in principle be computed through a polynomial, it is often advantageous and much more effective to use other algorithms.
The algorithm for the inverse can be given simply by referring to the quadratically convergent Newton algorithm for matrices [35] for computing inverses, i.e. we use Newton's method for (P ← ( − 1)) and apply it to the equation Hence the iteration function to be used in Algorithm 1 with starting vector v 0 := w is It is well known that the iteration converges quadratically [35], and according to what was explained in Subsection 2.1 this is true even with truncation.
The sign function is defined in Eq. (20). Many of the tasks explained in the previous Subsection 2.3 involve the sign function. Given the sign function, the characteristic, support, and level set functions are easily computed by the basic operations of the algebra. From the fact that each value in a tensor w is an eigenvalue in the Hadamard algebra, and from the isomorphy with the algebra of diagonal matrices diag(R, N ) ⊂ gl(R, N ) where each value on the diagonal is also obviously an eigenvalue, one sees that this is actually the same definition as the one for matrices [35], and is also the one which follows from general functional calculus in abstract algebras [60]. This means that sign(w ) is the application of the function f : t → sign(t) to w ∈ T , and again could in principle be computed with a polynomial. To compute it via the iterative Algorithm 1, i.e. (P ← sign), one uses the same algorithm [35] -the Roberts-Newton algorithm -as for matrices. It can be derived by applying Newton's method to the equation F (v ) := v v −  =  with starting value v 0 := w . This yields the iteration function Observe that this iteration function in a slightly more general form Φ ( √ w ) (v ) := 1 2 (v + v −1 w ) is the ancient Babylonian method [35] -in modern paralance Newton's method -to find the square root w 1/2 of w , as by inserting w =  in the Babylonian iteration Φ ( √ w ) , one obtains Eq. (29). This means that one is iterating to compute the square root of the unit element  with the specific starting value v 0 = w . We will come back to this point of view shortly.
It is known [35] that the iteration converges quadratically, and according to what was explained in Subsection 2.1, it does so even with truncation; observe that it needs the Hadamard inverse from the previous paragraph. Therefore it is not really practical as this would mean that it can only be applied to invertible w ∈ T = d =1 R M , but it also actually needs a new inverse v −1 i in each iteration. This last fact would computationally result in a nested iteration -using the algorithm from the previous paragraph in an inner iteration to compute the inverse v −1 i in each sweep of the outer iteration for the sign function -and such procedures are seldom computationally efficient.
It is thus simpler to use a device employed also for matrix sign computations [35], namely to replace the explicit inverse in Eq. (29) by one step of the iteration for the inverse in Eq. (28). The resulting Newton-Schulz algorithm [35] has the same starting point v 0 := w as before, but uses the iteration function in the Algorithm 1. It is now a simple iteration which still converges quadratically [35], even with truncation according to the explanations in Subsection 2.1. We refer to the discussion in [35] on how to get even faster algorithms using Padé approximations and other ways to accelerate the convergence through scaling, which is especially important in the initial stages.
Coming back to the somewhat curious idea of iterating for the square root of , one may observe that although the defining equation But this checks only how well the equation is satisfied, and not directly how accurate the iterate v i is. Regarding this latter issue, specifically for quadratic convergence, a natural criterion [35] to check for the accuracy of v i at step i is Further, referring to the discussion specifically for the sign function in [35], there are arguments to check for the exponents p = 0, 1, 2. For p = 0, this is the original general criterion Eq. (32), whereas the other values of p take specific consideration of the sign function.
Eigenvalue computations are involved in the first three tasks described in Subsection 2.3. One may use any algorithm developed for large scale matrices [24,57,58,65] which only uses the action of the matrix on a vector in the computation -this is the Hadamard product in our case -and we shall here only explain the main idea and a few variations. The simplest algorithm is power iteration (P ← pow-it), and the iteration map Φ pow-it is given in Algorithm 2, to be used in Algorithm 1. Assume that the datum wm > 0 of w with maximum absolute value or maximum modulus is indeed a (positive) maximum. Assume also that this occurs at one unique indexm, i.e. the datum λ 1 := wm is a simple eigenvalue, and denote the next by absolute value smaller element / eigenvalue by λ 2 . Assume the further eigenvalues ordered by decreasing absolute value.

Algorithm 2 One step power iteration
normalise output z to unit length. 5: Output z; As is well known [24,57], the generated sequence v i converges linearly to ±e (m) with a contraction factor of q = |λ 2 /λ 1 |, and with what was said in the explanations in Subsection 2.1, it still converges with truncation and will stagnate in the vicinity of ±e (m) . Also, as in the algorithm the input v i has unit length, the length of u -given by the inverse γ −1 of the scaling factor in line 3 -converges linearly to the desired λ 1 = wm with the same rate. Algorithmically, the eigenvalue approximation λ 1 ≈ γ −1 is a side-effect of Algorithm 2.
This can be immediately enhanced through the additional computation of the Rayleigh quotient (RQ) between line 2 and line 3 of Algorithm 2: w (v i ) = u|v i T -no need to divide by the length of the unit vector v i . The Rayleigh quotient is stationary at an eigenvalue and has thus a quadratic convergence, it usually represents a much better approximation to λ 1 = wm than γ −1 ; see Algorithm 3.
Stopping criteria for eigenvalues could be done similarly to normal convergent processes, but is is possible to have a posteriori error estimates which are actual bounds specifically for eigenvalues and -vectors, e.g. see [44,24,57]. The simplest seems to be the so-called Krylov-Bogolyubov bound. With an approximate eigenvalue µ and approximate eigenvector x of L w it is given by [44] min The right-hand side of Eq. (35) is minimised by µ = w (x). A short computation [44] shows that for the substitutions x ← v i and µ ← w (v i ) the right-hand side of Eq. (35) becomes Observe that this is not an a posteriori error estimate, but an actual bound. It is the best possible with the available information [44], and it can be very easily computed inside the iteration to control the possible termination. Similar bounds exist for the approximate eigenvector v i . Inserting these considerations into Algorithm 2 gives Algorithm 3 for the task power iteration with RQ computation (P ← pow-RQ) for the iteration map Φ pow-RQ to be used in Algorithm 1, which has the more accurate RQ eigenvalue approximation λ 1 ≈ w and error bound ε λ = ( u|u T − u|v i 2 T ) 1/2 included as side-effects.

Algorithm 3
One step power iteration with Rayleigh quotient (RQ): Φ pow-RQ inverse length of u. 7: z ← γ · u; normalise output z to unit length. 8: Output z; Once the determination of the eigenvalue is accurate enough such that one has an interval which contains only one eigenvalue, one may use the even tighter Temple-Kato bounds [44].  [27], which will be discussed later.

Starting vectors and other enhancements
Another technique in eigenvalue computations is deflation [24,57]: when one eigenvalue λ 1 = λm of largest modulus and corresponding eigenvector e (m) have been located, one may want to compute a further locationm where there is an element wm = λm of equal magnitude |wm| = |λ 1 | = |λm|. For this one can use deflation, and one could change line 2 in Algorithm 2 or Algorithm 3 to u ← w ( − e (m) ) v i to obtain a new iteration map for deflation. The factor ( − e (m) ) projects into the orthogonal complement of the invariant subspace span{e (m) }, and thus all the eigenvalues are left unchanged, except for λm, which is mapped to zero. But in this special case -where we know the form of all invariant subspaces -it is even simpler to keep the same iteration map and to choose as starting vector v 0 := ( − e (m) )/ (N − 1). It has a zero at positionm and is thus in the invariant subspace (span{e (m) }) ⊥ , and all iterates will stay in that subspace. Iterating in Algorithm 1 with that starting vector either with the iteration map in Algorithm 2 or in Algorithm 3 would give us either another eigenvalue λm of equal magnitude, or the eigenvalue with second largest magnitude λ 2 , and the corresponding eigenvector. In this manner all desired eigenvalues can be computed via deflation.
Some other acceleration techniques should be mentioned briefly: The convergence speed in power iteration is controlled by the ratio |λ 2 /λ 1 |, where λ 2 is the next smallest eigenvalue in absolute size. Sometimes this ratio can be very close to unity. If instead with a single vector v i one iterates with a whole block v (1) i , . . . , v (j) i of j mutually orthogonal vectors -effectively a subspace -the convergence speed changes to |λ 1+j /λ 1 |. One has to restore orthogonality though after each iteration sweep [57] of this block-or subspace iteration. This kind of technique also makes it possible to compute multiple eigenvalues, i.e. when the maximum occurs at several places.
Other well-known methods for eigenvalues with even faster convergence build on Krylov subspaces [57], here we have symmetric matrices and thus one would use the Lanczos method. These procedures rely on orthogonalisation though, and this may be problematic when combined with truncation. And certainly the Krylov subspace can be combined with the block iteration idea to give block-or subspace-Lanczos methods [45].
"Exponentiating" the power iteration can be easily achieved via the clever idea of [27]. Recall that with the starting vector v 0 := / √ N , power iteration would compute v 1 ∝ w 1 = w 1 , a scaled version of w , corresponding to the action of w 1 on . Starting actually with v 0 ∝ w and changing line 2 in Algorithm 2 or Algorithm 3 to u := v 2 i , one has v 1 ∝ w 2 = w 2 , i.e. the action of w 2 on . In the next iteration one has v 2 ∝ w 4 , and in iteration i one has v i ∝ w 2 i . Reformulating Algorithm 3 for this new task of "exponentiated" power iteration (P ← exp-pow) for the new iteration function Φ exp-pow results in Algorithm 4, to be used in Algorithm 1 with starting vector v 0 := w / w T . We also need to keep the auxiliary vector y = w 2 . Algorithm 4 One step "exponentiated" power iteration with RQ: Φ exp-pow Rayleigh quotient (RQ) λ 1 ≈ 1 4: 2 ← y |u T ; RQ of w 2 . inverse length of u. 7: z ← γ · u; normalise output z to unit length. 8: Output z; Thus the eigenvalue approximation is through the Rayleigh quotient 1 = w |u T = w |v 2 i T = w v i |v i T in line 3 of Algorithm 4. As y = w 2 , on line 4 the quantity 2 = y |u T = w 2 |v 2 i T = w 2 v i |v i T computes the RQ of w 2 . The "eigenvector" v i ∝ w 2 i converges with the rate |λ 2 /λ 1 | 2 i , and hence one achieves exponential convergence. This acceleration technique can be used in all eigenvalue computations described here.
Transformed eigenvalue computations are a well-known device to find e.g. the maximum of the function y = f (w ). After computing f (w ), one proceeds as before, but with y instead of w . Some functions f (t) are simple enough so that the computation of f (w ) is combined with the iteration. Recall that this transforms all eigenvalues λ w to λ y = f (λ w ) according to spectral calculus, leaving the eigenvectors unchanged.
Shifting and inverse shifting are functions needed to access other than the point of largest modulus in the spectrum of L w resp. w . These functions are the shift f : t → t + β and the inverse shift f : t → (t − ρ) −1 .
The simple function of shifting f (t) = t+β is needed if the actual element of maximum modulus is a minimum of w ; then this is a maximum of −w . In the iteration this will be picked up by the Rayleigh-quotient being negative: w < 0. Also note that for this case -where the maximum has smaller absolute value than the minimum -one may shift the tensor by a value of β = − w ≈ −λ 1 toŵ := w + β ·  so that the maximum ofŵ has larger absolute value than the minimum ofŵ . One then iterates withŵ . After determining the maximum ofŵ , we may subtract the number β > 0 again, to obtain the maximum of w . The index of the maximum ofŵ is of course the same as the one for the maximum of w . Completely analogous manipulations can be performed to find a minimum which is not of maximum modulus. These and similar techniques are well known from eigenvalue calculations of large / sparse symmetric matrices [24,57,58,65].
The inverse shift f (t) = (t − ρ) −1 is needed to access intermediate points in the spectrum σ(w ). Here one finds the data in w closest to ρ, as they have the maximum modulus under the transformation. Thus one computes first y := (w − ρ · ) −1 approximately through Algorithm 1 with the iteration function Eq. (28) -withw := w − ρ ·  instead of w -followed by the eigenvector computation with iteration function as in Algorithm 4, but with starting vector v 0 := y / y T .

Tensor formats
In this section we review definitions and properties of frequently used tensor formats. Many such formats are used in quantum physics under the name tensor networks, see [64,59,22,50,4,3]. We only look at the canonical polyadic (CP), the Tucker, and the Tensor Train (TT) format. The CP [36] and Tucker [62] formats have been well known for a long time and are therefore very popular. The TT format was originally developed in quantum physics and chemistry as "matrix product states" (MPS), see [64] and references therein, and rediscovered in [55,56] as tensor train.
A new class of tensor formats, which we do not consider here, is the hierarchical tensor (HT) format. It was introduced in [30], and further considered in [28].
We note that the sets of low-rank tensors of fixed rank in the CP format are not closed for d > 2, whereas in the Tucker, TT, or HT tensor formats these sets are closed. Therefore the minimisation problem Eq. (39) for the best approximation in a low-rank format has a solution, which can be computed by, for example, an appropriately modified Newton method. The computations in the TT and HT formats are based on the singular value decomposition (SVD) [53,28], and in the Tucker format on the higher order SVD (HOSVD) [63]. In all tensor formats the tensor rank doubles for addition and squares for the Hadamard product. To avoid unnecessary and harmful rank growth, the rank is usually truncated, e.g. by ALS-like or other optimisation algorithms [14,21]. A tensor format is described by a parameter vector space P = × d ν=1 P ν , where P ν = R dν , and a multilinear map U : P → T into the tensor space T := d =1 R M . For practical implementations of high dimensional problems we need to distinguish between a tensor w ∈ T and its tensor format representation P ∈ P, where w = U(P ). There are many possibilities to define tensor formats. Here, we consider the canonical (CP), the Tucker, and the tensor train (TT) format. In the following we briefly repeat definitions and properties of these formats [31]. We also note that the reader can invent his own tensor format, which especially well fits to his needs.
The CP format is cheap, it is simpler than the Tucker or TT format, but, compared to others, there are no reliable algorithms to compute CP decompositions for d > 2 [31,39]. The Tucker format has stable algorithms [38], but the storage and complexity costs are O(d r n + r d ), i.e. they grow exponentially with d. The TT format is a bit more complicated, but does not have this disadvantage. CP and Tucker rank-structured tensor formats have been applied in chemometrics and in signal processing [61,7].

The canonical polyadic tensor format
The canonical representation of multivariate functions [36] was introduced in 1927. A Figure 1: Schema of the CP tensor decomposition of a 3D tensor. schema of the CP tensor format for d = 3 is shown in Figure 1. The full tensor w ∈ R n 1 ×n 2 ×n 3 is shown on the left side of the equality sign, and its decomposition on the right. The lines denote the vectors w i1 , w i2 , and w i3 for i = 1, . . . , r, respectively.

Definition 3.1 (Canonical Polyadic (CP) Tensor Format). The canonical polyadic tensor format in T for variable r is defined by the multilinear mapping
P CP,r P := (w (ν) i We call the sum of elementary tensors w := U CP,r (P ) a tensor represented in the canonical tensor format with r terms. The system of vectors P = (w (ν) i ∈ R Mν : 1 ≤ i ≤ r, 1 ≤ ν ≤ d) is a representation system of w with representation rank r. One may think of P as of a vector valued r × d matrix with the vector w (ν) i ∈ R Mν at index position (i, ν).
The storage requirement for w = U CP,r (P ) is r × d ν=1 M ν , and in the simple case M 1 = · · · = M d = n it is O(r d n).

Basic operations with the canonical format
We denote the set of all tensors w ∈ T of rank r by T r . The set T r is a cone, i.e. w ∈ T r implies α · w ∈ T r for α ∈ R, and at the same time T r is not a vector space as for w 1 , w 2 ∈ T r one has in general w 1 + w 2 / ∈ T r , but w 1 + w 2 ∈ T 2r [21]. A complete description of fundamental operations in the canonical tensor format and a their numerical cost can be found in [31]. For recent algorithms in the canonical tensor format we refer to [14,15,16,17].
Multiplication by a scalar α ∈ R could be done for a w = U CP,r (P ) as in Eq. (37) by multiplying all of the vectors {w (ν) j , i = 1, . . . , r} for any ν by α. But to spread the effect equally and keep the vectors balanced in size, we recommend to define α ν := d |α| for all ν > 1, and α 1 := sign(α) d |α|. Then, with a computational cost of O(r n d), and without changing the rank, The sum of two tensors in the CP format w = u + v can be written as follows The result generally has rank r u + r v , and as this operation requires only concatenation of memory it has only a computing cost of O(1).
The Hadamard product w = u v can be written as follows The new rank is generally r u × r v , and the computational cost is O(r u r v n d) arithmetic operations.
The Euclidean inner product is computed as follows: The computational cost of the inner product is O(r u r v n d).
We note that in all operations the numerical cost grows only linearly with d, but the representation rank of the resulting tensors may increase. Therefore, a rank truncation procedure is needed, which approximates a given tensor represented in the canonical format with lower rank tensors up to a given accuracy.

Rank truncation in the CP format
Let w be a tensor of rank R. Truncating w to a new rank r < R is a fundamental problem [14]. This problem can be formulated as follows: find a w * with rank r such that ∀u with rank r : w − w * ≤ w − u .
Typical methods to solve this problem are the ALS-method and the Gauss-Newtonmethod. The ALS method may show slow convergence, see [14,21] and references therein. The Gauss-Newton method for d ≥ 3 requires some additional assumptions, and also may not show any convergence at all [14].
It is known, see for instance pp. 91-92 in [39], that the class of rank-r CP tensors is a non-closed set in the corresponding tensor product space for d > 2. Therefore there is no w * as above, and one may look for an ε-solution, i.e. minimising within a deviation of ε from the infimum: Definition 3.2 (Approximation Problem). For a given tensor w in CP format with rank R and ε > 0 we are looking for minimal r ε < R and a tensor w * of rank r ε in CP format, such that: This problem is discussed in [14,17].

The Tucker tensor format
The Tucker tensor format was introduced in [62]. A schema for the Tucker tensor format Figure 2: Schema of the Tucker decomposition of a 3D tensor.
for d = 3 is shown in Figure 2. The full tensor w ∈ R n 1 ×n 2 ×n 3 is shown on the left side of the equality sign, and its decomposition on the right. The small cube on the right side in the Tucker format denotes a core-tensorŵ of size r 1 × r 2 × r 3 , and the rectangles denote matrices of column vectors U 1 := [w 1 , . . . , w (1) 1 , . . . , w (2) r 2 ] ∈ R n 2 ×r 2 , U 3 := [w 1 , . . . , w (3) r 3 ] ∈ R n 3 ×r 3 , In this way one has where the symbol × k denotes the contraction of the k-th tensor index ofŵ with the first index of the matrix U T k in Figure 2.  Originally, it was applied for tensor decompositions of multidimensional arrays in chemometrics. The Tucker tensor format provides a stable algorithm for decomposition of full-size tensors. The higher order singular value decomposition (HOSVD) and the Tucker ALS algorithm for orthogonal Tucker approximation of higher order tensors were introduced in [8]. The storage cost for the Tucker tensor is bounded by O(d r n + r d ), with r := max r .

Basic operations with the Tucker format
Let us assume w = U T,r (P ) in Tucker format as in Eq. (38).
Multiplication by a scalar α can simply be done by multiplying the core tensorŵ : ) operations, and no change in rank and storage.

Addition of two Tucker tensors is computed as follows:
Similarly, the vectors w (i) as follows for i = 1, . . . , d and µ i = 1, . . . , ρ i + r i : For details see [31,39]. The computational cost is again O(1) since only more memory needs to be allocated, but the rank generally increases to the sum of the individual ranks of u and v .
The Hadamard product of two tensors given in Tucker format is computed as follows: Note [37] that the new Tucker core has size r 2d , i.e. the ranks get multiplied. The storage cost of the Hadamard product is O(d r 2 n + r 2d ).
The Euclidean inner product is computed as follows: The overall computational complexity of the inner product is O(d n r 2 + r 2d ).

Rank truncation in the Tucker format
To truncate the Tucker tensor rank from 2r or r 2 back to r the high order SVD (HOSVD) algorithm together with ALS iterations are applied [8]. As the set of Tucker tensors of fixed rank is closed, the best approximation exist, and there is an algorithm with quadratic convergence in the energy norm. One may compare this to the matrix case, where the SVD yields a best low-rank approximation. The truncated tensor w * resulting from the HOSVD could be sub-optimal, but nevertheless the following inequality holds For some improved estimates see [63]. The most time-consuming part of the Tucker algorithm is the HOSVD procedure. Namely, the computation of the initial guess using the SVD of the matrix unfolding of the original tensor [8]. The numerical cost of Tucker decomposition for full size tensors is O(n d+1 ).

The Tensor Train format
The tensor train (TT) format is described in [53,51,31,39]. As already noted, it was originally developed in quantum chemistry as "matrix product states" (MPS), see [64] and references therein, and rediscovered later [55,56]. For a motivation, let us start with a well-known example [1].

Example 3.4. Consider the d-dimensional
Laplacian operator discretised with standard finite differences over a uniform tensor grid with n degrees of freedom in each direction. It has the Kronecker (canonical) rank-d representation: I ⊗ · · · ⊗ I +I ⊗ A ⊗ I · · · ⊗ I + · · · + I ⊗ · · · ⊗ I ⊗ A ∈ R n d ×n d with A = tridiag{−1, 2, −1} ∈ R n×n , and I being the n × n identity. However, the same operator in the TT format is explicitly representable with all TT ranks equal to 2 for any dimension [1], where the "strong Kronecker" product operation " " is defined as a regular matrix product for the first level of TT cores, i.e. the block matrices, and the inner blocks are multiplied by means of the Kronecker or tensor product; for example In the element-wise matrix product notation, the Laplace operator reads Thus, the d-dimensional Laplacian operator is separable with TT ranks equal to 2, and each TT component (core) is defined by a one-dimensional Laplacian. Definition 3.5 (TT-Format, TT-Representation, TT-Ranks). The TT-tensor format is for variable TT-representation ranks r = (r 0 , . . . , r d ) ∈ N d+1 -with r 0 = r d = 1 and under the assumption that d > 2 -defined by the following multilinear mapping We call w := (w i 1 ...i d ) = U TT,r (P ) a tensor represented in the train tensor format. Note that the TT-cores W (ν) may be viewed as a vector valued r ν−1 × r ν matrix with the vector w (ν)

The representation in components is then
Alternatively, each TT-core W (ν) may be seen as a vector of r ν−1 × r ν matrices W which explains the name matrix product state. Observe that the first matrix is always a row vector as r 0 = 1, and the last matrix is always a column vector as r d = 1. The matrix components W (ν) iν of the TT-cores are also called "carriages" or "waggons" with "wheels" i ν at the bottom, coupled to the next "carriage" or "waggon" via the matrix product. This explains the tensor train name. If one notes more carefully W (ν) ∈ P r ν−1 ×rν (42) can be written more concisely as where U × k V is a contraction of the k-th index of U with the -th index of V , where one often writes just × k for × 1 k . Thus in Eq. (43) the contractions leave the indices from the R Mν untouched, so that the tensor w is formed.
Each TT-core (or block) W (ν) is defined by r ν−1 ×r ν ×M ν numbers. Assuming n = M ν for all ν = 1, . . . , d, the total number of entries scales as O(d n r 2 ), which is tractable as long as r = max{r k } is moderate.
A pictorial representation of the schema for the TT tensor format is shown in Figure 3. It shows d connected waggons with one wheel. The waggons denote the TT-cores, and Figure 3: Schema of the TT tensor decomposition. The waggons denote the TT cores and each wheel denotes the index i ν . Each waggon is connected with neighbours by indices j ν−1 and j ν . each wheel denotes the index i ν . The waggons for ν = 2, . . . , (d − 1) are connected with their neighbours by two indices j ν−1 and j ν . The first and the last waggons are connected by only one index, namely j 1 and j d−1 respectively. Since by the convention in Definition 3.5 above, r 0 = r d = 1, the indices j 0 and j d run from 1 to 1, i.e. are purely formal.
The waggon or carriage 2 -W (2) -is a tensor of degree 3, described by three indices j 1 (the left hitch), i 2 (the wheel), and j 2 (the right hitch). Multiplication of the third TT-core with the second and forth cores means the tensor contraction by the indices j 1 and j 2 . If we perform tensor contraction of all TT-cores over the indices j 1 , . . . , j d−1 , and disregard the purely formal constant indices j 0 and j d , then the indices j 0 , . . . , j d -the hitches -will disappear, and only the indices i 1 , . . . , i d -the wheels -will be left.

Basic operations with the TT format
We follow to the work of Oseledets [51] and list the major properties of the TT-tensor format.

Addition of two TT-tensors
Assume two tensors u and v are given in the TT-tensor format as in Eq. (42) iν . The sum w = u + v is given by the new cores W and the first and the last cores will be W (1) As only storage may have to be concatenated, the computational cost is O(1), but as the carriages resp. TT-cores grow, the final rank will generally be the sum of the ranks.
The Hadamard product w = u v in the TT format is computed as follows. Assume two tensors u and v are given in the TT tensor format as in Eq. (42), i.e. (u i 1 ..
The tensor w has also the TT-tensor format, namely with the new cores where ⊗ K is the Kronecker product of two matrices [31]. The rank of W (ν) = (W (ν) iν ) is the product of the ranks of the TT-cores U (ν) and V (ν) .
The Euclidean inner product of two tensors in the TT-format as in Eq. (40) with ranks r u and r v can be computed as follows: The computational complexity is O(d n r 4 ), and can be reduced further [54].

Rank truncation in the TT format
The rank truncation opearation is based on the SVD algorithm and requires O(d n r 3 ) operations [25]. The TT-rounding algorithm (p. 2305 in [56]) is based on QR decomposition and costs O(d n r 3 ). Corollary 2.4 in [56] states that for a given tensor w and rank bounds r k , the best approximation to w in the Frobenius norm with TT-ranks bounded by r k always exist (denote it by w * ), and the TT-approximation u computed by the TT-SVD algorithm (p. 2301 in [56]) is quasi-optimal: In [42] the authors suggested a new re-compression randomised algorithm for Tucker and TT tensor formats.

Numerical examples
The algorithms for post-processing high-dimensional data, or large volumes of data are shown on a few illustrative examples. Obviously, in any such application the actual computing times and even the possibility of executing such algorithms depend not only on that the huge amounts of data can be compressed to reasonable size, but also on the possibility of executing the algebraic operations with a reasonable speed on the compressed data. Additionally, as the compression may deteriorate in the course of the computation, one must assume that the intermediate results can be compressed to reasonable size, and that this re-compression can be done efficiently. All this depends very much on the data, and currently there are only very few general results to predict when this might or might not occur. Some results in this direction are already contained in [20,12,13] for the low-rank TT-format which was also used in our computations.

Finding the value and location of the maximum element
Already in Table 1 in Section 1 one could see an example from [14], illustrating the kinds of savings which are possible when compared to the method of inspecting every element.
Here is another very similar example, again the solution to the standard finite-difference discretisation of the same d-dimensional Poisson equation as detailed in Subsection 1.2 with n = 100 discretisation points in every direction. This is the special stationary case κ(ω, x) ≡ 1 of Eq. (2)   Here a 3 GHz CPU was used, so assume -unrealistically -that in each cycle one element could be inspected. Then the third column in Table 2 gives the required number of years to do so on such a machine, whereas the first columns give the dimension and the second column the full number n d of elements of the solution. Obviously this task could be parallelised, but no amount of conceivable computing power in form of more CPUs could beat the exponential growth of N = n d . Again full storage and the naïve algorithm of inspecting every element are -as far as we know -beyond any possibility in this universe. We used low-rank TT-format storage, and the Algorithm 3. In no case did it need more than 20 iterations. The actual computing times on a 3 GHz CPU are shown in the fourth column, and they behave like O(n d).

Finding level sets
We show two examples, the first is the same as in the previous Subsection 4.1, wheres the second one is a solution to a stochastic partial differential equation (SPDE).
Level set of the solution to a high-dimensional Poisson equation is the first example. The data vector is the same one as used in Subsection 4.1. But now we want to find the indices of the elements which are in the set S = ]− ∞, 0.25[, and according to Subsection 2.3 one has to compute the characteristic function χ S as given in Eq. (22).
For this one needs the sign function given in Eq. (20), computed with the Newton-Schulz iteration function in Eq. (30). The results are shown in Table 3. The first three columns are naturally the same as in Table 2, and the fourth column are the actual computing times in seconds on a 3 GHz CPU, which behave like O(n d).
Level sets of the solution to a stochastic partial differential equation (SPDE) is the next example. This example is described in a bit more detail, see also [20].  Table 3: Computing times (4th column) on 3 GHz CPU to compute χ S .
Two scientific software libraries sglib and TensorCalculus [66,18] were used. With the sglib procedures we discretise the SPDE Eq. (45), and with the TensorCalculus we solve the obtained tensor equation and compute the quantities of interests.
Consider the following stationary case of Eq. (2) in Subsection 1.1 a diffusion equation with uncertain diffusion coefficient κ(ω, x): Here G is the two-dimensional L-shaped domain [−1, 1] 2 \ [0, 1] 2 , and a triangular mesh with 557 mesh points was used for the spatial piece-wise linear finite element discretisation, and including the boundary values, the discrete spatial solution needs M d = 557 storage locations. After a stochastic discretisation with the Karhunen-Loève and Polynomial Chaos Expansion, and applying the stochastic Galerkin method as in [43,46,12,13,20,19], we can compute the solution u(ω, x) in a chosen tensor representation.
The random field κ(ω, x) was taken to have a shifted log-normal distribution for, i.e. log(κ(ω, x) − 1.1) has a normal distribution with parameters {µ = 0.5, σ 2 = 1.0}. The isotropic and homogeneous covariance function is of Gaussian type with covariance lengths x = y = 0.3.
Ten Karhunen-Loève terms were taken to represent the mean zero random part of the field κ(ω, x), and for the polynomial chaos expansion (PCE) multivariate Hermite polynomials of maximum second degree were taken.
For the right-hand side f (ω, x) a β-distribution {4, 2} was chosen for the random field. The covariance function is also of Gaussian type with covariance lengths x = y = 0.6. Again ten terms in the Karhunen-Loève expansion for the zero mean random part of f (ω, x) with a maximum second degree polynomial chaos were taken.
The total stochastic dimension of the solution u(ω, x) is 20 -ten random variables from the ten terms of the Karhunen-Loève expansion of the coefficient κ(ω, x), plus ten terms of the Karhunen-Loève expansion of the right hand side f (ω, x) -i.e. the multiindex α will consist of 20 indices (α = (α 1 , ..., α 20 )), and the solution will be represented as a function of 20 independent random variables. Together with the extra dimension from the spatial discretisation -which is lumped into one -one has a discrete version of a function u (ω, x)  For the numerical solution, a low-rank tensor solver was used [48,12,13,20,19], and the solution was represented in the CP-format Eq. (37) with a rank of r = 231. Thus  Finally, we computed χ S (u) for S =] − ∞, b[ by computing sign(b u ∞  − u) for b ∈ {0.2, 0.4, 0.6, 0.8}, which is shown in the first column of Table 4. The final representation ranks of χ S (u) resp. sign(b u ∞  − u) are given in the second column. In this numerical example, the ranks are in all cases smaller then 13. The sign-function was computed with the iteration from Algorithm 1. The iteration function is given in Eq. (30), and the iteration terminated after k max steps, shown in column four. The maximal representation rank of the iterates is documented in the third column. The error  − u kmax u kmax / (b u ∞  − u) is given in the last column. Each computationeach row of the table -took less than 10 min on a 3 GHz CPU.

Conclusion
Large volumes of data are used resp. generated in an increasing number of instances in science, technology, and business; some of which were mentioned in Subsection 1.1. Often this data is in the form of a sample of a high dimensional function, or it can be reshaped in this way. It may well happen that it is not possible to store such data in its entirety, so that it has to be compressed in some way. And it is usually not enough to compress these large data sets, efficient numerical algorithms are required to post-process such data. All these post-processing tasks may be trivial in low dimensions and for small data sets, but they become non-trivial for high-dimensional data, say d > 5, which is in some compressed format. We have introduced some numerical methods which may partially close this gap.
Here we have formulated algorithms for such post-processing, which only use the structure of an abstract commutative algebra. The Hadamard algebra is an example of such an algebra, and it can be defined for all real valued data sets. The Hadamard product of two tensors turned out to be a key tool in approximating the variance, pointwise inverse, level sets, frequency, and the maximal element of huge multi-dimensional data sets.
In our example application, the compression is low-rank tensor approximation, and we showed how to perform the algebra operations in some example low-rank tensor formats. We assumed that the initial data set has low-rank tensor representation, and showed that this low-rank property can be preserved during the whole computing process.
The choice of the tensor format is not crucial; it could be, for example, the canonical polyadic format, the Tucker format, the tensor train format, or any other. The only requirement is that it should be possible to compute the operations of the Hadamard algebra in a reasonable (linear in n) time, and the existence of stable rank truncation algorithms to truncate possibly large intermediate ranks. The algorithms needed were among others the multiplicative inverse and the sign function, and this could be fashioned after the algorithms for matrix algebras. Finally, the algorithms are demonstrated on some high dimensional data which comes from the solution of high dimensional or parametric resp. stochastic elliptic partial differential equations.