Complexity Estimation of Cubical Tensor Represented through 3D Frequency-Ordered Hierarchical KLT

: In this work is introduced one new hierarchical decomposition for cubical tensor of size 2 n , based on the well-known orthogonal transforms Principal Component Analysis and Karhunen–Loeve Transform. The decomposition is called 3D Frequency-Ordered Hierarchical KLT (3D-FOHKLT). It is separable, and its calculation is based on the one-dimensional Frequency-Ordered Hierarchical KLT (1D-FOHKLT) applied on a sequence of matrices. The transform matrix is the product of n sparse matrices, symmetrical at the point of their main diagonal. In particular, for the case in which the angles which deﬁne the transform coe ﬃ cients for the couples of matrices in each hierarchical level of 1D-FOHKLT are equal to π / 4, the transform coincides with this of the frequency-ordered 1D Walsh–Hadamard. Compared to the hierarchical decompositions of Tucker (H-Tucker) and the Tensor-Train (TT), the o ﬀ ered approach does not ensure full decorrelation between its components, but is close to the maximum. On the other hand, the evaluation of the computational complexity (CC) of the new decomposition proves that it is lower than that of the above-mentioned similar approaches. In correspondence with the comparison results for H-Tucker and TT, the CC decreases fast together with the increase of the hierarchical levels’ number, n. An additional advantage of 3D-FOHKLT is that it is based on the use of operations of low complexity, while the similar famous decompositions need large numbers of iterations to achieve the coveted accuracy.


Introduction
The well-known tensor decompositions Canonical Polyadic Decomposition (CPD), Higher-Order Singular Value Decomposition (HOSVD) [1][2][3], Tensor Train (TT) Decomposition [4], Hierarchical Tucker decomposition (H-Tucker) [5] and some of their modifications [6,7] are based on the calculation of the eigenvalues and eigenvestors of the decomposed tensor. The basic attribute of this group is that they are optimal regarding the minimization of the mean square approximation error derived from the "truncation" of the low-energy components. For the calculation of the retained components, we use iterative methods which need relatively large numbers of mathematical operations to achieve the requested accuracy (Power method [8], Jacoby method [9], etc.).
The hierarchical tensor decompositions based on the H-Tucker algorithm are generalized in publications [7] and [10]. In [7], Vasilescu and Kim introduced a compositional hierarchical tensor factorization that disentangles the hierarchical causal structure of object image formation, but the computational complexity (CC) of the method is not evaluated. In [10], Zniyed et al. offered the TT-based hierarchical decomposition of high-order tensors, for the calculation of which is used the Tensor-Train Hierarchical SVD (TT-HSVD). This approach permits parallel processing, which accelerates the calculations. Unlike the TT-SVD algorithm, the TT-HSVD is based on applying SVDs to matrices of smaller dimensions. As a result, the CC of TT-SVD is O(R 2 N 3 ) operations higher than that of TT-HSVD (here, R is the rank of the hypercubical tensor of size N).
As an alternative, in this work are presented new hierarchical 3D tensor decompositions, based on the famous statistical orthogonal transforms Principal Component Analysis (PCA) and the Karhunen-Loeve Transform (KLT) [11,12]. The offered decompositions are close to the optimal (which ensure full decorrelation of their components), but do not need iterations and have low computational complexity.
In Section 2, we present the methods for one-dimensional hierarchical adaptive transform based on the PCA (1D-HAPCA) and the Hierarchical Adaptive Karhunen-Loeve Transform (1D-HAKLT). In Section 3, on the basis of these 1D transforms we give the details of the cubical tensor decomposition through separable 3D Frequency-Ordered Hierarchical KLT (3D-FOHKLT), and Section 4 contains the related algorithm. Section 5 presents the comparative analysis of the computational complexity of the new approaches with respect to the well-known hierarchical Tucker decomposition (H-Tucker) and the TT decomposition. Section 6 is the Conclusion.

One-Dimensional Hierarchical Adaptive PCA/KLT Transform for a Sequence of Matrices
The detailed description of the algorithm for the one-dimensional adaptive PCA/KLT transform (1D-APCA/AKLT) for a couple of matrices is given below. This is the basis of the hierarchical algorithm for processing the sequence of matrices which represent the corresponding tensor.
Throughout this paper, scalars are represented as x, vectors and matrices as x and X correspondingly, and tensors as X. It is also assumed here that all tensor elements (entries) are nonnegative real numbers x(i,j,k).

One-Dimensional Hierarchical Adaptive KLT for a Sequence of Matrices
The relations above, used for the calculation of 1D-APCA [12], are the kernel for the processing of a sequence of matrices (for example, a sequence of moving images), and do not need iterative calculations. The procedure is executed in the following way. First, the sequence of matrices X k , each of size N × N for k = 1, 2, . . . , N is divided into N/2 couples (N = 2 n ). For each couple of matrices with elements x 1s (p,q) and x 2s (p,q), when q = 1, 2, . . . , N/2, in the level p = 1, 2, . . . , n, is executed the one-dimensional adaptive PCA (1D-APCA), in accordance with Equations (1)- (5). For m 1 (p, q) = m 2 (p, q) = 0 (valid for the KLT), the direct 1D-AKLT for the couple of matrices q in the level p of the 1D-HAKLT is defined by the relations y 1s (p, q) = x 1s (p, q) cos θ p,q + x 2s (p, q) sin θ p,q y 2s (p, q) = x 1s (p, q) sin θ p,q − x 2s (p, q) cos θ p,q for s = 1, 2, . . . , N 2 , where s is the sequential number of two elements of the same position in the couple of matrices, q, 1s (p, q));k 2 (p, q) = E(x 2 2s (p, q));k 3 (p, q) = E(x 1s (p, q)x 2s (p, q)).
In the equations, figures and text below are introduced the following symbols: cos θ p,q = c p,q and sin θ p,q = s p,q for p = 1, 2, 3 and q = 1, 2, 3, 4. The execution graph for 1D-HAKLT of three hierarchical levels (N = 8) is shown on Figure 1a. In every level p = 1, 2, 3, and for each couple of neighbor matrices q (fenced in red ellipses) 1D-APCA is executed. In result, the first transformed matrix calculated for each couple q has higher energy than the second matrix, and they are mutually decorrelated. To continue with the operations in the next level of 1D-HAKLT, these matrices are reordered, following the decrease of their energy, i.e., where P Y p,k = [y p,k (i, j)] 2 is the energy of matrices Y p,k with elements y p.k (i, j), obtained in result of the execution of 1D-HAKLT in levels p = 1, 2, 3, before the reordering. In result, the reordered matrices E r in the last level are decorrelated to the highest degree [12]. The reordering in each level p in accordance with Equation (9) needs significant computational resources. For the case when the connections between the outputs Yp,k for k = 1, 2, …, 8 in the level p and the inputs in the next level p + 1 are fixed, reordering is not performed, and the 1D-HAKLT turns into the one-dimensional Frequency-Ordered Hierarchical KLT (1D-FOHKLT), whose execution graph is shown in Figure 1b. In this case, only the output matrices in the last (third) level are reordered, which does not require mathematical operations. Besides, the energy distribution of matrices Yp,k is close to that defined by Equation (9). In particular, in accordance with Appendix A, if for p = 1, 2, 3 and q = 1, 2, 3, 4, then 1D-HAKLT coincides with the one-dimensional frequency-ordered Fast Walsh-Hadamard Transform (1D-FWHT) for N = 8 [14]. The reordering in each level p in accordance with Equation (9) needs significant computational resources. For the case when the connections between the outputs Y p,k for k = 1, 2, . . . , 8 in the level p and the inputs in the next level p + 1 are fixed, reordering is not performed, and the 1D-HAKLT turns into the one-dimensional Frequency-Ordered Hierarchical KLT (1D-FOHKLT), whose execution graph is shown in Figure 1b. In this case, only the output matrices in the last (third) level are reordered, which does not require mathematical operations. Besides, the energy distribution of matrices Y p,k is close to that defined by Equation (9). In particular, in accordance with Appendix A, if θ p,q = π/4 for p = 1, 2, 3 and q = 1, 2, 3, 4, then 1D-HAKLT coincides with the one-dimensional frequency-ordered Fast Walsh-Hadamard Transform (1D-FWHT) for N = 8 [14].
Here, x s and e s are respectively the input and the output column-vector of N components correspondingly, and the column-vectors y p,s for p = 2, 3, . . . , n − 1 are the intermediate results.
The relation between vectors e s and x s , defined on the basis of Equations (10)- (12), is given below: where is the matrix for the fast n-level 1D-HAKLT. The sparse transform matrices G p (2 n ) of size 2 n × 2 n have a symmetrical block-diagonal structure towards the main diagonal, and are defined by the relations where T p,q (2) = c p,q s p,q s p,q −c p,q , and ⊕ denotes the direct sum of matrices. The permutation matrices P p (2 n ) of size 2 n × 2 n for level p = 1, 2, . . . , n of 1D-HAKLT are defined on the basis of Equation (9) for k = 0, 1, 2, . . . , 2 n−1 .
The algorithm 1D-FOHKLT is depicted by a relation similar to Equation (13), but the reordering is executed only in the last level, n; where G p (2 n )] is the matrix for 1D-FOHKLT, P n (2 n ) is the permutation matrix P n (2 n ) for the last level n, and n p=1 G p (2 n is the product of n sparse transform matrices G p (2 n ) for p = 1, 2, 3, . . . , n. Each matrix G p (2 n ) is defined as follows: In the level n are rearranged the components of the column-vectors y n,s = G n (2 n )y n−1,s , and, respectively, the components of matrices Y n,k for k = 0, 1, . . . , N − 1. For this is used the permutation matrix P n (2 n ). From the components of the column-vectors e s = P n (2 n )y n,s are obtained the matrices E 0 , E 1 , . . . , E N−1 . The frequency-ordered matrices E r are calculated in accordance with the relation between their sequential number r and the corresponding sequential number k, for the matrices Y n,k . The relation which defines the matrix P n (2 n ) is as follows [14]: • the binary code k n−1 , k n−2 , . . . , k 0 of the sequential decimal number k = 0, 1, . . . , 2 n−1 of the component Y n,k is arranged inversely (i.e., k 0 , k 1 , . . . , k n−1 ), as follows: • the so-obtained code g n−1 , g n−2 , . . . , g 0 is transformed from Gray code into binary code r n−1 , r n−2 , . . . , r 0 , in accordance with the operations Here, ⊕ is the symbol for the operation "exclusive OR". The decimal number r = An example for the direct 1D-FOHKLT of three hierarchical levels is given in Appendix A.

Decomposition of Cubical Tensor of Size N × N × N through 3D-FOHKLT
The decomposition 3D-FOHKLT for a cubical tensor X of size 8 × 8 × 8, executed in three consecutive stages, is shown on Figure 2.
In the level n are rearranged the components of the column-vectors is as follows [14]: k ,k ,...,k − ), as follows: • the so-obtained code Here, ⊕ is the symbol for the operation "exclusive OR". The decimal number An example for the direct 1D-FOHKLT of three hierarchical levels is given in Appendix A.

Decomposition of Cubical Tensor of Size N × N × N through 3D-FOHKLT
The decomposition 3D-FOHKLT for a cubical tensor X of size 8×8×8, executed in three consecutive stages, is shown on Figure 2. In the first stage, the tensor X is divided into eight horizontal slices (mode-1). Then, in accordance with Figure 1, on each couple of matrices is applied 1D-FOHKLT. In result, the sequence of eight In the first stage, the tensor X is divided into eight horizontal slices (mode-1). Then, in accordance with Figure 1, on each couple of matrices is applied 1D-FOHKLT. In result, the sequence of eight transformed matrices is united into one new tensor, E. In the second stage, the input tensor E is divided into eight lateral slices (mode-2), then the couples of matrices are transformed, and after that united into the tensor F . In the third stage, the input tensor F is divided into eight frontal slices (mode-3). After applying the 1D-AKLT on each couple of matrices which build the tensor F , the final result is obtained-the general spectral tensor S, with elements s(u, v, l). The initial tensor X is restored from the tensor S through summation of all "weighted" basic tensors K u,v,l of size 8 × 8 × 8, in accordance with the relation Here, s(u, v, l) are the weight coefficients of the basic tensor, K u,v,l . Each basic tensor is represented as the outer product of three vectors: where o is the symbol for the outer product of the two column-vectors (x • y = x.y T ) and k l (t) = [k t l,1 , k t l,2 , k t l,3 , k t l,4 , k t l,5 , k t l,6 , k t l,7 , k t l,8 ] T .
The general number of the basic tensors in this case is 8 3 . For example, if u = 3, v = 2 and l = 6, the components of these vectors for each stage are defined by the relations below: The vectors k u (t),k v (t), k l (t) are orthonormalized, i.e., for their components are valid the relations: 1 for t 1 = t 2 , where u 1 , u 2 = 0, 1, . . . , 7 and t 1 , t 2 = 1, 2, 3.
The decomposition 3D-HAKLT, shown in Figure 2, could be easily generalized for the cases when N = 16, 32, . . . , 2 n through uniting the computational principles from Figures 1 and 2 for decompositions of high numbers of levels (n = 4, 5, . . . ). From Equation (29), it follows that the rows of the matrices T t 1 FOHKLT (8) and T t 2 FOHKLT (8) for t 1 t 2 are not correlated. This conclusion could also be generalized for the rows of the matrices T t 1 FOHKLT (2 n ) and T t 2 FOHKLT (2 n ). Hence, after the execution of each stage of 3D-FOHKLT, the so-calculated N slices in the corresponding direction are highly decorrelated [12,15].
In the general case, the spectral coefficients s(u,v,l) are calculated through direct 3D-KLT, in accordance with the relation Each third-order tensor X of size N × N × N is represented as the weighted sum of N 3 3D-KL functions obtained by using the tensors K u,v,l , each of size N × N × N: Then, the tensor decomposition based on the 3D-FOHKLT is defined by the relation Hence, in this case the tensor X of size N × N × N is decomposed in a way similar to that of the Tucker decomposition [2]; Here, d r 1 ,r 2 ,r 3 are the eigen values of the tensor X, a r 1 , b r 2 , c r 3 are the corresponding eigen vectors, and R is the tensor rank. The rank of the three-way tensor X of size I × J × K is bounded from above in accordance with the inequality, as follows: max{I, J, K} ≤ R ≤ min{I × J, I × K, J × K} [3]. If I = J = K = N, then N ≤ R ≤ N 2 . The difference between the decompositions from Equations (33) and (34) is the number of their components, which can reach up to N 3 and N 2 , respectively (for the maximum value of the rank, R). In the case that the decomposition components number as defined by Equation (33) is not reduced, the 3D-HAKLT is reversible, i.e., the coefficients s(u,v,l) of the tensor S and its basic tensors K u,v,l are sufficient to restore the tensor X.

Algorithm 3D-FOHKLT for Cubical Tensor Representation
The representation of the tensor X of size N × N × N, through 3D-FOHKLT based on the spectral tensor S, is executed in correspondence with Figure 2. Each stage of the 3D transform is calculated for N = 8 through 1D-FOHKLT, as shown in Figure 1b.
In result is obtained the spectral tensor S, which comprises n layers of coefficients. These coefficients participate as component weights in the decomposition represented by Equation (33).

Algorithm 1: Tensor representation based on the 3D-FOHKLT
Input: Third-order tensor X of size N × N × N (N = 2 n ) and with elements x(i,j,k) Output: Spectral tensor S of n layers, whose elements s(u,v,l) are the coefficients of the 3D-FOНKLT. 1 begin 2 Divide the tensor X into horisontal slices and compose N/2 couples of matrices X 1 p,q and X 2 p,q of size N × N, for q = 1, 2, . . . , N/2; 3 for each couple of matrices X 1 p,q and X 2 p,q when q = 1, 2, . . . , N/2 and p = 1, 2, . . . , n, do 4 Calculate the couple of matrices Y 1 p,q and Y 2 p,q obtained through 1D-FOНKLT in accordance with Equations (7) and (8); 5 Calculate the matrices P n (2 n ) and Y n,k for k = 0, 1, . . . , N−1 and p = n in accordance with Equations (18)-(20); 6 Define the matrices E r for r = 0, 1, . . . , N−1 in the level p = n by rearranging the matrices Y n,k for k = 0, 1, . . . , N−1 on the basis of the vector transform e s = P n (2 n )y n,s , for s = 1, 2, . . . , N 2 ; 7 Reshape all rearranged matrices E r into the corresponding tensor E of size N × N × N. 8 end 9 Divide the tensor E into lateral slices and compose N/2 couples of intermediate matrices E 1 p,q and E 2 p,q of size N × N, for q = 1, 2, . . . , N/2; 10 for each couple of matrices E 1 p,q and E 2 p,q when q = 1, 2, . . . , N/2 and p = 1, 2, . . . , n, do 11 Calculate the couple of matrices F 1 p,q and F 2 p,q transformed through 1D-FOНKLT in accordance with Equations (7) and (8); 12 Calculate the matrices P n (2 n ) and F n,k for k = 1, 2, . . . , N and p = n in accordance with Equations (18)-(20); 13 Define the matrices F r for r = 0, 1, . . . , N−1 in the level p = n by rearranging the matrices F n,k for k = 0, 1, . . . , N−1 on the basis of the vector transform F s = P n (2 n )f n,s , for s = 1, 2, . . . , N 2 ; 14 Reshape all rearranged matrices F r in level p = n into the corresponding tensor F of size N × N × N 15 end 16 Divide the tensor F into frontal slices and compose N/2 couples of matrices F 1 p,q and F 2 p,q of size N × N, for q = 1, 2, . . . , N/2; 17 for each couple of matrices F 1 p,q and F 2 p,q when q = 1, 2, . . . , N/2 and p = 1, 2, . . . , n, do 18 Calculate the couple of matrices S 1 p,q and S 2 p,q transformed through 1D-FOНKLT in accordance with Equations (7) and (8); 19 Calculate the matrices P n (2 n ) and S n,k for k = 1, 2, . . . , N and p = n in accordance with Equations (18)-(20); 20 Define the matrices S r for r = 0, 1, . . . , N−1 in the level p = n by rearranging the matrices S n,k for k = 0, 1, . . . , N−1 on the basis of the vector transform S s = P n (2 n )s n,s , for s = 1, 2, . . . , N 2 ; 21 Reshape all rearranged matrices S r in the level p = n into the spectral tensor S of size N × N × N 22 end 23 Arrange the coefficients s(u,v,l) of the spectral cubical tensor S layer by layer in accordance with their spatial frequencies (u,v,l) -from the lowest (0,0,0), to the highest (N−1, N−1, N−1), for u + v + l = const. 24 end Here, the angle θ p,q is defined by Equation (2), on the basis of the covariance matrix elements K(p, q) = k 1 (p, q) k 3 (p, q) k 3 (p, q) k 2 (p, q) for the couple of matrices X 1 p,q and X 2 p,q (in levels p = 1, 2, 3 for q = 1, 2, 3, 4).
Unlike 1D-FOHKLT, the decorrelation of the cubical tensor elements through 3D-FOHKLT is executed in three mutually orthogonal directions, which is its main advantage compared to 1D-FOHKLT.
The algorithm for 3D-HAHKLT is similar to Algorithm 1, given above. The only difference is that the reordering of the matrices in each consecutive level is executed in accordance with Equation (9).

Comparative Evaluation of the Computational Complexity of 3D-HAKLT and 3D-FOHKLT
The computational complexity is evaluated on the basis of the needed operations "addition" and "multiplication".

Comparison of the Computational Complexity of 3D-HAKLT and 3D-FOHKLT to H-Tucker and TT
The CC of 3D-HAKLT for the transformation of the cubical tensor X of size 2 n , evaluated through the needed number of operations, is O 3H (24 × 2 3n n). In accordance with [5], the CC of the H-Tucker decomposition is O HT ((d − 1)R 3 + d × 2 n × R), where R is the rank, 2 n is the size, and d = 3 is the order of the tensor X. For R = 2 n (the minimum value in the range 2 n ≤ R ≤ 2 2n ), the number of the operations is O HT (2 4n+1 + 3 × 2 3n ) = O HT (2 3n (2 n+1 + 3)). In this case, the relative CCs of 3D-HAKLT and 3D-FOHKLT compared to that of the H-Tucker is defined by the relations: The CC of the TT decomposition [4] for the cubical tensor X of size 2 n is O TT (2 n dR 3 ), i.e., O TT (3 × 2 4n ). Correspondingly, the relative CC of 3D-HAKLT and 3D-FOHKLT towards the TT decomposition is defined by the relations: The values of functions ψ i (n) for i = 1, 2, 3, 4 and n = 2, 3, . . . , 8 are given in Table 1, and their graphics are shown in Figure 3a,b. From Equations (47)-(50), it follows that the functions ψ i (n) for i = 1, 2, 3, 4 grow proportionally to 2 n /n. The obtained results show that for small values of n (in the range from 2 to 5), the decompositions 3D-HAKLT and 3D-FOHKLT have low efficiencies because their CCs are equal or higher than those of H-Tucker and TT.
However, for practical applications (image processing), these values are quite different (for example, n = 8 corresponds to tensor image of size 256 × 256 × 256, which is very small), and in this case the relative CC of the offered decompositions is much lower than those of H-Tucker and TT. An additional advantage is that unlike H-Tucker and TT, the proposed algorithms do not need iterative calculations.  However, for practical applications (image processing), these values are quite different (for example, n = 8 corresponds to tensor image of size 256 × 256 × 256, which is very small), and in this case the relative CC of the offered decompositions is much lower than those of H-Tucker and TT. An additional advantage is that unlike H-Tucker and TT, the proposed algorithms do not need iterative calculations.

Conclusions
In this work are presented new cubical tensor decompositions based on the 3D orthogonal hierarchical transforms HAKLT and FOHKLT. The comparative analysis for the relative CC of the new decompositions compared to those of H-Tucker and TT shows their advantages for tensors of size 256 and more. From these results it follows that 3D-FOHKLT has minimum CC, and is suitable for real-time applications. One more advantage of the offered algorithms is that they do not need iterative calculations, and instead use hierarchical calculations in the three mutually orthogonal directions of the tensor. Unlike the well-known 3D decompositions based on the N-dimensional eigen vectors of the tensor, the kernel of the offered hierarchical transforms is the AKLT of size 2 × 2, used repeatedly in each level. The presented algorithms 3D-HAKLT and 3D-FOHKLT could also be generalized for tensors with three different dimensions

Conclusions
In this work are presented new cubical tensor decompositions based on the 3D orthogonal hierarchical transforms HAKLT and FOHKLT. The comparative analysis for the relative CC of the new decompositions compared to those of H-Tucker and TT shows their advantages for tensors of size 256 and more. From these results it follows that 3D-FOHKLT has minimum CC, and is suitable for real-time applications. One more advantage of the offered algorithms is that they do not need iterative calculations, and instead use hierarchical calculations in the three mutually orthogonal directions of the tensor. Unlike the well-known 3D decompositions based on the N-dimensional eigen vectors of the tensor, the kernel of the offered hierarchical transforms is the AKLT of size 2 × 2, used repeatedly in each level. The presented algorithms 3D-HAKLT and 3D-FOHKLT could also be generalized for tensors with three different dimensions 2 n 1 × 2 n 2 × 2 n 3 for n 1 n 2 n 3 . The choice of the offered hierarchical 3D decompositions is made depending on the requirements and limitations of their CC imposed by the corresponding application area.
The future investigations of the 3D-HAKLT and 3D-FOHKLT will be aimed at the evaluation of their characteristics compared to the famous tensor decompositions, in order to define the best settings and to outline the most efficient applications in the 3D compression, filtration, analysis, search and recognition of multidimensional visual information, deep learning, etc.

(A9)
From this, it follows that k 1 .k 2 = 0 only if the conditions s 2,1 = c 2,1 and s 2,3 = c 2,3 are satisfied. In particular, if for a given input tensor X in the three stages of 3D-FOHKLT the angles θ p.q for p = 1, 2, 3 and q = 1, 2, 3, 4 are equal or close to π/4, then in correspondence with (A6) the matrix T π/4 FOHKLT (8) coincides with the frequency-ordered matrix of the Wash-Hadamard Transform T WHT (8) [14,16], i.e.,: . (A10) This property of the matrix T WHT (8) shows that when it is used, is got a frequency-ordered spectrum which corresponds to that obtained with the Fast Fourier Transform (FFT). Hence, the FOHKLT spectrum is frequency-ordered. In the cases when the angles θ p,q are equal or close to 0 or π/2, the matrix T FOHKLT (8) is transformed into the following two matrices: (A12) These matrices are not transform matrices: T 0 FOHKLT (8) is used to rearrange the components of the input vector and to change the signs in half of them, and the matrix T π/2 FOHKLT (8) only rearranges the input vector components.
The description given above, could be easily extended for N = 2 n on the basis of Equations (10)-(20).