1 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 realvalued 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 , 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
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
with an appropriate , reflecting the amount of data. The vector space operations are clear, as is the canonical Euclidean product, and for the commutative multiplication consider the pointwise or Hadamard product, i.e. the componentwise 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 realvalued function on a highdimensional 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 lowrank 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 lowrank 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 recompression or truncation.
1.1 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
(1) 
where is in some Hilbert space and is some parameter dependent operator; in particular could be some parameterdependent differential operator, for example
(2) 
where is a domain, is the time window of interest, is a parameterised random tensor field dependent on a random parameter in some probability space with probability measure , and one may take [43, 46].
For each one may thus seek for solutions in , where we have set . On the other hand, assume that for fixed
we are looking at the random variable
which we for simplicity assume to have finite variance, i.e. . Taking into account the parametric dependence, we are hence looking for a function which is defined on ; we are thus looking for a solution in . This applies equally well to the abstract eq:XXV, as on the right hand side the operator depends on , so will the solution to eq:XXV, and it will lie in a similar tensor product.Observe further that if the probability measure on in eq:XXV or eq:XXVI is a product measure on , the space can be split further , with . Thus in total the solution to both eq:XXV and eq:XXVI may be viewed as an element of a high order or degree tensor space [43, 46, 48]
(3) 
If we think of samples of at space points , time instances and parameter values , with a multiindex with , the samples of the solution
(4) 
form a high order or degree tensor [31].
In many chemical applications one has , and is few hundreds or even thousands [6]. Some other highdimensional problems from chemistry and physics like HartreeFock, Schrödinger, or Masterequations, and lowrank 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 highdimensional 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 multiparametric equation, e.g. weather research and forecasting, and climate models [23]. Also OilGas companies daily collect sensor data from multiple sources. Another source for huge volumes of data are highenergy particle accelerators like CERN [5].
For the sake of simplicity, we will treat all arguments of the function in eq:XXaX equally, denoting them in that case by . So just consider a realvalued function and its samples
(5) 
which form a tensor of order or degree representing a total of values, and where for the sake of convenience multiindices have been introduced.
As one may see, a tensor can be simply defined as a highorder matrix or multiindex array, where multiindices are used instead of indices. As an example, assume that we have data points. Now a vector is just a tensor of first order, but one may reshape the data, say into a matrix — a tensor of 2nd order. Reshaping further, one may take — a tensor of 4th order — as well as other combinations, such as — an array of size (16 times) — which is a tensor of 16th order; and finally all the way to — 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 , where , and consider its singular value decomposition (SVD) , where is the diagonal matrix of singular values , assumed arranged by decreasing value, and , collect the right and left singular vectors. Often there is a number , such that for some small one has for all . Then one can formulate a compressed or truncated version
where , , , and , . This reduced SVD is a special case of eq:CPI. For the sake of a concrete example, assume that . One may observe that the required storage has been reduced from to for the vectors .
The set of possible objects we want to work with will be denoted by . 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 pointwise addition and multiplication by scalars, such functions clearly form a vector space, and if one includes pointwise 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 . 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 matrices .
In our concrete examples we work directly with the data represented as tensors . Even without the identification with , it is obvious that this is naturally a vector space. As also , it is clear that can be equipped with the canonical Euclidean inner product from . The associative and commutative algebra product comes from the identification with the function algebra , i.e. the pointwise product, which is also known as the Hadamard product [31]. Hence it is clear that for data 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 , so that the tensor in eq:XXbX represents values. As already mentioned, our example of data compression is based on lowrank approximations to elements in . Although lowrank tensor data formats and techniques are almost unavoidable if working with large highdimensional 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 hence has terms, a compressed representation — some lowrank versions of which will be discussed in more detail in tensorrep — will have significantly fewer terms. For example, the CPdecomposition (canonical polyadic) representation, truncated to terms,
(6) 
If the rank is reasonably small compared to independent of and , then we have an approximation with much less storage, which also only depends on the product of and and is not exponential in the dimension . But now the maximum can not be easily looked up; to get a particular element, the expression eq:CPI has to be evaluated for that index combination. In the following we shall assume that we work with approximations such as in eq:CPI which need much less storage. One has to make sure then that the algebraic operations of the Hadamard algebra structure on can be performed efficiently in an approximative manner, so that they have much lower complexity than the obvious , although the single elements of which are usually needed for pointwise operations — and postprocessing, see postprocintro — are not directly available. Thus the motivating factors for applying compression, and in particular lowrank tensor techniques, include the following:

The storage cost is reduced, depending on the tensor format, from to , or to , where .

The lowrank tensor approximation is relatively new, but already a wellstudied technique with free software libraries available. Various lowrank 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 lowrank format. The basic fast Fourier transform on
would have complexity . These lowrank techniques can either be combined with the fast Fourier transform giving a complexity of , 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 lowrank tensor decomposition;

with sampling, it requires an axesparallel mesh;

although many functions have a lowrank 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.
1.2 Postprocessing 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:XXIX is something which allows efficient calculations to be performed on a sample of the solution like eq:XXaX or eq:XXbX with the help of expressions such as eq:CPI 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:XXV or eq:XXVI as function of 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 ,

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 , exponential in , which may be unacceptable for high
. 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 postprocessing algorithms, in the form of finding the maximum, was first presented in [14] for lowrank tensor approximations. Some further postprocessing tasks for such lowrank 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:XXV in xmplmotiv, or more specifically eq:XXVI in the stationary case with , a dimensional Poisson equation:
and righthandside
Assume further that this is solved numerically by a standard finitedifference method with gridpoints in each direction, so that the solution vector has data points in , but can also naturally be viewed as a tensor of degree in . So for full storage, one needs storage locations, and if one is looking for the maximum by inspecting each element, one would have to inspect elements.
We mention in passing that of course actually solving the discrete equation is a challenge for , but that is a different story, and for this we refer to [14] and the references therein.
# loc’s.:  years [a]  actual [14] time [s]  

inspect.  of Algorithm 2  
0.16  
0.42  
1.16  
2.58  
4.97  
8.56 
Now assume for the sake of simplicity that one can inspect elements per second — on an ideal 2 GHz CPU with one inspection per cycle. Then for the times needed to find the maximum for the full dataset are shown in tab:find_max 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 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 and grow exponentially with , whereas for a lowrank representation of not only can the data be stored on a modest laptop, but for the simple Algorithm 2 — to be explained later in auxil — the computing times in the fourth column are in terms of seconds and behave like .
1.3 State of the art
This is not a discussion of the state of the art regarding general tensor formats and their lowrank 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 lowrank approximations [26]. In the following, we concentrate on the question of postprocessing 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 pointwise inverse and the sign function were introduced, as well as the ranktruncation 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 lowlying 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
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 ranktruncated iterative eigensolvers. The authors were able to reduce the computational cost from to . They demonstrated the introduced algorithms by solving largescale spectral problems from quantum chemistry: the Schrödinger, the Hartree–Fock, and the Kohn–Sham equations in electronic structure calculations.1.4 Outline of the paper
In the following algos the necessary material for abstract algebras is quickly reviewed. Then the algorithms and functions on the algebra used to compute the postprocessing tasks outlined in postprocintro 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 tensorrep we present some concrete examples of compression of highdimensional data , once it is identified with a tensor with . For such tensors of degree we discuss the canonical polyadic (CP) representation, the Tucker representation, and the tensortrain (TT) representation, and the compression in terms of lowrank 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 lowrank representation. As the compression level may deteriorate, i.e. the rank of the approximation may grow, it is important that one is able to recompress. Pointers to such methods in the literature are included as well.
The numexp contains a few numerical examples and a discussion of their results used to illustrate the algorithms of algos and the representations from
tensorrep. These examples are solutions of elliptic highdimensional 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 conceptual —
reshape to consider it as an element of a tensor space. The conclusion is contained in concl.2 Algorithms for postprocessing
Given a dataset 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 postprocintro we want to consider are finding

the indices or values of maxima, minima of ,

the indices or values of maxima, minima of some function of ,

the index or value of the datum of closest to a given number,

level sets, i.e. (the indices of) all ,

the number of indices such that ,

the mean (sum of all items) and variance (sum of all items squared),

some auxiliary functions, such as the algebraic inverse of — in our case the Hadamard inverse — and others like the function which will be seen to be needed for the above tasks in postproctasks.
The above tasks and auxiliary functions will be computed by finding an iteration (mapping) for each postprocessing task , such that the fixedpoint of — i.e. — is either the sought solution for the task, or computes one of the auxiliary functions.
2.1 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 suboptimal, which means that the result has to be compressed again. Thus, when performing the algebraic operations for a fixedpoint iteration, the compression would either get worse and worse in each iteration, or one has to recompress or truncate the result again to a good compression level. Therefore a recompression 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 lowrank representation of tensors explained in tensorrep acts as compression, and the rank may increase through the algebraic operations, and hence we will use truncation to low rank with error [31, 2] to avoid the problem of a deteriorating compression level.
In other words, with a general compressed representation of the computation will be a truncated or perturbed iteration [34, 48]. If we denote the general compression mapping by — meaning compression with an accuracy — the iteration map is changed from to .
The general structure of the iterative algorithms for a postprocessing task 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 known that

if the iteration by is superlinearly convergent, the truncated iteration will still converge superlinearly, but finally stagnate in an neighbourhood of the fixed point [34]. One could loosely say that the superlinear convergence of is stronger than the truncation by .

if the iteration by is linearly convergent with contraction factor , the truncated iteration will still converge linearly, but finally stagnate in an neighbourhood of [48]. Again, one could loosely say that iteration by and truncation by 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 reconstruction accuracy and taking into account the possible influence due to the convergence behaviour of the iterator .
2.2 Preliminaries and basic algebraic operations
As already pointed out, we assume in general that the set is a vector space. Our example of the space of tensors of interest was already introduced in xmplmotiv. 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 tensorrep. Furthermore, as as vector spaces, we can carry the canonical Euclidean inner product on to , which for is denoted by
This makes into a Euclidean or Hilbert space. We assume that the computation of this inner product is feasible when both and are in a compressed representation; again this will be demonstrated in tensorrep for the lowrank representations we shall consider as examples.
On this space 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 point or componentwise multiplication on . It is the usual pointwise multiplication of functions, in this case the functions . For , define the Hadamard product as
(7) 
It is a bilinear operation (linear in each entry), and it is commutative. Thus this product makes into an associative and commutative algebra. The symbol will be used both for the abstract algebra product on , as well as for the Hadamard product on . As is any such algebra, for any it holds that .
Especially, as a function of , the product is a linear map on . This is the familiar canonical representation of an associative algebra as linear maps on itself, i.e. in this case in a commutative subalgebra of the algebra of all linear maps with concatenation as product.
It easy to see that this Hadamard algebra has a unit element for the product,
(8) 
a tensor with all entries equal to unity — this makes into a unital algebra. Observe that by defining for all the allones vectors , the Hadamard unit has rank one according to eq:CPI: . Obviously one has , the identity in .
Having defined a unit or neutral element for multiplication, we say that has a multiplicative inverse iff there is an element, denoted by , such that
(9) 
Obviously not all elements have an inverse, e.g. the zero element is never invertible. In any case, not every has a Hadamard inverse, for this it is necessary that all entries are nonzero, and then . Note that and as it should be, and it is easily seen that , the inverse of in .
Elements of the algebra which can be written as a square are called positive, they form a convex cone ; note that obviously the zero tensor and the multiplicative unit are positive, and that if an invertible is positive, so is its inverse . As usual, this defines an order relation , which implies that for any positive one has . Further one may observe that .
A certain compatibility of the inner product and the algebra product is needed, as we require that
(10) 
i.e. that the action of the algebra product is selfadjoint, or, in other words, that the maps are self adjoint. This condition is satisfied for the Hadamard algebra.
The inner product with the unit element defines, as function of , a positive linear functional , as for one has
It is a kind of “trace” or unnormalised state functional on the algebra, as in the Hadamard algebra. In particular
(11) 
Conversely, if the algebra comes equipped with such an unnormalised positive state functional, the eq:stateinnpr may be taken as the definition of an inner product, which automatically satisfies eq:selfadjmult.
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 lowrank formats in tensorrep), and that we can compute the multiplicative inverse also in compressed representation — again at least approximately.
As in any unital algebra, if for the element fails to be invertible — this means that is not invertible in — we shall say that is in the spectrum of ; here — as we are in finite dimensional spaces — it is an eigenvalue, i.e. there is an eigenvector such that . From this discussion it is immediate that in the Hadamard algebra case, if is an eigenvalue of or , there must be a multiindex such that , and hence in the real Hadamard algebra one has that always . Thus the spectrum of resp. is . Observe that according to eq:selfadjmult, each representing map from above is selfadjoint, which means that also in general all the spectra are real — .
Specifically, the Euclidean Hadamard algebra which we have constructed is obviously isomorphic — as a Euclidean algebra — to with the canonical Euclidean inner product when equipped with the point or elementwise Hadamard product . Let us denote this isomorphism by ; it implies some ordering of the terms of each .
More enlightening and less obvious may be the unital algebra isomorphism with the commutative subalgebra of diagonal matrices in the full matrix algebra with the usual matrix multiplication. Let us denote this isomorphism by . If is the vector containing all the elements of the tensor , then is the corresponding diagonal matrix, i.e. . Hence, by choosing the canonical Euclidean basis in and its image by as a basis in , the canonical representation is itself represented by the matrix .
Let be another tensor, with associated and , then from the definition of the isomorphism :
(12) 
which contains the important equality
(13) 
from where one sees that
(14) 
which helps in understanding the following algorithms. Note that this shows that the algebra 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 of an element 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 pointwise 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 of instead of with the “values of ”. These two notions only coincide for a function algebra.
2.3 Postprocessing Tasks
Here the different postprocessing 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 auxil.
Finding the maximum or minimum
of is the first task we consider. Observe, that the element of maximum modulus of the tensor is also equal to the norm .
Finding the maximum means finding the index
Comments
There are no comments yet.