# Adaptive Cross Tubal Tensor Approximation

Salman Ahmadi-Asl  Anh Huy Phan  Andrzej Cichocki  Anastasia Sozykina  Zaher
Al Aghbari  Jun Wang  and Ivan Oseledets

###### Abstract

In this paper, we propose a new adaptive cross algorithm for computing a low
tubal rank approximation of third-order tensors, with less memory and lower
computational complexity than the truncated tensor SVD (t-SVD). This makes it
applicable for decomposing large-scale tensors. We conduct numerical
experiments on synthetic and real-world datasets to confirm the efficiency and
feasibility of the proposed algorithm. The simulation results show more than
one order of magnitude acceleration in the computation of low tubal rank
(t-SVD) for large-scale tensors. An application to pedestrian attribute
recognition is also presented.

###### keywords:

Cross tensor approximation, tensor SVD, tubal product

###### MSC:

15A69 , 46N40 , 15A23

††journal: Journal

[label1]organization=Center for Artificial Intelligence Technology, Skolkovo
Institute of Science and Technology, Moscow, Russia<EMAIL_ADDRESS>[label2]organization=Systems Research Institute of Polish Academy of Science,
Warsaw, Poland [label3]organization=Department of Computer Science, University
of Sharjah, Sharjah, 27272, UAE

##  1 Introduction

Tensors are high-dimensional generalizations of matrices and vectors. Contrary
to the rank of matrices, the rank of tensors is not well understood and has to
be defined and determined. Different types of tensor decompositions can have
different rank definitions such as Tensor Train (TT) [1], Tucker decomposition
[2] and its special case, i.e. Higher Order SVD (HOSVD) [3], CANDECOMP/PARAFAC
decomposition (CPD) [4, 5], Block Term decomposition [6], Tensor Train/Tensor
Ring (TT-TR) decomposition [1, 7, 8], tubal SVD (t-SVD) [9]. The t-SVD
factorizes a tensor into three tensors, two orthogonal tensors and one
f-diagonal tensor (to be discussed in Section 3). Like the SVD for matrices,
the truncation version of the t-SVD provides the best tubal rank approximation
for every unitary invariant tensor norm. The t-SVD has been successfully
applied in deep learning [10, 11], tensor completion [12, 13], image
reconstruction [14] and tensor compression [15].

Decomposing big data tensors into the t-SVD format is a challenging task,
especially when the data is extremely massive and we can not view the entire
data tensor. The cross, skeleton, or CUR approximation is a useful paradigm
widely used for fast low-rank matrix approximation. Achieving a higher
compression ratio, and problems with data interpretation are other motivations
to use the cross approximation methods. The main feature of the cross
algorithms that makes them effective for managing very large-scale data
tensors is their ability to use less memory and have lower computational
complexity. When it comes to the higher compression capacity, for instance,
the cross matrix approximation provides sparse factor matrices, whereas the
SVD of sparse matrices fails to do so, resulting in a more compact data
structure. It is also known that the cross approximations can provide more
interpretable approximations, we refer to [16] for more details.

Due to the mentioned motivations, the cross matrix approximation [17, 18] has
been generalized to different types of tensor decompositions such as the TT-
Cross [19], Cross-3D [20], FSTD [21], and tubal Cross [22]. The cross matrix
approximation is generalized to the tensor case based on the tubal product
(t-product) in [22] where some individual lateral and horizontal slices are
selected and based on them a low tubal rank approximation is computed. The
main drawback of this approach is its dependency on the tubal rank estimation,
which may be a difficult task in real-world applications. To tackle this
problem, we propose to generalize the adaptive cross matrix approximation [23,
24, 25] to tensors based on the t-product. The idea is to select one actual
lateral slice and one actual horizontal slice at each iteration and adaptively
check the tubal rank of the tensor.

The generalization of the adaptive cross matrix approximation to tensors based
on the t-product is an interesting problem and in this paper, we discuss how
to perform it properly. The novelties done in this work include:

  * 1.

A new adaptive tubal tensor approximation algorithm, which estimates the tubal
rank and compute the low tubal rank approximation. The proposed algorithm does
not need to use the whole data tensor and at each iteration works only on a
part of the horizontal and lateral slices. This facilitates handling large-
scale tensors.

  * 2.

Presenting an application in pedestrian attribute recognition.

The rest of the paper is structured as follows. The basic definitions are
given in Section 2. The t-SVD model is introduced in Section 3. The cross
matrix approximation and its adaptive version are discussed in Section 4.
Section 5, shows how to generalize the adaptive cross approximation to the
tensor case based on the t-product. We compare the computational complexity of
the algorithms in Section 6. The experimental results are presented in Section
7 and Section 8 concludes the paper and presents potential future directions.

##  2 Preliminaries

The key notations and concepts used in the rest of the paper are introduced in
this section. A tensor, a matrix and a vector are denoted by an underlined
bold capital case letter, a bold capital case letter and a bold lower case
letter, respectively. Slices are subtensors generated with fixed all but two
modes. Our work is for real-valued third-order tensors but generalization to
complex higher order tensor is also straightforward. For a third-order tensor,
$\underline{{\bf X}},$ the three types of slices $\underline{{\bf
X}}(:,:,k),\,\underline{{\bf X}}(:,j,:),\,\underline{{\bf X}}(i,:,:)$ are
called frontal, lateral and horizontal slices. For a third-order tensor
$\underline{{\bf X}},$ the three types fibers $\underline{{\bf
X}}(:j,k),\,\underline{{\bf X}}(i,:,k),\,\underline{{\bf X}}(i,j,:)$ are
called columns, rows and tubes. The notation “${\rm conj}$” means the complex
conjugate of all elements (complex numbers) of a matrix. The notations ${\bf
X}_{(:,-j)}$ and ${\bf X}_{(-i,:)}$ are used to denote new sub-matrices of the
matrix ${\bf X}$ with the $j$-column and the $i$-th row removed. The Frobenius
norm of tensors/matrices is denoted by $\|.\|_{F}$ and the Euclidean norm of a
vector is shown by $\|.\|_{2}$. The notation $|.|$ stands for the absolute
value of a real number. We need the subsequent definitions to introduce the
tensor SVD (t-SVD) model.

######  Definition 1.

(t-product) Let $\underline{\mathbf{X}}\in\mathbb{R}^{I_{1}\times I_{2}\times
I_{3}}$ and $\underline{\mathbf{Y}}\in\mathbb{R}^{I_{2}\times I_{4}\times
I_{3}}$, the t-product
$\underline{\mathbf{X}}*\underline{\mathbf{Y}}\in\mathbb{R}^{I_{1}\times
I_{4}\times I_{3}}$ is defined as follows

$\underline{\mathbf{C}}=\underline{\mathbf{X}}*\underline{\mathbf{Y}}={\rm
fold}\left({{\rm circ}\left(\underline{\mathbf{X}}\right){\rm
unfold}\left(\underline{\mathbf{Y}}\right)}\right),$ (1)

where

${\rm
circ}\left(\underline{\mathbf{X}}\right)=\begin{bmatrix}{\mathbf{X}}^{(1)}&{\mathbf{X}}^{(I_{3})}&\cdots&{\mathbf{X}}^{(2)}\\\
{\mathbf{X}}^{(2)}&{\mathbf{X}}^{(1)}&\cdots&{\mathbf{X}}^{(3)}\\\
\vdots&\vdots&\ddots&\vdots\\\
{\mathbf{X}}^{(I_{3})}&{\mathbf{X}}^{(I_{3}-1)}&\cdots&{\mathbf{X}}^{(1)}\end{bmatrix},$

and

${\rm unfold}(\underline{\mathbf{Y}})=\begin{bmatrix}{\mathbf{Y}}^{(1)}\\\
{\mathbf{Y}}^{(2)}\\\ \vdots\\\ {\mathbf{Y}}^{(I_{3})}\end{bmatrix},\hskip
14.22636pt\underline{\mathbf{Y}}={\rm fold}\left({\rm
unfold}\left(\underline{\mathbf{Y}}\right)\right).$

Here, ${\bf X}^{(i)}=\underline{\bf X}(:,:,i)$ and ${\bf
Y}^{(i)}=\underline{\bf Y}(:,:,i)$ for $i=1,2,\ldots,I_{3}$.

We denote by $\widehat{\underline{\bf X}}$, the Fourier transform of
$\underline{\bf X}$ along its third mode, which can be computed as
$\widehat{\underline{\bf X}}={\rm fft}(\underline{\bf X},[],3)$. It is known
that the block circulant matrix, ${\rm circ}(\underline{\bf
X})\in\mathbb{R}^{I_{1}I_{3}\times I_{2}I_{3}}$, can be block diagonalized,
i.e.,

$\displaystyle({\bf F}_{I_{3}}\otimes{\bf I}_{I_{1}})\,{\rm
circ}\,(\underline{\bf X})({\bf F}^{-1}_{I_{3}}\otimes{\bf
I}_{I_{2}})=\widehat{{\bf X}},$ (2)

where ${\bf F}_{I_{3}}\in\mathbb{R}^{I_{3}\times I_{3}}$ is the discrete
Fourier transform matrix and $({\bf F}_{I_{3}}\otimes{\bf
I}_{I_{1}})/\sqrt{I_{3}}$ is a unitary matrix. Here, the block diagonal matrix
$\widehat{{\bf X}}$ is

$\displaystyle\widehat{{\bf X}}=\begin{bmatrix}\widehat{\underline{\bf
X}}(:,:,1)&&&\\\ &\widehat{\underline{\bf X}}(:,:,2)&&\\\ &&\ddots&\\\
&&&\widehat{\underline{\bf X}}(:,:,I_{3})\\\ \end{bmatrix},$ (3)

and we have the following important properties [26, 27]

$\displaystyle\widehat{\underline{\bf X}}(:,:,1)$ $\displaystyle\in$
$\displaystyle\mathbb{R}^{I_{1}\times I_{2}},$ (4) $\displaystyle{\rm
conj}(\widehat{\underline{\bf X}}(:,:,i))$ $\displaystyle=$
$\displaystyle\widehat{\underline{\bf X}}(:,:,I_{3}-i+2),$ (5)

for $i=2,\ldots,\lceil\frac{I_{3}+1}{2}\rceil+1$. The t-product can be
equivalently performed in the Fourier domain. Indeed, let $\underline{{\bf
C}}=\underline{\bf X}*\underline{\bf Y}$, then from the definition of the
t-product and the fact that the block circulant matrix can be block
diagonalized, we have

$\displaystyle{\rm unfold}(\underline{{\bf C}})$ $\displaystyle=$
$\displaystyle{\rm circ}(\underline{\bf X})\,{\rm unfold(\underline{\bf Y})}$
$\displaystyle=$ $\displaystyle({\bf F}_{I_{3}}^{-1}\otimes{\bf
I}_{I_{1}})(({\bf F}_{I_{3}}\otimes{\bf I}_{I_{1}})\,{\rm
circ}\,(\underline{\bf X})({\bf F}^{-1}_{I_{3}}\otimes{\bf I}_{I_{2}}))$
$\displaystyle(({\bf F}^{-1}_{I_{3}}\otimes{\bf I}_{I_{2}})\,{\rm
unfold}(\underline{\bf Y}))$ $\displaystyle=$ $\displaystyle({\bf
F}_{I_{3}}\otimes{\bf I}_{I_{1}})\,\widehat{{\bf X}}\,{\rm
unfold}(\widehat{\underline{\bf Y}}),$

where $\widehat{\underline{\bf Y}}={\rm fft}(\underline{\bf Y},[],3)$. If we
multiply both sides of (2) from the left-hand side with $({\bf
F}_{I_{3}}\otimes{\bf I}_{I_{1}})$, we get ${\rm
unfold}(\widehat{\underline{\bf C}})=\widehat{{\bf X}}\,{\rm
unfold}(\widehat{\underline{\bf Y}})$, where $\widehat{\underline{\bf C}}={\rm
fft}(\underline{\bf C},[],3)$. This means that
$\widehat{\underline{\mathbf{C}}}\left({:,:,i}\right)=\widehat{\underline{\mathbf{X}}}\left({:,:,i}\right)\,\widehat{\underline{\mathbf{Y}}}\left({:,:,i}\right)$.
So, it suffices to transform two given tensors into the Fourier domain and
multiply their frontal slices. Then, the resulting tensor in the Fourier
domain returned back to the original space via the inverse FFT. Note that due
to the equations in (4)-(5), half of the computations are reduced. This
procedure is summarized in Algorithm 1.

######  Definition 2.

(Transpose) The transpose of a tensor
$\underline{\mathbf{X}}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$ is
denoted by $\underline{\mathbf{X}}^{T}\in\mathbb{R}^{I_{2}\times I_{1}\times
I_{3}}$ produced by applying the transpose to all frontal slices of the tensor
$\underline{\mathbf{X}}$ and reversing the order of the second untill the last
transposed frontal slices.

######  Definition 3.

(Identity tensor) Identity tensor
$\underline{\mathbf{I}}\in\mathbb{R}^{I_{1}\times I_{1}\times I_{3}}$ is a
tensor whose first frontal slice is an identity matrix of size $I_{1}\times
I_{1}$ and all other frontal slices are zero. It is easy to show
$\underline{{\bf I}}*\underline{{\bf X}}=\underline{{\bf X}}$ and
$\underline{{\bf X}}*\underline{{\bf I}}=\underline{{\bf X}}$ for all tensors
of conforming sizes.

######  Definition 4.

(Orthogonal tensor) A tensor $\underline{\mathbf{X}}\in\mathbb{R}^{I_{1}\times
I_{1}\times I_{3}}$ is orthogonal (under t-product operator) if
${\underline{\mathbf{X}}^{T}}*\underline{\mathbf{X}}=\underline{\mathbf{X}}*{\underline{\mathbf{X}}^{T}}=\underline{\mathbf{I}}$.

######  Definition 5.

(f-diagonal tensor) If all frontal slices of a tensor are diagonal then the
tensor is called an f-diagonal tensor.

######  Definition 6.

(Inverse of a tensor) The inverse of a tensor $\underline{\bf
X}\in\mathbb{R}^{I_{1}\times I_{1}\times I_{3}}$ is denoted by $\underline{\bf
X}^{-1}\in\mathbb{R}^{I_{1}\times I_{1}\times I_{3}}$ is a unique tensor
satisfying $\underline{\bf X}*\underline{\bf X}^{-1}=\underline{\bf
X}^{-1}*\underline{\bf X}=\underline{\bf I},$ where $\underline{\bf
I}\in\mathbb{R}^{I_{1}\times I_{1}\times I_{3}}$ is the identity tensor. The
inverse of a tensor can also be computed in the Fourier domain and described
in Algorithm 2. The MATLAB command “inv” in Line 3 computes the inverse of a
matrix. The Moore–Penrose (MP) inverse of a tensor $\underline{\bf
X}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$ is denoted by $\underline{\bf
X}^{{\dagger}}\in\mathbb{R}^{I_{2}\times I_{1}\times I_{3}}$ and can be
computed by Algorithm 2 where “inv” is replaced with the MATLAB function
“pinv”. Here, “pinv” stands for the MP inverse of a matrix.

Input : Two data tensors
$\underline{\mathbf{X}}\in{\mathbb{R}^{{I_{1}}\times{I_{2}}\times{I_{3}}}},\,\,\underline{\mathbf{Y}}\in{\mathbb{R}^{{I_{2}}\times{I_{4}}\times{I_{3}}}}$

Output : t-product
$\underline{\mathbf{C}}=\underline{\mathbf{X}}*\underline{\mathbf{Y}}\in\mathbb{R}^{I_{1}\times
I_{4}\times I_{3}}$

1 $\widehat{\underline{\mathbf{X}}}={\rm
fft}\left({\underline{\mathbf{X}},[],3}\right)$;

2 $\widehat{\underline{\mathbf{Y}}}={\rm
fft}\left({\underline{\mathbf{Y}},[],3}\right)$;

3 for  _$i=1,2,\ldots,\lceil\frac{I_{3}+1}{2}\rceil$_ do

4
$\widehat{\underline{\mathbf{C}}}\left({:,:,i}\right)=\widehat{\underline{\mathbf{X}}}\left({:,:,i}\right)\,\widehat{\underline{\mathbf{Y}}}\left({:,:,i}\right)$;

5

6 end for

7for  _$i=\lceil\frac{I_{3}+1}{2}\rceil+1,\ldots,I_{3}$_ do

8 $\widehat{\underline{\mathbf{C}}}\left({:,:,i}\right)={\rm
conj}(\widehat{\underline{\mathbf{C}}}\left({:,:,I_{3}-i+2}\right))$;

9 end for

$\underline{\mathbf{C}}={\rm
ifft}\left({\widehat{\underline{\mathbf{C}}},[],3}\right)$;

Algorithm 1 Fast t-product of two tensors [9, 27]

Input : The data tensor
$\underline{\mathbf{X}}\in{\mathbb{R}^{{I_{1}}\times{I_{1}}\times{I_{3}}}}$

Output : Tensor inverse $\underline{\mathbf{X}}^{-1}\in\mathbb{R}^{I_{1}\times
I_{1}\times I_{3}}$

1 $\widehat{\underline{\mathbf{X}}}={\rm
fft}\left({\underline{\mathbf{X}},[],3}\right)$;

2 for  _$i=1,2,\ldots,\lceil\frac{I_{3}+1}{2}\rceil$_ do

3 $\widehat{\underline{\mathbf{C}}}\left({:,:,i}\right)={\rm
inv}\,(\widehat{{\bf X}}(:,:,i))$;

4

5 end for

6for  _$i=\lceil\frac{I_{3}+1}{2}\rceil+1,\ldots,I_{3}$_ do

7 $\widehat{\underline{\mathbf{C}}}\left({:,:,i}\right)={\rm
conj}(\widehat{\underline{\mathbf{C}}}\left({:,:,I_{3}-i+2}\right))$;

8 end for

$\underline{\mathbf{X}}^{-1}={\rm
ifft}\left({\widehat{\underline{\mathbf{C}}},[],3}\right)$;

Algorithm 2 Fast inverse computation of the tensor $\underline{\bf X}$

##  3 Tensor SVD (t-SVD)

The tensor SVD (t-SVD) represents a tensor as the t-product of three tensors.
The first and last tensors are orthogonal, while the middle tensor is an
f-diagonal tensor. To be more precise, let $\underline{\bf
X}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}$, then the t-SVD of the tensor
$\underline{\bf X},$ is $\underline{\bf X}=\underline{\bf U}*\underline{\bf
S}*{\underline{\bf V}^{T}},$ where $\underline{\bf
U}\in\mathbb{R}^{I_{1}\times R\times I_{3}}$, and $\underline{\bf
V}\in\mathbb{R}^{I_{2}\times R\times I_{3}}$ are orthogonal tensors and the
tensor $\underline{\bf S}\in\mathbb{R}^{R\times R\times I_{3}}$ is f-diagonal
[9, 28], see Figure 1 for an illustration on the t-SVD and its truncated
version. Note that Algorithm 3 only needs the truncated SVD of the
$\lceil\frac{I_{3}+1}{2}\rceil$ first frontal slices. The generalization of
the t-SVD to tensors of order higher than three is done in [29].

Figure 1: a) Tensor SVD (t-SVD) of the tensor $\underline{\bf X}$, b) The
truncated t-SVD of the tensor $\underline{\bf X}$ for the tubal rank $R$ [30].

Other types of classical matrix decompositions such as QR and LU
decompositions can be generalized based on the t-product in straightforward
ways.

The computational complexity of Algorithm 3 is dominated by the FFT of all
tubes of an input tensor and also the truncated SVD of the frontal slices in
the Fourier domain. In the literature, some algorithms have been developed to
accelerate these computations. For example, using the idea of randomization,
we can replace the classical truncated SVD with more efficient and faster
approaches such as the randomized SVD [31, 16] or cross matrix approximation.
Although this idea can somehow solve the mentioned computation difficulty of
Algorithm 3, still we need to access all elements of the underlying data
tensor. For very big data tensors where viewing the data tensor even once is
very prohibitive, it is required to develop algorithms that only use a part of
the data tensor at each iteration. In this paper, we follow this idea and
propose an efficient algorithm for the computation of the t-SVD, which uses
only a part of lateral and horizontal slices of a tensor at each iteration.
This significantly accelerates the computations and in some of our
simulations, we have achieved almost two orders of magnitude acceleration,
which shows the performance of the proposed algorithm. To the best of our
knowledge, this is the first adaptive cross algorithm developed for the
computation of the t-SVD.

Figure 2: a) One stage of the adaptive cross matrix approximation for rank
reduction. The corresponding column and row of the residual matrix ${\bf Z},$
with the same indices as the selected column and row of the original data
matrix ${\bf X}$ become zeros, i.e. ${\rm rank}({\bf X}-{\bf Y})={\rm
rank}({\bf X})-1.$ The rank one matrix approximation ${\bf Y}$ interpolates
${\bf X}$ at the selected column and row. b) One stage of the cross tubal
approximation for the tubal rank reduction. The corresponding lateral and
horizontal slices of the residual tensor $\underline{\bf Z},$ with the same
indices as the selected lateral and horizontal slices of the original data
tensor $\underline{\bf X}$ become zeros, ${\rm rank}(\underline{\bf
X}-\underline{\bf Y})={\rm rank}(\underline{\bf X})-1$. The tubal rank one
tensor approximation $\underline{\bf Y}$ interpolates $\underline{\bf X}$ at
the selected lateral and horizontal slices.

Input : The data tensor
$\underline{\mathbf{X}}\in{\mathbb{R}^{{I_{1}}\times{I_{2}}\times{I_{3}}}}$
and a target tubal rank $R$

Output : The truncated t-SVD of the tensor $\underline{\bf X}$

1 $\widehat{\underline{\mathbf{X}}}={\rm
fft}\left({\underline{\mathbf{X}},[],3}\right)$;

2 for  _$i=1,2,\ldots,\lceil\frac{I_{3}+1}{2}\rceil$_ do

3
$[\widehat{\underline{\mathbf{U}}}\left({:,:,i}\right),\widehat{\underline{\bf
S}}(:,:,i),\widehat{{\bf V}}(:,:,i)]={\rm svds}\,(\widehat{{\bf
X}}(:,:,i),R)$;

4

5 end for

6for  _$i=\lceil\frac{I_{3}+1}{2}\rceil+1,\ldots,I_{3}$_ do

7 $\widehat{\underline{\mathbf{U}}}\left({:,:,i}\right)={\rm
conj}(\widehat{\underline{\mathbf{U}}}\left({:,:,I_{3}-i+2}\right))$;

8
$\widehat{\underline{\mathbf{S}}}\left({:,:,i}\right)=\widehat{\underline{\mathbf{S}}}\left({:,:,I_{3}-i+2}\right)$;

9 $\widehat{\underline{\mathbf{V}}}\left({:,:,i}\right)={\rm
conj}(\widehat{\underline{\mathbf{V}}}\left({:,:,I_{3}-i+2}\right))$;

10 end for

$\underline{\mathbf{U}}_{R}={\rm
ifft}\left({\widehat{\underline{\mathbf{U}}},[],3}\right)$;
$\underline{\mathbf{S}}_{R}={\rm
ifft}\left({\widehat{\underline{\mathbf{S}}},[],3}\right)$;
$\underline{\mathbf{V}}_{R}={\rm
ifft}\left({\widehat{\underline{\mathbf{V}}},[],3}\right)$

Algorithm 3 The truncated t-SVD decomposition of the tensor $\underline{\bf
X}$

##  4 Matrix cross approximation and its adaptive version

The cross matrix approximation was first proposed in [17] for fast low-rank
approximation of matrices. It provides a low-rank matrix approximation based
on some actual columns and rows of the original matrix. It has been shown that
a cross approximation with the maximum volume of the intersection matrix leads
to close to optimal approximation [18]. The adaptive cross approximation or
Cross2D algorithm [32, 33, 23, 24, 25] sequentially selects a column and a row
of the original data matrix and based on them, computes a rank-1 matrix scaled
by the intersection element, as stated in the following theorem, which indeed
is the Gaussian elimination process.

######  Theorem 1.

(Rank-1 deflation) Let ${\bf X}\in\mathbb{R}^{I_{1}\times I_{2}}$ be a given
matrix, and we select the $i$-th row and the $j$-th column with the nonzero
intersection element ${\bf X}(i,j)$. Then the following residual matrix

${\bf Y}={\bf X}-\frac{1}{{\bf X}(i,j)}{\bf X}(:,j){\bf X}(i,:),$

vanishes at the $i$-th row and $j$-th column, so ${\rm rank}({\bf Y})={\rm
rank}({\bf X})-1$.

###### Proof.

It is obvious that the $j$-th column and $i$-th row of ${\bf Y}$ are zeros

$\displaystyle{\bf Y}(:,j)={\bf X}(:,j)-\frac{1}{{\bf X}(i,j)}{\bf X}(:,j){\bf
X}(i,j)=0,$ (7) $\displaystyle{\bf Y}(i,:)={\bf X}(i,:)-\frac{1}{{\bf
X}(i,j)}{\bf X}(i,j){\bf X}(i,:)=0.$ (8)

To prove Theorem 1, we consider two cases:

  1. 1.

If ${\bf X}$ is of full-rank, that is, either ${\bf X}(:,j)\notin\,{\rm
range}\,({\bf X}_{(:,-j)})$ or ${\bf X}(i,:)\notin\,{\rm range}\,({\bf
X}_{(-i,:)})$, then it is straightforward that ${\bf Y}$ has smaller rank than
${\bf X}$.

  2. 2.

Otherwise, we consider the case ${\bf X}$ is rank-deficient and ${\bf X}(:,j)$
and ${\bf X}(i,:)$ are non zero vectors that ${\bf X}(i,:)\in{\rm
range}\,({\bf X}_{(-i,:)})$ and ${\bf X}(:,j)\in{\rm range}\,({\bf
X}_{(:,-j)})$.

Without loss of generality, we assume that ${\bf X}(:,j)$ and ${\bf X}(i,:)$
are the last column and the last row of the matrix ${\bf X}$, i.e.
$i=I_{1},\,j=I_{2}$, (see illustration in Figure 3) and rank(${\bf X}$)
$<\min(I_{1},I_{2})$.

Since ${\bf X}(i,:)\in{\rm range}\,({\bf X}_{(-i,:)})$ and ${\bf
X}(:,j)\in{\rm range}\,({\bf X}_{(:,-j)})$, there exist linear combinations
such that

$\displaystyle{\bf X}(:,j)=\begin{bmatrix}{\bf a}\\\
{c}\end{bmatrix}=\begin{bmatrix}{\bf Z}\\\ {\bf
b}^{T}\end{bmatrix}{\boldsymbol{\alpha}},\quad\quad{\bf
X}(i,:)^{T}=\begin{bmatrix}{\bf b}\\\ {c}\end{bmatrix}=\begin{bmatrix}{\bf
Z}^{T}\\\ {\bf a}^{T}\end{bmatrix}{\boldsymbol{\beta}},$ (9)

where ${\boldsymbol{\alpha}}\neq 0,\,{\boldsymbol{\beta}}\neq 0$. This gives
$c={\boldsymbol{\beta}}^{T}{\bf Z}{\boldsymbol{\alpha}}$.

Figure 3: Partitioning the matrix $\bf X$ for the proof of Theorem 1.

The rank-1 matrix deflation yields

$\displaystyle{\bf Y}={\bf X}-\frac{1}{c}\begin{bmatrix}{\bf a}\\\
c\end{bmatrix}\begin{bmatrix}{\bf b}^{T}&c\end{bmatrix}=\begin{bmatrix}{\bf
Z}-\frac{1}{c}{\bf a}{\bf b}^{T}&{\bf 0}\\\ {\bf 0}&0\end{bmatrix}.$ (10)

where the top-left submatrix of ${\bf Y}$ has rank-1 reduction from $\bf Z$

$\displaystyle{\bf W}={\bf Z}-\frac{1}{c}{\bf a}{\bf b}^{T}={\bf
Z}-\frac{1}{c}{\bf Z}{\boldsymbol{\alpha}}\,{\boldsymbol{\beta}}^{T}{\bf Z}.$
(11)

We next substitute ${\bf Z}$ by its singular value decomposition ${\bf Z}={\bf
U}{\bf S}{\bf V}^{T}$, where $\bf S$ is a diagonal matrix of positive singular
values of $\bf Z$ and consider

$\displaystyle{\bf W}={\bf U}\left(\underbrace{{\bf S}-\frac{1}{c}{\bf S}({\bf
V}^{T}{\bf{\boldsymbol{\alpha}}})({\boldsymbol{\beta}}^{T}\bf U){\bf S}}_{\bf
K}\right){\bf V}^{T}.$ (12)

Assume ${\bf d}={\bf V}^{T}{\boldsymbol{\alpha}}$ and ${\bf e}={\bf
U}^{T}{\boldsymbol{\beta}}$, then from (12), we have

$\displaystyle{\bf K}={\bf S}-\frac{{\bf S}{\bf d}{\bf e}^{T}{\bf S}}{{\bf
e}^{T}{\bf S}{\bf d}}.$ (13)

It is straightforward to see that

$\displaystyle{\bf S}^{-1/2}{\bf K}{\bf S}^{-1/2}={\bf I}-\frac{{\bf u}{\bf
v}^{T}}{{\bf v}^{T}{\bf u}},$ (14)

where ${\bf u}={\bf S}^{1/2}{\bf d}$ and ${\bf v}={\bf S}^{1/2}{\bf e}$. Now,
it is readily seen that ${\bf K}$ has a zero singular value, i.e., its rank is
reduced from the rank of $\bf S$ by 1. Since the multiplication with
orthogonal matrices does not change the matrix rank, the proof of the theorem
is completed.

######  Remark 2.

An alternative proof for Theorem 1 adopts the following fact proved in [34].
Given ${\bf X}\in\mathbb{R}^{m\times n}$ and assume ${\bf
U}\in\mathbb{R}^{m\times k},\,{\bf R}\in\mathbb{R}^{k\times k}$, and ${\bf
V}\in\mathbb{R}^{n\times k}$. Then

$\displaystyle{\rm rank}({\bf X}-{\bf U}{\bf R}^{-1}{\bf V}^{T})={\rm
rank}({\bf X})-{\rm rank}({\bf U}{\bf R}^{-1}{\bf V}^{T}),$ (15)

if and only if there exist ${\bf A}\in\mathbb{R}^{n\times k}$ and ${\bf
B}\in\mathbb{R}^{m\times k}$ such that ${\bf U}={\bf X}{\bf A}$, ${\bf V}={\bf
X}^{T}{\bf B}$, and ${\bf R}={\bf B}^{T}{\bf X}{\bf A}$. If we define ${\bf
A}={\bf e}_{j}\in\mathbb{R}^{n},\,\,{\bf B}={\bf e}_{i}\in\mathbb{R}^{m}$ as
the $j$-th and $i$-th standard unit vectors111For the standard unit vector
$e_{i}$, the $i$th element is 1 and the rest are zero.. Then, we have ${\bf
U}={\bf X}{\bf A}={\bf X}(:,j),\,\,{\bf V}={\bf X}^{T}{\bf B}={\bf
X}(i,:)^{T}$ and ${\bf R}^{-1}=({\bf e}_{i}^{T}{\bf X}{\bf
e}_{j})^{-1}=\frac{1}{{\bf X}(i,j)}$. Now 15, demonstrates that

$\displaystyle{\rm rank}({\bf X}-\frac{1}{c}{\bf X}(:,j){\bf X}(i,:))$
$\displaystyle=$ $\displaystyle{\rm rank}({\bf X})-{\rm rank}(\frac{1}{c}{\bf
X}(:,j){\bf X}(i,:))$ $\displaystyle=$ $\displaystyle{\rm rank}({\bf X})-1.$

So, this completes the proof.

∎

In view of Theorem 1, we see that the corresponding column and row of the
residual matrix ${\bf Y},$ with the same indices as the selected column/row of
the original data matrix ${\bf X}$ become zero, which means that this
approximation interpolates the original data matrix at the mentioned indices
and reduces its rank by one order, see Figure 2 (a) for a graphical
illustration of this approach. This procedure is repeated by selecting a new
column and a new row of the residual matrix, so we can sequentially reduce the
matrix rank and interpolate the original matrix at the new columns/rows. The
adaptive cross matrix approximation method is summarized in Algorithm 4.
Clearly, the breakdown can happen in Algorithm 4 if the denominator ${\bf
u}_{k}(i_{k}),$ becomes zero. If we deal with such a case, we should select a
new index $i_{k}$ for which the ${\bf u}_{k}(i_{k})$ is not zero. For example
one can select randomly a new index, which has not been chosen previously.

Input : A data matrix ${\mathbf{X}}\in{\mathbb{R}^{{I_{1}}\times{I_{2}}}}$, an
approximation error $\epsilon$

Output : Low rank matrix approximation ${\mathbf{X}}={{\bf U}}{{\bf V}}$

1 ${\bf U}=0,\,{\bf V}=0,\,\mu=0,\,r_{0}=0,$ $j_{1}$ is a random column index
for  _$k=1,2,\ldots,{\rm min}(I_{1},I_{2})$_ do

2 ${\bf u}_{k}={\bf E}_{k-1}(:,j_{k})={\bf X}(:,j_{k})-{\bf U}{\bf
V}(:,j_{k})$;

3 $i_{k}={\rm\arg\max}_{i}|{\bf u}_{k}(i)|$;

4 ${\bf u}_{k}\leftarrow{\bf u}_{k}/{\bf u}_{k}(i_{k})$;

5 ${\bf v}_{k}^{T}={\bf E}_{k-1}(i_{k},:)={\bf X}(i_{k},:)-{\bf
U}(i_{k},:){\bf V}$;

6 $j_{k+1}=\arg\max_{j}|{\bf v}_{k}(j)|$;

7 $\rho^{2}=\|{\bf u}_{k}\|_{2}^{2}\|{\bf v}_{k}\|_{2}^{2}$

8 $\mu^{2}\leftarrow\mu^{2}+\rho^{2}+2\sum_{j=1}^{k-1}{\bf V}(j,:){\bf
v}_{k}{\bf u}_{k}^{T}{\bf U}(:,j)$;

9 ${\bf U}\leftarrow[{\bf U},{\bf u}_{k}],\,{\bf V}\leftarrow[{\bf V};{\bf
v}^{T}_{k}];$

10 $\,r_{k}=r_{k-1}+1$;

11 if  _$\rho <\epsilon\mu$_ then

12 Break

13 end if

14

15 end for

Algorithm 4 Adaptive cross approximation algorithm (ACA) [32, 33, 23, 24, 25]

##  5 Proposed adaptive tensor cross approximation based on the t-product

In this section, we show how to generalize the adaptive cross matrix
approximation to the tensor case based on the t-product. Compared to the
matrix case, instead of a column and a row, here we select a lateral slice and
a horizontal slice at each iteration but the important question is how to use
the intersection tube for scaling the corresponding tubal rank-1 tensor so
that the tubal rank of the residual tensor is reduced one in order. We found
that the inverse of the intersection tube should be used and this is proved in
Theorem 3.

######  Theorem 3.

(Tubal rank-1 deflation) Let $\underline{\bf X}\in\mathbb{R}^{I_{1}\times
I_{2}\times I_{3}}$ be a given data tensor with sampled lateral and horizontal
slices as $\underline{\bf X}(:,j,:)$ and $\underline{\bf X}(i,:,:)$ with the
nonzero intersection tube ${\bf X}(i,j,:)$. Then the residual tensor

$\underline{\bf Y}=\underline{\bf X}-\underline{\bf X}(:,j,:)*(\underline{\bf
X}(i,j,:))^{-1}*\underline{\bf X}(i,:,:),$ (16)

vanishes at its $i$-th lateral slice and $j$-th horizontal slice and ${\rm
rank}(\underline{\bf Y})={\rm rank}(\underline{\bf X})-1$.

Input : A data tensor $\underline{\bf
X}\in{\mathbb{R}^{{I_{1}}\times{I_{2}\times I_{3}}}}$, an approximation error
$\epsilon$

Output : Low tubal rank tensor approximation $\underline{\bf
X}={\underline{\bf U}}*{\underline{\bf V}}$

1 $\underline{\bf U}=0,\,\underline{\bf V}=0,\,\mu=0,\,r_{0}=0,$ $j_{1}$ is a
random lateral slice index for  _$k=1,2,\ldots,{\rm min}(I_{1},I_{2})$_ do

2 $\underline{\bf u}_{k}=\underline{\bf E}_{k-1}(:,j_{k},:)=\underline{\bf
X}(:,j_{k},:)-\underline{\bf U}*\underline{\bf V}(:,j_{k},:)$;

3 $i_{k}={\rm\arg\max}_{i}\,\,||\underline{\bf u}_{k}(i,j_{k},:)||_{2}^{2}$;

4 ${\underline{\bf u}}_{k}\leftarrow{\underline{\bf u}}_{k}*({\underline{\bf
u}}_{k}(i_{k},j_{k},:))^{-1}$;

5 ${\underline{\bf v}}_{k}^{T}=\underline{\bf
E}_{k-1}(i_{k},:,:)=\underline{\bf X}(i_{k},:,:)-\underline{\bf
U}(i_{k},:,:)*\underline{\bf V}$;

6 $j_{k+1}=\arg\max_{j}\,\,||{\underline{\bf v}}_{k}(i_{k},j,:)||_{2}^{2}$;

7 $\rho^{2}=\|{\underline{\bf u}}_{k}\|_{F}^{2}\|{\underline{\bf
v}}_{k}\|_{F}^{2}$;

8 $\mu^{2}\leftarrow\mu^{2}+\rho^{2}+2\|\sum_{j=1}^{k-1}{\underline{\bf
v}}(j,:,:)*{\underline{\bf v}}_{k}*{\underline{\bf u}}_{k}^{T}*{\underline{\bf
u}}_{k}(:,j,:)\|_{2}^{2}$;

9 $\underline{\bf U}\leftarrow[\underline{\bf U},{\underline{\bf
u}}_{k}],\,\underline{\bf V}\leftarrow[\underline{\bf V};{\underline{\bf
v}}^{T}_{k}];$

10 $r_{k}=r_{k-1}+1$;

11 if  _$\rho <\epsilon\mu$_ then

12 break

13 end if

14

15 end for

Algorithm 5 Proposed adaptive cross tubal tensor approximation algorithm
(ACTA)

###### Proof.

To prove Theorem 3, we show that the $i$-th row and the $j$-th column of each
frontal slice of the residual tensor $\underline{\bf Y}$ are zero. To do this,
let us consider the $k$-th frontal slice of the residual tensor
$\underline{\bf Y}$ as $\underline{\bf Y}(:,;,k)$. In the Fourier domain, it
can be represented as

$\displaystyle\widehat{\underline{\bf Y}}(:,:,k)=\widehat{\underline{\bf
X}}(:,:,k)-\frac{1}{\widehat{\underline{\bf X}}(i,j,k)}\widehat{\underline{\bf
X}}(:,j,k)\widehat{\underline{\bf X}}(i,:,k).$ (17)

In view of Theorem 1, this means that the $i$-th row and the $j$-th column of
the the $k$-th frontal slice of the tensor $\underline{\bf Y}$ in the Fourier
domain are zero and its rank is one order lower than the rank of the matrix
$\widehat{\bf X}(:,:,k)$. So the $i$-th row and the $j$-th column of the all
frontal slices $\underline{\bf X}(:,:,k),\,k=1,2,\ldots,I_{3}$ equal to zero.
This clearly completes the proof. ∎

It is not difficult to see that for second-order tensors (matrices), Equation
(16) is reduced to the classical matrix cross approximation. At each
iteration, we select a lateral slice and a horizontal slice and perform the
scaling using the pseudoinverse of the intersection tube. The corresponding
scaled tubal rank-1 tensor reduces the tubal rank of the underlying data
tensor by one order. It is interesting to note that similar to the matrix case
where after each iteration the corresponding selected column and row in the
residual matrix become zeros, here the corresponding lateral and horizontal
slices of the residual tensor vanish. So, naturally, this approximation
interpolates the original tensor at the mentioned slices. This procedure can
proceed with the residual tensor to reduce the tubal rank sequentially. The
generalized adaptive cross tubal approximation method is outlined in Algorithm
5. In Line 10 of Algorithm 5, the new horizontal and lateral slicers are
concatenated along the second and first modes, respectively. The relative
error accuracy is used for the stopping criterion as $\rho\leq\epsilon\mu$
according to

$\displaystyle\rho=\|\underline{\bf u}_{k}*\underline{\bf v}^{T}_{k}\|_{F}$
$\displaystyle\approx$ $\displaystyle\|\underline{\bf X}-\underline{\bf
U}*\underline{\bf V}\|_{F},$ $\displaystyle\mu=\|\underline{\bf
U}*\underline{\bf V}\|_{F}$ $\displaystyle\approx$
$\displaystyle\|\underline{\bf X}\|_{F},$

where $\underline{\bf U}=[\underline{\bf u}_{1},\ldots,\underline{\bf
u}_{k-1}]$ and $\underline{\bf V}=[\underline{\bf v}_{1};\ldots;\underline{\bf
v}_{k-1}]$. It is necessary to enforce $i_{k}\neq i_{1},i_{2},\ldots,i_{k-1}$
and $j_{k}\neq j_{1},j_{2},\ldots,j_{k-1},$ which means that each iteration
should produce new indices (different from the others). We also remark that if
the size of a frontal/lateral slice is big, one can compress it using the
classical cross methods, similar to [20] where the cross approximation is used
in two stages for the computation of the Tucker decomposition. Although our
results so far are for third-order tensors, clearly they can be
straightforwardly generalized to tensors of a higher order than three,
according to [29].

##  6 Computational complexity

The adaptive cross tensor algorithm is efficient as it only works on a lateral
slice and a horizontal slice at each iteration. The computational complexity
of Algorithm 5 is $\mathcal{O}((I+J)K\log(K))$. The computational complexity
of the truncated t-SVD for a tensor of the size $I\times J\times K$ is
$\mathcal{O}(IJK{\rm log}(K))+\mathcal{O}(IJK{\rm min}(I,J))$. Besides, the
truncated t-SVD needs to access and process the whole data tensor while the
proposed algorithm works only on a part of the lateral slice and horizontal
slices at each iteration. So, it is clearly seen that the proposed Algorithm 5
requires much less memory and computational operation than the t-SVD
algorithm. This makes it applicable for decomposing large-scale tensors.

##  7 Experimental Results

We have used Matlab and some functions of the toolbox

https://github.com/canyilu/Tensor-tensor-product-toolbox

to implement the proposed algorithm using a laptop computer with 2.60 GHz
Intel(R) Core(TM) i7-5600U processor and 8GB memory. We have used two metrics,
relative error and Peak signal-to-noise ratio (PSNR) to compare the efficiency
of the proposed algorithm with the baselines. The relative error is defined as
follows

${\rm Relative\,\,error}=\frac{\|\underline{\bf X}-\underline{\bf
U}*\underline{\bf S}*\underline{\bf V}^{T}\|_{F}}{\|\underline{\bf X}\|_{F}}.$

The PSNR is also defined as

${\rm{PSNR=10lo}}{{\rm{g}}_{{\rm{10}}}}\left({{{255}^{2}}/{\rm{MSE}}}\right),$

where ${\rm{MSE}}=\left\|{\underline{\bf X}-\widehat{\underline{\bf
X}}}\right\|_{F}^{2}/{\rm{num}}\left(\underline{\bf X}\right).$ Note “num”
denotes the number of parameters of a given data tensor. We mainly consider
three examples. In the first example, we examine the algorithms using low-rank
random data tensors. In the second example, we have used the functional based
tensors. In the last example, we used the images as real-world data with
application to the image completion problem.

######  Example 1.

In this example we consider a random data tensor ${\bf
X}\in\mathbb{R}^{N\times N\times N}$ with exact tubal rank $R=30$ for
$N=100,200,\ldots,600$. To generate such a tensor, we considered two standard
Gaussian tensors and orthonormalize them. Let us denote these orthogonal parts
by $\underline{\bf U}\in\mathbb{R}^{N\times R\times N}$ and $\underline{\bf
V}\in\mathbb{R}^{N\times R\times N}$. Then, we generate a tensor
$\underline{\mathbf{S}}\in\mathbb{R}^{R\times R\times N}$ with only $R$
nonzero diagonal tubes $\underline{\mathbf{S}}(i,i,:),\,i=1,2,\ldots,R$ whose
elements are also standard Gaussian and build the tensor $\underline{\bf
X}=\underline{\bf U}*\underline{\mathbf{S}}*\underline{\bf V}^{T},$ which is
used in our simulations. Assume that $R=30$ in our simulations, and set
$\epsilon=10^{-8}$ in Algorithm 5. Then, we apply the proposed algorithm to
find the tubal rank and the corresponding low tubal rank approximation. We
consider 100 Monte Carlo experiments and report the mean of our results
(accuracy and running time). In all our experiments, the proposed approach
retrieved the true tubal rank successfully, and this convinced us that it
works well for finding the tubal rank of a tensor. Then we used the truncated
t-SVD and the randomized t-SVD [22] to compute low tubal rank approximations
of the underlying data tensor. The running time of the proposed algorithm, the
truncated t-SVD, the randomized t-SVD are compared in Figure 4 (right). The
numerical results show almost two orders of magnitude speed-up of the proposed
approach compared with the truncated t-SVD algorithm, while it is also faster
than the randomized t-SVD. The accuracy comparison of the algorithms is also
presented in Figure 4 (left). This illustrates that the proposed algorithm can
provide acceptable results in much less time than the truncated t-SVD
algorithm and the randomized t-SVD.

Figure 4: The running time and accuracy comparisons of the proposed algorithm
and the truncated t-SVD for Example 1.

Figure 5: The running time comparison of the truncated t-SVD and the proposed
algorithm for Case I (upper left), Case II (upper right) and Case III (bottom)
for Example 2. Figure 6: The original image, the available with 70% missing
pixels (randomly) and the reconstructed images using the truncated t-SVD and
the proposed approach for Example 3.

######  Example 2.

In this example, we apply Algorithm 5 to compute low tubal rank approximations
of function based tensors. To do so, we consider the following case studies:

  * 1.

Case study I: $\quad\underline{\bf
X}(i,j,k)=\frac{1}{(i^{2}+j^{2}+k^{2})^{1/2}};$

  * 2.

Case study II: $\quad\underline{\bf X}(i,j,k)={\sin{(i+j+k)}}+\tanh(i+j+k);$

  * 3.

Case study III: $\quad\underline{\bf
X}(i,j,k)=\frac{1}{(i^{5}+j^{5}+k^{5})^{1/5}};$

where $1\leq i,j,k\leq n$ for $n=100,200,\ldots,600.$ It is not difficult to
see that these tensors have low tubal ranks. The numerical tubal rank for case
studies I, II and II for a tensor of size $100\times 100\times 100$ were 25, 5
and 43, respectively. However, for larger sizes the numerical tubal rank may
be slightly changed. We applied the proposed approach to the mentioned data
tensors in a similar way as for Example 1, to find the numerical tubal rank
and the corresponding low tubal rank approximation. Here, for the case studies
I, II and III, the proposed algorithm for $\epsilon=10^{(-8)}$ gave tubal
ranks 24, 5 and 42, respectively, which is very close to the true numerical
tubal ranks. Then for these tubal ranks, we applied the truncated t-SVD and
the randomized t-SVD to compute a low tubal rank approximation. The running
time and relative errors of the solutions of the algorithms are compared in
Figure 5 and Table 2, respectively. In view of Figure 5, the performance of
the proposed algorithm compared with the truncated t-SVD and the randomized
t-SVD is visible. These two experiments verified that the proposed approach is
applicable for large-scale tensors because it works only on a small a part of
the data tensor at each iteration while the classical approaches, e.g. the
truncated t-SVD deals with the whole data tensor. The results in Table 2 also
show that the proposed algorithm provides an approximation with almost the
same accuracy as the truncated t-SVD, which is known to be the best
approximation in the lease-squares sense222For any unitary invariant tensor
norm. for the low tubal rank approximation [9].

Table 1: Relative errors of results obtained by the truncated t-SVD and the
proposed Algorithm for Example 2. Case study I
---
Methods N | 100 | 200 | 300 | 400 | 500 | 600
Truncated t-SVD [9] | 1.6e-14 | 8.01e-13 | 6.9e-12 | 2.4e-11 | 5.9e-11 | 1.1e-10
Randomized t-SVD [22] | 4.9e-14 | 8.4e-12 | 4.4e-11 | 1.5e-10 | 3.7e-10 | 7.7e-10
Proposed algorithm | 3.1e-14 | 1.01e-12 | 3.5e-11 | 2.95e-10 | 3.2e-10 | 9.7e-10
Case study II
Truncated t-SVD [9] | 1.0e-15 | 1.6e-15 | 1.3e-15 | 1.7e-15 | 2.09e-15 | 2.5e-15
Randomized t-SVD [22] | 1.2e-15 | 1.5e-15 | 1.6e-15 | 2.4e-15 | 2.06e-15 | 2.3e-15
Proposed algorithm | 1.3e-15 | 1.5e-15 | 1.7e-15 | 1.8e-15 | 2.06e-15 | 2.3e-15
Case study III
Truncated t-SVD [9] | 2.7e-14 | 1.3e-11 | 1.3e-10 | 5.04e-10 | 1.1e-09 | 2.2e-09
Randomized t-SVD [22] | 1.9e-13 | 1.2e-10 | 1.01e-09 | 3.5e-09 | 6.5e-09 | 1.3e-08
Proposed algorithm | 5.3e-14 | 7.7e-11 | 4.7e-10 | 7.1e-09 | 3.2e-09 | 2.1e-08

######  Example 3.

Application to tensor completion. In this example, we show the application of
the proposed adaptive algorithm for the task of tensor completion. To this
end, we consider the benchmark images “Peppers”, “Lena” and “House” that of
size $256\times 256\times 3$ depicted in Figure 6 (left) and remove 70% of
their pixels randomly shown in Figure 6 (middle). We use the Peak signal-to-
noise ratio (PSNR) to compare the performance of the proposed algorithm with
the benchmark algorithm. The tensor decomposition formulation (18) for the
tensor completion problem is written as follows

$\begin{array}[]{cc}\displaystyle\min_{\underline{\bf X}}&{\|{{\bf
P}_{\underline{\bf\Omega}}}({\underline{\bf X}})-{{\bf
P}_{\underline{\bf\Omega}}}({\underline{\bf M}})\|^{2}_{F}},\\\
\textrm{s.t.}&{\rm rank}(\underline{\bf X})=R,\\\ \end{array}$ (18)

where the unknown tensor $\underline{\bf X}$ to be determined, and we assume
that it has low tensor rank representation, $\underline{\bf M}$ is the
original data tensor and ${\underline{\bf\Omega}}$ is the index of known
pixels. The projector ${\bf P}_{\underline{\bf\Omega}}$ is defined as follows

${\bf P}_{\underline{\bf\Omega}}({\bf X})=\left\\{\begin{array}[]{ll}{\bf
X}({\bf i})&{\bf i}\in\underline{\bf\Omega},\\\ 0&{\bf
i}\notin\underline{\bf\Omega},\end{array}\right.$

where ${\bf i}=(i_{1},i_{2},\ldots,i_{N})$ is an arbitrary multi-index with
$1\leq i_{n}\leq I_{n},\,\,n=1,2,\ldots,N$. Here, different kinds of tensor
ranks and associated tensor decompositions can be considered in the
formulation (18). The solution to the minimization problem (18) can be
approximated by the following iterative procedure

$\underline{\mathbf{Y}}^{(n)}=\mathcal{L}(\underline{\mathbf{X}}^{(n)}),$ (19)
$\underline{\mathbf{X}}^{(n+1)}=\underline{\mathbf{\Omega}}\oast\underline{\mathbf{Y}}^{(n)}+(\underline{\mathbf{1}}-\underline{\mathbf{\Omega}})\oast\underline{\mathbf{Y}}^{(n)},$
(20)

as described in [35] to complete the unknown pixels where $\mathcal{L}$ is an
operator, which computes a low-rank tensor approximation of the data tensor
$\underline{\mathbf{X}}^{(n)}$, $\underline{\mathbf{1}}$ is a tensor whose all
components are equal to one and $\oast$ is the Hadamard (elementwise) product.
For the low-rank computations in the first step (19), we apply the proposed
Algorithm 5 with a given number of iterations and not a given tolerance to
find the lateral and horizontal slice indices and compute the approximation
$\underline{\bf C}*\underline{\bf U}*\underline{\bf R}$ where $\underline{\bf
U}=\underline{\bf C}^{{\dagger}}*\underline{\bf X}*\underline{\bf
R}^{{\dagger}}$ and $\underline{\bf C},\,\underline{\bf R}$ are the sampled
lateral and horizontal slices, respectively. The tubal rank $R=70$, was used
in our computations. Beside applying the proposed algorithm, we also used the
truncated t-SVD and the randomized t-SVD in our computations. The
reconstructed images using the proposed approach and the truncated t-SVD are
displayed in Figure 6 (bottom). The running time required to compute these
reconstructions and also their PSNR are reported in Table 2. The results in
Table 2 and Figure 6, clearly illustrate that the proposed adaptive algorithm
provides comparable results in much less running time. This clearly shows the
feasibility and efficiency of the proposed algorithm for fast tensor
completion task.

Table 2: Running time (second ) and PSNR of the constructed images using the truncated t-SVD and the proposed Algorithm for Example 3, Time (s) and PSNR (dB). | Peppers | Lena | House
---|---|---|---
| Time | PSNR | Time | PSNR | Time | PSNR
Truncated t-SVD [9] | 55.03 | 27.74 | 45.2 | 27.36 | 61.94 | 28.05
Randomized t-SVD [22] | 10.12 | 27.55 | 11.23 | 27.40 | 10.83 | 29.40
Proposed algorithm | 7.60 | 27.55 | 8.34 | 27.46 | 7.27 | 29.59

######  Example 4.

Application to PEdesTrian Attribute Recognition task. In this experiment, we
show an application of the proposed method for the Pedestrian Attributes
Recognition (PAR) task [36]. We consider the PEdesTrian Attribute dataset
(PETA) dataset [37], which was widely used in the literature for the PAR
problem. It includes 8705 different persons in the 19000 pedestrian images,
which have 65 attributes (61 binary and 4 multi-class). Although the images
are not all the same size, we resize them in this experiment to $256\times
128\times 3$. We only take into account $1012$ images and construct a fourth-
order tensor with the size $256\times 128\times 3\times 1012$ and reshape it
into a third-order tensor with the size $256\times 384\times 1012$. With an
error bound of $\epsilon=0.1$, we applied the suggested approach to the above
dataset to determine the relevant tubal-rank. The truncated t-SVD of the
underlying dataset was then computed for this tubal-rank. The reconstructed
images that were obtained by them for two random samples using the proposed
algorithm and compared to the truncated t-SVD are shown in Figure 7. Here we
achieved $\times 5$ speed-up compared to the truncated t-SVD. The results
unequivocally show that the suggested technique can produce similar results in
less computational time. Additionally, the Attribute-specific Localization
(ASL) model [38], an effective deep neural network (DNN), was taken into
account as it had provided cutting-edge results for the PAR problem. In order
to create a lightweight model with fewer parameters and complexity [39], we
first compressed the underlying convolution layers in the ASL model using the
Error Preserving Correction-CPD [40] and the SVD. As our test datasets, we
also compressed $30\%$ of the PETA dataset’s images using Algorithm 5 (for
$\epsilon=0.1,0.2,$ and $0.3$) and compared the ASL model’s performance in
identifying pedestrian features in the compressed and original images. Table 3
displays the experiment’s outcomes. We see that the light-weight model’s
accuracy for both the original and compressed images is quite similar. It
should be noted that the model was not trained on compressed photos, though
one may do so to improve recognition accuracy. As a result, this concept may
be applied to Internet of Things (IoT) applications where tremendous amounts
of data in various shapes and formats are generated (for example, image
sensors embedded in mobile cameras produce enormous amounts of data in the
form of higher-resolution photographs and videos). Here, it is fundamental to
install compacted DNNs and portable DL models on the edge of the IoT network,
along with having fast data communication for real-time applications
(denoising, defogging, deblurring, segmentation, target detection, and
recognition). Using the suggested method, we may use the compressed form of
the data in these applications.

Figure 7: Comparing the original images and their compressed forms by the
proposed algorithm and the truncated t-SVD for Example 4. Table 3: Comparing
the running times (second) and relative errors achieved by the proposed
algorithm and the Truncated t-SVD [9] with original accuracy 0.8887 for
Example 4. $\epsilon=0.3$
---
Algorithms | Running Time (s) | Recognition accuracy
Truncated t-SVD [9] | 236.45 | 85.86%
Proposed algorithm | 45.56 | 84.36%
$\epsilon=0.2$
Truncated t-SVD [9] | 196.45 | 87.21%
Proposed algorithm | 35.97 | 86.16%
$\epsilon=0.1$
Truncated t-SVD [9] | 150.32 | 88.42%
Proposed algorithm | 27.12 | 87.21%

##  8 Conclusion and future works

In this work, we proposed an adaptive tubal tensor approximation algorithm for
the computation of the tensor SVD. The proposed algorithm can estimate the
tubal rank of a tensor and provide the corresponding low tubal rank
approximation. The experimental results verified the feasibility of the
proposed algorithm. Our future work will be developing a blocked version of
the proposed adaptive tubal tensor algorithm. The block version can be further
improved using the parallel hierarchical strategy [41] and we will investigate
this in future works. In the matrix case, it is known that the maximum volume
(maxvol) algorithm as a matrix cross approximation method provides close to
optimal low-rank approximations. Generalization of the maxvol approach from
the matrix case to tensors based on the t-product is our ongoing research
work.

##  9 Acknowledgement

The authors would like to thank the editor and two reviewer reviewers for
their constructive comments, which have greatly improved the quality of the
paper. The work was partially supported by the Ministry of Education and
Science of the Russian Federation (grant 075.10.2021.068).

##  10 Conflict of Interest Statement

The authors declare that they have no conflict of interest with anything.

## References

  * [1] I. V. Oseledets, Tensor-train decomposition, SIAM Journal on Scientific Computing 33 (5) (2011) 2295–2317.
  * [2] L. R. Tucker, et al., The extension of factor analysis to three-dimensional matrices, Contributions to mathematical psychology 110119 (1964).
  * [3] L. De Lathauwer, B. De Moor, J. Vandewalle, A multilinear singular value decomposition, SIAM journal on Matrix Analysis and Applications 21 (4) (2000) 1253–1278.
  * [4] F. L. Hitchcock, Multiple invariants and generalized rank of a p-way matrix or tensor, Journal of Mathematics and Physics 7 (1-4) (1928) 39–79.
  * [5] F. L. Hitchcock, The expression of a tensor or a polyadic as a sum of products, Journal of Mathematics and Physics 6 (1-4) (1927) 164–189.
  * [6] L. De Lathauwer, Decompositions of a higher-order tensor in block terms—Part II: Definitions and uniqueness, SIAM Journal on Matrix Analysis and Applications 30 (3) (2008) 1033–1066.
  * [7] Q. Zhao, G. Zhou, S. Xie, L. Zhang, A. Cichocki, Tensor ring decomposition, arXiv preprint arXiv:1606.05535 (2016).
  * [8] M. Espig, K. K. Naraparaju, J. Schneider, A note on tensor chain approximation, Computing and Visualization in Science 15 (6) (2012) 331–344.
  * [9] M. E. Kilmer, C. D. Martin, Factorization strategies for third-order tensors, Linear Algebra and its Applications 435 (3) (2011) 641–658.
  * [10] E. Newman, L. Horesh, H. Avron, M. Kilmer, Stable tensor neural networks for rapid deep learning, arXiv preprint arXiv:1811.06569 (2018).
  * [11] E. Newman, A step in the right dimension: Tensor algebra and applications, Ph.D. thesis, Tufts University (2019).
  * [12] Z. Zhang, S. Aeron, Exact tensor completion using t-svd, IEEE Transactions on Signal Processing 65 (6) (2016) 1511–1526.
  * [13] Z. Zhang, G. Ely, S. Aeron, N. Hao, M. Kilmer, Novel methods for multilinear data completion and de-noising based on tensor-svd, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2014, pp. 3842–3849.
  * [14] S. Soltani, M. E. Kilmer, P. C. Hansen, A tensor-based dictionary learning approach to tomographic image reconstruction, BIT Numerical Mathematics 56 (4) (2016) 1425–1454.
  * [15] M. Kilmer, L. Horesh, H. Avron, E. Newman, Tensor-tensor products for optimal representation and compression, arXiv preprint arXiv:2001.00046 (2019).
  * [16] M. W. Mahoney, et al., Randomized algorithms for matrices and data, Foundations and Trends® in Machine Learning 3 (2) (2011) 123–224.
  * [17] S. A. Goreinov, E. E. Tyrtyshnikov, N. L. Zamarashkin, A theory of pseudoskeleton approximations, Linear algebra and its applications 261 (1-3) (1997) 1–21.
  * [18] S. A. Goreinov, I. V. Oseledets, D. V. Savostyanov, E. E. Tyrtyshnikov, N. L. Zamarashkin, How to find a good submatrix, in: Matrix Methods: Theory, Algorithms And Applications: Dedicated to the Memory of Gene Golub, World Scientific, 2010, pp. 247–256.
  * [19] I. Oseledets, E. Tyrtyshnikov, TT-cross approximation for multidimensional arrays, Linear Algebra and its Applications 432 (1) (2010) 70–88.
  * [20] I. V. Oseledets, D. Savostianov, E. E. Tyrtyshnikov, Tucker dimensionality reduction of three-dimensional arrays in linear time, SIAM Journal on Matrix Analysis and Applications 30 (3) (2008) 939–956.
  * [21] C. F. Caiafa, A. Cichocki, Generalizing the column–row matrix decomposition to multi-way arrays, Linear Algebra and its Applications 433 (3) (2010) 557–573.
  * [22] D. A. Tarzanagh, G. Michailidis, Fast randomized algorithms for t-product based tensor operations and decompositions with applications to imaging data, SIAM Journal on Imaging Sciences 11 (4) (2018) 2629–2664.
  * [23] M. Bebendorf, Approximation of boundary element matrices, Numerische Mathematik 86 (4) (2000) 565–589.
  * [24] M. Bebendorf, R. Grzhibovskis, Accelerating galerkin bem for linear elasticity using adaptive cross approximation, Mathematical Methods in the Applied Sciences 29 (14) (2006) 1721–1747.
  * [25] K. Zhao, M. N. Vouvakis, J.-F. Lee, The adaptive cross approximation algorithm for accelerated method of moments computations of emc problems, IEEE transactions on electromagnetic compatibility 47 (4) (2005) 763–773.
  * [26] O. Rojo, H. Rojo, Some results on symmetric circulant matrices and on symmetric centrosymmetric matrices, Linear algebra and its applications 392 (2004) 211–233.
  * [27] C. Lu, J. Feng, Y. Chen, W. Liu, Z. Lin, S. Yan, Tensor robust principal component analysis with a new tensor nuclear norm, IEEE transactions on pattern analysis and machine intelligence 42 (4) (2019) 925–938.
  * [28] M. E. Kilmer, K. Braman, N. Hao, R. C. Hoover, Third-order tensors as operators on matrices: A theoretical and computational framework with applications in imaging, SIAM Journal on Matrix Analysis and Applications 34 (1) (2013) 148–172.
  * [29] C. D. Martin, R. Shafer, B. LaRue, An order-p tensor factorization with applications in imaging, SIAM Journal on Scientific Computing 35 (1) (2013) A474–A490.
  * [30] S. Ahmadi-Asl, An efficient randomized fixed-precision algorithm for tensor singular value decomposition, Communications on Applied Mathematics and Computation (2022) 1–20.
  * [31] N. Halko, P.-G. Martinsson, J. A. Tropp, Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions, SIAM review 53 (2) (2011) 217–288.
  * [32] D. Savostyanov, Polilinear approximation of matrices and integral equations, Ph. D. dissertation, Dept. Math., INM RAS, Moscow, Russia (2006).
  * [33] E. Tyrtyshnikov, Incomplete cross approximation in the mosaic-skeleton method, Computing 64 (4) (2000) 367–380.
  * [34] M. T. Chu, R. E. Funderlic, G. H. Golub, A rank–one reduction formula and its applications to matrix factorizations, SIAM review 37 (4) (1995) 512–530.
  * [35] S. Ahmadi-Asl, M. G. Asante-Mensah, A. Cichocki, A.-H. Phan, I. Oseledets, J. Wang, Cross tensor approximation for image and video completion, arXiv preprint arXiv:2207.06072 (2022).
  * [36] X. Wang, S. Zheng, R. Yang, A. Zheng, Z. Chen, J. Tang, B. Luo, Pedestrian attribute recognition: A survey, Pattern Recognition 121 (2022) 108220.
  * [37] Y. Deng, P. Luo, C. C. Loy, X. Tang, Pedestrian attribute recognition at far distance, in: Proceedings of the 22nd ACM international conference on Multimedia, 2014, pp. 789–792.
  * [38] C. Tang, L. Sheng, Z. Zhang, X. Hu, Improving pedestrian attribute recognition with weakly-supervised multi-scale attribute-specific localization, in: Proceedings of the IEEE/CVF International Conference on Computer Vision, 2019, pp. 4997–5006.
  * [39] A. Jha, D. Ermilov, K. Sobolev, A. H. Phan, S. Ahmadi-Asl, N. Ahmed, I. N. Junejo, Z. AL Aghbari, T. M. S. B. Shamsa, A. M. Khedr, A. Cichoki, Pedestrian attribute recognition using lightweight attribute specific localization, Submitted (2023).
  * [40] A.-H. Phan, P. Tichavskỳ, A. Cichocki, Error preserving correction: A method for cp decomposition at a target error bound, IEEE Transactions on Signal Processing 67 (5) (2018) 1175–1190.
  * [41] Y. Liu, W. Sid-Lakhdar, E. Rebrova, P. Ghysels, X. S. Li, A parallel hierarchical blocked adaptive cross approximation algorithm, The International Journal of High Performance Computing Applications 34 (4) (2020) 394–408.