A New Two-Sided Sketching Algorithm for Large-Scale Tensor Decomposition Based on Discrete Cosine Transformation
Zhiguang Cheng Gaohang Yu Xiaohao Cai Liqun Qi


Understanding Privacy Risks of Embeddings Induced by Large Language Models
    Zhihao Zhu, Ninglu Shao, Defu Lian, Chenwang Wu, Zheng Liu, Yi Yang, Enhong Chen
    April 29, 2024
====================================================================================




Large tensors are frequently encountered in various fields such as computer vision, scientific simulations, sensor networks, and data mining. However, these tensors are often too large for convenient processing, transfer, or storage. Fortunately, they typically exhibit a low-rank structure that can be leveraged through tensor decomposition. Despite this, performing large-scale tensor decomposition can be time-consuming. Sketching is a useful technique to reduce the dimensionality of the data. In this study, we introduce a novel two-sided sketching method based on the t-product decomposition and the discrete cosine transformation. We conduct a thorough theoretical analysis to assess the approximation error of the proposed method. Specifically, we enhance the algorithm with power iteration to achieve more precise approximate solutions. Extensive numerical experiments and comparisons on low-rank approximation of color images and grayscale videos illustrate the efficiency and effectiveness of the proposed approach in terms of both CPU time and approximation accuracy.











12pt Key words. Large-scale tensor, tensor decomposition, sketching, t-product, discrete cosine transformation, power iteration.















§ INTRODUCTION



Multidimensional arrays, known as tensors, are often used to represent real-world high-dimensional data, such as videos <cit.>, hyperspectral images <cit.>, multilinear signals <cit.>, and communication networks <cit.>. In most cases, these tensor data usually have a low-rank structure and can be approximated by tensor decomposition. Nevertheless, computing the tensor decomposition of these large-scale data is usually computationally demanding, and thus finding an accurate approximation of large-scale data with great efficiency plays a key role in tensor data analysis. Sketching is a useful technique for data compression, utilizing random projections or sampling to imitate the original data. Although the sketching technique may slightly reduce the accuracy of the approximation, it can significantly reduce the computational and storage complexity <cit.>. As a result, sketching is commonly used in low-rank tensor approximation, and many researchers have proposed various tensor sketching algorithms <cit.>. They have also been successfully applied to a variety of tasks, such as Kronecker product regression, polynomial approximation, and the construction of deep convolutional neural networks <cit.>, etc.

Recently, extensive research has been carried out on the application of sketching algorithms for low-rank matrix approximation <cit.>. Woodruff et al. <cit.> examined the numerical linear algebra algorithms of linear sketching techniques and identified their limitations. Tropp et al. <cit.> developed the two-sided matrix sketching algorithm, which can maintain the structural properties of the input matrix and generate a low-rank matrix approximation with a given rank. Tropp et al. <cit.> also proposed a new two-sided matrix sketching algorithm to construct a low-rank approximation matrix from streaming data. These matrix sketching algorithms are very effective in reducing storage and computational costs when computing low-rank approximations of large-scale matrices. Subsequently, many researchers have applied matrix sketching algorithms to tensor decomposition and developed various low-rank tensor approximation algorithms. The following are some low-rank tensor approximation algorithms based on different decompositions <cit.>.

Based on the CANDECOMP/PARAFAC (CP) decomposition, the robust tensor power method based on tensor sketch (TS-RTPM) can quickly mine the potential features of tensor, but in some cases, its approximation performance is limited, Cao et al. <cit.> propose a data-driven framework called TS-RTPM-Net, which improves the estimation accuracy of TS-RTPM by jointly training the TS value matrices with the RTPM initial matrices. Li et al. <cit.> introduced a randomized algorithm for CP tensor decomposition in least squares regression, aiming to achieve reduced dimensionality and sparsity in the randomized linear mapping. Wang et al. <cit.> developed innovative techniques for performing randomized tensor contractions using FFTs, avoiding the explicit formation of tensors. Numerous studies have also been conducted based on Tucker decomposition <cit.>. For instance, Che et al. <cit.> devised an adaptive randomized approach for approximating Tucker decomposition. Malik et al. <cit.> proposed two randomized algorithms for low-rank Tucker decomposition, which entail a single pass of the input tensor by integrating sketching. Ravishankar et al. <cit.> introduced the hybrid Tucker TensorSketch vector quantization (HTTSVQ) algorithm for dynamic light fields. Sun et al. <cit.> developed a randomized method for Tucker decomposition, which can provide a satisfactory approximation of the original tensor data through the incorporation of two-sided sketching. Minster et al. <cit.> devised randomized adaptations of the THOSVD and STHOSVD algorithms. Dong et al. <cit.> present
two practical randomized algorithms for low-rank Tucker approximation of
large tensors based on sketching and power scheme, with a rigorous errorbound
analysis.

The work based on tensor-train (TT) decomposition can be found e.g. in <cit.>. In particular, Che et al. <cit.> designed an adaptive random algorithm to calculate the tensor column approximation. Hur et al. <cit.> introduced a sketching algorithm to construct a TT representation of a probability density from its samples, which can avoid the curse of dimensionality and sample complexities of the recovery problem. 

Qi and Yu <cit.> proposed a tensor sketching method based on the t-product, which can quickly obtain a low tubal rank tensor approximation. As pointed out by Kernfeld et al. <cit.>, the t-product has a disadvantage in that, for real tensors, implementation of the t-product and factorizations using the t-product require intermediate complex arithmetic, which, even taking advantage of complex symmetry in the Fourier domain, is more expensive than real arithmetic. 



In this paper, based on the ⋆_L-product decomposition <cit.>,

we investigate some new two-sided sketching algorithms for low tubal rank tensor approximation. The main contributions of this paper are threefold. 
Firstly, we propose a new one-pass sketching algorithm (i.e., Algorithm 4), which can significantly improve the computational efficiency of T-sketch and rt_SVD, for low tubal rank approximation. Secondly, we first establish a low tubal rank tensor approximation model based on ⋆_L-product factorization, extending the two-sided matrix sketching algorithm proposed by Tropp et al. <cit.>. A rigorous theoretical analysis is conducted to evaluate the approximation error of the proposed Algorithm 5.  Thirdly, the proposed Algorithm 6 can deliver a more accurate tensor approximation tensor than simpler randomized algorithms by combining the subspace power iteration. 



The following paper is organized as follows. Section 2 introduces some common notations and preliminaries. Section 3 presents two new two-sided sketching algorithms based on the transformed domain (i.e., Algorithm 4  and Algorithm 6). Section 4 provides strict theoretical guarantees for the approximation error of the proposed algorithms. Section 5 presents numerical experiments and comparisons demonstrating the efficiency and effectiveness of the proposed algorithms. We conclude in Section 6.




§ NOTATIONS AND PRELIMINARIES

In this paper, matrices and tensors are represented by capital letters (e.g. A, B,…) and curly letters (e.g. 𝒜,ℬ,…), respectively. We can use the Matlab command A' represent conjugate transpose of the matrix A. ℝ and ℂ represent the real number space and the complex number space, respectively.

For matrix A∈ℂ^n_1× n_2, its (i,j)-th element is represented by a_i,j. For  third-order tensor 𝒜∈ℂ^n_1× n_2 × n_3, its (i,j,k)-th element is represented by a_i,j,k. 
The Matlab notations 𝒜(i,:,:), 𝒜(:,i,:) and 𝒜(:,:, i) are used to represent the i-th horizontal, lateral and facial slices of 𝒜, respectively. 
The facial slice 𝒜(:,:,i) is also represented by 𝒜^(i). The Frobenius norm of a tensor is defined as the square root of the sum of the squares of its elements:

    𝒜 _F:=𝒜(:) _2=√(⟨𝒜,𝒜⟩)=√(∑_ijk|a_ijk|^2) .

𝒜^H and 𝒜^† represent the conjugate transpose and pseudo-inverse of 𝒜, respectively.

This paper focuses on the low tubal rank tensor approximation meeting desired accuracy in an efficient manner.
For tensor 𝒜∈ℝ^n_1× n_2× n_3, the mathematical model for finding the low-rank approximation 𝒜̂ of 𝒜 can be expressed as

    𝒜-𝒜̂ _F = min_ rank_L(ℬ)≤ l𝒜-ℬ _F,

where l, l≪min{n_1,n_2}, is the target rank, L (by abuse of notation) represents an arbitrary invertible linear transform, and rank_L(ℬ) denotes the transformed tubal rank of tensor ℬ.




 §.§ Transformed Tensor SVD



We below briefly recall the transformed tensor SVD of third-order tensors; more details can be found in <cit.>.

Let F_n_3 and C_n_3 represent the discrete Fourier transform matrix and the DCT matrix, respectively. For any third-order tensor 𝒜∈ℝ^n_1× n_2× n_3, let 𝒜̅_L represent a third-order tensor obtained via multiplied by L (an arbitrary invertible linear transform) on all tubes along the third-dimension of 𝒜, i.e., 𝒜̅_L(i,j,:)= L(𝒜(i,j,:)), i=1,…,n_1,j=1,…,n_2. Here we write 𝒜̅_L=L[𝒜]. Moreover, one can get 𝒜 from 𝒜̅_L by using L^H along the third-dimension of 𝒜̅_L, i.e., 𝒜=L^H[𝒜̅_L].  The ⋆_L-product is defined in Definition <ref> below.


	<cit.> For any two tensors 𝒳∈ℂ^n_1× n_2× n_3 and 𝒴∈ℂ^n_2× n_4× n_3, and an arbitrary invertible linear transform L, the ⋆_L-product of 𝒳 and 𝒴 is a tensor 𝒵∈ℂ^n_1× n_4× n_3 given by
	
    𝒵=𝒳⋆_L𝒴=L^H[fold(Block(𝒳̅_L)Block(𝒴̅_L))],

	where fold(Block(𝒳̅_L))=𝒳̅_L and X̅_L=Block(𝒳̅_L)=[   𝒳̅_L^(1)                                 ;              𝒳̅_L^(2)                      ;                                ⋱;                                  𝒳̅_L^(n_3) ]. 



	The Kronecker product of matrices A∈ℂ^m× l and B∈ℂ^p× r is
	
    A ⊗ B=[ a_11B a_12B     … a_1lB; a_21B a_22B     ⋱ a_2lB;     ⋮     ⋮     ⋱     ⋮; a_m1B a_m2B     … a_mlB ].



The t-product <cit.> of 𝒜∈ℂ^n_1× n_2× n_3 and ℬ∈ℂ^n_2× n_4× n_3 is a tensor 𝒞∈ℂ^n_1× n_4× n_3 given by

    𝒞=𝒜*ℬ=Fold_vec(Circ(𝒜)×Vec(ℬ)),


where

    Vec(ℬ)=[   ℬ^(1);   ℬ^(2);       ⋮; ℬ^(n_3) ]∈ℝ^n_2n_3× n_4,   Fold_vec(Vec(ℬ))=ℬ,



    Circ(𝒜)=[     𝒜^(1)   𝒜^(n_3) 𝒜^(n_3-1)         …     𝒜^(2);     𝒜^(2)     𝒜^(1)   𝒜^(n_3)         …     𝒜^(3);         ⋮         ⋮         ⋮         ⋱         ⋮;   𝒜^(n_3) 𝒜^(n_3-1) 𝒜^(n_3-2)         …     𝒜^(1) ].



The t-product in Eq. (<ref>) can be seen as a special case of Definition <ref>. Recall that the block circulant matrix Circ(𝒜) can be
diagonalized by the fast Fourier transform matrix F_n_3 and the block diagonal matrices are the frontal slices of 𝒜̅_F_n_3 ,
i.e.,

    (F_n_3⊗ I_n_1)×Circ(𝒜)× (F^H_n_3⊗ I_n_2)=Block(𝒜̅_F_n_3)=A̅_F_n_3.

It follows that

    𝒜*ℬ=Fold_vec(Circ(𝒜)×Vec(ℬ)) 
        =    Fold_vec((F^H_n_3⊗ I_n_1)×Block(𝒜̅_F_n_3) × (F_n_3⊗ I_n_2)×Vec(ℬ))
        =    fold((F^H_n_3⊗ I_n_1)×Block(𝒜̅_F_n_3) ×Block(ℬ̅_F_n_3))
        =    F^H_n_3[fold(Block(𝒜̅_F_n_3)Block(ℬ̅_F_n_3))]
        =    F^H_n_3[fold(A̅_F_n_3B̅_F_n_3)]
        =    𝒜⋆_F_n_3ℬ .


The definitions of the conjugate transpose of tensor, the identity tensor, the unitary tensor, the invertible tensor, the diagonal tensor, and the inner product between
two tensors related to the ⋆_L-product are given as follows.

	
  * <cit.> The conjugate transpose of 𝒜∈ℂ^n_1× n_2× n_3 with respect to L is the tensor 𝒜^H∈ℂ^n_1× n_2× n_3 obtained by 𝒜^H=L^H[fold(Block(𝒜̅_L)^H)]=L^H[fold(A̅^H_L)].
	
	
  * <cit.> The identity tensor ℐ_L∈ℂ^n× n× n_3 (with respect to L) is defined to be a
	tensor such that ℐ_L=L^H[ℐ], where ℐ with each frontal slice being the n× n identity matrix.
	
	
  * <cit.> A tensor 𝒬∈ℂ^n× n× n_3 is unitary with respect to ⋆_L-product if it satisfies 𝒬^H⋆_L𝒬=𝒬⋆_L𝒬^H=ℐ_L, where ℐ_L is the identity tensor.
	
	
  * <cit.> For tensors 𝒜∈ℂ^n× n× n_3,ℬ∈ℂ^n× n× n_3, if 𝒜⋆_Lℬ=ℬ⋆_L𝒜=ℐ_L, then tensor ℬ is the invertible tensor under the ⋆_L-product of tensor 𝒜.
	
  * <cit.> A tensor is a diagonal tensor if each frontal slice of the tensor is a diagonal matrix. For a third-order tensor, if all of its frontal slices are upper or lower triangles, then the tensor is called f-upper or f-lower.




	<cit.> Suppose that 𝒜∈ℝ^m× k × p and ℬ∈ℝ^k× n × p are two arbitrary tensors. Let 𝒵=𝒳⋆_L𝒴. Then the following properties hold:
	
	
		
  (1) 𝒜 _F^2=1/pA̅_F_n_3 _F^2=1/pA̅_L _F^2=1/p∑_i=1^p𝒜̅^(i)_L _F^2.
		
  (2) 𝒵=𝒳⋆_L𝒴 is equivalent to Block(𝒵̅_L)=Block(𝒜̅_L)Block(ℬ̅_L)(i.e., Z̅_L=A̅_LB̅_L).
		
	





	




	(Gaussian random tensor) <cit.> A tensor 𝒢∈ℝ^m× n× p
	is called a Gaussian random tensor if the elements of 𝒢^(1) satisfy the standard normal distribution
	(i.e., Gaussian with mean zero and variance one) and the other frontal slices are all
	zero.
	


Based on the above
definitions, we have the following transformed tensor SVD with respect to L.
<cit.>
	For any 𝒜∈ℂ^n_1× n_2× n_3, the transformed tensor SVD is given by
	
	
    𝒜=𝒰⋆_L𝒮⋆_L𝒱^H,

	where 𝒰∈ℂ^n_1× n_1× n_3 and 𝒱∈ℂ^n_2× n_2× n_3 are unitary tensors with respect to the ⋆_L-product, and 𝒮∈ℂ^n_1× n_2× n_3 is diagonal.



The transformed tubal rank, denoted as
rank_L(𝒜), is defined as the number of nonzero singular tubes of 𝒮, where 𝒮 comes from the transformed tensor SVD of 𝒜=𝒰⋆_L𝒮⋆_L𝒱^H, i.e.,

    rank_L(𝒜)=#{i:𝒮(i,i,:) 0},

where # denotes the cardinality of a set.

The transformed tensor SVD could be implemented efficiently by the SVDs of the frontal slices
in the transformed domain. We also refer the readers to <cit.> for more details about the
computation of the transformed tensor SVD.

For the Kernfeld-Kilmer transformed tensor SVD (i.e., Eq. (<ref>)), by <cit.>, we have

    ∑_k=1^p𝒮(1,1,k)^2≥∑_k=1^p𝒮(2,2,k)^2≥…≥∑_k=1^p𝒮(min{m,n}, min{m,n},k)^2.


	(Transformed Tensor Singular Values) 
	Suppose 𝒜∈ℝ^m× n× p with a Kernfeld-Kilmer transformed tensor SVD such that Eq. (<ref>) is satisfied. The i-th largest transformed tensor singular value of 𝒜 is defined as
	
	
    σ_i=√(∑_k=1^p𝒮(i,i,k)^2),    for i=1,2,…, min{m,n}.

	













Similarly to the definition of the matrix tail energy, the tail energy of tensor is defined below.



	(Tail Energy) For tensor 𝒜∈ℝ^m× n× p, the i-th largest transformed tensor singular value is σ_i, i=1,…, min{m,n}. Then the j-th tail energy of 𝒜 is defined as
	
	
    τ_j^2(𝒜):=min_ rank_L(ℬ)<j𝒜-ℬ_F^2=∑_i≥ jσ_i^2(𝒜) .
	
	


According to the above definition and using Lemma <ref> and the linearity, we can obtain

    τ_j^2(𝒜)
    	    =    ∑_i≥ jσ_i^2(𝒜)=∑_i≥ j𝒜_i-𝒜_i-1_F^2
        =    ∑_i≥ j1/p∑_k=1^p𝒜̅_(i)_L-𝒜̅_(i-1)_L_F^2
        =    1/p∑_i=1^pτ_j^2(𝒜̅_L^(i)),

where 𝒜̅_(i)_L represents a third-order tensor obtained via being multiplied by L on all tubes along the third-dimension of 𝒜_(i). We then have the following proposition.
 
	Suppose 𝒜∈ℝ^m× n× p and 𝒜̅_L represents a third-order tensor obtained via being multiplied by L on all tubes along the third-dimension of 𝒜. Let j be a positive integer satisfying j≤min{m,n}. Then
	
	
    τ_j^2(𝒜)=1/p∑_i=1^pτ_j^2(𝒜̅_L^(i)) .
	











 §.§ Random Projection


This subsection presents three matrix sketching techniques, i.e., Gaussian projection, subsampled randomized Hadamard transform (SRHT), and count sketch <cit.>.



  §.§.§ Gaussian Projection

The n × s Gaussian random projection matrix S is a matrix formed by S = G/√(s), where
each entry of G is sampled i.i.d. from 𝒩 (0, 1). For matrix A∈ℝ^m × n, the time complexity of Gaussian projection is 𝒪(mns); see Algorithm 1 for the Gaussian projection process.


	
	
	
	
	
	
	

  §.§.§ Subsampled Randomized Hadamard Transform

	The SRHT matrix is defined by S = DH_nP/√(sn), where
	
	
		
  * D ∈ℝ^n × n
		is a diagonal matrix with diagonal entries sampled uniformly from {+1, 1}.
		
  * H_n∈ℝ^n × n
		is defined recursively by 
    H_n=[  H_n/2  H_n/2;  H_n/2 -H_n/2 ]   and   H_2=[ +1 +1; +1 -1 ] .

		∀ y ∈ℝ^n, the matrix vector product y^HH_n can be performed in 𝒪(n log n) by the fast Walsh-Hadamard transform algorithm in a divide-and-conquer fashion.
		
  * P ∈ℝ^n × s samples s from the n columns.
	
	
	Algorithm 2 is not introduced in the main text here.
	
	
	
	

  §.§.§  Count Sketch

	
	There are two different ways to implement count sketch in “map-reduce fashion” and “streaming fashion”. These two approaches are equivalent, and thus we only discuss the streaming fashion below. The streaming fashion has two steps. Initially, the m × s matrix C is set to all zero. Then, for each column of A, the sign is randomly flipped with a probability of 0.5 and added to a randomly selected column of C. This is outlined in Algorithm 3.
	
	The streaming fashion keeps the sketch C in memory and scans the data A in a single pass. If A does not fit in memory, this approach is more efficient than the map-reduce fashion as it scans the columns in sequence. Additionally, if A is a sparse matrix, randomly accessing the entries may not be efficient, and then it is better to access the columns sequentially.
	
	
	
	
	
	It is noticed that the count sketch does not explicitly form the sketching
	matrix S. In fact, S is such a matrix that each row has only one nonzero entry. The time complexity of Count sketch is 𝒪(nnz(A)) (nnz(A) represent number of nonzeros of the matrix A).
	
	
	

 §.§ Tensor Sketching Operator

	
	Using the three matrix sketching operators introduced in the previous section, three corresponding tensor sketching operators can be generated. As for the new two-sided sketching algorithm, we need to generate four random linear dimension reduction maps, i.e.,
	
    Υ∈ℝ^k× m× p,   Ω∈ℝ^k× n × p,   Φ∈ℝ^s× m× p,    and Ψ∈ℝ^s× n × p .

	
	Different ways of generating tensor sketch operators are shown below.
	
		
  * Gaussian tensor sketching operator.
		Set

		Υ=zeros(k,m,p);   Υ(:,:,1)=(𝒜(:,:,1),k)';

		Ω=zeros(k,n,p); Ω(:,:,1)=(𝒜(:,:,1)',k)';
Φ=zeros(s,m,p);  Φ(:,:,1)=(𝒜(:,:,1),s)';

		Ψ=zeros(s,n,p);  Ψ(:,:,1)=(𝒜(:,:,1)',s)'.
 
		Then Υ, Ω, Φ and Ψ are said to be Gaussian tensor sketching operators. 
		
		
		
  * SRHT tensor sketching operator.
		Set
 
		Υ=zeros(k,m,p); Υ(:,:,1)=(𝒜(:,:,1),k)';

		Ω=zeros(k,n,p);  
		Ω(:,:,1)=(𝒜(:,:,1)',k)';
Φ=zeros(s,m,p); Φ(:,:,1)=(𝒜(:,:,1),s)';

		Ψ=zeros(s,n,p); Ψ(:,:,1)=(𝒜(:,:,1)',s)'. 

		Then Υ, Ω, Φ and Ψ are said to be SRHT tensor sketching operators.
		
		
  *  Count tensor sketching operator.
		For i = 1,2,…, p, set

		Υ=zeros(k,m,p);  Υ(:,:,i)=(𝒜(:,:,1), k)';

		Ω=zeros(k,n,p);  Ω(:,:,i)=(𝒜(:,:,1)', k)';
Φ=zeros(s,m,p);  Φ(:,:,i)=(𝒜(:,:,1), s)';

		Ψ=zeros(s,n,p);  Ψ(:,:,i)=(𝒜(:,:,1)', s)'.

		Then Υ, Ω, Φ and Ψ are said to be count tensor sketching operators.
	





§ PROPOSED THE NEW TWO-SIDED SKETCHING ALGORITHMS BASED ON THE TRANSFORMED DOMAIN


In this section, we present the framework of our new two-sided sketching algorithm for low tubal rank tensor approximation. 

Given the input tensor 𝒜∈ℝ^m× n× p and the objective tubal rank k, using the appropriate tensor sketching operators in Eq. (<ref>), we may realize the randomized sketches (𝒳,𝒴,𝒵)  such as


    𝒳:=Υ⋆_L𝒜∈ℝ^k× n × p,  𝒴:=𝒜⋆_LΩ^H∈ℝ^m× k × p;


    𝒵:=Φ⋆_L𝒜⋆_LΨ^H∈ℝ^s× s × p.


The first two tensor sketches 𝒳 and 𝒴 capture the co-range and range of 𝒜, respectively. The core sketch 𝒵 contains new information that improves our estimates of the transformed tensor singular values and the transformed tensor singular vectors of 𝒜 and is responsible for the superior performance of the algorithms.

Once the sketches (𝒳,𝒴,𝒵) of the input tensor 𝒜 are obtained, we can find the low-rank approximation 𝒜̂ by following a three-step process below.

1. Form an L-orthogonal-triangular factorization

    𝒳^H:=𝒫⋆_Lℛ_1,   𝒴:=𝒬⋆_Lℛ_2,

where 𝒫∈ℝ^n× k × p and 𝒬∈ℝ^m× k × p
are  partially orthogonal tensors, and ℛ_1∈ℝ^k× k × p and ℛ_2∈ℝ^k× k × p
are f-upper
triangular tensors, in the sense of the ⋆_L-product operation.

2. Solve the least-squares problem based on ⋆_L-product

    min_ C1/2𝒢⋆_L𝒞⋆_Lℳ^H-ℬ_F^2 .

The above objective function (<ref>) can be reformulated as

    min_𝒞̅_L1/2p𝒢̅_L𝒞̅_Lℳ̅^H_L-ℬ̅_L_F^2 .


Then there is 𝒞̅^(i)_L=(𝒢̅^(i)_L)^†ℬ̅^(i)_L((ℳ̅ ^(i)_L)^†)^H. According to Lemma <ref>, the least-squares solution of problem (<ref>) is 𝒞 = 𝒢^†⋆_Lℬ⋆_ L(ℳ^†)^H. Thus we have

    𝒞:=(Φ⋆_L𝒬)^†⋆_L𝒵⋆_L ((Ψ⋆_L𝒫)^†)^H∈ℝ^k× k × p .


3. Construct the transformed tensor tubal rank k approximation

    𝒜̂:=𝒬⋆_L𝒞⋆_L𝒫^H .


The storage cost for the sketches (𝒳,𝒴,𝒵) is p(nk+mk+s^2) floating point numbers. The storage complexity of LTRP-SKETCH Algorithms (i.e., Algorithm 4)

for the original data 𝒜∈ℝ^m× n× p is 𝒪(mkp+nkp+kkp) floating point numbers. The algorithm pseudo-code is given in Algorithm 4.



The original data 𝒜∈ℝ^m× n× p modulo-3 is expanded into a matrix W∈ℝ^p× mn, 

and then performing SVD on W to obtain the unitary transformation matrix U. The U Transformed Domain skethching algorithm (U-TRP-SKETCH) is based on a linear transformation L that uses a discrete cosine transformation matrix C_n_3. When four random linear dimensionality reduction mappings in Eq. (<ref>) are chosen as the count tensor sketching operator, the Gaussian tensor sketching operator, and the SRHT tensor sketching operator, the algorithm is referred to as the U Count-Sketch, the U Gaussian-Sketch, and the U SRHT-Sketch algorithms, respectively. 

Analogously, the DFT Transformed Domain skethching algorithm (DFT-TRP-SKETCH) is based on a linear transformation L that uses a fast Fourier transform matrix F_n_3. When four random linear dimensionality reduction mappings in Eq. (<ref>) are chosen as the count tensor sketching operator, the Gaussian tensor sketching operator, and the SRHT tensor sketching operator, the algorithm is referred to as the DFT Count-Sketch, the DFT Gaussian-Sketch, and the DFT SRHT-Sketch algorithms, respectively.



 §.§ Principle of the new two-sided Sketching Algorithms



    𝒜≈𝒬⋆_L(𝒬^H⋆_L𝒜⋆_L𝒫)⋆_L𝒫^H.

The core tensor 𝒬^H⋆_L𝒜⋆_L𝒫 cannot be calculated directly from the linear sketch since 𝒫 and 𝒬 are functions of 𝒜. Using the approximation in Eq. (<ref>), the core sketch 𝒵 estimating the core tensor can be achieved by

    𝒵 = Φ⋆_L𝒜⋆_LΨ^H≈ (Φ⋆_L𝒬)⋆_L(𝒬^H⋆_L𝒜⋆_L𝒫)⋆_L(𝒫^H⋆_LΨ^H).

Transferring the external matrix to the left-hand side, the core approximation 𝒞 defined in Eq. (<ref>) is found to satisfy

    𝒞=(Φ⋆_L𝒬)^†⋆_L𝒵⋆_L ((Ψ⋆_L𝒫)^†)^H≈𝒬^H⋆_L𝒜⋆_L𝒫.

Given Eq. (<ref>) and (<ref>), we have

    𝒜≈𝒬⋆_L(𝒬^H⋆_L𝒜⋆_L𝒫)⋆_L𝒫^H≈𝒬⋆_L𝒞⋆_L𝒫^H = 𝒜̂.

The error in the last relation depends on the error in the best transformed tubal rank k approximation of 𝒜.










	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	




 §.§ The New Two-Sided Sketching Algorithm Based on the DCT Transformed Domain





The DCT transformed domain sketching algorithm 
(DCT-TRP-SKETCH) is based on a linear transformation L that uses a DCT matrix C_n_3. When four random linear dimensionality reduction mappings  in Eq. (<ref>) are chosen as the count tensor sketching operator, the Gaussian tensor sketching operator, and the SRHT tensor sketching operator, the algorithm is referred to as the DCT count-sketch, the DCT Gaussian-sketch, and the DCT SRHT-sketch algorithms, respectively.








 §.§ Improvement with the Power Iteration Technique


As shown in <cit.>, the power iteration technique is useful to improve sketching algorithms for low rank matrix approximation. Here, we combine the power iteration technique with the DCT Gaussian-sketch algorithm, in which we exploit the third order tensor 𝒵=(𝒜⋆_L𝒜^H)^q ⋆_L𝒜 (q is a nonnegative integer) instead of the original tensor 𝒜, and the DCT Gaussian-Sketch algorithm is applied to the new tensor 𝒵. According to the transformed tensor SVD, 𝒜 = 𝒰⋆_L𝒮⋆_L𝒱^H, we have 𝒵=𝒰⋆_L (𝒮)^2q+1⋆_L𝒱^H. Therefore, the transformed tensor singular values of 𝒵 have a faster decay rate. This can improve the solution obtained by the DCT Gaussian-sketch algorithm. The DCT Gaussian-sketch algorithm can be equipped with the power iteration technique when the transformed tensor singular values do not decay fast. We summarize the resulting scheme in Algorithm 6. 









§ ERROR BOUND OF THE DCT GAUSSIAN-SKETCH ALGORITHM


In this section, taking the DCT Gaussian-Sketch algorithm as an example, we derive the error bound of the new two-sided sketching algorithm.


A fact about random tensors are firstly given below.


	
	Assume t > q+1. Let 𝒢_1∈ℝ^t× q × p and 𝒢_2∈ℝ^t× l × p be Gaussian random tensors. For any tensor ℬ with conforming dimensions,
	
    𝔼∥𝒢_1^†⋆_L𝒢_2⋆_Lℬ∥_F^2 = q/t-q-1∥ℬ∥_F^2 .






 §.§ Results from Randomized Linear Algebra


	
	Fix 𝒜∈ℝ^m× n× p and let ϱ < k be a natural number.
	Then the tensor 𝒬∈ℝ^m× k× p calculated by Eq. (<ref>) satisfies
	
    𝔼_Ω𝒜-𝒬⋆_L𝒬^H⋆_L𝒜_F^2≤ (1+f(ϱ,k))τ_ϱ+1^2(𝒜⋆_L𝒜^H)  .

	An analogous result holds for the tensor 𝒫∈ℝ^n× k× p
	computed by Eq. (<ref>), i.e.,
	
    𝔼_Υ𝒜-𝒜⋆_L𝒫⋆_L𝒫^H_F^2≤ (1+f(ϱ,k))τ_ϱ+1^2(𝒜⋆_L𝒜^H)  ,

	where f(ϱ,k) = ϱ/(k-ϱ-1).

proof  By Lemma <ref> and the linearity of the expectation, we have
	
    𝔼_Ω𝒜 - 𝒬_L⋆_L𝒬^H_L⋆_L𝒜_F^2
    		    =    1  p(∑_i=1^p 𝔼𝒜̅^(i)_L - 𝒬̅^(i)_L(𝒬̅^(i)_L)^H𝒜̅^(i)_L_F^2).

	By <cit.>, we have
	
    𝔼𝒜̅^(i)_F_n_3-𝒬̅^(i)_F_n_3(𝒬̅^(i)_F_n_3)^H𝒜̅^(i)_F_n_3_F^2≤ (1+f(ϱ,k))·τ_ϱ+1^2(𝒜̅^(i)_F_n_3(𝒜̅^(i)_F_n_3)^H)  .

	Similarly, we know from <cit.> that multiplying by a unitary matrix will not change the Frobenius norm, i.e., the Frobenius norm keeps the unitary matrix unchanged. Thus we have
	
    𝔼𝒜̅^(i)_L-𝒬̅^(i)_L(𝒬̅^(i)_L)^H𝒜̅^(i)_L_F^2≤ (1+f(ϱ,k))·τ_ϱ+1^2(𝒜̅^(i)_L(𝒜̅^(i)_L)^H)

	and
	
    𝔼_Ω𝒜 - 𝒬_L⋆_L𝒬^H_L⋆_L𝒜_F^2
    		   ≤1/p(1+f(ϱ,k))·(∑_i=1^pτ_ϱ+1^2(𝒜̅^(i)_L(𝒜̅^(i)_L)^H))
        = (1+f(ϱ,k))·τ_ϱ+1^2(𝒜⋆_L𝒜^H)  .

	Similarly, there is 
	
    𝔼_Υ𝒜-𝒜⋆_L𝒫⋆_L𝒫^H_F^2≤ (1+ϱ/k-ϱ-1)τ_ϱ+1^2(𝒜⋆_L𝒜^H)  ,

	which completes the proof.



Let 𝒜̂ be the transformed tensor tubal rank k approximation of 𝒜 obtained by the DCT Gaussian-sketch algorithm. We now split the error 𝒜-𝒜̂ _F^2 into two parts by the proposition below.
 
	Let 𝒜̂ be the transformed tensor tubal rank k approximation of 𝒜∈ℝ^m× n× p obtained by the DCT Gaussian-sketch algorithm, with 𝒬, 𝒞 and 𝒫 being the intermediate tensors obtained by the DCT Gaussian-Sketch algorithm satisfying 𝒬^H⋆_L𝒬=ℐ_kkp,𝒫⋆_L𝒫^H=ℐ_nnp. Then
	
    𝒜-𝒜̂ _F^2=𝒜-𝒬⋆_L𝒬^H⋆_L𝒜⋆_L𝒫⋆_L𝒫^H_F^2+𝒞-𝒬^H⋆_L𝒜⋆_L𝒫 _F^2 .
		

proof  Since 𝒜̂=𝒬⋆_L𝒞⋆_L𝒫^H, we have
	
    𝒜 - 𝒜̂_F^2 
        =    1/p∑_i=1^p𝒜̅_Φ-𝒜̅̂̅_Φ_F^2
        =    1/p∑_i=1^p(𝒜̅_Φ-𝒬̅_Φ𝒬̅^H_Φ𝒜̅_Φ𝒫̅_Φ𝒫̅^H_F^2+
    		𝒞̅_Φ-𝒬̅^H_Φ𝒜̅_Φ𝒫̅_Φ_F^2) 
        =    𝒜-𝒬⋆_L𝒬^H⋆_L𝒜⋆_L𝒫⋆_L𝒫^H_F^2+𝒞-𝒬^H⋆_L𝒜⋆_L𝒫 _F^2 .

	The first and last equations are known by Lemma <ref>, and the second equation is known by <cit.> Appendix A.

























	For any natural number ϱ < k-1, it holds that
	
    𝔼_Υ𝔼_Ω𝒜-𝒬⋆_L𝒬^H⋆_L𝒜⋆_L𝒫⋆_L𝒫^H_F^2≤(1+2ϱ/k-ϱ-1)τ_ϱ+1^2(𝒜⋆_L𝒜^H)  .

	
	
	
	
	
	
	

proof  We have
	
    𝔼_Υ𝔼_Ω𝒜-𝒬⋆_L𝒬^H⋆_L𝒜⋆_L𝒫⋆_L𝒫^H_F^2
        =    1/p∑_i=1^p𝔼_Υ𝔼_Ω𝒜̅_Φ-𝒬̅_Φ𝒬̅^H_Φ𝒜̅_Φ𝒫̅_Φ𝒫̅^H_F^2
       ≤   (1+2ϱ/k-ϱ-1)τ_ϱ+1^2(𝒜⋆_L𝒜^H)  .

	
	The first equation is known by Lemma <ref>, and the inequality is known by <cit.> Appendix A.






 §.§ Decomposition of the Core Tensor Approximation Error



The first step in the argument is to obtain a formula for the error in the approximation (𝒞-𝒬^H⋆_L𝒜⋆_L𝒫).
The core tensor 𝒞∈ℝ^k× k × p
is defined in Eq. (<ref>), and the
orthonormal tensors 𝒫∈ℝ^n× k × p and 𝒬∈ℝ^m× k × p are constructed in
Eq. (<ref>).

Let us introduce tensors whose ranges are complementary to those of 𝒫 and 𝒬, i.e.,


    𝒫_⊥⋆_L𝒫_⊥^H    = ℐ_mmp-𝒫⋆_L𝒫^H,𝒫_⊥∈ℝ^n× (n-k)× p, 
    𝒬_⊥⋆_L𝒬_⊥^H    = ℐ_mmp-𝒬⋆_L𝒬^H,𝒬_⊥∈ℝ^m× (m-k)× p.

The columns of 𝒫_⊥ and 𝒬_⊥ are orthonormal, separately. Next, we introduce the subtensors, i.e.,


    Φ_1:=Φ⋆_L𝒬∈ℝ^s× k× p,  Φ_2:=Φ⋆_L𝒬_⊥∈ℝ^s× (m-k)× p, 
    Ψ_1^H:=𝒫^H⋆_LΨ^H∈ℝ^k× s× p,  Ψ_2^H:=𝒫^H_⊥⋆_LΨ^H∈ℝ^(n-k)× s× p.

With these notations at hand, we can state and prove the lemma below.


	(Decomposition of the Core Tensor Approximation) Assume the tensor transformed tubal rank of Φ_1 and Ψ_1 is k and s, respectively. Then
	
	
	
	
    𝒞-𝒬^H⋆_L𝒜⋆_L𝒫
    
    		=     Φ_1^†⋆_LΦ_2⋆_L (𝒬^H_⊥⋆_L𝒜⋆_L𝒫)+(𝒬^H⋆_L𝒜⋆_L𝒫_⊥)⋆_LΨ_2^†⋆_L ( Ψ^†_1)^H
        +Φ_1^†⋆_LΦ_2⋆_L(𝒬^H_⊥⋆_L𝒜⋆_L𝒫_⊥)⋆_LΨ_2^†⋆_L ( Ψ^†_1)^H.


proof  Adding and subtracting terms, we write the core sketch 𝒵 as
	
    𝒵=     Φ⋆_L𝒜⋆_LΨ^H
    
    		=     Φ⋆_L (𝒜-𝒬⋆_L𝒬^H⋆_L𝒜⋆_L𝒫⋆_L𝒫^H)⋆_LΨ^H
        + ( Φ⋆_L𝒬)⋆_L (𝒬^H⋆_L𝒜⋆_L𝒫)⋆_L (𝒫^H⋆_LΨ^H).

	Using Eq. (<ref>), we identify the tensors Φ_1 and Ψ_1. After left-multiplying by Φ^†_1 and right-multiplying by (Ψ^†_1)^H, we have
	
    𝒞    =  Φ^†_1⋆_L𝒵⋆_L (Ψ^†_1)^H
        =  Φ^†_1⋆_LΦ⋆_L (𝒜-𝒬⋆_L𝒬^H⋆_L𝒜⋆_L𝒫⋆_L𝒫^H)⋆_LΨ^H⋆_L (Ψ^†_1)^H
           +𝒬^H⋆_L𝒜⋆_L𝒫.
  
	
	We have identified the core tensor 𝒞 defined in (<ref>). Moving the term 𝒬^H⋆_L𝒜⋆_L𝒫 to the
	left-hand side to isolate the approximation error, we have
	
	
	
    Φ^†_1⋆_LΦ    = Φ^†_1⋆_LΦ⋆_L𝒬⋆_L𝒬^H+Φ^†_1⋆_LΦ⋆_L𝒬_⊥⋆_L𝒬^H_⊥
        = 𝒬^H+Φ^†_1⋆_LΦ_2⋆_L𝒬^H_⊥ .

	
	Likewise,
	
    Ψ^H⋆_L (Ψ^†_1)^H
    
    		=     𝒫⋆_L𝒫^H⋆_LΨ^H⋆_L (Ψ^†_1)^H+𝒫_⊥⋆_L𝒫_⊥^H⋆_LΨ^H⋆_L (Ψ^†_1)^H
    
    		=     𝒫+𝒫_⊥⋆_LΨ_2^H⋆_L(Ψ^†_1)^H.

	
	Combining the last three displays, we have
	
    𝒞-𝒬^H⋆_L𝒜⋆_L𝒫
    
    		=     (𝒬^H+Φ^†_1⋆_LΦ_2⋆_L𝒬_⊥^H)⋆_L (𝒜-𝒬⋆_L𝒬^H⋆_L𝒜⋆_L𝒫⋆_L𝒫^H)⋆_L (𝒫
       +𝒫_⊥⋆_LΨ_2^H⋆_L(Ψ^†_1)^H).

	Expanding the expression and using the orthogonality relations 𝒬^H⋆_L𝒬=ℐ_kkp, 𝒬_⊥^H⋆_L𝒬=𝒪, 𝒫^H⋆_L𝒫=ℐ_kkp, and 𝒫_⊥^H⋆_L𝒫=𝒪, we complete the proof.
	





 §.§ Probabilistic Analysis of the Core Tensor
	

We can then study
the probabilistic behavior of the error (𝒞-𝒬^H⋆_L𝒜⋆_L𝒫),  conditional on 𝒬 and 𝒫.




	(Probabilistic Analysis of the Core Tensor) Assume that the dimension reduction tensors Φ and Ψ are Gaussian linear sketching operators. When s ≥ k, it holds that
	
	
    𝔼_Φ,Ψ[ 𝒞-𝒬^H⋆_L𝒜⋆_L𝒫 ]=𝒪.

	
	When s > k + 1, the error can be expressed as
	
	
    𝔼_Φ,Ψ𝒞-𝒬^T⋆_L𝒜⋆_L𝒫 _F^2
    
    		=     k/s-k-1𝒜-𝒬⋆_L𝒬^H⋆_L𝒜⋆_L𝒫⋆_L𝒫^H _F^2
        +k(2k+1-s)/(s-k-1)^2𝒬_⊥^H⋆_L𝒜⋆_L𝒫_⊥ _F^2 .

	
	In particular, when s < 2k + 1, the last term is nonnegative; and  when s ≥ 2k + 1, the last term is nonpositive.
	
	
proof  Since Φ is a Gaussian linear sketching operator, the orthogonal subtensors Φ_1 and Φ_2 are also	Gaussian linear sketching operators because of the marginal property of the normal distribution. Likewise, Ψ_1 and Ψ_2 are Gaussian linear sketching operators. Provided that s ≥ k, the tensor transformed tubal rank
	of Φ_1 and Ψ_1 is k.
	Thus,
	
    𝔼_Φ,Ψ[ 𝒞-𝒬^H⋆_L𝒜⋆_L𝒫 ]
    
    		=     𝔼_Φ_1𝔼_Φ_2[ Φ_1^†⋆_LΦ_2⋆_L (𝒬_⊥^H⋆_L𝒜⋆_L𝒫) ]
        +𝔼_Ψ_1𝔼_Ψ_2[ (𝒬^H⋆_L𝒜⋆_L𝒫_⊥)⋆_LΨ_2^H⋆_L(Ψ_1^†)^H ]
        + 𝔼[ Φ_1^†⋆_LΦ_2⋆_L (𝒬_⊥^H⋆_L𝒜⋆_L𝒫_⊥)⋆_LΨ_2^H⋆_L( Ψ_1^†)^H]  .

	
	We have used the decomposition of the approximation error from Lemma <ref>. Then
	we invoke independence to write the expectations as iterated expectations. Since Φ_2
	and Ψ_2 have to mean zero. This formula makes it clear that the approximation error has to
	mean zero.
	
	To study the fluctuations, applying the independence and zero-mean property of Φ_2
	and Ψ_2, we have
	
	
    𝔼_Φ,Ψ𝒞-𝒬^H⋆_L𝒜⋆_L𝒫_F^2
    
    		=     𝔼_ΦΦ_1^†⋆_LΦ_2⋆_L (𝒬_⊥^H⋆_L𝒜⋆_L𝒫) _F^2
        +𝔼_Ψ (𝒬^H⋆_L𝒜⋆_L𝒫_⊥)⋆_LΨ_2^H⋆_L (Ψ_1^†)^H_F^2
        + 𝔼_Φ𝔼_ΨΦ_1^†⋆_LΦ_2⋆_L (𝒬_⊥^H⋆_L𝒜⋆_L𝒫_⊥)⋆_LΨ_2^H⋆_L (Ψ_1^†)^H_F^2.

	Invoking Proposition <ref> yields
	
    𝔼_Φ,Ψ𝒞-𝒬^H⋆_L𝒜⋆_L𝒫_F^2
    
    		=     k/s-k-1[𝒬_⊥^H⋆_L𝒜⋆_L𝒫_F^2+𝒬^H⋆_L𝒜⋆_L𝒫_⊥_F^2
        + k/s-k-1𝒬_⊥^H⋆_L𝒜⋆_L𝒫_⊥_F^2]
    
    		=     k/s-k-1[𝒬_⊥^H⋆_L𝒜⋆_L𝒫_F^2+𝒬^H⋆_L𝒜⋆_L𝒫_⊥_F^2+𝒬_⊥^H⋆_L𝒜⋆_L𝒫_⊥_F^2
        + 2k+1-s/s-k-1𝒬_⊥^H⋆_L𝒜⋆_L𝒫_⊥_F^2]  .

	Using the Pythagorean Theorem to combine the terms on the above equation, 
	
	we have
	
	
    𝔼_Φ,Ψ𝒞-𝒬^H⋆_L𝒜⋆_L𝒫 _F^2
    
    		=     k/s-k-1𝒜-𝒬⋆_L𝒬^H⋆_L𝒜⋆_L𝒫⋆_L𝒫^H _F^2
        +k(2k+1-s)/(s-k-1)^2𝒬_⊥^H⋆_L𝒜⋆_L𝒫_⊥ _F^2 .

	This completes the proof. 



	Assume that the sketch parameters satisfy s≥ 2k+1 and 𝒜̂ is the transformed tensor tubal rank-k approximation of 𝒜 defined by (<ref>). Then
	
	
    𝔼𝒜-𝒜̂ _F^2≤ (1+f(k,s))(1+2ϱ/k-ϱ-1)τ_ϱ+1^2(𝒜⋆_L𝒜^H),

	where ϱ is a natural number less than k-1, f(ϱ,k):=ϱ/(k-ϱ-1), and the tail energy τ^2_ϱ+1
	is defined by Definition <ref>.
	

proof  We have
	
    𝔼𝒜-𝒜̂ _F^2
    
    		=     𝔼_Υ𝔼_Ω𝒜-𝒬⋆_L𝒬^H⋆_L𝒜*𝒫⋆_L𝒫^H _F^2+𝔼𝒞-𝒬^H⋆_L𝒜⋆_L𝒫 _F^2
    
    		=     (1+f(k,s))𝔼_Υ𝔼_Ω𝒜-𝒬⋆_L𝒬^H⋆_L𝒜⋆_L𝒫⋆_L𝒫^H _F^2
        + k(2k+1-s)/(s-k-1)^2𝔼𝒬_⊥^H⋆_L𝒜⋆_L𝒫_⊥ _F^2
    ≤    (1+f(k,s))(1+2ϱ/k-ϱ-1)τ_ϱ+1^2(𝒜⋆_L𝒜^H) 
        +k(2k+1-s)/(s-k-1)^2𝔼𝒬_⊥^H⋆_L𝒜⋆_L𝒫_⊥ _F^2
    ≤    (1+f(k,s))(1+2ϱ/k-ϱ-1)τ_ϱ+1^2(𝒜⋆_L𝒜^H).

	The first equation is known by Proposition <ref>, the second is known by Lemma <ref>, the first inequality is known by Theorem <ref>, and the last inequality is because we require s≥ 2k+1, the missing item k(2k+1-s)/(s-k-1)^2𝔼𝒬_⊥^H⋆_L𝒜⋆_L𝒫_⊥ _F^2
	
	is negative. This completes the proof.
	






§ NUMERICAL EXPERIMENTS
	

This section presents numerical experiments to verify the performance of the proposed the new two-sided sketching algorithm, i.e., L-TRP-SKETCH. The following relative error ϵ_ err and the peak signal-to-noise ratio (PSNR) ρ_ psnr are used as metrics of the low-rank approximation to the input tensor data, i.e.,

    ϵ_ err :=||𝒜-𝒜̂_ out||_F^2/||𝒜||_F^2,
    	  PSNR:=10log_10n_1n_2n_3||𝒜||_∞^2/||𝒜-𝒜̂_out||_F^2,

where 𝒜 and 𝒜̂_ out are the original tensor and the low-rank approximation, respectively. 

The sketch size	parameters are set as s = 2k+1.

Extensive numerical experiments on synthetic and real-world data are conducted and compared with other algorithms. For example, the two-sided sketching algorithm <cit.> based on the t-product with the power iteration technique, the truncated t-SVD algorithm <cit.>, and the
rt-SVD algorithm <cit.> (i.e., a one-sided randomized algorithm based on t-SVD).






 §.§ Synthetic Experiment

We firstly conduct numerical tests on some synthetic input tensors 𝒜∈ℝ^10^3× 10^3× 10 with decaying spectrum.















Polynomial decay:	These tensors are f-diagonal tensors. Considering their j-th frontal slices with the form 
    𝒜^(j)=(1,…,1_min(r,j),
    2^-p,3^-p,4^-p,…,(n-min(r,j)+1)^-p) ∈ℝ^n× n,

we study three examples, i.e., PolyDecaySlow (p = 0.5), PolyDecayMed (p = 1), and PolyDecayFast (p = 2).







Figures <ref>–<ref> show examples regarding the relative error, PSNR and CPU time of the T-Sketch algorithm (Algorithm 2 in <cit.>), T-Sketch algorithm PI (q=1) <cit.>, truncated t-SVD algorithm <cit.>, rt-SVD algorithm (Algorithm 6 in <cit.>) and our L-TRP-SKETCH algorithm, as the size of the sketch parameter k varies.







From Figures <ref>–<ref>, we can see that the best one among the new two-sided Gaussian sketching algorithms based on the transformed domain (referred to the U Gaussian-Sketch, DCT Gaussian-Sketch, and DFT Gaussian-Sketch algorithms) is the DCT Gaussian-Sketch algorithm. Thus, the new two-sided Gaussian sketching algorithms based on the transformed domain chose the DCT transform to work best. 

In terms of CPU time, the DCT Gaussian-Sketch method is the fastest. In particular, as shown in Figures <ref>–<ref>, the accuracy of the DCT Gaussian-Sketch method is better than the rt-SVD method and T-Sketch method for input tensor with PolyDecayMed decay spectrum and PolyDecayFast decay spectrum. In this case, with less storage and manipulation, the DCT Gaussian-Sketch method achieves better accuracy for the low-rank approximation. The DCT Gaussian-Sketch method is second only to the truncated t-SVD method in terms of accuracy but with faster speed.







From Figures <ref>–<ref>, we see that the DCT Gaussian-Sketch algorithm has the best accuracy among the new two-sided sketching algorithms based on the DCT transformed domain. The new two-sided sketching algorithms based on the DCT transformed domain select the Gaussian tensor sketching operator with the best results. In terms of CPU time, the DCT Gaussian-Sketch method is the fastest. In particular, as shown in Figures <ref>–<ref>, the accuracy of the DCT Gaussian-Sketch method is better than the rt-SVD method and T-Sketch method for input tensors with PolyDecayMed decay spectrum and PolyDecayFast decay spectrum. In this case, with less storage and manipulation, the DCT Gaussian-Sketch method achieves better accuracy for the low-rank approximation. The DCT Gaussian-Sketch method is second only to the truncated t-SVD method in terms of accuracy but is the fastest and far faster than the truncated t-SVD method.








Figures <ref>–<ref> demonstrate the superiority
of the proposed sketching algorithms, compared with the DCT Gaussian-Sketch (q=0) method and rt-SVD method, the DCT Gaussian-Sketch PI (q=1) method achieves smaller relative error and higher PSNR. In terms of CPU time, the DCT Gaussian-Sketch method is the fastest. In particular, as shown in Figures <ref>–<ref>, the accuracy of the DCT Gaussian-Sketch method is better than the rt-SVD method and T-Sketch method for input tensors with PolyDecayMed decay spectrum and PolyDecayFast decay spectrum. In this case, with less storage and manipulation, the DCT Gaussian-Sketch method achieves better accuracy for the low-rank approximation. The DCT Gaussian-Sketch PI (q=1) method is far superior to the T-Sketch PI (q=1) method in terms of time and accuracy. The DCT Gaussian-Sketch PI (q=1) method is second only to the truncated t-SVD method in terms of accuracy but is much faster in speed.





 §.§ Real-world Data

We now conduct experiments on real-world data, i.e., color images and grayscale videos.	




  §.§.§ Color Images

Using the new two-sided sketching algorithm based on the DCT transformed domain low-rank approximation, we compress three color images, i.e., HDU picture[https://www.hdu.edu.cn/landscape], with size of 1200×1800×3, LONDON picture with size of 4775×7155×3, and HDU 2D code with size of 800×800×3. As shown in Figures <ref>, <ref> and <ref>, tensor 𝒜 has slow decaying spectrum (q=0), while it has faster decaying spectrum for 𝒵=(𝒜⋆_L𝒜^⊤)^q ⋆_L𝒜 with q=1.






























From Figures <ref>–<ref>, Figures <ref>–<ref> and Figures <ref>–<ref>, we can see that the new two-sided sketching algorithm based on the DCT transformed domain achieves a better low-rank approximation as the sketch parameter size k increases. In terms of CPU time, the DCT Gaussian-Sketch (q=0) method is the fastest. In particular, as shown in Figures <ref> and <ref>, the accuracy of the DCT Gaussian-Sketch method is better than the rt-SVD method and T-Sketch method. In this case, with less storage and manipulation, the DCT Gaussian-Sketch method achieves better accuracy for the low-rank approximation. The DCT Gaussian-Sketch PI (q=1) method is far superior to the T-Sketch PI (q=1) method in terms of CPU time and accuracy. The DCT Gaussian-Sketch PI (q=1) method is second only to the truncated t_SVD method in terms of accuracy but is much faster in speed.




  §.§.§ Grayscale Video
	
We finally evaluate the proposed sketching algorithms on the widely used YUV Video Sequences[http://trace.eas.asu.edu/yuv/index.html]. Taking the `hall monitor' video as an example and using the first 30 frames, a three order tensor with size of 144×176×30 is then formed for this test. As shown in Figures <ref>–<ref>, similar performance is obtained. In terms of CPU time, the DCT Gaussian-Sketch (q=0) method is the fastest. The DCT Gaussian-Sketch PI (q=1) method is far superior to the T-Sketch PI (q=1) method in terms of CUP time and accuracy. Moreover, the DCT Gaussian-Sketch PI (q=1) method is second only to the truncated t-SVD method in terms of accuracy but is much faster in speed.











§ CONCLUSION


In this paper, we introduced a low tubal rank tensor approximation model based on ⋆_L-product factorization. It is an extension of the new two-sided matrix sketching algorithm <cit.>. We proposed a novel two-sided sketching algorithm in the transformed domain, called DCT Gaussian-sketch algorithm,  combining the DCT transformation with Gaussian sketching operators. 
We also presented the DCT Gaussian-sketch algorithm with power iteration, which can achieve more accurate approximate solutions. A rigorous theoretical analysis of the approximation error of the proposed DCT Gaussian-Sketch algorithm was also provided. Thorough experiments on low-rank approximation of synthetic data, color images and grayscale videos demonstrate the effectiveness and efficiency of the proposed algorithm in terms of computation time and low-rank approximation effects.







9
T.-X17 T.-X. Jiang, T.-Z. Huang, X.-L. Zhao, L.-J. Deng, and Y. Wang, “A
novel tensor-based video rain streaks removal approach via utilizing
discriminatively intrinsic priors,” in Proceedings of the 2017 IEEE
Conference on Computer Vision and Pattern Recognition. IEEE, 2017,
pp. 4057–4066.

T.-K07 T.-K. Kim, S.-F. Wong, and R. Cipolla, “Tensor canonical correlation
analysis for action classification,” in Proceedings of the 2007 IEEE
Conference on Computer Vision and Pattern Recognition. IEEE, 2007,
pp. 1–8.

Du16 B. Du, M. Zhang, L. Zhang, R. Hu, and D. Tao, “Pltd: Patch-based low-rank tensor decomposition for hyperspectral images,” IEEE Transactions
on Multimedia, vol. 19, no. 1, pp. 67–79, 2016.

Renard08 N. Renard, S. Bourennane, and J. Blanc-Talon, “Denoising and dimensionality reduction using multilinear tools for hyperspectral images,”
IEEE Geoscience and Remote Sensing Letters, vol. 5, no. 2, pp. 138–142, 2008.

Wu22 Wu W H, Huang T Z, Zhao X L, et al. Hyperspectral Image Denoising via Tensor Low-Rank Prior and Unsupervised Deep Spatial–Spectral Prior[J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 1-14.

Wang22 Wang J L, Huang T Z, Zhao X L, et al. CoNoT: Coupled Nonlinear Transform-Based Low-Rank Tensor Representation for Multidimensional Image Completion[J]. IEEE Transactions on Neural Networks and Learning Systems, 2022.

Lathauwer97 L. De Lathauwer, “Signal processing based on multilinear algebra,”
Ph.D. dissertation, Katholieke Universiteit Leuven Leuven, 1997.

Comon02 P. Comon, “Tensor decompositions,” Mathematics in Signal Processing
V, pp. 1–24, 2002.

Papalexakis14 E. Papalexakis, K. Pelechrinis, and C. Faloutsos, “Spotting misbehaviors
in location-based social networks using tensors,” in Proceedings of the
2014 23rd International Conference on World Wide Web, 2014, pp. 551–552.

Nakatsuji17 M. Nakatsuji, Q. Zhang, X. Lu, B. Makni, and J. A. Hendler, “Semantic
social network analysis by cross-domain tensor factorization,” IEEE
Transactions on Computational Social Systems, vol. 4, no. 4, pp. 207–217, 2017.

Signoretto14 M. Signoretto, Q. T. Dinh, L. De Lathauwer, and J. A. Suykens. Learning with tensors: a framework based on convex optimization and spectral regularization. Mach. Learn., 94(3):303–351, 2014.

T19 X.-T. Li, X.-L. Zhao, T.-X. Jiang, Y.-B. Zheng, T.-Y. Ji, and T.-Z. Huang. Low-rank tensor
completion via combined non-local self-similarity and low-rank regularization. Neurocomputing,
367:1–12, 2019.

Liu13 J. Liu, P. Musialski, P. Wonka, and J. Ye. Tensor completion for estimating missing values in
visual data. IEEE Trans. Pattern Anal. Mach. Intell., 35(1):208–220, 2013.

Zhang19 X. Zhang. A nonconvex relaxation approach to low-rank tensor completion. IEEE Trans. Neural
Netw. Learn. Syst., 30(6):1659–1671, 2019.

Liu22 Liu Y, Liu J, Long Z, et al. Tensor Sketch[M]//Tensor Computation for Data Analysis. Springer, Cham, 2022: 299-321.

Cao23 Cao X, Zhang X, Zhu C, et al. TS-RTPM-Net: Data-driven Tensor Sketching for Efficient CP Decomposition[J]. IEEE Transactions on Big Data, 2023.

Li17 Li X, Haupt J, Woodruff D. Near optimal sketching of low-rank tensor regression. Advances in Neural Information Processing Systems, 2017, 30.

WangY15 Wang Y, Tung H Y, Smola A J, et al. Fast and guaranteed tensor decomposition via sketching. Advances in neural information processing systems, 2015, 28.

Che19 Che M, Wei Y. Randomized algorithms for the approximations of Tucker and the tensor train decompositions. Adv. in Comput.Math, 45(2019), pp.395-428.

Che21 Che M , Wei Y , Yan H . Randomized algorithms for the low multilinear rank approximations of tensors. Journal of Computational and Applied Mathematics, 2021, 390(2):113380.

Che21J Che M , Wei Y , Yan H . An Efficient Randomized Algorithm for Computing the Approximate Tucker Decomposition. Journal of Scientific Computing, 2021.

Malik18 Malik O A, Becker S. Low-rank tucker decomposition of large tensors using tensorsketch[J]. Advances in neural information processing systems, 2018, 31.

Ravishankar21 Ravishankar J, Sharma M, Khaidem S. A Novel Compression Scheme Based on Hybrid Tucker-Vector Quantization Via Tensor Sketching for Dynamic Light Fields Acquired Through Coded Aperture Camera[C]//2021 International Conference on 3D Immersion (IC3D). IEEE, 2021: 1-8.

Sun20 Sun Y, Guo Y, Luo C, et al. Low-rank tucker approximation of a tensor from streaming data. SIAM Journal on Mathematics of Data Science, 2020, 2(4): 1123-1150.

Hur22 Hur Y, Hoskins J G, Lindsey M, et al. Generative modeling via tensor train sketching. arXiv preprint arXiv:2202.11788, 2022.

Qi21 L. Qi and G. Yu, T-singular values and T-Sketching for third order tensors, arXiv:2013.00976, 2021.

Minster20 Minster R , Saibaba A K , Kilmer M E . Randomized Algorithms for Low-Rank Tensor Decompositions in the Tucker Format. SIAM Journal on Mathematics of Data Science, 2020, 2(1):189-215.

Dong23Dong W, Yu G, Qi L, et al. Practical Sketching Algorithms for Low-Rank Tucker Approximation of Large Tensors[J]. Journal of Scientific Computing, 2023, 95(2): 52.

Diao19 Diao, H., Jayaram, R., Song, Z., Sun, W., & Woodruff, D. (2019). Optimal sketching for kronecker product regression and low rank approximation. Advances in neural information processing systems, 32.

Han20 Han, I., Avron, H., & Shin, J. (2020, November). Polynomial tensor sketch for element-wise function of low-rank matrix. In International Conference on Machine Learning (pp. 3984-3993). PMLR.

Kasiviswanathan18 Kasiviswanathan, S. P., Narodytska, N., & Jin, H. (2018, July). Network Approximation using Tensor Sketching. In IJCAI (pp. 2319-2325).

Ma20 Ma, J., Zhang, Q., Ho, J. C., & Xiong, L. (2020, September). Spatio-temporal tensor sketching via adaptive sampling. In Joint European Conference on Machine Learning and Knowledge Discovery in Databases (pp. 490-506). Springer, Cham.

Shi19 Shi Y, Anandkumar A. Higher-order count sketch: Dimensionality reduction that retains efficient tensor operations[J]. arXiv preprint arXiv:1901.11261, 2019.

Yang18 Yang B, Zamzam A, Sidiropoulos N D. Parasketch: Parallel tensor factorization via sketching[C]//Proceedings of the 2018 SIAM International Conference on Data Mining. Society for Industrial and Applied Mathematics, 2018: 396-404.

Carmon19 Carmon Y, Duchi J C, Aaron S, et al. A rank-1 sketch for matrix multiplicative weights[C]//Conference on Learning Theory. PMLR, 2019: 589-623.

David14 David P. Woodruff. Sketching as a tool for numerical linear algebra. Foundations and Trends in Theoretical
Computer Science, 10(12):1-157, 2014. ISSN 1551-305X. doi: 10.1561/0400000060.

Tropp17 J.A. Tropp, A. Yurtsever, M. Udell and V. Cevher, Practical sketching algorithms for low-rank matrix approximation, SIAM Journal on Matrix Analysis
and Applications 38 (2017) 1454-1485.

Tropp19 J. A. Tropp, A. Yurtsever, M. Udell, and V. Cevher, Streaming low-rank matrix approximation
with an application to scientific simulation, SIAM J. Sci. Comput., 41 (2019), pp. A2430–A2463,
https://doi.org/10.1137/18M1201068.

Kressner22 Kressner D, Vandereycken B, Voorhaar R. Streaming tensor train approximation[J]. arXiv preprint arXiv:2208.02600, 2022.

Che2022 Che M ,  Wei Y ,  Xu Y . Randomized algorithms for the computation of multilinear rank- (μ _1,μ _2,μ _3) approximations. Journal of Global Optimization, 2022:1-31.

Kernfeld15 Kernfeld E, Kilmer M, Aeron S. Tensor-tensor products with invertible linear transforms. Linear Algebra and its Applications, 2015, 485: 545-570.

Martucci94 S. A. Martucci, Symmetric convolution and the discrete sine and cosine transforms,
IEEE Transactions on Signal Processing 42 (5) (1994) 1038–1051.

V.95 V. Sánchez, P. García, A. M. Peinado, J. C. Segura, A. J. Rubio, Diagonalizing
properties of the discrete cosine transforms, IEEE Transactions on Signal Processing
43 (11) (1995) 2631–2641.

Kailath05 T. Kailath, V. Olshevsky, Displacement structure approach to discrete-trigonometric-transform based preconditioners of G. Strang type and of T. Chan type, SIAM Journal
on Matrix Analysis and Applications 26 (3) (2005) 706–734.

Strang99 G. Strang, The discrete cosine transform, SIAM Review 41 (1) (1999) 135–147.

Song20 Song G, Ng M K, Zhang X. Robust tensor completion using transformed tensor singular value decomposition[J]. Numerical Linear Algebra with Applications, 2020, 27(3): e2299.

Kilmer11  
M. Kilmer and C.D. Martin, Factorization strategies for third-order tensors, Linear Algebra and Its Applications 435 (2011) 641-658.

Zhang18 J. Zhang, A.K. Saibaba, M.E. Kilmer and S. Aeron, “A randomized tensor singular value decomposition based on the t-product”, Numerical Linear Algebra with Applications 25 (2018) e2179.

Kilmer08 M. Kilmer, C.D. Martin and L. Perrone, A third-order generalization of the matrix svd as a product of third-order tensors, Tech. Report TR-2008-4 Tufts University, Computer Science Department, 2008.

WangS15  Wang S. A practical guide to randomized matrix computations with MATLAB implementations. arXiv preprint arXiv:1505.07570, 2015.

Halko11 N. Halko, P.G. Martinsson and J.A. Tropp, Find structure with randomness:
Probablistic algorithms for constructing approximate matrix decompositions,
SIAM Review 53 (2011) 217-288.







