Hierarchical Cubical Tensor Decomposition through Low Complexity Orthogonal Transforms

: In this work, new approaches are proposed for the 3D decomposition of a cubical tensor of size N × N × N for N = 2 n through hierarchical deterministic orthogonal transforms with low computational complexity, whose kernels are based on the Walsh-Hadamard Transform (WHT) and the Complex Hadamard Transform (CHT). On the basis of the symmetrical properties of the real and complex Walsh-Hadamard matrices are developed fast computational algorithms whose computational complexity is compared with that of the famous deterministic transforms: the 3D Fast Fourier Transform, the 3D Discrete Wavelet Transform and the statistical Hierarchical Tucker decomposition. The comparison results show the lower computational complexity of the o ﬀ ered algorithms. Additionally, they ensure the high energy concentration of the original tensor into a small number of coe ﬃ cients of the so calculated transformed spectrum tensor. The main advantage of the proposed algorithms is the reduction of the needed calculations due to the low number of hierarchical levels compared to the signiﬁcant number of iterations needed to achieve the required decomposition accuracy based on the statistical methods. The choice of the 3D hierarchical decomposition is deﬁned by the requirements and limitations related to the corresponding application area. decomposition of 3D cubical tensor; hierarchical 3D Fast Walsh-Hadamard Transform (3D-FWHT); 3D Hierarchical Fast Complex Walsh-Hadamard Transform


Introduction
The famous tensor decompositions-Canonical Polyadic Decomposition (CPD), Higher-Order Singular Value Decomposition (HOSVD) [1][2][3], Tensor Trains Decomposition (TTD) [4], and Hierarchical Tucker decomposition (H Tucker) [5] -and their modifications [6] are based on the calculation of the eigen values and eigen vectors of the decomposed tensor. Their basic advantage is that they are optimum with respect to the Mean Square Error (MSE) of the approximation in the case of the truncation of the low-energy decomposition components. For the calculation of the retained components are used various iterative methods (the power method [7], the Jacoby method [8], etc.), which require relatively high numbers of calculations to achieve the needed approximation accuracy.
As an alternative, in this work are presented new hierarchical approaches for 3D tensor decomposition, based on well-known deterministic orthogonal transforms: the Walsh-Hadamard Transform (WHT) [9][10][11][12] and the Complex Hadamard Transform (CHT) [13]. The offered decompositions are not optimum with respect to MSE minimization, but due to the lack of iterations, they have low Computational Complexity (CC), and as a result, they are suitable for the fast processing of the 3D tensors obtained, for example, from single or sequences of 2D correlated images. In the first case, each single 2D image is divided into N 2 blocks of size N × N for N = 2 n , from which is obtained a sequence of N cubical tensors, each of size N × N × N. In the second case, from each sequence of N 2D images of size N × N is obtained a single cubical tensor of size N × N × N. Examples for the transformation of a single 2D image and of a sequence of 2D video frames into a cubical (3D) tensor are shown on Figure 1.
Symmetry 2020, 12, x FOR PEER REVIEW 2 of 17 images. In the first case, each single 2D image is divided into N 2 blocks of size N×N for N=2 n , from which is obtained a sequence of N cubical tensors, each of size N×N×N. In the second case, from each sequence of N 2D images of size N×N is obtained a single cubical tensor of size N×N×N. Examples for the transformation of a single 2D image and of a sequence of 2D video frames into a cubical (3D) tensor are shown on Figure 1. The presented approach for fast hierarchical cubical tensor decomposition is applicable not only for images but also for various kinds of multidimensional signals (medical, seismic, spectrometric, etc.) and big data analysis. In the next sections, Sections 2-5, are presented the proposed algorithms for 1D Hierarchical Fast Walsh-Hadamard Transform (1D-FWHT), 1D Hierarchical Fast Complex Hadamard Transform (1D-FCHT), and the corresponding 3D Hierarchical Fast real and complex transforms (3D-FWHT and 3D-FCHT). In Section 6 is shown a comparative analysis of the CC of the presented algorithms for hierarchical orthogonal tensor transforms, with respect to the famous deterministic and statistical transforms: the 3D Fast Fourier Transform (3D-FFT), 3D Discrete Wavelet Transform (3D-DWT) and Hierarchical Tucker Decomposition (H-Tucker). The last section, Section 7, contains the conclusions.
In Appendices A.1 and A.2 is given the factorization of the matrices for the n-level one-dimensional fast hierarchical transforms 1D-FWHT and 1D-FCHT.

One-Dimensional Hierarchical Fast Walsh-Hadamard Transform
The following symbols are introduced for tensors, matrices and vectors, respectively: X for a tensor, X for a matrix and x  for a vector.
The forward and inverse 1D Walsh-Hadamard Transform (1D-WHT) (frequency-ordered) is represented in scalar form by the equations below [9][10][11] where N=2 n and x(i) and s(u) are, respectively, the N-dimensional discrete signal and its spectrum. The discrete Wash functions are defined by the following relations [10,12]: The presented approach for fast hierarchical cubical tensor decomposition is applicable not only for images but also for various kinds of multidimensional signals (medical, seismic, spectrometric, etc.) and big data analysis. In the next sections, Sections 2-5, are presented the proposed algorithms for 1D Hierarchical Fast Walsh-Hadamard Transform (1D-FWHT), 1D Hierarchical Fast Complex Hadamard Transform (1D-FCHT), and the corresponding 3D Hierarchical Fast real and complex transforms (3D-FWHT and 3D-FCHT). In Section 6 is shown a comparative analysis of the CC of the presented algorithms for hierarchical orthogonal tensor transforms, with respect to the famous deterministic and statistical transforms: the 3D Fast Fourier Transform (3D-FFT), 3D Discrete Wavelet Transform (3D-DWT) and Hierarchical Tucker Decomposition (H-Tucker). The last section, Section 7, contains the conclusions.
In Appendices A.1 and A.2 is given the factorization of the matrices for the n-level one-dimensional fast hierarchical transforms 1D-FWHT and 1D-FCHT.

One-Dimensional Hierarchical Fast Walsh-Hadamard Transform
The following symbols are introduced for tensors, matrices and vectors, respectively: X for a tensor, X for a matrix and → x for a vector. The forward and inverse 1D Walsh-Hadamard Transform (1D-WHT) (frequency-ordered) is represented in scalar form by the equations below [9][10][11]: where N = 2 n and x(i) and s(u) are, respectively, the N-dimensional discrete signal and its spectrum. The discrete Wash functions are defined by the following relations [10,12]: ; q r (i) = i n−r ⊕ i n−r−1 for r = 1, 2, . . . , n − 1; q 0 (i) = i n−1 ; i = n−1 r=0 i r 2 r , u = n−1 r=0 u r 2 r . (2) Here, the operation "exclusive OR" is represented by the symbol ⊕.
The algorithm for the one-dimensional hierarchical fast Walsh-Hadamard transform (1D-FWHT) for the cubical tensor X of size N × N × N when N = 2 n (a sequence of matrices X k , each of size N × N for k = 1, 2, . . . ,N) is presented for N = 8 (i.e., for n = 3 levels). The general case of the transform when N = 2 n and n >3, is given in Appendix A.1.
The execution of the direct and inverse 1D-WHT is grounded on the basic operation "butterfly" for N = 2, in accordance with the following relations [9,11]: The application of the operation "butterfly" to the elements of the couples of matrices in each hierarchical level is represented by the equations below.
In Level 1 of the 1D-FWHT: The matrix G 1 (8) of size 8 × 8 is used to execute the direct 1D-WHT when N = 2 for each couple of neighbor matrices in the sequence X k for k = 1,2, . . . ,8. These matrices are the components of the matrix-column X = [X 1 , X 2 , X 3 , X 4 , X 5 , X 6 , X 7 , X 8 ] T . From Equation (3), it follows that the matrix G 1 (8) could be represented in the following way: Here, I(4) is the identity matrix of size 4 × 4, the symbol ⊗ stands for the Kroneker product of the matrices [3], and Y = [Y 1 , Y 2 , Y 3 , Y 4 , Y 5 , Y 6 , Y 7 , Y 8 ] T corresponds to the transformed matrix-column whose components are the matrices Y k , for k = 1,2, . . . ,8.
After the rearrangement 1 is obtained: where In the Level 2 of the 1D-FWHT is obtained: After the rearrangement 2 is obtained: where In the Level 3 of the 1D-FWHT is obtained: where After the rearrangement 3 is obtained: where P 3 (8) = P 2 (8); E = [E 1 , E 2 , E 3 , E 4 , E 5 , E 6 , E 7 , E 8 ] T -output matrix-column with components The relation between the input and output matrix-column X and E, correspondingly, defined by the relations in Equations (4)-(13), is: where H w (8) is a frequency-ordered Hadamard matrix, defined as follows: Symmetry 2020, 12, 864

of 17
The rearrangement of the intermediate matrix components in each consecutive transform level must satisfy the requirement for the frequency-ordering of the output matrices [9,14]. The chosen method for matrix H w (8) factorization permits the operation "butterfly" in the first two levels to be applied to the elements of the neighbor couples of matrices X k and X k+1 (for k = 1,2, . . . ,8) from the sequence, which builds the cubical tensor X of size 8 × 8 × 8. In this way, the decorrelation efficiency for the elements of each couple of transformed matrices is enhanced, which results in higher energy concentration in the low-frequency components of the output matrix-column, E. In the result of the 1D-FWHT is obtained the cubical tensor E of size 8 × 8 × 8, represented by the sequence of matrices E l , each of size 8 × 8, when l = 1,2, . . . ,8.
The graph of the 3-level 1D-FWHT algorithm for the tensor X of size 8 × 8 × 8 is shown in Figure 2. Here, the basic operation "butterfly" in the first two levels is framed in red. An example 1D-FWHT is given below for the processing of a sequence of matrices Xk for k=1,2,..,8, in a particular case: the transform of one element ) k , m , i ( x of the same spatial position (i,m) in each matrix, k: Let: In Level 2 of the 1D-HFWHT: After the rearrangement 2: An example 1D-FWHT is given below for the processing of a sequence of matrices X k for k = 1, 2, . . . , 8, in a particular case: the transform of one element x(i, m, k) of the same spatial position (i,m) in each matrix, k: Let: Then, in Level 1 of the 1D-HFWHT: In Level 2 of the 1D-HFWHT: After the rearrangement 2: In Level 3 of the 1D-HFWHT: After the rearrangement 3: Here, y(i, m, k), z(i, m, k), d(i, m, k) and e(i, m, k) are, correspondingly, the elements of the sequences of matrices Y k , Z k , D k and E k for k = 1,2, . . . ,8, calculated in levels 1, 2 and 3, and for the graph output.
In accordance with the generalized equation of Parseval [15] is obtained: Here, x(i.m, k) and e(u, v, l) are, respectively, the elements of the input tensor X and of the spectral tensor, E. Then, for the example above, for N = 8, from Equation (16) In Appendix A.1 is given, in detail, the frequency-ordered hierarchical n-level 1D-FWHT.

One-Dimensional Hierarchical Fast Complex Hadamard Transform
The direct and inverse 1D Complex Hadamard Transform (1D-CHT) (frequency-ordered) is represented in a scalar form by the equations below [14]: where N = 2 n ≥4 (the minimum possible value is n = 2); x(i) and s(u) are, correspondingly, the N-dimensional discrete signal and its complex spectrum; j = √ −1; and s Re (u), and s Im (u) are the real and the imaginary parts of the spectrum s(u). The coefficients of the 1D-CHT matrix C(N) of size N × N are represented by the equation: where i, u = 0,1, . . . ,2 n −1. The sign function h(i,u) is: Symmetry 2020, 12, 864 Here, • is an operator, which represents the integer part of the result, obtained after the division. The even coefficients of the 1D-CHT are real, and the odd ones are complex-conjugated.
In Figure 3 is shown the calculation graph of the direct 3-level fast hierarchical 1D-CHT (n = 3) for N = 8, built in accordance with Figure 2. The graph for the 1D-FCHT when N=8 is created on the basis of Equation (20), to which are added four "butterfly" operations with weight coefficients (+1,−1) in the last, third level. In the result of the use of the frequency-ordered 1D-FCHT [14] is obtained the following system of equations, which represent the relations between the sequences of the input and output matrices, respectively Xk and Ek, for k=1,2,..,8 and N=8: Unlike Equation (14), part of the equations in the system above contain the complex variable The matrix Φk permits the phase modification of its elements, which could be used (for The graph for the 1D-FCHT when N = 8 is created on the basis of Equation (20), to which are added four "butterfly" operations with weight coefficients (+1,−1) in the last, third level. In the result of the use of the frequency-ordered 1D-FCHT [14] is obtained the following system of equations, which represent the relations between the sequences of the input and output matrices, respectively X k and E k , for k = 1, 2, . . . , 8 and N = 8: Unlike Equation (14), part of the equations in the system above contain the complex variable j = √ −1, and they represent complex matrices (i.e., these are relations E 3 , E 4 , E 5 and E 6 ). The remaining four equations (E 1 , E 2 , E 7 and E 8 ) represent real matrices. The components of each complex matrix E k = A k + jB k (k = 3,4,5,6) are A k = p=1,3,5,7 α p X p and B k = p=2,4,6,8 β p X p ; α p and β p are sign functions with values of +1 or −1, defined by Equation (21). The elements e(i,m.k) of the complex matrix E k are defined as follows: where M k (i, m) = e(i,m,k) 2 Re + e(i,m,k) 2 Im is a module; and Φ k (i, m) = arctg[e(i,m,k) Im /e(i,m,k) Re ], a phase of the element e(i,m,k).
The matrix Φ k permits the phase modification of its elements, which could be used (for example) for the resistant digital watermarking of multidimensional signals and images [14].
In Appendix A.2 is presented in detail the factorization of the matrix for an n-level hierarchical frequency-ordered 1D-FCHT, based on the "butterfly" operation (Equation (20)) used for the 1D-CHT execution, if N = 4.

Hierarchical Cubical Tensor Decomposition through the 3D-FWHT
At first, here is described the decomposition of the cubical tensor X of size 8 × 8 × 8. In this case, the decomposition could be represented in the simplified form shown in Figure 4.  The decomposition based on the 3D-FWHT comprises three stages. In the first stage, the tensor is divided into horizontal slices (mode-1), after which, on each couple of matrices, is applied the 1D-FWHT, and from the so obtained eight slices is restored the tensor E. In a similar way are executed the operations in the second and third stages. Unlike in the first, in the second stage, the input tensor E is divided into eight lateral slices (mode-2). After their transform through the 1D-FWHT, from the obtained eight slices is restored the tensor F. In the third stage, the tensor F is divided into eight frontal slices (mode-3), from which, after the 1D-FWHT, is obtained the spectrum tensor S with elements ) u,v,l ( s . From this tensor, after the reverse execution of the three stages of the inverse 1D-FWHT, is restored the initial tensor X with elements x(i,j,k). In the result of the decomposition, it is represented as a sum of weighted "basic" tensors u,v,l W of size 8×8×8 with elements wal(u,v,l), which are 3D Walsh functions: The decomposition based on the 3D-FWHT comprises three stages. In the first stage, the tensor is divided into horizontal slices (mode-1), after which, on each couple of matrices, is applied the 1D-FWHT, and from the so obtained eight slices is restored the tensor E. In a similar way are executed the operations in the second and third stages. Unlike in the first, in the second stage, the input tensor E is divided into eight lateral slices (mode-2). After their transform through the 1D-FWHT, from the obtained eight slices is restored the tensor F. In the third stage, the tensor F is divided into eight frontal slices (mode-3), from which, after the 1D-FWHT, is obtained the spectrum tensor S with elements s(u, v, l). From this tensor, after the reverse execution of the three stages of the inverse 1D-FWHT, is restored the initial tensor X with elements x(i,j,k). In the result of the decomposition, it is represented as a sum of weighted "basic" tensors W u,v,l of size 8 × 8 × 8 with elements wal(u,v,l), which are 3D Walsh functions: The 3D functions wal(u,v,l) are divisible [16,17], and they could be represented as the product of three 1D Walsh functions: wal(u).wal(v).wal(l). The coefficients s(u,v,l) for u = 0,1, v = 0,1 and l = 0,1 in the first layer correspond to the lowest spatial frequencies of the spectrum tensor, S. In Figure 5a is shown an example for a spectrum tensor S of size 4 × 4 × 4, and in Figure 5b, the 8 "basic" tensors, which correspond to coefficients s(u,v,l) from the initial layer of the tensor S [18]. The 3D-FWHT algorithm, shown in Figure 4, could be generalized for the case when N=2 n by creating an extended n-level computational graph in correspondence with the transform matrix factorization explained in Appendix A.1. The number of stages needed for the execution of the 1D-FWHT (in accordance with Figure 4) is three, not only for N=8 but also for the higher values of N=16, 32,… In the general case, the spectrum coefficients s(u,v,l) are calculated through the direct 3D-WHT, in accordance with the relation [12,18]: The 3D inverse WHT (3D-IWHT) is defined by the relation: Here, the discrete one-dimensional Walsh-Hadamard (WH) functions of N th order wal(i,u), wal(m,v), wal(k,l) for i,m,k = 0,1,..,N−1, are defined by the following relations: Each "basic" tensor v,l u, W with frequency (u,v,l) could be represented through the outer x(i, j, k )wal(i, u)wal(j, v)wal(k, l) for u, v, l = 0, 1, . . . , N − 1.
The 3D inverse WHT (3D-IWHT) is defined by the relation: x(i, m, k) = where q r (i), q r (m), q r (k) and u r , v r , l r are defined by analogy with Equation (2)  Each "basic" tensor W u,v,l with frequency (u,v,l) could be represented through the outer product of the vectors where the vectors, which represent the tensor W u,v,l , are defined by the relations below: Then, the tensor decomposition based on the 3D-WHT, is defined by the relation: As follows from the relations in Equations (25) and (32), the 3D-WHT is reversible, i.e., X could be restored from the tensor S through the 3D-IWTH. The decomposition in Equation (32) of the tensor X of size N × N × N corresponds to the Tucker decomposition [1,3,5]: Here, g r 1 ,r 2 ,r 3 are the entries of the core tensor G of size R 1 × R 2 × R 3 , and R 1 , R 2 , R 3 , the multilinear rank of the tensor X. For the case, when it is a cube of size N × N × N, the tensor G is diagonal and of size R×R×R; g r (for r = 1,2, . . . ,R) are its eigen values, and → a r 1 , → b r 2 , → c r 3 , the eigen vectors. In this case, the Tucker decomposition is transformed into CPD [1], and is represented by the relation below: Here, R is the rank of the tensor X, whose value is limited in the range N ≤ R ≤ N 2 [19]. The difference between the decompositions, represented by Equations (33) and (32), is in the numbers of their components; in the first case, it is ≤N 2 , and in the second, it is ≤N 3 . This shows the higher energy concentration in the first components of the H-Tucker decomposition (respectively, CPD) compared to that of the 3D-WHT.

Hierarchical Cubical Tensor Decomposition through the 3D-FCHT
The comparison of the relations in Equations (14), (15) and (21) shows that the 1D-FWHT and 1D-FCHT are executed in accordance with algorithms with similar structure (Figures 2 and 3). From this, it follows that for the decomposition of the tensor X of size 8 × 8 × 8 through 3D-FCHT could be used the three-stage scheme shown in Figure 4. The difference is that in the three decomposition stages, the 1D-FWHT algorithm from Figure 2 must be replaced by this for the 1D-FCHT, shown in Figure 3. After the triple execution of the algorithm is obtained the spectrum tensor S (shown on Figure 4) with elements s(u, v, l), half of which are complex. The 3D-FCHT algorithm, shown in Figure 3, could be generalized for the case when N = 2 n by creating an extended computational graph of n levels. The number of stages needed for the execution of the 1D-FCHT by analogy with the 3D-FWHT is also three (not for N = 8 only but also for the next values of N = 16, 32, . . . ).
In the general case, the decomposition of the tensor X of size N × N × N in accordance with Equation (27) is represented as the sum of the weighted "basic" tensors: where The vectors which define the tensor C u,v,l are defined by the relations below: In the above relations, the sign functions h(i,u), h(m,v) and h(k,l) for i,m,k = 0, 1, . . . , N−1 are defined by Equation (19). The spectrum coefficients s(u,v,l) in Equation (35)  (40)

Comparative CC Evaluation for the 3D-FWHT and 3D-FCWHT Algorithm
In this section is evaluated the CC of the decomposition algorithms for a tensor of size N × N × N, when N = 2 n .

CC of the Algorithm Based on the Real Orthogonal Transform, 3D-FWHT
The number of additions, A W (1), needed for the execution of the 1D-WHT based on the "butterfly" operation for all N 2 elements of a couple of matrices, both of size N × N, is: For all N/2 couples of matrices in one stage of the hierarchical 3D-FWHT of n levels, this number is: A 1H (n) = (N/2)2N 2 lg 2 N = 2 3n n.
Then, for the three stages of the hierarchical 3D-FWHT is obtained: The needed number of multiplications for the execution of the direct 3D-FWHT is M 3H (n) = 0. Then, the CC of the 3D-FWHT evaluated through the total number of operations O 3H (n) is defined by the relation: 6.2. CC of the Algorithm Based on the Complex Orthogonal Transform, 3D-FCHT The number of additions for all N 2 elements of N/2 couples of matrices of size N × N in one stage of the hierarchical n-level 1D-FCHT is: In the relation above is taken into account that in the first level of the 1D-FCHT are added N real numbers; in the second level, there are N/2 real numbers, and in the remaining (n-2) levels, N/4 real and 3N/4 complex numbers. Additionally, the summarizing of two complex numbers is executed as a sum of two couples of real numbers.
Then, for the three stages of the 3D Fast CHT (3D-FCHT), the number of needed additions is: The number of complex multiplications needed for the direct 3D-FCHT is M 3C (n) = 0 (the multiplication of j by a real number is not an arithmetical operation but only needs larger memory in which to save the so obtained complex number). Then, the CC of the 3D-FCHT, evaluated on the basis of the global number of needed operations O 3C (n), is defined by the relation: (47)

CC of the Algorithm Based on the Complex Orthogonal Transform, 3D-DFT
The divisibility of the 3D fast Fourier transform (3D-DFT) permits the tensor decomposition to be executed through applying the 1D-DFT on all fibers (vectors) obtained in the result of the tensor transform, called vectorization [1,3]. The vectors are oriented in three directions: x (mode-1), y (mode-2) and z (mode-3). The calculation of the 1D-DFT is based on the Cooley-Tukey algorithm for Fast Fourier Transform (FFT) [15]. In the result of the calculation of the radix-2 FFT for an N-dimensional vector, the number of operations (respectively, CC) is reduced from N 2 to Nlog 2 N = Nn. Then, the CC of the 3D-FFT is O(N c log 2 N c ), where N c = N 1 × N 2 × N 3 . In the case that N 1 = N 2 = N 3 = N = 2 n , for a tensor of N 3 voxels, the number of the complex additions is A 3F (n) = N 3 log 2 N 3 = 3N 3 n, and that of the complex multiplications, M 3F (n) = (N/2) 3 log 2 N 3 = (3/8)N 3 n [20]. The multiplication of two complex numbers is equivalent to four multiplications of real numbers, and one complex addition, to two additions of real numbers. Then, for the CC of the 3D-FFT is obtained: O 3F (n) = A 3F (n) + M 3F (n) = 6.75N 3 n + 0.75N 3 n = 7.5 × 2 3n n. (48)

CC of the Algorithm Based on the 3D Discrete Wavelet Transform, 3D-DWT and H-Tucker
The global number of additions and multiplications needed to calculate the cubical tensor of size N × N × N through the 3D Discrete Wavelet Transform (3D-DWT) of n levels with orthogonal filters (3,5), are correspondingly [20]: Then, the CC of the 3D-DWT is: The H-Tucker decomposition [5] of the cubical tensor of minimum rank R = 2 n (as defined in [19]), size N = 2 n and order d = 3 requires O HT (3 × 2 3n + 2 × 2 4n ) operations. For the TT-decomposition [4], the CC for the same tensor is O TT (3 × 2 4n ), i.e., it is approximately 1.5 times higher than that for the H-Tucker. This is why the H-Tucker transform was selected for the CC comparison with the analyzed 3D deterministic orthogonal transforms.

CC of 3D-FWHT and 3D-FCHT, Compared to 3D-FFT, 3D-DWT and H-Tucker
The relative acceleration ψ(n) of the calculations needed for the execution of the 3D-FWHT, 3D-FCHT and 3D-FFT compared to that of the 3D-DWT and H-Tucker is given in detail below: 3D-FWHT, 3D-FCHT and 3D-FFT compared to 3D-DWT: In Figure 6a are shown the graphics of functions ψ i (n) for i = 1,3,5, which define the relative CC of the 3D-FWHT, 3D-FCHT and 3D-FFT towards the 3D-DWT, and in Figure 6b, the functions ψ i (n) for i = 2,4,6, which define the relative CC for same transforms towards the H-Tucker. The graphics in Figure 6a show that, together with the growth of the hierarchical levels n = log 2 N, the relative CC of the 3D-FWHT, 3D-FCHT and 3D-FFT decreases inversely proportionally to n towards the 3D-DWT, while towards H-Tucker, it grows proportionally to 2 n /n. As follows from the graphics on Figure 6b, the CC of the 3D-FCHT is the lowest compared to that of the 3D-FWHT, H-Tucker and 3D-FFT only for levels n = 2,3, while for the next levels (n = 4, 5, . . . , 10), the CC of the 3D-FWHT is the lowest. For the same range of n (from 4 to 10), the value of the function ψ 2 (n) changes form 2.6 to 67.5, and that of ψ 4 (n), from 2.5 to 44. The comparison results permit the choosing of the number of hierarchical decomposition levels n for which the developed new algorithms 3D-FWHT and 3D-FCHT are more efficient than 3D-FFT, 3D-DWT and H-Tucker. graphics on Figure 6 b, the CC of the 3D-FCHT is the lowest compared to that of the 3D-FWHT, H-Tucker and 3D-FFT only for levels n=2,3, while for the next levels (n=4,5,..,10), the CC of the 3D-FWHT is the lowest. For the same range of n (from 4 to 10), the value of the function ) n ( ψ 2 changes form 2.6 to 67.5, and that of

Conclusions
In this work is presented one new approach for the decomposition of the 3D tensors of size N × N × N, based on the deterministic orthogonal transforms 3D-FWHT and 3D-FCHT. The results of the comparative analysis for the CC of these decompositions outline their advantages over 3D-FFT, 3D-DWT and H-Tucker. Another advantage is that no iterative calculations are needed, and instead are used hierarchical calculations only, of low CC. The comparison results show that after the 1D-transform of the couples of tensor slices (matrices) in three mutually orthogonal directions, the energy is concentrated into a small number of low-frequency spectrum coefficients. The disadvantage of the method is that unlike the famous 3D decompositions, which are based on the tensor eigen vectors, the offered deterministic orthogonal transforms do not ensure minimum MSE, due to low-energy decomposition component truncation. The choice of the most suitable 3D deterministic hierarchical tensor decomposition depends on the corresponding application area.
The method could be also used for parallelepiped tensors. In such a case, the original tensor is divided into cubical sub-tensors, and on each should be applied the presented decomposition. The incomplete sub-tensors at the main tensor boundaries should be supplemented through 3D interpolation based on the closest boundary data. Besides, instead of WHT or CHT could be used orthogonal transforms whose matrices are of size different from 2 n . Future investigations will be aimed at the modeling of the offered algorithms, to permit the full evaluation of their characteristics. Additionally, the algorithm for tensor decomposition will be further developed to process tensors of any size (not cubical only). As a result, the most efficient applications will be detected for various areas such as data compression and filtration; the analysis, recognition and searching of visual information; deep learning, digital image watermarking, etc.

Conclusions
In this work is presented one new approach for the decomposition of the 3D tensors of size N × N × N, based on the deterministic orthogonal transforms 3D-FWHT and 3D-FCHT. The results of the comparative analysis for the CC of these decompositions outline their advantages over 3D-FFT, 3D-DWT and H-Tucker. Another advantage is that no iterative calculations are needed, and instead are used hierarchical calculations only, of low CC. The comparison results show that after the 1D-transform of the couples of tensor slices (matrices) in three mutually orthogonal directions, the energy is concentrated into a small number of low-frequency spectrum coefficients. The disadvantage of the method is that unlike the famous 3D decompositions, which are based on the tensor eigen vectors, the offered deterministic orthogonal transforms do not ensure minimum MSE, due to low-energy decomposition component truncation. The choice of the most suitable 3D deterministic hierarchical tensor decomposition depends on the corresponding application area.
The method could be also used for parallelepiped tensors. In such a case, the original tensor is divided into cubical sub-tensors, and on each should be applied the presented decomposition. The incomplete sub-tensors at the main tensor boundaries should be supplemented through 3D interpolation based on the closest boundary data. Besides, instead of WHT or CHT could be used orthogonal transforms whose matrices are of size different from 2 n . Future investigations will be aimed at the modeling of the offered algorithms, to permit the full evaluation of their characteristics. Additionally, the algorithm for tensor decomposition will be further developed to process tensors of any size (not cubical only). As a result, the most efficient applications will be detected for various areas such as data compression and filtration; the analysis, recognition and searching of visual information; deep learning, digital image watermarking, etc.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Multi-Level Representation of One-Dimensional Hierarchical Hadamard Transforms
The mathematical descriptions of the n-level 1D-FWHT and 1D-FCHT are given below in correspondence with the transformation graphs from Figures 2 and 3, for N = 8. The factorization of the matrices used for the n-level transforms permits their execution based on matrix processors, which will additionally accelerate all calculations.