On identification of self-similar characteristics using the Tensor Train decomposition method with application to channel turbulence flow
Thomas von Larcher Rupert Klein




Application of the Tensor Train decomposition to channel turbulence flow





Th. von Larcher Freie Universität Berlin, Institute of Mathematics, Arnimallee 9, 14195 Berlin, Germany
Tel.: +49-30-838-56093
 larcher@math.fu-berlin.de R. KleinFreie Universität Berlin, Institute of Mathematics, Arnimallee 6, 14195 Berlin, Germany


 

On identification of self-similar characteristics using the Tensor Train decomposition method with application to channel turbulence flow
    Thomas von Larcher Rupert Klein
    
=========================================================================================================================================




A study on the application of the Tensor Train decomposition method to 3D direct numerical simulation data of channel turbulence flow is presented. 
The approach is validated with respect to compression rate and storage requirement. In tests with synthetic data, it is found that grid-aligned self-similar patterns are well captured, and also the application to non grid-aligned self-similarity yields satisfying results. It is observed that the shape of the input Tensor significantly affects the compression rate. 
Applied to data of channel turbulent flow, the Tensor Train format allows for surprisingly high compression rates whilst ensuring low relative errors.






§ INTRODUCTION



Multidimensional data sets, data of dimension 3 or higher, require massive storage capacity that strongly depends on the (Tensor) dimension, d, and on the number of entities per dimension, n. The datasize or storage requirement scales with 𝒪 (n^d), n=max_i{n_i}. It is that curse of dimensionality that makes it difficult to handle higher-order Tensors or big data in an appropriate manner. Tensor product decomposition methods, <cit.>, were originally developed to yield low-rank, i.e. data-sparse, representations or approximations of high-dimensional data in mathematical issues. It has been shown that those methods are as good as approximations by classical functions, e.g., like polynomials and wavelets, and that they allow very compact representations of large-data sets. Novel developments focus on hierarchical Tensor formats as e.g. the Tree-Tucker format, <cit.>, and the Tensor Train format, <cit.>, <cit.>. Nowadays, those hierarchical methods are successfully applied in e.g., physics or chemistry, where they are used in many body problems and quantum states.

Here, we apply  the Tensor Train decomposition to data of a 3D direct numerical simulation (DNS) of a turbulence channel flow. We aim at capturing self-similar structures that might be hidden in the data as explained in the next paragraph. However, as those multiscale flow structures in highly irregular flows are not commonly aligned with the underlying grid but are translated, stretched, and rotated, we, at first, use synthetic data to evaluate the suitability of the method to generally detect self-similarity. 

In wall-bounded turbulent flows characteristic coherent patterns, heterogeneous structures limited in time appearing irregularly at varying positions, are observed in distinct spatial regions <cit.>, <cit.>. A number of experimental and numerical studies were able to describe some (quasi-) coherent features rather accurately, e.g. quasi-streamwise vortical structures are identified as a form of quasi-coherent structures <cit.>, <cit.>.  A primary focus of research is the near-wall region, i.e. y^+<40 (y^+ as wall coordinate), where low-speed streaks are observed moving away from the wall, and, consequently, a flow towards the wall is required for continuity reasons. Indeed, regions of such wall-directed rapidly moving fluid, called sweeps, are observed, e.g. <cit.>, and it has been discussed that these patterns as well as the low-speed streaks positively contribute to the production of turbulent energy, e.g. <cit.>, <cit.>. 

Today, in large eddy simulation (LES) it is still subject of intense research to model the unresolved so-called small scales. One branch follows the idea to reconstruct the fluctuations of the unresolved scales. For example, <cit.> uses a stochastic approach of an adaptive deconvolution model for the subgrid scale closure; <cit.> suggests a variational finite element formulation; <cit.> develops a variational multiscale method where adaption controls the influence of an eddy viscosity model. Self-similarity, however, is not resolved per-se in multiscale models. For this purpose, wavelet decomposition methods have been applied to data of turbulent flows, e.g. <cit.>, and they were recently used to split the coherent signal from the incoherent part, e.g. <cit.>, <cit.>. Similarity models, e.g. <cit.>, and fractal models, e.g. <cit.>, are promising tools as they extrapolate self-similarity to the small scales.

Our study is concerned with the question of whether Tensor decomposition methods can support the development of improved understanding and quantitative characterisation of multiscale behavior of turbulent flows, cf. e.g.<cit.>. Recent high-resolution numerical simulations obviously confirm established theoretic approaches that turbulent flows contain hierarchies of self-similar structures like vortex tubes or sheet-like flow patterns, e.g. <cit.>. Provided that our tests yield promising results, those quantitative features could be helpful in developing a LES closure approach based on and extending the idea of fractal or dynamic SGS models, e.g. <cit.>, <cit.>. Therefore, if proved positively, a long-term goal would be the construction of a self-consistent closure for LES of turbulent flows that explicitly exploits the Tensor decomposition approach's capability of capturing self-similar structures. Our approach is automatically linked with the following questions: (i) Can real data from multiscale dynamics be approximated or represented by the Tensor Train decomposition technique and how compact are the resulting storage schemes, i.e what compression rate can be achieved at which level of accuracy? (ii) Does the Tensor Train approximated data retain the dynamics? (iii) Is the Tensor Train format suitable for detecting cascades of scales in real data and in turbulence data in particular? In this article, we present results of evaluating the Tensor Train format especially to (i). Promising results would encourage us for future work.



§ THE TENSOR TRAIN FORMAT


In this section, we introduce the Tensor Train format, that is the Tensor decomposition method of choice for our study. A number of Tensor formats were developed in the past, e.g. the r-term approximation or so-called canonical format and the Tucker format <cit.>,and we refer to e.g. <cit.> and <cit.> for surveys of low-rank approximation techniques and specifically Tensor approximation formats. Previous studies show that the Tensor Train format allows for a  much higher attainable compression rate compared to the low rank approximation formats just mentioned.

The Tensor Train format, however, is based on the Tucker decomposition. In the Tucker format a Tensor is decomposed into a product of a core Tensor and factor matrices. The factor matrices can be treated as the principal components of the Tensor and the core Tensor implies the level of interaction between the factor matrices. The Tucker format uses the advantages of the higher order singular value decomposition, hereafter HOSVD for short, with a closed set of low rank Tensors. The main disadvantage of the Tucker format is its high storage requirement as the core Tensor takes r^d elements, r as rank of the cores, resulting in a total storage requirement of 𝒪(dnr+r^d). The Tensor Train format on the other hand is designed such that the storage requirement scales linearly with order d: 𝒪(dnr^2), cf. <cit.>. 

The matrix SVD decomposes a matrix into a product of matrices, which then represents the original matrix: A=U S V^T, where A∈ℝ^m× n is the original matrix, U∈ℝ^m× m and V∈ℝ^n× n are orthonormal matrices and S∈ℝ^m× n is a diagonal matrix with the singular values, σ_n, of A on its diagonal. Note that σ_1≥σ_2≥…≥σ_n≥ 0 and the number of non zero singular values is equal to the rank of A: rank(A)=r. If the spectrum of the singular values contains only a few large entries and a sharp cut-off to the remaining tail, considerable data compression rates can be realized by truncating the matrices to the significant part of the spectrum. E.g., an approximation of the original matrix, A, with the first three singular values reads A≈σ_1u_1v_1^T + σ_2u_2v_2^T + σ_3u_3v_3^T. Moreover, a truncation to rank r results in a storage requirement of mr+r+nr instead of mn for the original matrix A.

The key idea of Tensor Train decomposition is to separate a high-dimensional Tensor into its component Tensors (modes). The modes are then Tensors of order 3 by construction. 
Principally, the representation of a d-dimensional Tensor X in the Tensor Train format reads

    X(n_1,…,n_d) = ∑_k_1=1^r_1…∑_k_d-1=1^r_d-1U_1(n_1,k_1) U_2(k_1,n_2,k_2) … U_d(k_d-1,n_d).

The component Tensors U_i are obtained by the step-by-step application of the matrix SVD. The procedure is as follows: in the first step, the d-dimensional Tensor X is reshaped into a 2-dimensional matrix A_1∈ℝ^n_1× n_2… n_d to which a SVD is applied, A_1=U_1S_1V_1^T, and the first mode U_1∈ℝ^n_1× r_1 is obtained. The remaining matrices of the SVD are contracted to one matrix A_2=(S_1V_1^T)∈ℝ^r_1× n_2… n_d which is used for the second step. In the second step, A_2 is reshaped, at first, so that A_2∈ℝ^r_1n_2× n_3… n_d. Then, a SVD is applied, again, and the second mode U_2∈ℝ^r_2× n_3 is obtained. By following this procedure successively, the remaining modes U_3,…,U_d are obtained. Note, that we need to apply the SVD d-times in total. As an example, figure <ref> shows a graphical representation of the resulting Tensor Train network for an input Tensor of order 4.

The Tensor Train decomposition reveals some important topics: (a) in the Tensor Train network, the different modes U_i are linked by the ranks r_i (as sketched in figure <ref>). (b) the parameters r_i determine the approximation quality. In case of Tensor approximation, as opposed to its representation, the r_i are limited to some maximum rank, r, and a lower value of r usually results in lower storage requirement for the approximated Tensor but also results in decreased approximation quality. The ranks r_k, therefore, are also called compression ranks or TT-ranks (TT as Tensor Train). This is because, in the sense of the procedure described previously, the value of r determines the number of rows to be kept for the next step and sets the cut-off of the matrices A_i. The hope is, that the omitted rows will be of small norm and the truncated component Tensors then will represent the significant part of the input Tensor. However, the optimal value of r is not known a priori and it is not guaranteed that the omitted rows are indeed of small norm. That implies consequences for the application of the Tensor Train decomposition to real data as we will see later.
	




§ RESULTS


Here, we present the results of the Tensor Train approximation applied to data. We use the open source Tensor library Xerus developed at the Technical University Berlin, Germany, <cit.>. The relative error of the approximation w.r.t. the original Tensor is measured in the Frobenius norm that, for a d-dimensional Tensor X, reads

    ||X|| =√(∑^N_1_n_1=1∑^N_2_n_2=1⋯∑^N_d_n_d=1 x^2_n_1⋯ n_d),
 
and the relative error in the Frobenius norm then reads

    e = ||(U-X)||/||X||,

with U as approximation of X. 



 §.§ Principal operation of the Tensor Train approximation

We compute a data series of a function, f(x), that involves an overlay of sine-functions of different periods (figure <ref>). The function reads

    f(x)=3.125sin(x)sin(2x)sin(4x),     x∈[0,2π]

Here, f(x) is sampled with 2^14=16384 digits. By construction, f(x) involves similarity but not self-similarity in the proper sense.

Formally, the data series is a Tensor of order 1, i.e. a vector, and X(n_1)∈ℝ with n_1=1,…,2^14. Without loss of generality, we can reshape the 1d-Tensor into a Tensor of order 3, that is T(n_1,n_2,n_3)∈ℝ with T[2,2,4096]. Note, that the choice of the shape of the Tensor T is initially arbitrary but has a significant impact on the result as we will see below.

Applying the Tensor Train decomposition to T, it turns out that the input Tensor is approximated exactly, i.e. represented, for TT-rank ≥ 2, and we find r_1=1 and r_2=2. Figure <ref>a shows a sketch of the approximated Tensor written in the Tensor Train format. The modes U_1 and U_2 read

    U_1=1/√(2)[ -1  1 ]    U_2=1/√(2)[ -1  1;  1  1 ],

and the basis is given by the components of mode U_3∈ℝ^2×4096 (see figure <ref>b). The storage of T in the Tensor Train format requires 8100 entries, which is less than half the storage requirement in the original Tensor notation. Note that, in the limit r_k=1 ∀ k, the Tensor approximation yields a relative error e=0.33 with a storage requirement of 4102 entries.

Now, we perform a test on the sensitivity of the Tensor Train approximation on the shape of the input Tensor. Suppose, we reshape the given data series, (<ref>), into a Tensor T_2 of order 4 where T_2[2,2,2,2048]. Written in the Tensor Train format, TT-rank 4 is needed as minimum TT-rank to represent T_2 (instead of TT-rank 2 for T in the former example). Note, that TT-rank 2 approximates the input Tensor with a relative error e≈0.05. Table <ref> shows the relative error and the storage requirement for TT-rank 1 to 4.   









 §.§ Detection of self-similar structures

Now, we extent the first example and test the Tensor Train format to detect self-similar structures hidden in data series. For this purpose, we compute a data series that consists of a 2^4=16 digit-sequence of triangles. Thus, it is in line with the power of two sequence. The 16 digit sequence is not only repeated but is also scaled-up in every loop, see figure <ref>. The data series comprises 2^10=1024 digits in total. As shown previously, the impact of the shape of the input Tensor on the resulting storage requirement can be significant. Therefore, we reshape the data series in a Tensor of order 9 where each of the 9 dimensions have 2 entries except the second to last one that has 4 entries, i.e. T_3[2,2,2,2,2,2,2,4,2]. Thus, the self-similar structure in the data series is maintained and, furthermore, remains grid-aligned and is not disturbed by splitting of the Tensor. We find TT-rank 2 as minimum rank to represent the input Tensor in the Tensor Train format. The storage requirement, 66 entries (≈6.4 % of the original Tensor notation), is very low. Variations of the input Tensor's shape results in very low storage requirements, too, but it turns out that the shape of T_3 demands the lowest storage requirement and the lowest TT-rank. Table <ref> shows results for varieties of the input Tensor entries.




 

Finally, we consider the properties of the Tensor Train format in case the input Tensor involves self-similar structures not aligned with the grid, that is, the self-similar patterns are not in line with the power of two sequence. For that purpose we compute a top-hat function where the squares are scaled up, see figure <ref>a. The data series contains 2^14=16384 entries in total. To test the effect of noisy data, we perform two additional tests and add random noise of different amplitude, here amplitude factor 1 and 10, to the data series. We reshape the data series into a Tensor of order 14, thus each dimension has 2 entries, T_4[2,2,2,2,2,2,2,2,2,2,2,2,2,2]. With this arbitrary choice, we intentionally neglect our a priori knowledge on the character of the data series. This is reasonable as we can not expect detailed knowledge of the data structure in case of real data.

Figure <ref>b shows the storage requirement (hereafter denoted as datasize) over the relative error for both, the data series without noise and with noise amplitude factor 1 and 10, for a number of TT-ranks. In case of no noise Tensor T_4 is represented with TT-rank 5 where relative error is about 10^-15, and storage in the Tensor Train format requires about 500 entries thus ≈3.0 % of the original datasize. Limiting the TT-rank to 2 (4) yields approximation error of about 0.2 (0.02) %, and storage of the approximated Tensor in the Tensor Train format requires about 100 (350) entries, i.e. 0.6 (2.1) % of the datasize in original Tensor notation. Approximation of the noisy data series for TT-rank 5 to about 80 shows a drastic increase in the storage requirement whilst only small changes in the relative error. Obviously, this is due to the small-scale noise the approximation of which naturally requires many more data entries.

Both data series, with noise amplitude 1 and with noise amplitude 10, are approximated exactly at TT-rank 128 at which the relative error drops to 10^-15. In both cases, the storage requirement of the representations in the Tensor Train format (43690 entries for both cases) substantially exceeds the storage requirement of the source data in the original Tensor format (16384 entries). For both noisy data series, we find the original datasize at TT-rank 45 where e≈0.03 (noise factor 10) and e≈0.003 (noise factor 1), resp. Furthermore, for a given datasize, the relative error strongly depends on the noise amplitude as long as the noise is not well approximated.

Interestingly, the relative error agrees well with that of the no-noise case up to TT-rank 4, TT-rank 3 for noise amplitude 10, indicating that the self-similar top-hat structure, large-scaled relative to the small scale noise, but not the noise itself is represented by low TT-ranks. Figure <ref>c shows the euclidian norm of the difference between the original Tensor and the approximated Tensor. At low TT-ranks, the euclidian norm agrees well in all cases underlining the previous statement that large-scale structures are resolved at low TT-ranks. At large TT-ranks, the graphs of both noisy data series show a significant drop at TT-rank 127 specifying the exact approximation of the input Tensor.  

Figure <ref>d confirms that picture since the self-similar top-hat structure but not the noisy part (small-scaled w.r.t. the top-hat structure) is approximated well at TT-rank 4 (except the small-scale top-hats at the beginning of the data series). With TT-rank 2 as maximum permitted TT-rank, the approximation of the top-hat structures is less accurate but the general trend is obtained, at least. That result is of course linked with the mathematics of the algorithm discussed previously since the given TT-rank is linked with the number of rows, i.e. number of singular values, to be kept in the Tensor Train decomposition procedure (cf. section <ref>). Low rank, therefore, means a loss of information particularly about small-scales, that is usually hidden in higher eigenvalues of the matrices. In other words, large-scale structures are already being represented whilst small-scale patterns are still coarsely approximated.

In summary, the tests with different synthetic data series show that approximation in the Tensor Train format is efficient and promising in detecting self-similar patterns in particular if those structures are grid-aligned. Even if the self-similar structures are not grid-aligned, the approximation in the Tensor Train format yields satisfying results particularly for data series without noise. Furthermore, we found that the minimal TT-rank to represent a data series in the Tensor Train format is sensitive to the shape of the input Tensor. The latter statement has been also discussed by <cit.>.





 §.§ Channel turbulence flow

We use data of 3D direct numerical simulation (DNS) - by using a pseudo-spectral Fourier-Chebyshev method- of an incompressible, isothermal plane channel flow, figure <ref>a,  computed by <cit.>, <cit.>, to which we refer for an in-depth description of data generation. 

Data are generated at Re_τ=590, where Re_τ is the friction-based Reynolds number. The channel geometric factor h has been set to h=1.
The original grid spatial resolution is 600× 385× 600 in (x,y,z). Here, we make not use of data of the original grid but of the so-called fine grid that is 600× 352× 600 in (x,y,z), and we have applied a post-processing procedure to obtain a constant grid increment throughout the y-direction instead of the polynomial distribution of the original grid.  We refer to <cit.> for a detailed description of the post-processing algorithms that convert the original data to the equidistant fine grid data.

In wall units the grid increment of the fine grid is Δ y^+=3.37. Compared to the signatures of the original grid data, e.g. fig. 3 and 4 in <cit.>, no quantitative loss of information is recognized due to the fine grid post-processing procedure. However, the resolution of the near-all region is significantly reduced. In particular, the viscous sublayer (y^+<5) is not resolved anymore.

The scalar q-criterion proposed by <cit.>, <cit.> is a measure to demarcate vortex structures within turbulent flows. The q-criterion reads

    q = 1/2 (|| A||^2 - || S||^2),

where || S|| = [tr(SS^t)]^1/2, || A|| = [tr(AA^t)]^1/2, with S and A as the symmetric and anti-symmetric  part of the velocity gradient Tensor ∇ U, resp. Figure <ref>d shows an iso-surface map of the q-criterion for a given snapshot. Here, vortex tubes of different size and shape, rotated and stretched, can be identified indicating the highly turbulent flow.

Here, we apply the Tensor Train decomposition to velocity data (fig. <ref>b), the datasize is 380160000≈ 10^9 entries per snapshot. 


Since the tests with synthetic data series already have revealed the sensitivity to the shape of the input Tensor, see the previous subsection, we consider two different approaches. Firstly, the input Tensor retains the original data structure, that is, it is a 4th-order Tensor T[600,352,600,3] where the last dimension contains the components of velocity. This approach is referred to as TT-approximation hereafter. On the other hand, we decompose every dimension into its prime factors, resulting in a Tensor of order 19

    T[n_1,…,n_19]=T[2,2,2,3,5,5, 2,2,2,2,2,11, 2,2,2,3,5,5, 3].

 The latter approach is referred to as QTT-approximation hereafter as it is similar to the Quantics Tensor Train (QTT) approach, <cit.>, <cit.>, inasmuch as the binary representation is realized as well as possible.

Figure <ref> shows the results of the Tensor Train decomposition for both approaches. Here, the diagram shows storage requirement over relative error for various TT-ranks.  Obviously, the QTT-approximation yields less relative error with significantly less storage requirement. For example, given a relative error of 0.085 the QTT-approximation demands a storage requirement of about 18000 entries but TT-approximation demands about 60000 entries. Induced by the different shapes of the input Tensors the necessary TT-rank is much larger in the QTT-approach.

Figures <ref> for the TT-approximation and figure <ref> for the QTT-approximation displays 2D-slices of the approximated data for selected TT-ranks, and also a 2D-slice of the original data is shown for comparison. Figure <ref> displays the magnitude of velocity computed from the components after TT-approximation and figure <ref> shows the streamwise component of the velocity. In addition, 1D cross-sections along x-coordinate are given in figure <ref>. Generally, images of lower TT-ranks show a smooth reflection of the turbulent flow, and with increasing TT-rank small-scale patterns are mapped with enhanced resolution. Interestingly, in the TT-approximation, TT-rank 2 involves the box-shaped mean velocity profile typical for channel turbulence flows, figure <ref>e, cf. e.g. <cit.>. In that sense, Tensor Train decomposition restricted to TT-rank 2 apparently acts as an averaging process.

Now, we are interested in identifying those regions in the channel flow that show the most significant deviations from the original in the approximated data. In both cases, TT-approximation and QTT-approximation, reconstruction of the velocity field for higher TT-ranks, for example TT-rank 100 in figure <ref>c and <ref>b, already shows several details of small-scale features of the turbulent flow by eye, even if the relative error is significant. For the QTT-approximation, figure <ref>e for TT-rank 20 and <ref>f for TT-rank 100 shows the difference between the approximated and the original data. Obviously, significant deviations are located towards the near-wall region and the flow interior is mostly unaffected. This picture is becoming more and more visible the more details of the small-scale features are resolved, thus with increasing TT-rank.
 








§ CONCLUDING REMARKS


In this article, the Tensor Train decomposition technique, a specific branch of the family of Tensor product decomposition methods, has been applied to DNS data of channel turbulence flow. We showed that the Tensor Train format yields significant compression rates whilst ensuring low relative errors. In particular, the QTT-approximation, that in our definition is decomposing the entries of each dimension of the input Tensor into its prime factors, shows considerably better results than the approximation of the original data Tensor of order 4.

Applied to channel turbulence data, the Tensor Train format allows for surprisingly high compression rates. QTT-approximation at TT-rank 100 results in a compression factor of roughly 1000 and in a relative error of about 5 %. At TT-rank 2000, the relative error is about 0.2 % with a compression factor of about 5, that is a storage requirement of ≈10^8 (instead of ≈10^9 for the original data). However, it demonstrates that low-rank representation of those highly irregular and chaotic flows in the Tensor Train format can not be expected.

Our study is concerned with the detection of and quantitative characterization of self-similar patterns in turbulent flow data. Therefore, we, firstly, applied synthetic data to demonstrate the Tensor Train decomposition´s capability to capture self-similarities. The synthetic data tests with and without noise provide promising results. Grid aligned self-similarity is well captured and non grid-aligned self-similarity is approximated at low TT-ranks strengthening the suitability of the method. Moreover, with no noise, non grid-aligned self-similar structures are approximated exactly at low TT-ranks.


This research has been funded by Deutsche Forschungsgemeinschaft (DFG) through grant CRC 1114  'Scaling Cascades in Complex Systems', Project B04 'Multiscale Tensor decomposition methods for partial differential equations'. The authors thank Prof. Illia Horenko (CRC 1114 Mercator Fellow) as well as Prof. Reinhold Schneider and Prof. Harry Yserentant for rich discussions and for steady support. Thomas von Larcher thanks Sebastian Wolf and Benjamin Huber (both at TU Berlin, Germany) very much for developing the Tensor library xerus which has been used for data analysis, as well as for their round the clock support in the project. The data were generated and processed using resources of the North-German Supercomputing Alliance (HLRN), Germany, and of the Department of Mathematics and Computer Science, Freie Universität Berlin, Germany. The authors thank Alexander Kuhn and Christian Hege (both at Zuse Institute Berlin, Germany) for steady support in data processing and data visualisation. 


spmpsci 
 

