Tensor Approximation for Multidimensional and Multivariate Data

Tensor decomposition methods and multilinear algebra are powerful tools to cope with challenges around multidimensional and multivariate data in computer graphics, image processing and data visualization, in particular with respect to compact representation and processing of increasingly large-scale data sets. Initially proposed as an extension of the concept of matrix rank for 3 and more dimensions, tensor decomposition methods have found applications in a remarkably wide range of disciplines. We brieﬂy review the main concepts of tensor decompositions and their application to multidimensional visual data. Furthermore, we will include a ﬁrst outlook on porting these techniques to multivariate data such as vector and tensor ﬁelds.


Introduction
Data approximation is widely used in the fields of computer graphics and scientific visualizations.One way to achieve it is to decompose the data into a more compact and compressed representation.The general idea is to express a dataset by a set of bases, which are used to reconstruct the dataset to its approximation when needed (see Fig. 1).More specifically, such a representation usually consists of the actual bases and the corresponding coefficients describing the relationship between the original data and the actual bases.Typically, such compact representations consist of less data than the original dataset, capture the most significant features, and, moreover, describe the data in a format that is convenient for adaptive data loading and access.
A A decompose reconstruct

bases + cients compact data representation
Fig. 1 Compact data representation for a 3rd-order tensor A (a volume) by bases and coefficients that can be used to reconstruct the data to its approximation A Bases for compact data representation can be classified into two types: pre-defined and learned bases.Pre-defined bases comprise a given function or filter, which is applied to the dataset without any a-priori knowledge of the correlations in the dataset.In contrast, learned bases are generated from the dataset itself.Established examples of the former are the Fourier transform (FT) and the wavelet transform (WT).Wellknown examples of the latter are the principal component analysis (PCA), the singular value decomposition (SVD) and vector quantization.Using predefined bases is often computationally cheaper, but learned bases potentially remove more redundancy from a dataset.
Generally, PCA-like methods are able to extract the main data direction of the dataset and represent the data in another coordinate system such that it makes it easier for the user to find the major contributions within the dataset.To exploit this, PCAs higher-order extension-tensor approximation (TA)-can be used for multidimensional datasets.The first occurrence of TA was in [21].The idea of multiway analysis, however, is generally attributed to Catell in 1944 [11].It took a few decades until tensor approximations received more widespread attention, namely by several authors in the field of psychometrics [10,19,52].

Higher-Order Data Decompositions
The matrix SVD works on 2D data arrays and exploits the fact that the dataset can be represented with a few highly significant coefficients and corresponding reconstruction vectors based on the matrix rank reduction concept.The SVD, being a multilinear algebra tool, computes (a) a rank-R decomposition and (b) orthonormal row and column vector matrices.Unlike for 2D, the extension to higher-orders is not unique and these two main properties are captured by two different models that are both called tensor approximations : the Tucker model [14,15,25,49,52] preserves the orthonormal factor matrices while the CP model (from CANDECOMP [10] and PARAFAC [19]) preserves the rank-R decomposition.
Generally speaking, a tensor is a term for a higher-order generalization of a vector or a multidimensional array.In TA approaches, a multidimensional input dataset in array form, i.e., a tensor, is factorized into a sum of rank-one tensors or into a product of a core tensor (coefficients that describe the relationship to input data) and matrices (bases), i.e., one for each dimension.This factorization process is generally known as tensor decomposition, while the reverse process of the decomposition is the tensor reconstruction.
Tensor decompositions have been widely studied in other fields and were reviewed [13,25,34] and summarized [27,44].Since TA was emerging from various disciplines, it was developed under multiple names.The CP model was independently developed under the terms CANDECOMP and PARAFAC, therefore, it is sometimes referenced with a single name.The Tucker model takes its name from Tucker, who initially worked on the three-mode factor analysis (3MFA), which is sometimes referred to as the Tucker3 model, also called the three mode PCA (3MPCA) [27,28,49].Similarly the model was referenced to as n-mode PCA [23] since it is equivalent to applying PCA n times, each time along a different mode of the tensor.In [14] all these previous works were captured and written down as generalization of the SVD as multilinear singular SVD, which is usually termed as higher-order SVD or HOSVD.Furthermore, it was also called n-mode SVD in [54,55].
Given an input tensor, different decompositions capture different types of structures and result in varying numbers of coefficients.For a given accuracy, the number of CP ranks required to decompose a tensor is usually much larger to that of Tucker ranks.On the other hand, CP's storage cost grows only linearly with respect to its ranks, whereas that relationship becomes exponential in the case of Tucker.In sum, there is no silver bullet: CP is more suitable than Tucker for certain types of data, and vice versa.As a rule of thumb: • Dense tensors of moderate dimensionality n over continuous variables, for example n = 3 or n = 4 including spatial and temporal axes, can often be compressed more compactly via the Tucker model.• Tensors with categorical variables, sparse tensors, and tensors of higher dimensionality (say, n ≥ 5) often benefit more from the CP model.
The data sets addressed in this chapter fit in the first category, and so we restrict our experiments to the Tucker model.To further illustrate the usual advantage of Tucker over CP for spatial visual data, see Fig. 2, we compare the number of coefficients and root mean squared error (RMSE) obtained with CP vs. Tucker using different numbers of ranks for a 256 × 256 × 256 CT scan of a bonsai.
In scientific visualization, TA methods have first been introduced for interactive multiresolution and multiscale direct volume rendering [6][7][8][46][47][48].Additionally, their compact representation power has been exploited for 3D volume data compression [4,6] with notable advantages over other approaches at extreme compression ratios [2].In this work, we explore the multiscale feature expressiveness of TA methods for the first time on vector fields, i.e. multidimensional multivariate data .Hence we interpret the vector field as a 4D array, or as a 4th-order data tensor.However, the rank-reduction of TA is only applied to the three spatial dimension.

Motivation and Contributions
In the results section (Sect.9) we demonstrate that the feature sensitive approximation power of TA methods carries over from scalar to vector fields.These first results are promising and encourage the extension of tensor compression techniques to multivariate data fields, e.g. for compact storage and quick visualization at variable feature scales of large vector data in the fields of computational fluid dynamics or biomedicine.Moreover, based on the compression-domain data filtering and processing capabilities demonstrated in [5], we expect that important vector-field operators such as divergence or vorticity as well as other feature extraction operations can be analyzed and performed directly in the compressed TA format.

Singular Value Decomposition
The SVD is a widely used matrix factorization method to solve linear least-square problems.The SVD can be applied to any square or rectangular matrix A ∈ R M×N .Hence, the decomposition is always possible.The aim of the SVD is to produce a diagonalization of the input matrix A. Since the input matrix A is not symmetric, two bases (matrices) are needed to diagonalize A. Therefore, the SVD produces a matrix factorization into two orthogonal bases U ∈ R M×M and V ∈ R N ×N and a diagonal matrix ∈ R M×N , as expressed in Eq. (1) (matrix form) or Eq.(2) (summation form). (1) The bases U and V contain orthogonal unit length vectors u j and v j , respectively, and represent a r -dimensional column space (R M ) and a r -dimensional row space (R N ).Hence, the bases U and V are even orthonormal.The diagonal matrix contains the singular values σ r , where A singular value and a pair of singular vectors of a square or rectangular matrix A correspond to a non-negative scalar σ and two non-zero vectors u j and v j , respectively.The vectors u j are the left singular vectors, and the vectors v j are the right singular vectors (see Fig. 3).The number of non-zero singular values determines the rank R of the matrix A.
In applications truncated versions of the SVD are frequently desired.That is, only the first K singular values σ 1 . . .σ K and the corresponding K singular vectors u 1 . . .u K and v 1 . . .v K are used for the reconstruction.This approach is referred to as low-rank approximation of a truncated SVD.Basically, each weighted outer vectorproduct term σ j • u j • v j corresponds to a rank-one component (see also Fig. 4), and the SVD of matrices or images consequently represents a 2D data array eventually as a sum of such rank-one components.

Tensor Approximation Notation and Definitions
The notation taken here is mostly taken from De Lathauwer et al. [13,14], Smilde et al. [44], as well as Kolda and Bader [25], who follow the notation proposed by Kiers [24].To illustrate higher-order extensions we mostly make examples of order three.

General Notation
A tensor is a multidimensional array (or an N -way data array): a 0th-order tensor (tensor0) is a scalar, a 1st-order tensor (tensor1) is a vector, a 2nd-order tensor is a matrix, and a 3rd-order (tensor3) is a volume, see Fig. 5.We consistently use the letter 'a' to represent the data, following the notation of, e.g., [14,15,51,61,62].We use lower case letters for a scalar a, lower case boldface letters for a vector a in R I 1 , capital boldface letters for a matrix A in R I 1 ×I 2 , and calligraphic letters for a 3rd-order tensor A in R I 1 ×I 2 ×I 3 .A Fig. 5 A tensor is a multidimensional array: a 2nd-order tensor is a matrix A and a 3rd-order tensor is a volume A Fig. 6 Slices of a 3rd-order tensor A The order of a tensor is the number of data directions, also referred as ways or modes.Along a mode n, the index i n runs from 1 to I N .By using lower script indices for the modes, we can extend the index scheme to any order, i.e., I 1 , I 2 , I 3 , I 4 , . . . .The ith entry of a vector a is denoted by a i , an element (i 1 , i 2 ) of a matrix A is denoted by a i 1 i 2 , and an element (i 1 , i 2 , i 3 ) of a 3rd-order tensor A is denoted by The general term fibers is used as a generalization for vectors taken along different modes in a tensor.A fiber is defined by fixing every index but one.A matrix column is thus a mode-1 fiber and a matrix row is a mode-2 fiber.3rd-order tensors have column, row, and tube fibers, denoted by a i 1 , a i 2 , and a i 3 , respectively.Sometimes, fibers are also called mode-n vectors.
Slices are two-dimensional sub-sections of a tensor (e.g., one fixed index in a 3rdorder tensor).For a 3rd-order tensor A, there are, for example, frontal, horizontal, and lateral slices, denoted by A i 3 , A i 1 , and A i 2 , respectively as illustrated in Fig. 6.
For computations, a tensor is often reorganized into a matrix what we denote as tensor unfolding (sometimes called matricization).There are two main unfolding strategies, backward cyclic unfolding [14] and forward cyclic unfolding [24] as shown in Fig. 7.An unfolded tensor in matrix shape is denoted with a subscript in parentheses, e.g., A (n) .

Computing with Tensors
While many more operations on tensors exist, here we only outline the most common products used within the scope of this work.
• An N th-order tensor is defined as • The tensor product is denoted here by ⊗: however, other symbols are used in the literature, too.For rank-one tensors, the tensor product corresponds to the vector outer product (•) of N vectors b (n) ∈ R I n and results in an N th-order tensor A. The tensor product or vector outer product for a 3rd-order rank-one tensor is illustrated Fig. 7 Forward cyclic unfolding [24] of a 3rd-order tensor in Fig. 8: A = b (1) • b (2) • b (3) , where an element (i 1 , i 2 , i 3 ) of A is given by the products of their entries, i.e., Eq. ( 3).
• The n-mode product [14] multiplies a tensor by a matrix (or vector) in mode n.
The n-mode product of a tensor That is, element-wise we have Eq.( 4).
Each mode-n fiber is multiplied by the matrix C. The idea can also be expressed in terms of unfolded tensors (reorganization of a tensor into a matrix as described in Sect.3.1).
The n-mode product of a tensor with a matrix is related to a change of basis in the case when a tensor defines a multilinear operator [25].The n-mode product is the generalized operand to compute tensor times matrix (TTM) multiplications, and can best be illustrated using unfolded tensors as in Fig. 9.
Fig. 8 Three-way outer product for a 3rd-order rank-one tensor A = b (1) • b (2) • b (3)   Fig. 9 TTM multiplication C • B (n) , multiplying the unfolded 3rd-order tensor • The norm of a tensor A ∈ R I 1 ×I 2 ×•••×I N is defined analogously to the matrix Frobenius norm A F and is the square root of the sum squares of all its elements, i.e., Eq. (6).

Rank of a Tensor
In order to describe the definitions of the tensor rank, the definition for the matrix rank is recaptured.The matrix rank of a matrix A is defined over its column and row ranks, i.e., the column and row matrix rank of a matrix A is the maximal number of linearly independent columns and rows of A that can be chosen, respectively.For matrices, the column rank and the row rank are always equal and, a matrix rank is therefore simply denoted as rank(A).A tensor rank is defined similarly to the matrix rank, however, there are differences.In fact, the extension of the rank concept is not uniquely defined in higher orders and we review the definitions for the tensor ranks from [14] here.
• The n-rank of a tensor A, denoted by R n = rank n (A), is the dimension of the vector space spanned by mode-n vectors, where the mode-n vectors of A are the column vectors of the unfolding A (n) , and rank n (A) = rank(A (n) ).Unlike matrices, the different n-ranks of a tensor are not necessarily the same.
• A rank-one tensor is an N -way tensor A ∈ R I 1 ×I 2 ×•••×I N under the condition that it can be expressed as the outer product of N vectors, as in Eq. ( 7) (see also [12,30]).A rank-one tensor is also known under the term Kruskal tensor.
• The tensor rank R = rank(A) is the minimal number of rank-one tensors that yield A in a linear combination (see [12,14,25,30]).Except for the special case of matrices, the tensor rank is not necessarily equal to any of its n-ranks, but it always holds that R n ≤ R.

Tensor Decompositions
In tensor decompositions an input tensor the relationship/interactivity between A and the set of U (n) .Historically, as seen earlier, tensor decompositions are a higher-order extension of the matrix SVD.The nice properties of the matrix SVD, i.e., rank-R decomposition and orthonormal row-space vectors and column-space vectors, do not extend uniquely to higher orders.The rank-R decomposition can be achieved with the socalled CP model, while the orthonormal row and column vectors are preserved in the so-called Tucker model.An extensive review of the two models and further hybrid models can be found in [25].Here, we only outline the Tucker model that we apply in our experiments.

Tucker Model
The Tucker model is a widely used approach for tensor decompositions.As given in Eq. ( 8), any higher-order tensor is approximated by a product of a core tensor B ∈ R R 1 ×R 2 ×•••×R N and its factor matrices U (n) ∈ R I n ×R n , where the products × n denote the n-mode product as outlined in Sect.3.2.This decomposition can then again be reconstructed to its approximation A. The missing information of the input tensor A that cannot be captured by A is denoted with the error e.The Tucker decomposition is visualized for a 3rd-order tensor in Fig. 10.Equivalently, a Tucker decomposition can also be represented as a sum of rank-one tensors as in Eq. ( 9) and illustrated in Fig. 11.
Fig. 11 Tucker 3rd-order tensor as a sum of rank-one tensors: The column vectors u (n) r n of the factor matrices U (n) ∈ R I n ×R n are usually orthonormal and can be thought of as principal components R n in each mode n [25].The core tensor its factor matrices and is always of the same order as the input data.The core tensor is computed in general, as shown in Eq. (10), and for orthogonal factor matrices as in Eq. (11).The element-wise core tensor computation is denoted in Eq. (12).In other words, the core tensor coefficients b r 1 r 2 ...r N represent the relationship between the Tucker model and the original data. (10) . . .
The Tucker decomposition is not unique, which means that we can modify the core tensor B without affecting the model fit as long as we apply the same changes to the factor matrices (so-called core tensor rotations), for more details see [25].
Often, we are interested in compact models, which enable a compression of the input data.For example, after computing a (full) Tucker decomposition the core tensor B has the same size as the original input A and all the factor matrices are square.However, we are more interested in reduced-size, approximative Tucker decompositions, where B is an element of R R 1 ×R 2 ×R 3 with R 1 < I 1 , R 2 < I 2 and R 3 < I 3 .Using so-called rank-reduced tensor decompositions or truncated tensor decompositions one can directly obtain more compact decompositions.

Tensor Rank Reduction
As seen in Sect.3.3, the extension of the matrix rank concept to higher orders is not unique and we will mostly follow the rank-(R 1 , R 2 , . . ., R N ) tensor decomposition and reduced-rank approximation of the Tucker model here.

. , R N ) Approximations
A simple rank-one approximation is defined as A = λ • u (1) • u (2) • • • • u (N ) from the rank-one tensor (vector) product (•) of its N basis vectors u (n) ∈ R I n and a weight factor λ. Hence a tensor A could be approximated by a linear combination of many rank-one approximations as in Eq. ( 13).This approximation, also known as a CP model, is called a rank-R approximation.
Alternatively, if we allow all weighted tensor (vector) products u (1) i N of any arbitrary index combinations i 1 i 2 . . .i N , we end up with the Tucker model of Sect.4.1 and Eq. 12 where the weight factors for all index combinations form the core tensor B. Choosing R 1,...,N < I 1,...,N we end up with a rank-(R 1 , R 2 , . . ., R N ) approximation of A, which is given by a decomposition into a lower-rank tensor a given reduced rank space (Eq.( 14)).This rank-(R 1 , R 2 , . . ., R N ) approximation was previously introduced as the Tucker model.
In general, a rank-reduced approximation is sought such that the least-squares cost function of Eq. ( 15) is minimized.

A = arg min
Given that (R 1 , R 2 , . . ., R N ) are sufficiently smaller than the initial dimensions (I 1 , I 2 , . . ., I N ), the core coefficients B ∈ R R 1 ×R 2 ×•••×R N and the factor matrices U (n) ∈ R I n ×R n can lead to a compact approximation of A of the original tensor A. In particular, the multilinear rank-(R 1 , R 2 , . . ., R N ) is typically explicitly chosen to be smaller than the initial ranks in order to achieve a compression of the input data (see also [2,4,6]).

Truncated Tensor Decomposition
Similar to matrix SVD, tensor rank reduction can be used to generate lower-rank reconstructions A of the input A. The tensor rank parameters R n are chiefly responsible for the number of TA bases and coefficients that are used for the reconstruction and hence are responsible for the approximation level.In higher orders, the CP decomposition is not directly rank-reducible, however, the truncation of the Tucker decomposition is possible due to the all-orthogonality property of the core tensor.
For a 3rd-order tensor, all-orthogonality means that the different horizontal matrix slices of the core B (the first index i 1 is kept fixed, while the two other indices, i 2 and i 3 , are free) are mutually orthogonal with respect to the scalar product of matrices (i.e., the sum of the products of the corresponding entries vanishes).The same holds for the other slices with fixed indices i 2 and i 3 (see [14]).Therefore, given an initial sufficiently accurate rank-(R 1 , R 2 , R 3 ) Tucker model, we can progressively choose lower ranks K n < R n for reduced quality reconstructions.As indicated in Fig. 12, the ranks K n indicate how many factor matrix columns and corresponding core tensor entries are used for the reconstruction.
Note that the ordering of the coefficients in the Tucker core tensor B is not strictly decreasing in contrast to the decreasing singular values in the matrix SVD case.However, in practice it can be shown that progressive tensor rank reduction in the Tucker model works very well for adaptive reconstruction of the data at different accuracy levels.
Fig. 12 Illustration of a rank reduced Tucker tensor reconstruction: A reduced range of factor matrix columns with corresponding fewer core tensor entries reconstructs a lower quality approximation but at full resolution

Tucker Decomposition Algorithms
There are various strategies for how to compute and generate a tensor decomposition.The most popular and widely used group of algorithms belongs to the alternating least squares (ALS) algorithms, the other group of algorithms uses various Newton methods.The respective algorithms differ for the computation of the different tensor models, and we will mainly focus on the Tucker model in our review.
For the Tucker model, the first decomposition algorithms used were a simple higher-order SVD (HOSVD) (see [14]), the so-called Tucker1 [52], a three-mode SVD.However, the truncated decompositions of higher orders are not optimal in terms of best fit, which is measured by the Frobenius norm of the difference.Starting from a HOSVD algorithm, tensor approximation ALS algorithms [26,28] were developed, where one of the first Tucker ALS was the so-called TUCKALS [49].Later various improvements accelerated [1] or optimized the basic TUCKALS method.The higher-order orthogonal iteration (HOOI) algorithm [15] is an iterative algorithm that performs a better fit for a truncated HOSVD version.
Newton methods are also used for the Tucker decomposition or rank-(R 1 , R 2 , . . ., R N ) approximation.They typically start with a HOOI initialization and then converge faster to the final point.[16] developed a Newton-Grassman optimization approach, which takes much fewer iterations than the basic HOOI -even though one single iteration is more expensive due to the computation of the Hessian.While the HOOI is not guaranteed to converge, the Newton-Grassmann Tucker decomposition is guaranteed to converge to a stationary point.Another Newton method was proposed by [22], who developed a differential-geometric Newton algorithm with a fast quadratic convergence of the algorithm in a neighborhood of the solution.Since this method is not guaranteed to converge to a global maximum, they support the method by starting with an initial guess of several HOOI iterations, which increases the chances of converging to a solution.

Tensor Reconstruction
The tensor reconstruction from a reduced-rank tensor decomposition can be achieved in multiple ways.One alternative is a progressive reconstruction: Each entry in the core tensor B is considered as weight for the outer product between the corresponding column vectors in the factor matrices U (n) .This gives rise to Eq. ( 16) for the Tucker reconstruction.
r N (16) This reconstruction strategy corresponds to forming rank-one tensors and cumulatively summing them up.The accumulated weighted subtensors then form the approximation A of the original data A. In particular for the Tucker model, this is an expensive reconstruction strategy since it involves multiple for-loops to run over all the summations, for a total cost of O(R N • I N ) operations.

Element-Wise Reconstruction
A simple approach is to reconstruct each required element of the approximated dataset individually, which we call element-wise reconstruction.Each element a i 1 i 2 i 3 is reconstructed, as shown in Eq. ( 17) for the Tucker reconstruction.That is, all core coefficients multiplied with the corresponding coefficients in the factor matrices are summed up (weighted sum).
Element-wise reconstruction requires O(R N ) operations on average.It can be beneficial for applications where only a sparse set of reconstructed elements are needed.

Optimized Tucker Reconstruction
A third reconstruction approach-applying only to the Tucker reconstruction-is to build the n-mode products along every mode, which leads to a TTM multiplication for each mode, e.g., TTM1 along mode 1, (see also Eq. ( 5)).This is analogous to the Tucker model given by Eq. ( 14).The intermediate results are then multiplied along the next modes, e.g., TTM2 and finally TTM3.In Fig. 13 we visualize the TTM reconstruction, and the intermediate results B and B as well as the final approximation A, applied to a 3rd-order tensor using n-mode products.
Fig. 13 Forward cyclic TTM multiplications after [24] along the three modes (n-mode products) Since it exploits matrix multiplications, this optimized algorithm is much faster than progressive reconstruction (Eq.( 16)).Its cost is dominated by

Useful TA Properties for Scientific Visualization
As stated in the introduction, TA is the higher-order generalization of the matrix SVD, which can offer either properties of (a) rank-R decomposition or (b) orthonormal row-space and column-space vectors.In higher orders, the orthonormal row and column vectors are preserved in the Tucker model which thus supports progressive rank-reduced approximations and reconstructions.

Spatial Selectivity and Subsampling
The Tucker model (Sect.4.1) consists of one factor matrix per mode (data direction) The core tensor B is in effect a projection of the original data A onto the basis of the factor matrices U (n) .In case of a volume, the Tucker model has three modes, as illustrated in Fig. 10, and defines an approximation A = B × 1 U (1) × 2 U (2) × 3 U (3) of the original volume A (using n-mode products × n ).
The row and column axes of the factor matrices represent two different spaces: (1) the rows correspond to the spatial dimension in the corresponding mode, and (2) the columns to the approximation quality.These two properties can be exploited for multiresolution modeling (spatial selection and subsampling of rows) and multiscale approximation (rank reduction on the columns) (see also Fig. 14).In [5,47] we demonstrated how these features can be exploited for multiresolution and multiscale reconstruction as well as filtering of compressed volume data in the context of interactive visualization.

Approximation and Rank Reduction
As described earlier, the Tucker model defines a rank-(R 1 , R 2 , R 3 ) approximation, where a small R n corresponds to a low-rank approximation (with details removed) and a large R n corresponds to a more accurate approximation of the original.The highest rank R n for the initial Tucker decomposition has to be given explicitly.However, rank reductions can be applied after the initial decomposition (similar to the rank reduction in matrix SVD).Even though the core tensor coefficients are not guaranteed to be in decreasing order, as in matrix SVD, in practice it can be shown that progressive tensor rank reduction in the Tucker model works well for adaptive visualization of the data at different feature scales [2,[4][5][6]47].
As illustrated in Fig. 12 in Sect.5.2, the ranks indicate how many factor matrix columns and corresponding core tensor entries are used for a desired reconstruction.Thus, given a rank-(R 1 , R 2 , R 3 ) Tucker model, we can specifically or progressively choose lower ranks K n < R n for a reduced quality reconstruction, at the original spatial output resolution given by I n (or also subsampled).
Figure 15 shows the progressive rank reduction from an initial rank-(256, 256, 256) Tucker decomposition of an original 512 3 example volume.Shown are the visual results and the data reduction of the approximation at variable reduced-rank reconstructions.The numbers of coefficients include all core tensor and factor matrix entries that are used, e.g. a rank- (32,32,32) reconstruction corresponds to 32 3 + 3 • 512 • 32 = 81 920 coefficients.The data reduction ratio can be derived by dividing the number of coefficients by 512 3 , which for R = 32 results in using only 0.06% of the original amount of data.In particular, Fig. 15 demonstrates the power of lowrank tensor approximations that can be used for multiscale feature visualization or progressive image refinement in volume rendering.

Application to Multivariate Data
Encouraged by the data reduction power and the approximation quality of tensor approximations as shown in Fig. 16, we have extended the Tucker decomposition and reduced rank reconstruction to vector field data.Compared to scalar 3D volumes (e.g. from MRI or CT scans), a 3D vector field V (e.g. of velocities from a weather or fluid simulation) can be interpreted as a multivariate data field with D-dimensional vector-valued entries at each position in a I 1 × I 2 × I 3 grid.Thus we can interpret V as a 4th-order tensor V ∈ R I 1 ×I 2 ×I 3 ×D .Note that the TA rank-reduction (Sect.5) is only applied to the three spatial dimension as we are not interested in a dimensionality reduction of the vector-valued field value itself.

Dataset
For our first experiments we used a data set from the Johns Hopkins Turbulence Database http://turbulence.pha.jhu.edu/representing a direct numerical simulation of incompressible magnetohydrodynamic (MHD) equations (see also Fig. 17).The data set contains, among other output variables, 3 velocity components that we used, hence D = 3, covers a 3D grid of 256 3 cells (downsampled subset from the original 1024 3 ), hence I 1 = I 2 = I 3 = 256, and is the first time step from the output of the simulation.
The vector field therefore covers a 3D cube which we visualize using direct volume and streamline rendering techniques with color-coding of velocity, divergence, vorticity or error magnitudes in ParaView.Transparency is used to reduce clutter and opacity such as to focus the visualizations on the high-magnitude value regions.Note that this cubic vector field is dense and rendered over a black background, thus always appearing as a cube like object in the images.

Vector Field Magnitude and Angle
As can be seen in Fig. 17, the velocity magnitude and structure of the flow directions is very well maintained down to fairly course reconstructions using R = 32.Especially the (important) regions with larger velocity magnitudes are very well preserved with respect to their flow orientation as visible from the streamline visualization.The rendering applies an opacity transfer function setting which is almost linear, slight curve below the diagonal, to the vector field magnitude, highlighting the important high-velocity regions.
To see how many details are lost after rank reduction of the velocity field data set, we calculated and visualized the error information for both magnitude and angle deviation in Figs.18 and 19, respectively.For the magnitude, we normalized the error to be the percentage relative to the local vector field magnitude value.The rendering uses an opacity mapping which is almost linear, slight curve below the diagonal, to the vector field magnitude.This makes low magnitude and low opacity areas, and their errors, transparent or dark and highlights any errors in the more critical high-velocity regions.The completely dark images in Fig. 18e-h demonstrate that the magnitude error after rank reduction using Tucker decomposition is close to zero even with for R = 16, which is 1/16 of the original full rank of 256.The barely noticeable dark blue regions in Fig. 18a-d correspond to low errors in high-velocity regions, with errors in the range of less than 10% down to 10 −8 % (with 60% being white).
Compared to the almost negligible relative magnitude errors of the velocity field, however, the angular error of the same data set after rank reduction is more prominent as shown in Fig. 19.The maximum angular error π = 180 • is shown in red while small errors are shown in dark blue (0 • ).The rendering applies an opacity setting almost linear to the angle offset, hence in Fig. 19 large angular errors are highlighted.However, this does not relate to the local strength of the vector field since large angle differences in low-magnitude areas are not that important.Directly comparing Figs 17 and 19, which have the same viewing configuration, one can observe that large angular errors occur mostly in very low-velocity areas and may thus not be that relevant.In particular, since for very short vectors small errors can cause large angular changes.

Vorticity and Divergence
We also inspected the influence of low-rank tensor compression on the two most fundamental features of vector field data: divergence and vorticity.As can be seen in Fig. 20, the global structure of the magnitude from both the divergence field and the vorticity field are still very well preserved with rank R = 64, which is 1/4 of the original full rank of 256.The values are given in their absolute ranges of [0.36, −0.43] and [0.52, 10 −5 ] for the divergence and vorticity respectively.The rendering applies an opacity setting which is almost linear, slight curve below the diagonal, to the magnitude values, highlighting the important regions.
For the vorticity vector field, we also calculated the relative magnitude error as for the velocity vector field data shown in Fig. 18, thus as percentage of the vorticity magnitude, and similar conclusions can also be drawn here for the results shown in Fig. 21.All errors are very low, and barely noticeable for rank-reductions down to R = 16, which corresponds to 1/16 of the original data volume.
Additionally, we also computed the absolute angular error for the vorticity vector field, in this case applying an opacity proportional to the magnitude of the raw vector field, and the results are shown in Fig. 22.We can observe that the angular error almost everywhere is in the color range of dark to light blue which maps to errors from 0 up to 45 • in this plot, and thus the directional vorticity information seems to be well preserved.In the analysis conducted in this first study on vector fields, we have shown that tensor approximation methods are not only very useful for multidimensional scalar fields but can also be applied to multivariate data, thus extending to vector and possibly tensor fields.In general we can observe that for the important high-velocity vector field regions, the low-rank reconstructions maintain the important overall structures of the flow features, in particular also the vorticity.Furthermore, we note that the MHD simulation model in theory should be divergence-free and thus the numerical results report very low divergence numbers.Therefore, we are very satisfied with the result that the low-rank tensor reconstructions do not result in an unexpected and uncontrolled enlargement of these divergence values, keeping them within the numerical range of the simulation.61.Wang, H., Wu, Q., Shi, L., Yu, Y., Ahuja, N.: Out-of-core tensor approximation of multidimensional matrices of visual data.ACM Trans.Graph.The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

Fig. 3 Fig. 4
Fig. 3 Visualization of the summed form of the SVD as shown in Eq. (2)-illustrating the singular values with the corresponding left and right singular vector pairs

Fig. 16 a
Fig. 16 a A 512 3 turbulence volume of 512MB, b visually identical compression result using 51.2MB and c result after extreme compression down to only 1.71MB using the TTHRESH method [2].All the visualizations are generated using ParaView