remarkRemark hypothesisHypothesis claimClaim Variable T-Product and Tensor
CompletionLiqun Qi, Rui Yan, Ziyan Luo, Hong Yan, and Gaohang Yu
[][nocite]ex_supplement

# Variable T-Product and Zero-Padding Tensor Completion with
Applications††thanks: Submitted to the editors DATE. This work was supported
in part by National Natural Science Foundation of China (No.12071104), Natural
Science Foundation of Zhejiang Province (No.LD19A010002), Beijing Natural
Science Foundation (No.Z190002), Hong Kong Research Grants Council (Project
11204821), Hong Kong Innovation and Technology Commission(InnoHK Project
CIMDA) and City University of Hong Kong (Project 9610034).

Liqun Qi Department of Mathematics, Hangzhou Dianzi University, Hangzhou,
310018, China<EMAIL_ADDRESS>maghyu@hdu.edu.cn).
<EMAIL_ADDRESS>Rui Yan22footnotemark: 2 Ziyan Luo School of Mathematics
and Statistics, Beijing Jiaotong University, Beijing 100044, China ().
<EMAIL_ADDRESS>Hong Yan Department of Electrical Engineering, and
Center for Intelligent Multidimensional Data Analysis, City University of Hong
Kong, Kowloon, Hong Kong, China<EMAIL_ADDRESS>Gaohang
Yu22footnotemark: 2

###### Abstract

The T-product method based upon Discrete Fourier Transformation (DFT) has
found wide applications in engineering, in particular, in image processing. In
this paper, we propose variable T-product, and apply the Zero-Padding Discrete
Fourier Transformation (ZDFT), to convert third order tensor problems to the
variable Fourier domain. An additional positive integer parameter is
introduced, which is greater than the tubal dimension of the third order
tensor. When the additional parameter is equal to the tubal dimension, the
ZDFT reduces to DFT. Then, we propose a new tensor completion method based on
ZDFT and TV regularization, called VTCTF-TV. Extensive numerical experiment
results on visual data demonstrate that the superior performance of the
proposed method. In addition, the zero-padding is near the size of the
original tubal dimension, the resulting tensor completion can be performed
much better.

###### keywords:

Variable T-product, Zero-Padding Discrete Fourier Transform, Tensor Completion

46A32, 46B28, 47A80

##  1 Introduction

The Discrete Fourier Transformation (DFT) based the tensor SVD (T-SVD)
factorizations were proposed by Kilmer, Martin and others [14, 13], and have
found extensive applications for tensor completion in machine learning,
imaging processing, quantum information and other engineering fields [16, 7,
29, 33, 35, 34, 37, 40, 41, 42, 46, 38, 30, 32, 43]. In 2015, Knerfeld, Kilmer
and Aeron [15] further extended T-SVD to the so-called transformed T-SVD by
replacing the DFT with an arbitrary invertible linear transformation. Since
then, transformed T-SVD has attracted attention of many researcheers [11, 12,
28, 27, 31, 47, 18, 8, 44, 45].

A typical application of the T-product method is for color image processing,
in which the data tensor is an $m\times n\times p$ real third order tensor
${\mathcal{A}}=(a_{ijk})$. Here, $m\times n$ is the image size, and $p=3$ for
three color channels, red, green, and blue (RGB). One way to handle this
problem is to use a quaternion matrix to represent the third order tensor
${\mathcal{A}}$ [6, 17, 22, 26]. Another way to deal with this problem is to
regard ${\mathcal{A}}$ as a tubal matrix ${\mathcal{A}}=({\bf a}_{ij})$, where
the entries of the tubal matrix ${\mathcal{A}}$ are $p$-dimensional vectors.
For each pixel, a $p$-dimensional vector of data is acquired. Then, via the
DFT, the scalar product of entries in the classical matrix product is replaced
by a convolution-like product of the vector entries of the tubal matrix. In
this regard, the matrix SVD can be extended to T-SVD. This is the essence of
the T-product method.

Since the T-SVD method suffers from expensive computational cost, Zhou, Lu,
Lin and Zhang [46] proposed a tensor factorization model for the low rank
tensor completion problem (TCTF). Instead of analyzing the data in the
original image domain, they solve an optimization problem in the Fourier
domain. If the original data tensor is in $\mathbb{R}^{m\times n\times p}$,
then the Fourier domain is in $\mathbb{C}^{m\times n\times p}$.

However, in signal processing and image processing, direct application of the
DFT may not be enough. Zero-Padding [23, 9] is needed to obtain the accurate
result. Otherwise, artifacts may be produced. With the above motivation, in
this paper, we propose variable T-product, and apply the Zero-Padding Discrete
Fourier Transform (ZDFT) to the data tensor. Then we solve an optimization
problem in the variable Fourier domain in $\mathbb{C}^{m\times n\times v}$.
Here, $v$ is an additional positive integer parameter, with $v>p$. If $v=p$,
then the variable Fourier domain and the ZDFT reduce to the Fourier domain and
the DFT. The ZDFT matrix consists of the first $p$ columns of the $v\times v$
DFT matrix, and it maps an original $m\times n\times p$ tensor to the $m\times
n\times v$ variable Fourier domain. In theory, shorter signals $(v<2p-1)$
would produce artifacts, and longer zero-padding $(v>2p-1)$ should produce the
same result, but the result can be slightly different due to quantization of
the frequency variable in the Fourier domain. Therefore, when $v=2p-1$ or near
$2p$, the Zero-Padding takes effect, and the resulting tensor completion can
be performed much better. Numerical experiments show that the new method
performs better when $v$ is $2p-1$.

Throughout this paper, we assume that $m,n,p$ and $v$ are positive integers,
and $v>p$.

###  1.1 Contributions

Our contributions are summarized as below:

  * (1)

We propose variable T-product and apply the Zero-Padding Discrete Fourier
Transform (ZDFT), to convert third order tensor problems to the variable
Fourier domain $\mathbb{C}^{m\times n\times v}$.

  * (2)

We present a novel Zero-Padding tensor completion method, termed as VTCT F-TV,
for low-rank tensor completion problem, based on the ZDFT. Numerical
experiments verify the superior performance of the proposed method in image
and video inpainting when $v$ is or near $2p-1$.

  * (3)

We design an efficient proximal alternating minimization (PAM) algorithm to
solve the proposed tensor completion problem and establish its convergence
results.

##  2 Preliminaries

###  2.1 Tensor Notations and Definitions

An $M$th-order tensor $\mathcal{X}\\!\in\\!\mathbb{R}^{I_{1}\times\dots\times
I_{M}}$ is an $M$-way array consisting of entries $x_{i_{1}\cdots i_{m}}$ with
each $i_{j}$ varying among ${1,\cdots,I_{j}}$ for all $j=1,\cdots,M$. Vectors
and matrices are typical low order tensors with $M=1$ and $M=2$, respectively.
Define the inner product, the Frobenius norm and the $\infty$ norm for tensors
with symbols $\langle\cdot,\cdot\rangle$, $\|\cdot\|_{F}$,
$\|\cdot\|_{\infty}$ respectively, formulated by

$\langle\mathcal{X},\mathcal{Y}\rangle:=\sum_{i_{1}=1}^{I_{1}}\dots\sum_{i_{M}=1}^{I_{M}}x_{i_{1}\dots
i_{M}}y_{i_{1}\dots i_{M}},$
$\|\mathcal{X}\|_{F}:=\sqrt{\langle\mathcal{X},\mathcal{X}\rangle},~{}~{}\|\mathcal{X}\|_{\infty}:=\max\limits_{i_{1},\cdots,i_{M}}|x_{i_{1}\dots
i_{M}}|,$

for all $\mathcal{X},\mathcal{Y}\in\mathbb{R}^{I_{1}\times\dots\times I_{M}}$.

The Hadamard product of two vectors or matrices or tensors with the same
dimensions, is the component-wise product of these two vectors, or matrices,
or tensors.

In this paper, we deal with third order real tensors. A third order tensor
${\mathcal{A}}=(a_{ijk})\in\mathbb{R}^{m\times n\times p}$ has three modes:
row, column and tubal. Let $A^{(k)}={\mathcal{A}}(:,:,k)$ be the $k$th frontal
slice of ${\mathcal{A}}$ for $k=1,\cdots,p$. Given a positive integer $t\leq
p$, let ${\mathcal{A}}(:,:,[1:t])\in\mathbb{R}^{m\times n\times t}$ be the
subtensor of ${\mathcal{A}}$, generated by the first $t$ frontal slices. As
stated in the introduction, we may regard the third order tensor
${\mathcal{A}}=(a_{ijk})$ as a tubal matrix ${\mathcal{A}}=({\bf a}_{ij})$,
where ${\bf a}_{ij}={\mathcal{A}}(i,j,:)$ is a $p$-dimensional vector, also
called a tubal scalar [13].

We use calligraphic letters ${\mathcal{A}},\mathcal{B},\mathcal{C},\cdots$ to
denote tensors and tubal matrices, capital letters $A,B,C,\cdots$ to denote
matrices, small bold letters ${\bf a},{\bf b},{\bf c}$ to denote vectors and
tubal scalars, small letters $a,b,c.\cdots$ to denote scalars.

###  2.2 The DFT process

Let $\mathcal{C}\in\mathbb{R}^{m\times n\times p}$ be the data tensor. Let
$F=F_{p}$ be the $p\times p$ Discrete Fourier Transform (DFT) matrix with the
form

$F=\begin{bmatrix}1&1&\cdots&1&1\\\
1&\omega_{p}&\cdots&\omega_{p}^{p-2}&\omega_{p}^{p-1}\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\
1&\omega_{p}^{p-2}&\cdots&\omega_{p}^{(p-2)(p-2)}&\omega_{p}^{(p-2)(p-1)}\\\
1&\omega_{p}^{p-1}&\cdots&\omega_{p}^{(p-1)(p-2)}&\omega_{p}^{(p-1)(p-1)}\end{bmatrix},$

where

$\omega_{p}=e^{-{2\pi{\bf i}\over p}},$

and ${\bf i}=\sqrt{-1}$ is the imaginary unit. Then $F$ is invertible and
$F^{-1}={1\over p}F^{H}$. Apply $F$ to $\mathcal{C}$ along the third
dimension, we have $\bar{\mathcal{C}}(F)\in\mathbb{C}^{m\times n\times p}$
such that

(1) $\bar{\mathcal{C}}(F)(i,j,:)=F\mathcal{C}(i,j,:),$

for any $i$ and $j$. Then

(2) $\mathcal{C}(i,j,:)={1\over p}F^{H}\bar{\mathcal{C}}(F)(i,j,:),$

for any $i$ and $j$, and

(3) $\|\mathcal{C}\|_{F}^{2}={1\over p}\|\bar{\mathcal{C}}(F)\|_{F}^{2}.$

Define $\bar{C}(F)\in\mathbb{C}^{mp\times np}$ as a block diagonal matrix

(4) $\bar{C}(F)={\rm
diag}\left(\bar{C}(F)^{(1)},\bar{C}(F)^{(2)},\cdots,\bar{C}(F)^{(p)}\right).$

######  Definition 2.1.

For $\mathcal{C}\in\mathbb{R}^{m\times n\times p}$, its tensor tubal rank is
defined as

${\rm Rank}_{p}(\mathcal{C})=\max\\{{\rm Rank}(\bar{C}(F)^{(1)}),\cdots,{\rm
Rank}(\bar{C}(F)^{(p)})\\}.$

###  2.3 The TCTF Method

Since T-SVD suffers from expensive computational cost, Zhou, Lu, Lin and Zhang
[46] proposed a tensor factorization model for the low rank tensor completion
problem (TCTF). They work in the Fourier domain by solving the following
optimization problem:

(5) $\displaystyle\mathop{\min}\limits_{\hat{\cal X},\hat{\cal
Y},\mathcal{C}}$
$\displaystyle\frac{1}{2p}\sum_{k=1}^{p}\left\|{\hat{X}^{(k)}\hat{Y}^{(k)}-\bar{C}(F)^{(k)}}\right\|_{F}^{2}$
s.t. $\displaystyle{{\cal P}_{\Omega}}({\cal C}-{\cal G})=0$
$\displaystyle\hat{\mathcal{X}}\in\mathbb{C}^{m\times r\times
p},\hat{\mathcal{Y}}\in\mathbb{C}^{r\times n\times
p},\mathcal{C}\in\mathbb{R}^{m\times n\times p},$

where $\mathcal{C}$ and $\bar{\mathcal{C}}(F)$ are related by (1) and (2),
$\mathcal{G}=(g_{i_{1}i_{2}i_{3}})$ is the observed incomplete tensor data in
$\mathbb{R}^{m\times n\times p}$, $\Omega$ is the index set corresponding to
the observed entries of $\mathcal{G}$, and $P_{\Omega}$ is the linear operator
to extract known elements in the subset $\Omega$ and to fill the elements that
are not in $\Omega$ with zero values. Here, $r$ is the beforehand estimated
tubal rank and is usually much smaller than ${\min}\\{m,n\\}$.

##  3 Main Results

In Subsection 3.1, we provide the definition of variable T-product of two
third order tensors. In Subsection 3.2, we describe the ZDFT. We define the
variable tensor tubal rank in Subsection 3.3. A Zero-Padding tensor completion
method is proposed in Subsection 3.4.

###  3.1 Variable T-Product

We first define the variable product between two vectors in $\mathbb{R}^{p}$.
For ${\bf a}\in\mathbb{R}^{p}$ or $\mathbb{C}^{p}$, denote its $i$th component
as ${\bf a}(i)$.

######  Definition 3.1 (Variable product of two $p$-dimensional vectors).

For any ${\bf a},{\bf b}\in\mathbb{R}^{p}$ or $\mathbb{C}^{p}$, we have

(6) $\displaystyle({\bf a}\odot_{v}{\bf b})(k)$ (9) $\displaystyle=$
$\displaystyle\sum\left\\{{\bf a}(i){\bf b}(j):\begin{array}[]{l}i+j-k-1=0\
{\rm mod}(v),\\\ i,j=1,\cdots,p\end{array}\right\\},$

for $k=1,\cdots,p$. We call ${\bf a}\odot_{v}{\bf b}$ the variable product of
${\bf a}$ and ${\bf b}$.

One can see that when $v=2p-1$, ${\bf a}\odot_{v}{\bf b}$ is exactly the
$p$-truncated convolution of ${\bf a}$ and ${\bf b}$, since it consists of the
first $p$ components of the convolution. It is well-known that the convolution
of two $p$-dimensional vectors is a $(2p-1)$-dimensional vector. In order that
the variable T-product of two third order tensors can be defined, we need to
keep the variable product of two $p$-dimensional vectors in dimension $p$.
Thus, we cannot use convolution of two vectors to define the variable product
of two vectors directly here. But convolution is the best way to combine two
$p$-dimensional signals, and the variable product is the closest to the
convolution when $v=2p-1$. This explains why the tensor completion method
proposed later performs much better when $v=2p-1$.

Let ${\bf e}\in\mathbb{R}^{p}$ be defined by ${\bf e}(1)=1$, ${\bf e}(k)=0$,
for $k=2,\cdots,p$.

######  Proposition 3.2.

For any ${\bf a},{\bf b}\in\mathbb{R}^{p}$, we have

(10) ${\bf a}\odot_{v}{\bf b}={\bf b}\odot_{v}{\bf a},$

and

(11) ${\bf a}\odot_{v}{\bf e}={\bf e}\odot_{v}{\bf a}={\bf a}.$

######  Proof 3.3.

By (6), we have (10). For $k=1,\cdots,p$, we have

$\displaystyle({\bf a}\odot_{v}{\bf e})(k)$ $\displaystyle=$
$\displaystyle\sum\left\\{{\bf a}(i){\bf e}(j):\begin{array}[]{l}i+j-k-1=0\
{\rm mod}(v),\\\ i,j=1,\cdots,p\end{array}\right\\}$ $\displaystyle=$
$\displaystyle\sum\left\\{{\bf a}(i):i=k,i=1,\cdots,p\right\\}$
$\displaystyle=$ $\displaystyle{\bf a}(k).$

Thus, ${\bf a}\odot_{v}{\bf e}={\bf a}$. Similarly, we have ${\bf
e}\odot_{v}{\bf a}={\bf a}$. This proves (11).

Then $\mathbb{K}_{p}=(\mathbb{R}^{p},+,\odot_{v})$ is a commutative ring with
unity ${\bf e}$, with $+$ as the usual addition of vectors. As discussed in
[4], $\mathbb{K}_{p}$ is also a module. As in [13], we call an element ${\bf
a}\in\mathbb{K}_{p}$ a tubal scalar, and call $\mathbb{K}_{p}$ the
$p$-dimensional tubal scalar module. The zero tubal scalar ${\bf
0}\in\mathbb{K}_{p}$ has components ${\bf 0}(k)=0$ for $k=1,\cdots,p$.

For ${\bf a}\in\mathbb{K}_{p}$, define its modulus as

$|{\bf a}|=\sqrt{\sum_{k=1}^{p}{\bf a}(k)^{2}}.$

######  Definition 3.4.

Let ${\bf a}\in\mathbb{K}_{p}$. We say that ${\bf b}\in\mathbb{K}_{p}$ is the
transpose of ${\bf a}$, and denote ${\bf a}^{\top}={\bf b}$, if ${\bf
b}(1)={\bf a}(1)$, and ${\bf b}(k)={\bf a}(p+2-k)$ for $k=2,\cdots,p$.

For any ${\bf a}\in\mathbb{K}_{p}$, if ${\bf a}^{\top}={\bf a}$, then we say
that ${\bf a}$ is symmetric. The following proposition can be proved by
definition.

######  Proposition 3.5.

Let ${\bf a},{\bf b}\in\mathbb{K}_{p}$. Then $({\bf a}^{\top})^{\top}={\bf a}$
and $({\bf a}\odot_{v}{\bf b})^{\top}={\bf b}^{\top}\odot_{v}{\bf a}^{\top}$.

As stated early, a third order tensor
${\mathcal{A}}=(a_{ijk})\in\mathbb{R}^{m\times n\times p}$ can be regarded as
a tubal matrix ${\mathcal{A}}=({\bf a}_{ij})$, where ${\bf a}_{ij}(k)=a_{ijk}$
for $k=1,\cdots,p$.

######  Definition 3.6 (Variable T-product of third order tensors in the real
domain).

Suppose that ${\mathcal{A}}=({\bf a}_{il})\in\mathbb{R}^{m\times q\times p}$
and $\mathcal{B}=({\bf b}_{lj})\in\mathbb{R}^{q\times n\times p}$. The
variable T-product of ${\mathcal{A}}$ and $\mathcal{B}$, denoted as
$\mathcal{C}={\mathcal{A}}*_{v}\mathcal{B}$, is defined by $\mathcal{C}=({\bf
c}_{ij})$ with

(13) ${\bf c}_{ij}=\sum_{l=1}^{q}{\bf a}_{il}\odot_{v}{\bf b}_{lj}.$

Apparently, if $v=p$, then “$*_{v}$” reduces to the $T$-product “$*$”.

###  3.2 Zero-Padding Discrete Fourier Transform

Let $T=T_{vp}$ be the $v\times p$ Zero-Padding Discrete Fourier Transform
(ZDFT) matrix with the form

$T=\begin{bmatrix}1&1&\cdots&1&1\\\
1&\omega_{v}&\cdots&\omega_{v}^{p-2}&\omega_{v}^{p-1}\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\
1&\omega_{v}^{v-2}&\cdots&\omega_{v}^{(v-2)(p-2)}&\omega_{v}^{(v-2)(p-1)}\\\
1&\omega_{v}^{v-1}&\cdots&\omega_{v}^{(v-1)(p-2)}&\omega_{v}^{(v-1)(p-1)}\end{bmatrix},$

where

$\omega_{v}=e^{-{2\pi{\bf i}\over v}},$

and ${\bf i}=\sqrt{-1}$ is the imaginary unit. Then $T$ consists of the first
$p$ columns of the $v\times v$ DFT matrix $F_{v}$. Thus $T$ is of full column
rank $p$, and $T^{H}T=vI_{p}$, where $I_{p}$ is the $p\times p$ identity
matrix. For any ${\bf a}\in\mathbb{C}^{p}$, zero-padding ${\bf a}$ to obtain
$\tilde{{\bf a}}:=({\bf a};{\bf 0})\in\mathbb{C}^{v}$, we have

(14) $T{\bf a}=F_{v}\tilde{{\bf a}}.$

For this reason, $T$ is called the ZDFT matrix.

We may see that the variable T-product can be interpreted by the ZDFT.

For ${\bf a}\in\mathbb{C}^{p}$ and ${\bf c}\in\mathbb{C}^{v}$, define a
mapping $\phi:\mathbb{C}^{p}\to\mathbb{C}^{v}$ by $\phi({\bf a}):=T{\bf a}$,
and its conjugate transpose mapping $\phi^{H}:\mathbb{C}^{v}\to\mathbb{C}^{p}$
by $\phi^{H}({\bf c}):=T^{H}{\bf c}$. Then for any ${\bf a}\in\mathbb{C}^{p}$
and ${\bf c}\in\mathbb{C}^{v}$, we have

(15) $\phi({\bf a})(l)=\sum_{i=1}^{p}\omega^{(i-1)(l-1)}{\bf a}(i),$

for $l=1,\cdots,v$, and

(16) $\phi^{H}({\bf c})(k)=\sum_{l=1}^{v}\bar{\omega}^{(l-1)(k-1)}{\bf c}(l),$

for $k=1,\cdots,p$.

######  Proposition 3.7.

For any ${\bf a},{\bf b}\in\mathbb{R}^{p}$, we have

(17) ${\bf a}\odot_{v}{\bf b}={1\over v}\phi^{H}[\phi({\bf a})\circ\phi({\bf
b})],$

where $\circ$ represents the Hadamard product.

######  Proposition 3.8.

Suppose that ${\mathcal{A}}=({\bf a}_{il})\in\mathbb{R}^{m\times q\times p}$
and $\mathcal{B}=({\bf b}_{lj})\in\mathbb{R}^{q\times n\times p}$. Let
${\mathcal{A}}_{0}\in\mathbb{R}^{m\times q\times v}$ and
$\mathcal{B}_{0}\in\mathbb{R}^{q\times n\times v}$ be the zero-padding
counterparts of ${\mathcal{A}}$ and $\mathcal{B}$, with
${\mathcal{A}}_{0}(:,:,k)$ and $\mathcal{B}_{0}(:,:,k)$ all zero frontal
slices for any $k=p+1,\cdots,v$. Then

(18)
${\mathcal{A}}*_{v}\mathcal{B}=\left({\mathcal{A}}_{0}*\mathcal{B}_{0}\right)(:,:,[1:p]),$

where $*$ is the T-product.

###  3.3 Variable Tensor Tubal Rank and H-product

We now use $T$ to replace $F$ in Subsection 2.2. Again, let
$\mathcal{C}\in\mathbb{R}^{m\times n\times p}$ be the data tensor. Apply $T$
to $\mathcal{C}$ along the third dimension, we have
$\bar{\mathcal{C}}(T)\in\mathbb{C}^{m\times n\times v}$ such that

(19) $\bar{\mathcal{C}}(T)(i,j,:)=T\mathcal{C}(i,j,:),$

for any $i$ and $j$. Then

(20) $\mathcal{C}(i,j,:)={1\over v}T^{H}\bar{\mathcal{C}}(T)(i,j,:),$

for any $i$ and $j$, and

(21) $\|\mathcal{C}\|_{F}^{2}={1\over v}\|\bar{\mathcal{C}}(T)\|_{F}^{2}.$

Define $\bar{C}(T)\in\mathbb{C}^{mv\times nv}$ as a block diagonal matrix

(22) $\bar{C}(T)={\rm
diag}\left(\bar{C}(T)^{(1)},\bar{C}(T)^{(2)},\cdots,\bar{C}(T)^{(v)}\right).$

Clearly, the equations (19) and (20) define two mappings between $\mathcal{C}$
and $\bar{C}(T)$, respectively. For simplicity, let’s use
$\bar{C}(T)=\mbox{fft}_{p2v}(\mathcal{C},[],3)$ to denote the relation (19),
and $\mathcal{C}=\mbox{ifft}_{v2p}(\bar{C}(T),[],3)$ to denote the relation
(20).

######  Definition 3.9 (Variable tensor tubal rank).

For $\mathcal{C}\in\mathbb{R}^{m\times n\times p}$, its variable tensor tubal
rank is defined as

${\rm Rank}_{v}(\mathcal{C})=\max\\{{\rm Rank}(\bar{C}(T)^{(1)}),\cdots,{\rm
Rank}(\bar{C}(T)^{(v)})\\}.$

Obviously, if $v=p$, then ${\rm Rank}_{v}(\mathcal{C})$ is exactly the tubal
rank ${\rm Rank}_{t}(\mathcal{C})$.

We may define H-product for tensors in $\mathbb{C}^{m\times n\times v}$.
Again, we regard tensors in $\mathbb{C}^{m\times n\times v}$ as tubal
matrices.

######  Definition 3.10 (H-product for third order tensors in the variable
Fourier domain).

Suppose that ${\mathcal{A}}=({\bf a}_{il})\in\mathbb{C}^{m\times q\times v}$
and $\mathcal{B}=({\bf b}_{lj})\in\mathbb{C}^{q\times n\times v}$. The
Hadamard tensor product, or simply called the H-product of ${\mathcal{A}}$ and
$\mathcal{B}$, denoted as $\mathcal{C}={\mathcal{A}}*_{H}\mathcal{B}$, is
defined by $\mathcal{C}=({\bf c}_{ij})$ with

(23) ${\bf c}_{ij}=\sum_{l=1}^{q}{\bf a}_{il}\circ{\bf b}_{lj}.$

######  Theorem 3.11.

Suppose that $\mathcal{C}\in\mathbb{R}^{m\times n\times p}$. Then the variable
tensor tubal rank of $\mathcal{C}$ is not greater than $r$ if and only if
there are tensors $\bar{{\mathcal{X}}}\in\mathbb{C}^{m\times r\times v}$ and
$\bar{{\mathcal{Y}}}\in\mathbb{C}^{r\times n\times v}$ such that
$\bar{\mathcal{C}}(T)=\bar{{\mathcal{X}}}*_{H}\bar{{\mathcal{Y}}}$.
Furthermore, if there are ${\mathcal{A}}\in\mathbb{R}^{m\times r\times p}$ and
$\mathcal{B}\in\mathbb{R}^{r\times n\times p}$ such that
$\mathcal{C}={\mathcal{A}}*_{v}\mathcal{B}$, then the variable tensor tubal
rank of $\mathcal{C}$ is not greater than $r$.

###  3.4 A Zero-Padding Tensor Completion Method

Based on the Zero-Padding Discrete Fourier Transform, we propose the following
Zero-Padding Tensor Completion Method. In order to accurately recover the
missing data, spatio-temporal prior knowledge is introduced, and TV
regularization is applied to tensor completion problem by the application of
total variation smoothing prior method in spatio-temporal video [19]. In this
paper, the low-rank completion problem of third-order tensor based on variable
tubal rank and TV regularization is studied. The model is as follows:

(24)
$\displaystyle\min_{\mathcal{X},\mathcal{Y},\mathcal{C}}\frac{1}{2}\left\|\mathcal{X}\ast_{v}\mathcal{Y}-\mathcal{C}\right\|_{F}^{2}+\alpha_{1}\left\|\mathcal{D}_{1}\ast_{v}\mathcal{C}\right\|_{l_{1}}+\alpha_{2}\left\|\mathcal{C}\ast_{v}\mathcal{D}_{2}\right\|_{l_{1}}$
$\displaystyle s.t.\mathcal{C}_{\Omega}=\mathcal{G}_{\Omega}$
$\displaystyle\mathcal{X}\in\mathbb{R}^{m\times q\times
p},\mathcal{Y}\in\mathbb{R}^{q\times n\times
p},\mathcal{C}\in\mathbb{R}^{m\times n\times p},$

where $\mathcal{G}$ represents the observation tensor, $\Omega$ represents the
index set of observation elements, and
$\mathcal{C}_{\Omega}=\mathcal{G}_{\Omega}$ represents that the elements of
$\mathcal{C}$ and $\mathcal{G}$ in $\Omega$ are consistent, and

$\mathcal{D}_{1}\left(:,:,i\right)=\begin{cases}L_{m},i=1\\\
O_{m},i>1\end{cases},$
$\mathcal{D}_{2}\left(:,:,i\right)=\begin{cases}L_{n}^{\ast},i=1\\\
O_{n},i>1\end{cases},$ $L_{m}=\begin{pmatrix}0&0&0&\cdots&0&0\\\
-1&1&0&\cdots&0&0\\\ 0&-1&1&\cdots&0&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&0&\cdots&1&0\\\
0&0&0&\cdots&-1&1\end{pmatrix}\in\mathbb{R}^{m\times m},$

$\alpha_{1}$ and $\alpha_{2}$ are regularization parameters and
$\alpha_{1}\left\|\mathcal{D}_{1}\ast_{v}\mathcal{C}\right\|_{l_{1}}+\alpha_{2}\left\|\mathcal{C}\ast_{v}\mathcal{D}_{2}\right\|_{l_{1}}$
is TV regularization term. The above optimization problems can be transformed
into the following unconstrained optimization problems:

(25) $\min
f\left(\mathcal{X},\mathcal{Y},\mathcal{C}\right)=\frac{1}{2}\left\|\mathcal{X}\ast_{v}\mathcal{Y}-\mathcal{C}\right\|_{F}^{2}+\alpha_{1}\left\|\mathcal{D}_{1}\ast_{v}\mathcal{C}\right\|_{l_{1}}+\alpha_{2}\left\|\mathcal{C}\ast_{v}\mathcal{D}_{2}\right\|_{l_{1}}+\Phi\left(\mathcal{C}\right),$

where

$\Phi\left(\mathcal{C}\right)=\begin{cases}0,\mathcal{C}_{\Omega}=\mathcal{G}_{\Omega},\\\
\infty,otherwise.\end{cases}$

Optimization problem (25) is not a joint convex function about
$\left(\mathcal{X},\mathcal{Y},\mathcal{C}\right)$, but it is convex for each
variable $\mathcal{X}$, $\mathcal{Y}$ and $\mathcal{C}$. Because of the high
efficiency of Alternate Minimization (AM) algorithm, it is often used to solve
multivariate optimization problems. In order to improve the theoretical
convergence and numerical stability of the algorithm, this paper adds a
neighboring term to the sub-problem of AM algorithm, that is, the Proximal
Alternate Minimization (PAM) algorithm. Let
$f\left(\mathcal{X},\mathcal{Y},\mathcal{C}\right)$ be the objective function
of problem (25) and given the initial point
$\left(\mathcal{X}^{k},\mathcal{Y}^{k},\mathcal{C}^{k}\right)$, use PAM
algorithm to solve the framework, and update each variable alternately as
follows:

(26)
$\begin{cases}\mathcal{X}^{k+1}=\underset{\mathcal{X}}{argmin}f\left(\mathcal{X},\mathcal{Y}^{k},\mathcal{C}^{k}\right)+\frac{\rho_{1}}{2}\left\|\mathcal{X}-\mathcal{X}^{k}\right\|_{F}^{2}\\\
\mathcal{Y}^{k+1}=\underset{\mathcal{Y}}{argmin}f\left(\mathcal{X}^{k+1},\mathcal{Y},\mathcal{C}^{k}\right)+\frac{\rho_{2}}{2}\left\|\mathcal{Y}-\mathcal{Y}^{k}\right\|_{F}^{2}\\\
\mathcal{C}^{k+1}=\underset{\mathcal{C}}{argmin}f\left(\mathcal{X}^{k+1},\mathcal{Y}^{k+1},\mathcal{C}\right)+\frac{\rho_{3}}{2}\left\|\mathcal{C}-\mathcal{C}^{k}\right\|_{F}^{2},\end{cases}$

where $\rho_{1}$, $\rho_{2}$ and $\rho_{3}$ are given parameters, and
$\rho_{1},\rho_{2},\rho_{3>0}$. It can be seen that all subproblems are
strongly convex optimization problems, the existence and uniqueness of
solutions are guaranteed, and all of them have explicit solutions. The details
are as follows:
$\mathcal{X}$-subproblem:

(27)
$\mathcal{X}^{k+1}=\left(\rho_{1}\mathcal{X}^{k}+\mathcal{C}^{k}\ast_{v}\left(\mathcal{Y}^{k}\right)^{\ast}\right)\left(\mathcal{Y}^{k}\ast_{v}\left(\mathcal{Y}^{k}\right)^{\ast}+\rho_{1}I\right)^{\dagger},$

$\mathcal{Y}$-subproblem:

(28)
$\mathcal{Y}^{k+1}=\left(\left(\mathcal{X}^{k+1}\right)^{\ast}\ast_{v}\mathcal{X}^{k+1}+\rho_{2}I\right)^{\dagger}\left(\left(\mathcal{X}^{k+1}\right)^{\ast}\ast_{v}\mathcal{C}^{k}+\rho_{2}\mathcal{Y}^{k}\right),$

$\mathcal{C}$-subproblem:

$\mathcal{C}^{k+1}=\underset{\mathcal{C}}{argmin}f\left(\mathcal{X}^{k+1},\mathcal{Y}^{k+1},\mathcal{C}\right)+\frac{\rho_{3}}{2}\left\|\mathcal{C}-\mathcal{C}^{k}\right\|_{F}^{2}.$

The above problem is equivalent to the following equality constraint problem:

(29) $\displaystyle\mathcal{C}^{k+1}\\!=$
$\displaystyle\underset{\mathcal{C},\mathcal{Q}_{1},\mathcal{Q}_{2}}{argmin}\frac{1}{2}\left\|\mathcal{X}^{k+1}\\!\ast_{v}\mathcal{Y}^{k+1}\\!-\\!\mathcal{C}\right\|_{F}^{2}+\\!\alpha_{1}\left\|\mathcal{Q}_{1}\right\|_{l_{1}}+\\!\alpha_{2}\left\|\mathcal{Q}_{2}\right\|_{l_{1}}+\\!\frac{\rho_{3}}{2}\left\|\mathcal{C}\\!-\\!\mathcal{C}^{k}\right\|_{F}^{2}+\\!\Phi\left(\mathcal{C}\right)$
$\displaystyle
s.t.\mathcal{Q}_{1}=\mathcal{D}_{1}\ast_{v}\mathcal{C},\mathcal{Q}_{2}=\mathcal{C}\ast_{v}\mathcal{D}_{2}.$

Let $\mathcal{S}$ and $\mathcal{T}$ be Lagrangian multipliers of (29), then

(30)
$\mathcal{Q}_{1}=\underset{\mathcal{Q}_{1}}{argmin}\alpha_{1}\left\|\mathcal{Q}_{1}\right\|_{l_{1}}+\frac{\beta}{2}\left\|\mathcal{Q}_{1}-\mathcal{D}_{1}\ast_{v}\mathcal{C}+\frac{\mathcal{S}}{\beta}\right\|_{F}^{2},$
(31)
$\mathcal{Q}_{2}=\underset{\mathcal{Q}_{2}}{argmin}\alpha_{2}\left\|\mathcal{Q}_{2}\right\|_{l_{1}}+\frac{\mu}{2}\left\|\mathcal{Q}_{2}-\mathcal{C}\ast_{v}\mathcal{D}_{2}+\frac{\mathcal{T}}{\mu}\right\|_{F}^{2},$
(32) $\displaystyle\mathcal{C}^{k+1}=$
$\displaystyle\underset{\mathcal{C}}{argmin}\frac{1}{2}\left\|\mathcal{X}^{k+1}\ast_{v}\mathcal{Y}^{k+1}-\mathcal{C}\right\|_{F}^{2}+\frac{\beta}{2}\left\|\mathcal{Q}_{1}-\mathcal{D}_{1}\ast_{v}\mathcal{C}+\frac{\mathcal{S}}{\beta}\right\|$
$\displaystyle+\frac{\mu}{2}\left\|\mathcal{Q}_{2}-\mathcal{C}\ast_{v}\mathcal{D}_{2}+\frac{\mathcal{T}}{\mu}\right\|_{F}^{2}+\frac{\rho_{3}}{2}\left\|\mathcal{C}-\mathcal{C}^{k}\right\|+\Phi\left(\mathcal{C}\right),$
(33)
$\mathcal{S}=\mathcal{S}+\beta\left(\mathcal{Q}_{1}-\mathcal{D}_{1}\ast_{v}\mathcal{C}\right),$
(34)
$\mathcal{T}=\mathcal{T}+\mu\left(\mathcal{Q}_{2}-\mathcal{C}\ast_{v}\mathcal{D}_{2}\right).$

According to the soft threshold, equations (30) and (31) have the following
unique solutions:

(35)
$\mathcal{Q}_{1}\left(i,j,k\right)=T_{\frac{\alpha_{1}}{\beta}}\left(\mathcal{D}_{1}\ast_{v}\mathcal{C}-\frac{\mathcal{S}}{\beta}\right)\left(i,j,k\right),$
(36)
$\mathcal{Q}_{2}\left(i,j,k\right)=T_{\frac{\alpha_{2}}{\mu}}\left(\mathcal{C}\ast_{v}\mathcal{D}_{2}-\frac{\mathcal{T}}{\mu}\right)\left(i,j,k\right),$

where

$T_{\eta}\left(x\right)=\begin{cases}\left(\left|x\right|-\eta\right)sign\left(x\right),\left|x\right|>\eta\\\
0,\left|x\right|\leq\eta,\end{cases}$

where $x\in\mathbb{R}$, and formula (32) is equivalent to

(37) $\displaystyle\bar{C}^{k+1}=$
$\displaystyle\underset{\bar{C}}{argmin}\frac{1}{2v}\left\|\bar{X}^{k+1}\bar{Y}^{k+1}-\bar{C}\right\|_{F}^{2}+\frac{\beta}{2v}\left\|\bar{Q}_{1}-\bar{D}_{1}\bar{C}+\frac{\bar{S}}{\beta}\right\|_{F}^{2}$
$\displaystyle+\frac{\mu}{2v}\left\|\bar{Q}_{2}-\bar{C}\bar{D}_{2}+\frac{\bar{T}}{\mu}\right\|+\frac{\rho_{3}}{2v}\left\|\bar{C}-\bar{C}^{k}\right\|_{F}^{2}+\Phi\left(\bar{C}\right).$

The only solution of the above formula is actually the solution of the
following matrix equation:

(38)
$\bar{C}+\beta\bar{D}_{1}^{\ast}\bar{D}_{1}\bar{C}+\mu\bar{C}\bar{D}_{2}\bar{D}_{2}^{\ast}+\rho_{3}\bar{C}=R^{k},$

where

$R^{k}=\bar{X}^{k+1}\bar{Y}^{k+1}+\bar{D}_{1}^{\ast}\left(\beta\bar{Q}_{1}+\bar{S}\right)+\left(\mu\bar{Q}_{2}+\bar{T}\right)\bar{D}_{2}^{\ast}+\rho_{3}\bar{C}^{k}.$

The matrices in formula (38) all have block diagonal structures. From the
definitions of $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$, it can be seen that
the diagonal blocks of $\bar{D}_{1}$ and $\bar{D}_{2}$ are the same, so this
formula is equivalent to $v$ matrix equations with smaller sizes, namely

(39) $\bar{C}_{l}+\beta
H_{m}\bar{C}_{l}+\mu\bar{C}_{l}H_{n}+\rho_{3}\bar{C}_{l}=R_{l}^{k},l=1,2,\cdots,v,$

where $\bar{C}_{l}$ and $R_{l}^{k}$ represent the $l$-th diagonal blocks of
$\bar{C}$ and $R^{k}$, respectively, and

$H_{m}:=L_{m}^{T}L_{m}=\begin{pmatrix}1&-1&&&\\\ -1&2&-1&&\\\
&\ddots&\ddots&\ddots&\\\ &&-1&2&-1\\\
&&&-1&1\end{pmatrix}\in\mathbb{R}^{m\times m},H_{m}^{0}:=I_{m}.$

Figure 1: The difference between TCTF and V-TCTF

It is easy to verify that $H_{m}$ has the following orthogonal diagonalization
form:

$H_{m}=K_{m}\Lambda_{m}K_{m}^{T},$
$K_{m}=\sqrt{\frac{2}{m}}\left[\sqrt{\left(1+\delta_{j,1}\right)^{-1}}cos\left(\frac{\pi\left(2i-1\right)\left(j-1\right)}{2m}\right)\right]_{i,j=1}^{m},$
$\delta_{j,1}=\begin{cases}1,j=1\\\ 0,otherwise\end{cases},\\\
\Lambda_{m}=4\times
diag\left(sin^{2}\left(\frac{\left(i-1\right)\pi}{2m}\right)\right)_{i=1}^{m}.$

That is, $\Lambda_{m}$ is a nonnegative diagonal matrix of $m\times m$, so
equation (39) is equivalent to

(40) $\bar{C}_{l}+\beta
K_{m}\Lambda_{m}K_{m}^{T}\bar{C}_{l}+\mu\bar{C}_{l}K_{n}\Lambda_{n}K_{n}^{T}+\rho_{3}\bar{C}_{l}=R_{l}^{k},l=1,2,\cdots,v.$

Multiply $K_{m}^{T}$ by left and $K_{n}$ by right on both sides of equation
(40) respectively to obtain

(41)
$\left(1+\rho_{3}\right)\hat{C}_{l}+\beta\Lambda_{m}\hat{C}_{l}+\hat{C}_{l}\mu\Lambda_{n}=\hat{R}_{l}^{k},l=1,2,\cdots,v,$

where

$\hat{C}_{l}=K_{m}^{T}\bar{C}_{l}K_{n},\hat{R}_{l}^{k}=K_{m}^{T}R_{l}^{k}K_{n},$

and

(42)
$\hat{C}_{l}\left(i,j\right)=\frac{\hat{R}_{l}^{k}\left(i,j\right)}{1+\rho_{3}+\beta\Lambda_{m}\left(i\right)+\mu\Lambda_{n}\left(j\right)},\left(i,j\right)\in\Xi\left(m\right)\times\Xi\left(n\right).$

By using

$\bar{C}_{l}=K_{m}\hat{C}_{l}K_{n}^{T},$

then we obtain

(43)
$\mathcal{C}^{k+1}\left(i,j,k\right)=\begin{cases}ifft_{v2p}\left(\bar{\mathcal{C}},\left[\right],3\right),\left(i,j,k\right)\notin\Omega\\\
\mathcal{G}\left(i,j,k\right),\left(i,j,k\right)\in\Omega.\end{cases}$

For clarity, the algorithmic framework is presented in Algorithm 1.

When $\alpha_{1}=\alpha_{2}=0$, the model is simplified as:

(44) $\displaystyle\mathop{\min}\limits_{{\cal X},{\cal Y},{\mathcal{C}}}$
$\displaystyle\frac{1}{2}\left\|{\cal X}*_{v}{\cal
Y}-\mathcal{C}\right\|_{F}^{2}$ s.t.
$\displaystyle\mathcal{C}_{\Omega}=\mathcal{G}_{\Omega}$
$\displaystyle{\mathcal{X}}\in\mathbb{R}^{m\times q\times
p},{\mathcal{Y}}\in\mathbb{R}^{q\times n\times
p},\mathcal{C}\in\mathbb{R}^{m\times n\times p}.$

Using the same computational scheme proposed in [46] to solve the Zero-Padding
TCTF (44) in the variable Fourier domain, we term it as V-TCTF. In summary,
the difference between TCTF and our proposed V-TCTF is illustrated in Fig. 1.

Algorithm 1 The VTCTF-TV algorithm to solve (25)

Input: The tensor data $\mathcal{G}\in\mathbb{R}^{m\times n\times p}$, the
observed set $\Omega$, the parameter $v=2p-1$, the initialized rank
$r^{0}\in\mathbb{R}^{v}$, and $\varepsilon=1e-5$, iteration $N$=200.

Initialize: ${\cal X}^{0}$, ${\cal Y}^{0}$;

While $k\leq N$ and not converged do

1\. Fix $\mathcal{Y}^{k}$ and $\mathcal{C}^{k}$ to update $\mathcal{X}^{k+1}$
via (27).

2\. Fix $\mathcal{X}^{k+1}$ and $\mathcal{C}^{k}$ to update
$\mathcal{Y}^{k+1}$ via (28).

3\. Fix $\mathcal{X}^{k+1}$ and $\mathcal{Y}^{k+1}$ to update
$\mathcal{C}^{k+1}$ via (43).

4\. Check the termination
criterion:$\frac{\left\|\mathcal{C}^{k+1}-\mathcal{C}^{k}\right\|_{F}^{2}}{\left\|\mathcal{C}^{k+1}\right\|_{F}^{2}}\leq\varepsilon$.

end while

Output: $\mathcal{C}^{k+1}$.

Obviously, the objective function in model (24) is proper lower semi-
continuous, and only contains a simple projection constraint condition.
According to the algorithm design, the generated point sequence
$\left\\{\left(\mathcal{X}^{k},\mathcal{Y}^{k},\mathcal{C}^{k}\right)\right\\}$
meets the conditions H1, H2 and H3 in Section 2.3 of reference[1], so it can
be proved that it converges to the critical point of model (24). The concrete
convergence analysis is given in Appendix D.

##  4 Numerical Experiments

In this section, Algorithm 1 is used to restore color images and multispectral
images to evaluate its performance, and compared with the following four data
completion methods, namely TCTF[46], Tmac[36], TCTFTVT[19] and MTRTC[39]. In
order to quantitatively evaluate the image quality restored by each method,
PSNR[5], SSIM[21] and CPU time are used as numerical indicators, and PSNR and
SSIM are defined as follows:

$\begin{split}&PSNR=10\cdot
log_{10}\frac{mnp\left\|\mathcal{C}_{true}\right\|_{\infty}^{2}}{\left\|\mathcal{C}-\mathcal{C}_{true}\right\|_{F}^{2}},\\\
&SSIM=\frac{\left(2\mu_{\mathcal{C}}\mu_{\mathcal{C}_{true}}\right)\left(2\sigma_{\mathcal{C}\mathcal{C}_{true}}+c_{2}\right)}{\left(\mu_{\mathcal{C}}^{2}\mu_{\mathcal{C}_{true}}^{2}+c_{1}\right)\left(\sigma_{\mathcal{C}}^{2}+\sigma_{\mathcal{C}_{true}}^{2}+c_{2}\right)},\end{split}$

where $\mathcal{C}_{true}$ is the real tensor, $\mathcal{C}$ is the recovery
tensor, $\mu_{\mathcal{C}}$ and $\mu_{\mathcal{C}_{true}}$ are the mean values
of images $\mathcal{C}$ and $\mathcal{C}_{true}$ , $\sigma_{\mathcal{C}}$ and
$\sigma_{\mathcal{C}_{true}}$ are the standard deviations of images
$\mathcal{C}$ and $\mathcal{C}_{true}$ ,
$\sigma_{\mathcal{C}\mathcal{C}_{true}}$ is the covariance of images
$\mathcal{C}$ and $\mathcal{C}_{true}$ , $c_{1}$ and $c_{2}>0$ are constants.
The higher the $PSNR$ and $SSIM$ values, the better the image quality.

The parameters of the VTCTF-TV method proposed in this paper are set as
follows:

$\alpha_{1}=\alpha_{2}=\beta=\mu=\varepsilon=1e-5,\rho_{1}=\rho_{2}=\rho_{3}=5e-6,N=200.$

Set the initial rank to $r^{0}=\left[30,\cdots,30\right]\in\mathbb{R}^{v}$,
the experiment is implemented on Matlab R2016b under Windows 10, equipped with
a 3.00GHz CPU and a PC with 8GB of memory.

###  4.1 Gray video

In order to find the optimal value of the introduced additional parameter $v$,
we evaluate our method on the widely used YUV Video Sequences111
http://trace.eas.asu.edu/yuv/. Each sequence contains at least 150 frames. In
the experiments, we first test V-TCTF on the Hall Monitor video. The frame
size of this video is 144$\times$176 pixels. Due to the computational
limitation, we only use the first 20 frames of the sequences, thus p = 20
here. The sampling rate is set to 0.3, 0.5 and 0.7, and the initialized rank
is set to $r^{0}=\left[30,\cdots,30\right]\in\mathbb{R}^{v}$.

$\begin{array}[]{cc}\includegraphics[width=203.80193pt,height=108.405pt]{figure2.pdf}&\includegraphics[width=203.80193pt,height=108.405pt]{figure3.pdf}\\\
\includegraphics[width=203.80193pt,height=108.405pt]{figure4.pdf}&\includegraphics[width=203.80193pt,height=108.405pt]{figure5.pdf}\\\
\includegraphics[width=203.80193pt,height=108.405pt]{figure6.pdf}&\includegraphics[width=203.80193pt,height=108.405pt]{figure7.pdf}\\\
\end{array}$ Figure 2: The PSNR and SSIM value v.s. different v for V-TCTF
method on the hall video, from top to bottom, the sampling rate is 0.3, 0.5
and 0.7.

As shown in Fig. 2, we display the PSNR and SSIM comparison of the testing
video. We can see from the Fig. 2 that the resulting tensor completion can be
performed much better when $v=2p-1$ or near $2p$.

###  4.2 Color images

In this subsection, we conduct experiments on six popular color images (Lena,
Panda, Sailboat, Barbara, House and Pepper) that are widely used in the
literature. All images are of size 256 $\times$ 256 $\times$ 3\. We compared
the performance of our method with five other image completion methods,
including TCTF method, TMac, TCTFTVT method and MTRTC method.

Figure 3: Restored results of color images with sampling rate 0.7. From top to
bottom, Lena, Panda, Sailboat, Barbara, House and Pepper. From left to right:
the original data, the observed data, the recovered results by TCTF, Tmac,
TCTFTVT, MTRTC and VTCTF-TV, respectively.

In Table 1, by comparing the peak signal-to-noise ratio, structural similarity
and operation time at different sampling rates of 0.6, 0.7 and 0.8, it is easy
to see that the peak signal-to-noise ratio and structural similarity of VTCTF-
TV are higher than those of TCTF, TMac, TCTFTVT and MTRTC.

Color image | SR | TCTF | Tmac | TCTFTVT | MTRTC | VTCTF-TV
---|---|---|---|---|---|---
Lena | 0.6 | 28.40/0.83/4.10 | 29.90/0.86/9.00 | 30.18/0.90/7.63 | 29.47/0.84/3.31 |  30.60/0.93/2.78
0.7 | 29.78/0.86/3.68 | 30.79/0.88/8.73 | 31.52/0.92/7.83 | 30.77/0.86/3.10 |  32.47/0.96/3.85
0.8 | 31.68/0.91/3.23 | 31.69/0.90/8.49 | 33.35/0.95/7.80 | 33.15/0.91/2.76 |  34.87/0.97/3.07
Panda | 0.6 | 30.40/0.85/3.25 | 30.99/0.83/7.74 | 31.67/0.88/8.29 | 31.28/0.87/3.86 |  31.88/0.90/3.70
0.7 | 31.83/0.88/2.97 | 31.92/0.86/8.84 | 33.00/0.91/7.89 | 33.14/0.91/2.83 |  33.32/0.92/3.10
0.8 | 33.70/0.92/3.51 | 32.79/0.88/8.44 | 34.79/0.94/7.89 | 35.32/0.94/2.83 |  35.65/0.95/3.10
Sailboat | 0.6 | 25.73/0.79/3.24 | 27.32/0.81/7.66 | 27.48/0.85/7.28 | 26.97/0.83/3.83 |  27.76/0.90/14.98
0.7 | 27.16/0.84/3.75 | 28.21/0.84/9.68 | 28.86/0.89/7.22 | 28.68/0.87/3.45 |  29.81/0.94/4.10
0.8 | 29.06/0.89/3.68 | 29.13/0.87/8.46 | 30.74/0.92/7.60 | 30.92/0.91/2.79 |  32.17/0.96/2.93
Barbara | 0.6 | 27.85/0.84/3.42 | 29.41/0.86/8.09 |  29.73/0.89/7.45 | 28.16/0.86/3.41 | 29.70/0.91/14.81
0.7 | 29.27/0.87/3.80 | 30.26/0.88/8.53 | 31.10/0.92/7.98 | 29.37/0.88/3.16 |  32.01/0.94/3.78
0.8 | 31.13/0.91/3.38 | 31.14/0.90/8.59 | 32.96/0.95/8.04 | 31.57/0.91/2.63 |  34.37/0.97/3.02
House | 0.6 | 30.52/0.86/2.89 | 33.09/0.89/7.84 | 33.18/0.91/8.12 | 28.98/0.80/3.12 |  33.48/0.92/4.63
0.7 | 31.90/0.89/3.83 | 34.04/0.91/8.74 | 34.58/0.94/7.53 | 30.98/0.85/2.84 |  35.56/0.95/3.39
0.8 | 33.73/0.93/3.62 | 34.99/0.92/9.31 | 36.39/0.96/7.66 | 33.32/0.90/2.65 |  37.83/0.96/2.92
Pepper | 0.6 | 25.15/0.72/3.37 | 28.97/0.84/8.77 | 29.48/0.88/8.18 | 24.49/0.71/3.12 |  30.66/0.93/5.11
0.7 | 26.57/0.77/3.59 | 29.92/0.87/8.31 | 30.79/0.91/8.26 | 25.72/0.75/2.83 |  33.02/0.96/3.90
0.8 | 28.41/0.83/3.27 | 30.96/0.89/8.94 | 32.69/0.94/8.47 | 27.45/0.80/2.61 |  35.30/0.97/3.39

Table 1: The PSNR/SSIM/CPU time of different methods on color images “Lena”,
“Panda”, “Sailboat”, “Barbara”, “House” and “Pepper” with different sampling
rates.

Fig. 3 is a restored image of six color images in models TCTF, TMac, TCTFTVT,
MTRTC and VTCTF-TV with the sampling rate set to 0.7. It can be observed that
the recovery effect of VTCTF-TV is better than that of TCTF, TMac, TCTFTVT and
MTRTC.

###  4.3 Multispectral images

In this section, the performance of different methods on multi-spectral images
is tested, and three multi-spectral
images222https://www.cs.columbia.edu/CAVE/databases/multispectral/stuff/ are
numerically tested, and their peak signal-to-noise ratio, structural
similarity and CPU time are compared. The size of three multispectral images
is 512 $\times$512 $\times$31.

Figure 4: PSNR values of each band of the recovered MSI Pompoms .
Figure 5: SSIM values of each band of the recovered MSI Pompoms .
Figure 6: The one mode-3 tube of recovered MSIs by different methods with
SR=0.7. From top to bottom: Pompoms , Flowers and Stuffed toys, respectively.
From left to right: the original data, the observed data, the recovered
results by TCTF, Tmac, TCTFTVT, MTRTC and VTCTF-TV, respectively.

In Table 2, by comparing the average peak signal-to-noise ratio, average
structural similarity and CPU time at different sampling rates of 0.6, 0.7 and
0.8, it is easy to see that the peak signal-to-noise ratio and structural
similarity of VTCTF-TV are higher than those of TCTF, TMac, TCTFTVT and MTRTC.
Fig. 4 and Fig. 5 give the PSNR and SSIM values comparison of each band of the
MSI Pompoms recovered by all methods, respectively. The results on the other
two datasets are similar.

In addition, Fig. 6 displays one mode-3 tube of the recovered MSIs by
different methods with SR = 0.7, we clearly observe that VTCTF-TV yields the
closest mode-3 tubes in all cases.

Fig. 7 is the restored image of the 10th band of three multispectral images
when the sampling rate of models TCTF, TMac, TCTFTVT, MTRTC and VTCTF-TV is
set to 0.7. From left to right are the original picture, the observed image
and the experimental results under TCTF, TMac, TCTFTVT, MTRTC and VTCTF-TV
algorithms. From the visual comparison, it is clear that VTCTF-TV performs
best in preserving multispectral image edges and details. In a word, VTCTF-TV
has better recovery effect than TCTF, TMac, TCTFTVT and MTRTC.

Multispectral image | SR | TCTF | Tmac | TCTFTVT | MTRTC | VTCTF-TV
---|---|---|---|---|---|---
Pompoms | 0.6 | 28.21/0.74/95.81 | 36.80/0.94/100.94 | 36.72/0.94/175.47 | 34.53/0.90/85.09 |  41.31/0.97/148.09
0.7 | 29.47/0.77/123.47 | 37.42/0.95/122.71 | 38.04/0.95/178.49 | 36.84/0.94/75.43 |  43.08/0.98/96.18
0.8 | 31.06/0.82/105.76 | 38.09/0.95/113.28 | 39.86/0.97/178.73 | 36.97/0.94/67.87 |  44.80/0.99/74.49
Flowers | 0.6 | 31.89/0.87/100.76 | 36.78/0.94/106.15 | 36.53/0.95/215.37 | 36.92/0.94/192.47 |  39.17/0.97/212.02
0.7 | 32.78/0.89/139.82 | 37.41/0.95/151.71 | 37.91/0.97/247.35 | 37.54/0.95/140.43 |  41.10/0.98/108.00
0.8 | 34.11/0.92/108.01 | 38.09/0.96/115.71 | 39.76/0.98/177.81 | 39.02/0.96/146.64 |  43.50/0.99/86.34
Stuffed toys | 0.6 | 39.73/0.97/107.05 | 40.87/0.97/121.59 | 40.83/0.98/156.97 | 37.10/0.92/107.25 |  41.83/0.96/116.93
0.7 | 41.29/0.97/108.68 | 41.55/0.98/121.62 | 42.20/0.98/160.13 | 38.65/0.94/118.04 |  43.53/0.97/89.61
0.8 | 34.97/0.97/116.43 | 42.28/0.98/126.03 | 44.06/0.99/162.50 | 43.96/0.98/118.74 |  45.89/0.99/73.23

Table 2: The PSNR/SSIM/CPU time of different methods on multispectral images “
Pompoms ”, “ Flowers ” and “ Stuffed toys ” with different sampling rates.
Figure 7: Restored results on MSIs with 0.7 observed entries. From top to
bottom: The tenth band of MSI Pompoms , Flowers and Stuffed toys. From left to
right: the original data, the observed data, the recovered results by TCTF,
Tmac, TCTFTVT, MTRTC and VTCTF-TV, respectively.

##  5 Conclusion

In this paper, we propose variable T-product, and apply the Zero-Padding
Discrete Fourier Transformation (ZDFT), to convert third order tensor problems
to the variable Fourier domain, and a novel tensor completion model based on
ZDFT and TV regularization is proposed. In the model, tensor decomposition
based on variable T-product is used to describe the low rank of tensor, and TV
regularization is used to describe the smoothness of image data. Numerical
experiments on color images and multispectral images show the effectiveness of
the proposed algorithm.

##  Appendix A Proof of Proposition 3.7

######  Proof A.1.

Let ${\bf a},{\bf b}\in\mathbb{R}^{p}$, and ${\bf c}=\phi({\bf
a})\circ\phi({\bf b})$. For $k=1,\cdots,p$, we have

$\displaystyle\phi^{H}({\bf c})(k)$ $\displaystyle=$ $\displaystyle{1\over
v}\sum_{l=1}^{v}\bar{\omega}^{(l-1)(k-1)}{\bf c}(l)$ $\displaystyle=$
$\displaystyle{1\over
v}\sum_{l=1}^{v}\sum_{i=1}^{p}\sum_{j=1}^{p}\bar{\omega}^{(l-1)(k-1)}\omega^{(i-1)(l-1)}\omega^{(j-1)(l-1)}{\bf
a}(i){\bf b}(j)$ $\displaystyle=$ $\displaystyle{1\over
v}\sum_{i=1}^{p}\sum_{j=1}^{p}{\bf a}(i){\bf
b}(j)\left[\sum_{l=1}^{v}\omega^{(l-1)(i+j-k-1)}\right].$

Note that

$\sum_{l=1}^{v}\omega^{(l-1)(i+j-k-1)}=\left\\{\begin{aligned} v,&\ \ {\rm if\
}i+j-k-1=0,\ {\rm mod}(v)\\\ 0,&\ \ {\rm otherwise}.\end{aligned}\right.$

Thus,

$\phi^{H}({\bf c})(k)=\sum\left\\{{\bf a}(i){\bf
b}(j):\begin{array}[]{l}i+j-k-1=0\ {\rm mod}(v),\\\
i,j=1,\cdots,p\end{array}\right\\}.$

By (6), $\phi^{H}({\bf c})={\bf a}\odot_{v}{\bf b}$. This proves (17).

##  Appendix B Proof of Proposition 3.8

######  Proof B.1.

Denote $F_{v}=[T~{}N]$, where $N\in\mathbb{C}^{(v-p)\times v}$ consists of the
last $v-p$ columns of $F_{v}$. Note that $F_{v}^{-1}=\frac{1}{v}F^{H}$. Thus,

(45) $F_{v}^{-1}([1:p],:)=\frac{1}{v}T^{H}.$

Let $\mathcal{C}:={\mathcal{A}}*_{v}\mathcal{B}$ and ${\cal
D}:={\mathcal{A}}_{0}*\mathcal{B}_{0}$. It follows from Proposition 3.7 that
for any $i\in\\{1,\cdots,m\\}$ and any $j\in\\{1,\cdots,n\\}$,

$\displaystyle{\bf c}_{ij}$ $\displaystyle=$
$\displaystyle\frac{1}{v}\phi^{H}\left[\sum\limits_{l=1}^{q}\phi({\bf
a}_{il})\circ\phi({\bf b}_{lj})\right]$ $\displaystyle=$
$\displaystyle\frac{1}{v}T^{H}\left[\sum\limits_{l=1}^{q}T{\bf a}_{il}\circ
T{\bf b}_{lj}\right]$ $\displaystyle=$
$\displaystyle\frac{1}{v}T^{H}\left[\sum\limits_{l=1}^{q}F_{v}\tilde{{\bf
a}}_{il}\circ F_{v}\tilde{{\bf b}}_{lj}\right]$ $\displaystyle=$
$\displaystyle F_{v}^{-1}([1:p],:)\left[\sum\limits_{l=1}^{q}F_{v}({\bf
a}_{0})_{il}\circ F_{v}({\bf b}_{0})_{lj}\right]$ $\displaystyle=$
$\displaystyle\left(F_{v}^{-1}\left[\sum\limits_{l=1}^{q}F_{v}({\bf
a}_{0})_{il}\circ F_{v}({\bf b}_{0})_{lj}\right]\right)(1:p)$ $\displaystyle=$
$\displaystyle{\bf d}_{ij}(1:p),$

where the third equality is from (14) and the fourth one from (45). Thus,
${\cal C}={\cal D}(:,:,[1:p])$. This completes the proof.

##  Appendix C Proof of Theorem 3.11

######  Proof C.1.

The equivalence in the first part can be derived by using

$\bar{\mathcal{C}}(T)=\bar{{\cal X}}*_{H}\bar{{\cal
Y}}\Longleftrightarrow\bar{C}(T)^{(k)}=\bar{X}^{(k)}\bar{Y}^{(k)},k=1,\cdots,v$

from the definitional expression (23), and by the definition of the variable
tubal rank as stated in Definition 3.9. To show the “furthermore” part, let
${\mathcal{A}}_{0}$ and $\mathcal{B}_{0}$ be the zero-padding counterparts of
${\mathcal{A}}$ and $\mathcal{B}$, respectively, and ${\cal D}$ be the
T-product of ${\mathcal{A}}_{0}$ and $\mathcal{B}_{0}$, i.e., ${\cal
D}:={\mathcal{A}}_{0}*\mathcal{B}_{0}$. It is known from Proposition 3.8 that
Rank${}_{v}(\mathcal{C})\leq{\rm Rank}_{t}{\cal D}$. Additionally, by Lemma 2
in [46], we have ${\rm Rank}_{t}{\cal D}\leq\min\\{{\rm
Rank}_{t}({\mathcal{A}}_{0}),{\rm Rank}_{t}(\mathcal{B}_{0})\\}\leq r$. Thus,
Rank${}_{v}(\mathcal{C})\leq r$. This completes the whole proof.

##  Appendix D Convergence analysis for Algorithm 1

In this section, we will prove the global convergence of the proposed
algorithm. First and foremost, for positive integers $I_{1}$,$I_{2}$,$I_{3}$,
define bijections
$\mathbb{V}_{[I_{1},I_{2},I_{3}]}:\mathcal{A}\in\mathbb{R}^{I_{1}\times
I_{2}\times
I_{3}}\mapsto\mathbb{V}_{[I_{1},I_{2},I_{3}]}(\mathcal{A})\in\mathbb{R}^{I_{1}I_{2}I_{3}\times
1}$ and
$\mathbb{I}_{[I_{1},I_{2},I_{3}]}:\left(i_{1},i_{2},i_{3}\right)\in\Xi\left(I_{1}\right)\times\Xi\left(I_{2}\right)\times\Xi\left(I_{3}\right)\mapsto\mathbb{I}_{\left[I_{1},I_{2},I_{3}\right]}\left(i_{1},i_{2},i_{3}\right)\in\Xi\left(I_{1}I_{2}I_{3}\right)$
by:

$\mathbb{V}_{\left[I_{1},I_{2},I_{3}\right]}\left(\mathcal{A}\right):=\mathcal{A}\left(:\right),\\\
$
$\mathcal{A}\left(i_{1},i_{2},i_{3}\right)=\left[\mathcal{A}\left(:\right)\right]\left(\mathbb{I}_{\left[I_{1},I_{2},I_{3}\right]}\left(i_{1},i_{2},i_{3}\right)\right)\mbox{holding
for each }\mathcal{A}\in\mathbb{R}^{I_{1}\times I_{2}\times I_{3}}.$

Denote $d_{1}=mqp$ , $d_{2}=qnp$ ,$d_{3}=mnp$, $d=d_{1}+d_{2}+d_{3}$, and
denote $\mathbb{V}_{\left[m,q,p\right]}$, $\mathbb{V}_{\left[q,n,p\right]}$,
$\mathbb{V}_{\left[m,n,p\right]}$, $\mathbb{I}_{\left[m,n,p\right]}$ as
$\mathbb{V}_{1}$,$\mathbb{V}_{2}$,$\mathbb{V}_{3}$,$\mathbb{I}$, respectively.
For the sake of convenience, we will represent the variables
$\mathcal{X},\mathcal{Y}$, and $\mathcal{C}$ in tensor form as vectors, as
follows:

$x=V_{1}\left(\mathcal{X}\right),y=V_{2}\left(\mathcal{Y}\right),c=V_{3}\left(\mathcal{C}\right)$

Then, we rewrite the objective function (25) as

(46)
$\tilde{F}\left(x,y,c\right)=G\left(x,y,c\right)+H\left(c\right)+\Phi\left(c\right)$

where $x\in\mathbb{R}^{d_{1}},y\in\mathbb{R}^{d_{2}},c\in\mathbb{R}^{d_{3}}$

$\displaystyle\tilde{F}\left(x,y,c\right):=G\left(x,y,c\right)+H\left(c\right)+\Phi\left(c\right),$
$\displaystyle
G\left(x,y,c\right)=\frac{1}{2}\|\mathbb{V}_{1}^{-1}\left(x\right)*\mathbb{V}_{2}^{-1}\left(y\right)-\mathbb{V}_{3}^{-1}\left(c\right)\|_{\ell_{2}}^{2},$
$\displaystyle
H\left(c\right):=\alpha_{1}\|D_{1}\left(c\right)\|_{\ell_{1}}+\alpha_{2}\|D_{2}\left(c\right)\|_{\ell_{1}},$
$\displaystyle
D_{1}\left(c\right):=\mathcal{D}_{1}*_{v}\left[\mathbb{V}_{3}^{-1}\left(c\right)\right],$
$\displaystyle
D_{2}\left(c\right):=\left[\mathbb{V}_{3}^{-1}\left(c\right)\right]*_{v}\mathcal{D}_{2},$
$\displaystyle\tilde{S}:=\left\\{w\in\mathbb{R}^{d_{3}}|w\left(i\right)=\left[\mathbb{V}_{3}\left(\mathcal{G}\right)\right]\left(i\right)\mbox{
for
}i\in\mathbb{I}\left(\Omega\right),w\left(i\right)\in\left[0,1\right]\mbox{
for }i\notin\mathbb{I}\left(\Omega\right)\right\\},$
$\displaystyle\Phi\left(c\right)=\begin{cases}0,c\in\tilde{S}\\\
+\infty,c\notin\tilde{S}\end{cases}$

Because of the bilinear operation of ”$*_{v}$,” it is evident that both
$D_{1}:\mathbb{R}^{mnp\times 1}\rightarrow\mathbb{R}^{m\times n\times p}$ and
$D_{2}:\mathbb{R}^{mnp\times 1}\rightarrow\mathbb{R}^{m\times n\times p}$ are
linear operators. Moreover, it is easy to see that $\tilde{S}$ is a non-empty
closed set, which means that $\Phi\left(c\right)$ is a proper lower semi-
continuous (PLSC) function on $\mathbb{R}^{d_{3}}$ and
$\tilde{F}\left(\cdot\right)$ is a PLSC function on $\mathbb{R}^{d}$. (26) is
equivalent to

(47)
$\begin{cases}x^{k+1}\in\mathop{\arg\min}_{x}\tilde{F}\left(x,y^{k},c^{k}\right)+\frac{\rho_{1}}{2}\|x-x^{k}\|_{2}^{2},\\\
y^{k+1}\in\mathop{\arg\min}_{y}\tilde{F}\left(x^{k+1},y,c^{k}\right)+\frac{\rho_{2}}{2}\|y-y^{k}\|_{2}^{2},\\\
c^{k+1}\in\mathop{\arg\min}_{c}\tilde{F}\left(x^{k+1},y^{k+1},c\right)+\frac{\rho_{3}}{2}\|c-c^{k}\|_{2}^{2}\end{cases}$

To show the convergence of PAM algorithm, the following convergence theory is
needed.

######  Definition D.1.

(KŁ property[1])

  1. (a)

The function
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\left\\{+\infty\right\\}$ is said to
have the KŁ property at $\bar{x}\in dom\left(\partial f\right)$, if there
exist $\eta\in\left(0,+\infty\right]$, a neighbourhood U of $\bar{x}$ and a
continuous concave function
$\phi:\left[0,\eta\right)\rightarrow\left[0,+\infty\right)$, such that:

    1. (i)

$\phi\left(0\right)=0$

    2. (ii)

$\phi$ is first-order continuous on $\left(0,\eta\right)$ ,

    3. (iii)

$\phi^{\prime}$ is positive on $\left(0,\eta\right)$ ,

    4. (iv)

for each $x\in
U\cap\left[f\left(\bar{x}\right)<f<f\left(\bar{x}\right)+\eta\right]$, the KŁ
inequality holds:

$\phi^{\prime}\left(f\left(x\right)-f\left(\bar{x}\right)\right)dist\left(0,\partial
f\left(x\right)\right)\geq 1.$

  2. (b)

PLSC functions which satisfies the KŁ property at each point of
$dom\left(\partial f\right)$ are called KŁ functions, where the norm involved
in $dist\left(\cdot,\cdot\right)$ is $\|\cdot\|_{2}$ and the convention
$dist\left(0,\emptyset\right):=+\infty$.

######  Lemma D.2.

[1] Let $f:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\left\\{+\infty\right\\}$ be
a PLSC function. Let
$\left\\{x^{k}\right\\}_{k\in\mathbb{N}}\subset\mathbb{R}^{n}$ be a sequence
such that

  * H1

$\left(\mbox{Sufficient decrease condition}\right)$ For each $k\in\mathbb{N}$,
there exits $a\in\left(0,+\infty\right)$ such that
$f\left(x^{k+1}\right)+a\|x^{k+1}-x^{k}\|_{2}^{2}\leq f\left(x^{k}\right)$
hold

  * H2

$\left(\mbox{Relative error condition}\right)$ For each $k\in\mathbb{N}$,
there exits $w^{k+1}\in\partial f\left(x^{k+1}\right)$ and a constant
$b\in\left(0,+\infty\right)$ such that$\|w^{k+1}\|_{2}\leq
b\|x^{k+1}-x^{k}\|_{2}$ hold

  * H3

$\left(\mbox{Continuity condition}\right)$There exists a subsequence
$\left\\{x^{k_{j}}\right\\}_{j\in\mathbb{N}}$ and $\bar{x}\in\mathbb{R}^{n}$
such that

$x^{k_{j}}\rightarrow\bar{x}\mbox{ and }f\left(x^{k_{j}}\right)\rightarrow
f\left(\bar{x}\right),j\rightarrow\infty.$

If $f$ has the KŁ property at $\bar{x}$, then

  1. (i)

$x^{k}\rightarrow\bar{x}$

  2. (ii)

$\bar{x}$ is a critical point of $f$, i.e., $0\in\partial
f\left(\bar{x}\right)$;

  3. (iii)

the sequence $\left\\{x^{k}\right\\}_{k}\in\mathbb{N}$ has a finite length,
i.e.,

$\sum_{k=0}^{+\infty}\|x^{k+1}-x^{k}\|_{2}<+\infty.$

Second, we prove that the objective function $\tilde{F}$ in (46) and the
iterative sequence $(x^{k},y^{k},c^{k})_{k\in\mathbb{N}}$ generated by
iteration (47) satisfy the assumptions in LEMMA D.2. Thus we establish the
convergence of the proposed algorithm.

For convenience, we denote
$\hat{S}=\mathbb{R}^{d_{1}}\times\mathbb{R}^{d_{2}}\times\tilde{S}\subset\mathbb{R}^{d}$.
Next, we need to prove that $\tilde{F}$ satisfies the KŁ property at each
$\left(x,y,c\right)\in\hat{S}$.

######  Lemma D.3.

$\tilde{F}$ satisfies the KŁ property at each $\left(x,y,c\right)\in\hat{S}$.

######  Proof D.4.

For any $\left(x,y,c\right)\in\hat{S}$, $\tilde{F}$ can be expressed as
follows:

$\tilde{F}:=G\left(x,y,c\right)+H\left(c\right)$

Obviously, $\hat{S}$ is a semi-algebraic set([19] Proposition 2).

For convenience, we denote
$u:=\left(u_{1},u_{2},u_{3}\right),u_{1}:=x,u_{2}:=y,u_{3}:=c$ . $D_{1}$ and
$D_{2}$ are linear mappings between finite dimensional spaces, each element of
$D_{1}\left(c\right)$ and $D_{2}\left(c\right)$ can be regarded as a linear
polynomial of $u:=(x,y,c)$, so $H\left(c\right)$ is only a finite sum of the
composition of the absolute value function and linear polynomial of $u$.
Therefore, $H\left(c\right)$ is semi-algebraic on $\hat{S}$([19] Proposition
2,3). $\mathbb{V}_{1}^{-1}\left(x\right)$,
$\mathbb{V}_{2}^{-1}\left(y\right)$, and $\mathbb{V}_{3}^{-1}\left(c\right)$
are linear mappings between finite-dimensional spaces. Therefore, each element
of $\mathbb{V}_{i}^{-1}\left(u_{i}\right)\left(i=1,2,3\right)$ is a linear
polynomial of $u:=\left(u_{1},u_{2},u_{3}\right)$. In addition, according to
the definition of $\|\cdot\|_{\ell_{2}}$, we know that $G\left(u\right)$ is a
polynomial of $u$. Thus, $G$ is semi-algebraic on $\hat{S}$([19] Proposition
2).

In a word, $\tilde{F}$ is semi-algebraic on $\hat{S}$([19] Proposition 2). A
semi-algebraic function is a KŁ function (see[2, 3]). Hence, $\tilde{F}$
satisfies the KŁ property at each $\left(x,y,c\right)\in\hat{S}$.

######  Theorem D.5.

Assume that the sequence $u^{k}=(x^{k},y^{k},c^{k})$ generated by iteration
(47) is bounded, then it converges to a critical point of $\tilde{F}$.

######  Proof D.6.

As mentioned above, $\tilde{F}$ is a PLSC function on $\mathbb{R}^{d}$. From
(47), and we see that

$\displaystyle\begin{cases}\tilde{F}(x^{k+1},y^{k},c^{k})+\frac{\rho_{1}}{2}\|x^{k+1}-x^{k}\|_{2}^{2}\leq\tilde{F}(x^{k},y^{k},c^{k}),k\in\mathbb{N}\\\
\tilde{F}(x^{k+1},y^{k+1},c^{k})+\frac{\rho_{2}}{2}\|y^{k+1}-y^{k}\|_{2}^{2}\leq\tilde{F}(x^{k+1},y^{k},c^{k}),k\in\mathbb{N}\\\
\tilde{F}(x^{k+1},y^{k+1},c^{k+1})+\frac{\rho_{3}}{2}\|c^{k+1}-c^{k}\|_{2}^{2}\leq\tilde{F}(x^{k+1},y^{k+1},c^{k}),k\in\mathbb{N}.\end{cases}$

Combining the above three inequalities, we get

(48) $\tilde{F}(u^{k+1})+\frac{\rho_{\rm
min}}{2}\|u^{k+1}-u^{k}\|_{2}^{2}\leq\tilde{F}(u^{k}),k\in\mathbb{N}.$

where $\rho_{\rm min}=\min\\{\rho_{1},\rho_{2},\rho_{3}\\}$. Therefore,
H1(sufficient decrease condition) is satisfied, and $a=\rho_{\rm min}/2$.

Denote $u:=\left(u_{1},u_{2},u_{3}\right),u_{1}:=x,u_{2}:=y,u_{3}:=c$,
$M=\\{u^{k}|k\in\mathbb{N}\\}$, for $u=(u_{1},u_{2},u_{3})\in\mathbb{R}^{s}$,
$u_{i}\in\mathbb{R}^{d_{i}}(i=1,2,3)$. Because $G$ is a polynomial function,
so it is infinitely differentiable. From (47) and Proposition 1[19],

(49)
$\begin{cases}0\in\nabla_{u_{1}}G(u_{1}^{k+1},u_{2}^{k},u_{3}^{k})+\rho_{1}(u_{1}^{k+1}-u_{1}^{k}),k\in\mathbb{N}\\\
0\in\nabla_{u_{2}}G(u_{1}^{k+1},u_{2}^{k+1},u_{3}^{k})+\rho_{2}(u_{2}^{k+1}-u_{2}^{k}),k\in\mathbb{N}\\\
0\in\nabla_{u_{3}}G(u_{1}^{k+1},u_{2}^{k+1},u_{3}^{k+1})+\partial(\Phi+H)(u_{3}^{k+1})+\rho_{3}(u_{3}^{k+1}-u_{3}^{k}),k\in\mathbb{N}.\end{cases}$

(49) means the existence of $w_{1}^{k+1}\in\partial(\Phi+H)(u_{3}^{k+1})$, so
that

$-\tilde{w}^{k+1}=\tilde{u}^{k}+(\rho_{1}(u_{1}^{k+1}-u_{1}^{k}),\rho_{2}(u_{2}^{k+1}-u_{2}^{k}),\rho_{3}(u_{3}^{k+1}-u_{3}^{k})),k\in\mathbb{N},$

where

$\displaystyle\tilde{w}^{k+1}=(0,0,w_{1}^{k+1}),$
$\displaystyle\tilde{u}^{k}=(\nabla_{u_{1}}G(u_{1}^{k+1},u_{2}^{k},u_{3}^{k}),\nabla_{u_{2}}G(u_{1}^{k+1},u_{2}^{k+1},u_{3}^{k}),\nabla_{u_{3}}G(u_{1}^{k+1},u_{2}^{k+1},u_{3}^{k+1})).$

Let $w^{k+1}=\tilde{w}^{k+1}+\nabla G(u^{k+1})$, obviously
$w^{k+1}\in\partial\tilde{F}(u^{k+1})$ and

(50) $-w^{k+1}=\tilde{u}^{k}-\nabla
G(u^{k+1})+(\rho_{1}(u_{1}^{k+1}-u_{1}^{k}),\rho_{2}(u_{2}^{k+1}-u_{2}^{k}),\rho_{3}(u_{3}^{k+1}-u_{3}^{k})),k\in\mathbb{N},$

Define the coordinate projections by

$\Psi_{i}(u):=u_{i},i=1,2,3.$

Denote $\hat{M}=\Psi_{1}(u)\times\Psi_{2}(u)\times\Psi_{3}(u)$, since $M$ is
bounded, so is $\hat{M}\subset\mathbb{R}^{d}$. Since $G$ is a polynomial, it
is easy to prove that $\nabla G$ is Lipschitz continuous on any bounded subset
of $\mathbb{R}^{d}$. Therefore, for any $u,w\in\hat{M}$, there is a constant
$c>0$ such that

$\|\nabla G(u)-\nabla G(w)\|_{2}^{2}\leq c\|u-w\|_{2}^{2}$

hence,

$\displaystyle\|\tilde{u}^{k}-\nabla G(u^{k+1})\|_{2}=$
$\displaystyle\left[\|\nabla_{u_{1}}G(u_{1}^{k+1},u_{2}^{k},u_{3}^{k})-\nabla_{u_{1}}G(u^{k+1})\|_{2}^{2}\right.$
$\displaystyle\left.+\|\nabla_{u_{2}}G(u_{1}^{k+1},u_{2}^{k+1},u_{3}^{k})-\nabla_{u_{2}}G(u^{k+1})\|_{2}^{2}\right]^{\frac{1}{2}}$
$\displaystyle\leq$
$\displaystyle\left[c^{2}(\|u_{2}^{k}-u_{2}^{k+1}\|_{2}^{2}+\|u_{3}^{k}-u_{3}^{k+1}\|_{2}^{2})+c^{2}(\|u_{3}^{k}-u_{3}^{k+1}\|_{2}^{2})\right]^{\frac{1}{2}}$
$\displaystyle\leq$
$\displaystyle\sqrt{2}c\|u^{k+1}-u^{k}\|_{2}^{2},k\in\mathbb{N}$

according to (50), there are

$\displaystyle\|w^{k+1}\|_{2}$ $\displaystyle\leq\|\tilde{u}^{k}-\nabla
G(u^{(k+1)})\|_{2}+\rho_{\max}\|u^{k+1}-u^{k}\|_{2}$
$\displaystyle\leq(\sqrt{2}c+\rho_{\max})\|u^{k+1}-u^{k}\|_{2},k\in\mathbb{N}$

where $\rho_{\rm{max}}=\max\\{\rho_{1},\rho_{2},\rho_{3}\\}$. Therefore,
H2(Relative error condition) is satisfied, and $b=\sqrt{2}c+\rho_{\max}$.

In addition, because
$\left\\{u^{k}\right\\}_{k\in\mathbb{N}}\subset\mathbb{R}^{d}$ is bounded, it
is relatively compact, so there exists a subsequence
$\left\\{u^{k_{j}}\right\\}_{j\in\mathbb{N}}$ and $\bar{u}\in\mathbb{R}^{d}$
such that $u^{k_{j}}\rightarrow\bar{u}$, $j\rightarrow+\infty$. Since
$\left\\{u_{3}^{k}|k\in\mathbb{N}\right\\}\subset\tilde{S}$, so
$\left\\{u^{k}|k\in\mathbb{N}\right\\}\subset\hat{S}$ holds. Since $\hat{S}$
is closed, $\bar{u}\in\hat{S}$, and $\tilde{F}$ is continuous on $\hat{S}$.
Therefore,
$\tilde{F}\left(u^{k_{j}}\right)\rightarrow\tilde{F}\left(\bar{u}\right)$ ,
$j\rightarrow+\infty$. Hence, H3(Continuity condition) is satisfied.

Finally, according to Lemma D.3, $\tilde{F}$ satisfies the KŁ property at
$\bar{u}\in\hat{S}$. By Lemma D.2, the sequence $\\{u^{k}\\}_{k\in\mathbb{N}}$
converges to the critical point of $\tilde{F}$. This completes the proof.

## References

  * [1] H. Attouch, J. Bolte, and B. F. Svaiter, Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward–backward splitting, and regularized gauss–seidel methods, Mathematical Programming, 137 (2013), pp. 91–129.
  * [2] J. Bolte, A. Daniilidis, and A. Lewis, The łojasiewicz inequality for nonsmooth subanalytic functions with applications to subgradient dynamical systems, SIAM Journal on Optimization, 17 (2007), pp. 1205–1223.
  * [3] J. Bolte, A. Daniilidis, A. Lewis, and M. Shiota, Clarke subgradients of stratifiable functions, SIAM Journal on Optimization, 18 (2007), pp. 556–572.
  * [4] K. Braman, Third-order tensors as linear operators on a space of matrices, Linear Algebra and Its Applications, 433 (2010), pp. 1241–1253.
  * [5] D. Chen, Y. Chen, and D. Xue, Fractional-order total variation image denoising based on proximity algorithm, Applied Mathematics and Computation, 257 (2015), pp. 537–545.
  * [6] J. Chen and M. K. Ng, Color image inpainting via robust pure quaternion matrix completion: Error bound and weighted loss, SIAM Journal on Imaging Sciences, 15 (2022), pp. 1469–1498.
  * [7] Y. Chen, X. Xiao, and Y. Zhou, Multi-view subspace clustering via simultaneously learning the representation tensor and affinity matrix, Pattern Recognition, 106 (2020).
  * [8] Y.-J. Deng, H.-C. Li, S.-Q. Tan, J. Hou, Q. Du, and A. Plaza, t-linear tensor subspace learning for robust feature extraction of hyperspectral images, IEEE Transactions on Geoscience and Remote Sensing, 61 (2023), pp. 1–15.
  * [9] C. Donciu and M. Temneanu, An alternative method to zero-padded DFT, Measurement, 70 (2015), pp. 14–20.
  * [10] G. H. Golub and C. F. Van Loan, Matrix Computation 4th ed., The Johns Hopkins University Press, Baltimore, USA, 2013.
  * [11] J. Hou, F. Zhang, and J. Wang, One-bit tensor completion via transformed tensor singular value decomposition, Applied Mathematical Modelling, 95 (2021), pp. 760–782.
  * [12] T. X. Jiang, M. K. Ng, X. L. Zhao, and T. Z. Huang, Framelet representation of tensor nuclear norm for third-order tensor completion, IEEE Transaction on Image Processing, 29 (2020), pp. 7233–7244.
  * [13] M. E. Kilmer, K. Braman, N. Hao, and R. 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 (2013), pp. 148–172.
  * [14] M. E. Kilmer and C. D. Martin, Factorization strategies for third-order tensors, Linear Algebra and Its Applications, 435 (2011), pp. 641–658.
  * [15] E. Knerfeld, M. E. Kilmer, and S. Aeron, Tensor-tensor products with invertible linear transforms, Linear Algebra and Its Applications, 485 (2015), pp. 545–570.
  * [16] T. G. Kolda and B. W. Bader, Tensor decompositions and applications, SIAM review, 51 (2009), pp. 455–500.
  * [17] N. Le Bihan and S. Sangwine, Quaternion principal component analysis of color images, IEEE International Conference on Image Processing, 1 (2003), pp. 809–812.
  * [18] B.-Z. Li, X.-L. Zhao, T.-Y. Ji, X.-J. Zhang, and T.-Z. Huang, Nonlinear transform induced tensor nuclear norm for tensor completion, Journal of Scientific Computing, 92 (2022), p. 83.
  * [19] X.-L. Lin, M. K. Ng, and X.-L. Zhao, Tensor factorization with total variation and tikhonov regularization for low-rank tensor completion in imaging data, Journal of Mathematical Imaging and Vision, 62 (2020), pp. 900–918.
  * [20] C. Ling, J. Liu, C. Ouyang, and L. Qi, St-svd factorization and s-diagonal tensors, arXiv:2104.05329, (2021).
  * [21] R. W. Liu, L. Shi, W. Huang, J. Xu, S. C. H. Yu, and D. Wang, Generalized total variation-based mri rician denoising model with spatially adaptive regularization parameters, Magnetic Resonance Imaging, 32 (2014), pp. 702–720.
  * [22] J. Miao, K. I. Kou, and W. Liu, Low-rank quaternion tensor completion for recovering color videos and images, Pattern Recognition, 107 (2020).
  * [23] B. Muquet, Z. Wang, G. B. Giannakis, M. de Courville, and P. Duhamel, Cyclic prefixing or zero padding for wireless multicarrier transmissions?, IEEE Transactions on Communications, 5 (2002), pp. 2136–2148.
  * [24] L. Qi, C. Ling, J. Liu, and C. Ouyang, An orthogonal equivalence theorem for third order tensors, to appear in: Journal of Industrial and Management Optimization, (2021).
  * [25] L. Qi and Z. Luo, Tubal matrices, arXiv:2105.00793v3, (2021).
  * [26] L. Qi, Z. Luo, Q.-W. Wang, and X. Zhang, Quaternion matrix optimization: Motivation and analysis, Journal of Optimization Theory and Applications, 193 (2022), pp. 621–648.
  * [27] D. Qiu, M. Bai, M. K. Ng, and X. Zhang, Nonlocal robust tensor recovery with nonconvex regularization, Inverse Problems, 37 (2021), p. 035001.
  * [28] D. Qiu, M. Bai, M. K. Ng, and X. Zhang, Robust low transformed multi-rank tensor methods for image alignment, Journal of Scientific Computing, 87 (2021), p. 24.
  * [29] O. Semerci, N. Hao, M. E. Kilmer, and E. L. Miller, Tensor-based formulation and nuclear norm regularization for multienergy computed tomography, IEEE Transactions on Image Processing, 23 (2014), pp. 1678–1693.
  * [30] Q. Shi, Y.-M. Cheung, and J. Lou, Robust tensor svd and recovery with rank estimation, IEEE Transactions on Cybernetics, 52 (2021), pp. 10667–10682.
  * [31] G. Song, M. K. Ng, and X. Zhang, Robust tensor completion using transformed tensor singular value decomposition, Numerical Linear Algebra with Applications, 27 (2020).
  * [32] D. A. Tarzanagh and G. Michailidis, Fast randomized algorithms for t-product based tensor operations and decompositions with applications to imaging data, SIAM Journal on Imaging Sciences, 11 (2018), pp. 2629–2664.
  * [33] X. Wang, L. Gu, H. W. Lee, and G. Zhang, Quantumn context-aware recommendation systems based on tensor singular value decomposition, Quantum Information Processing, 20 (2021).
  * [34] X. Xiao, Y. Chen, Y.-J. Gong, and Y. Zhou, Prior knowledge regularized multiview self-representation and its applications, IEEE Transactions on neural networks and learning systems, 32 (2020), pp. 1325–1338.
  * [35] X. Xiao, Y. Chen, Y. J. Gong, and Y. Zhou, Low-rank reserving t-linear projection for robust image feature extraction, IEEE Transactions on Image Processing, 30 (2021), pp. 108–120.
  * [36] Y. Xu, R. Hao, W. Yin, and Z. Su, Parallel matrix factorization for low-rank tensor completion, arXiv preprint arXiv:1312.1254, (2013).
  * [37] L. Yang, Z. H. Huang, S. Hu, and J. Han, An iterative algorithm for third-order tensor multi-rank minimization, Computational Optimization and Applications, 63 (2016), pp. 169–202.
  * [38] Q. Yu and X. Zhang, T-product factorization based method for matrix and tensor completion problems, Computational Optimization and Applications, (2022), pp. 1–28.
  * [39] Q. Yu, X. Zhang, and Z.-H. Huang, Multi-tubal rank of third order tensor and related low rank tensor completion problem, arXiv preprint arXiv:2012.05065, (2020).
  * [40] 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).
  * [41] Z. Zhang and S. Aeron, Exact tensor completion using t-SVD, IEEE Transactions on Signal Processing, 65 (2017), pp. 1511–1526.
  * [42] Z. Zhang, G. Ely, S. Aeron, N. Hao, and M. E. 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, ser. CVPR ’14, IEEE, 2014, pp. 3842–3849.
  * [43] X. Zhao, M. Bai, D. Sun, and L. Zheng, Robust tensor completion: Equivalent surrogates, error bounds, and algorithms, SIAM Journal on Imaging Sciences, 15 (2022), pp. 625–669.
  * [44] Y.-B. Zheng, T.-Z. Huang, X.-L. Zhao, T.-X. Jiang, T.-Y. Ji, and T.-H. Ma, Tensor n-tubal rank and its convex relaxation for low-rank tensor recovery, Information Sciences, 532 (2020), pp. 170–189.
  * [45] Y.-B. Zheng, T.-Z. Huang, X.-L. Zhao, T.-X. Jiang, T.-H. Ma, and T.-Y. Ji, Mixed noise removal in hyperspectral image via low-fibered-rank regularization, IEEE Transactions on Geoscience and Remote Sensing, 58 (2019), pp. 734–749.
  * [46] P. Zhou, C. Lu, Z. Lin, and C. Zhang, Tensor factorization for low-rank tensor completion, IEEE Transactions on Image Processing, 27 (2018), pp. 1152–1163.
  * [47] Y. Zhou and Y. Cheung, Bayesian low-tubal-rank robust tensor factorization with multi-rank determination, IEEE Transactions on Pattern Analysis and Machine Intelligence, 43 (2021), pp. 62–76.


# Mapping the complete evolution of magnetic excitation in beam-plasma system
driven by an ultra-intense, femtosecond laser

Moniruzzaman Shaikh  Tata Institute of Fundamental Research, 1 Homi Bhabha
Road, Mumbai 400005, India  Amit D Lad  Tata Institute of Fundamental
Research, 1 Homi Bhabha Road, Mumbai 400005, India  Devshree Mandal  Institute
of Plasma Research, Gandhinagar, Bhat, Gandhinagar 382428 India  Kamalesh Jana
Tata Institute of Fundamental Research, 1 Homi Bhabha Road, Mumbai 400005,
India  Deep Sarkar  Tata Institute of Fundamental Research, 1 Homi Bhabha
Road, Mumbai 400005, India  Amita Das  Department of Physics, Indian Institute
of Technology Delhi, Hauz Khas, New Delhi 110016, India  G Ravindra Kumar
<EMAIL_ADDRESS>Tata Institute of Fundamental Research, 1 Homi Bhabha Road,
Mumbai 400005, India

###### Abstract

Plasmas are beset with instabilities of all types, hydrodynamic, magneto-
hydrodynamic, and electromagnetic. These instabilities are complex, occur over
a large range of temporal and spatial scales, are most often unmanageable, and
have seriously challenged our efforts at applications, even as they have shed
light on the understanding of the physics of plasmas in the laboratory and
astrophysical environments. A major reason for our limited success in their
containment is the lack of direct experimental information on their origins
and evolution, both temporal and spatial. In plasmas produced by high-
intensity, short, and ultrashort pulse lasers, our knowledge of the
instability stems from the (secondary) signals they generate e.g. scattering
of electromagnetic waves in the form of Raman or Brillouin scattering. Rarely,
if ever, has a direct measurement been made of the instantaneous evolution of
the instabilities in plasmas. In this paper, we present direct measurements of
the femtosecond evolution of the electromagnetic beam-driven instability that
arises from the interaction of forward and return currents in an ultrahigh-
intensity laser-produced plasma on a solid target. We track the evolution from
the initial linear stage to the later nonlinear stage, capturing the entire
instability evolution- all in one experiment, via femtosecond pump-probe
measurements of the giant (megagauss) magnetic field created during the
intense laser-matter interaction. Our experimental findings are excellently
supported by particle–in–cell (PIC) simulations. These first results indicate
heretofore unknown rich physics associated with beam-driven instabilities. In
particular, it has been shown that the emission of electromagnetic radiation
triggers the excitation of instability as negative energy modes feed from the
free energy source associated with the beam propagation. This has been known
in the context of gravitational interaction where the emission of
gravitational radiation is responsible for driving certain negative energy
modes in rotating black holes. In the electromagnetic context, this is the
first such demonstration. Femtosecond laser-produced plasmas can serve as
model templates for experimental instability studies across various branches
of physics.

##  I Introduction

The high energy density (HED) plasma created by ultrahigh intense, femtosecond
laser pulses is a test bed for the fundamental understanding of intense light-
matter interactions and a facilitator for tabletop creation of extreme states,
shocks, and giant magnetic fields [1, 2]. They readily enable table-top
simulations of astrophysics, novel particle acceleration schemes, ultra-bright
X-ray sources, and light-driven nuclear physics. While these applications have
aroused numerous, widespread investigations, several basic questions continue
to challenge us. An important and ubiquitous one is about the onset and
evolution of the instabilities in HED plasmas [3]. In particular, the charged
particle beam-plasma instabilities [4] have recently become very important due
to their role in the fast ignition of laser fusion and several astrophysical
processes.

While there have been numerous simulation studies of this class of
instabilities [5, 6, 7, 8], experimental information is rare and highly
wanted. More specifically, the dynamics, including the onset and early
evolution of instabilities in HED plasmas have not been captured in any
experimental study.

It is precisely this question that we address in this paper and present an
important advance. In a single experimental setting, we study the origin, the
linear growth phase, and the nonlinear stage of a crucial, widely prevalent
beam-plasma instability that is very relevant to the transport of relativistic
charged particle beams in supercritical, high-temperature plasmas, and the
cause of giant magnetic fields in HED plasmas [9, 10, 11, 12] and in
astrophysical systems. We accomplish this using a 25 femtosecond, multi-
terawatt table-top laser, and pump-probe spatio-temporal polarimetry [13] of a
HED plasma created on a solid target, achieving micron-scale spatial and
femtosecond temporal resolution simultaneously. Our novel technique involved
examining the evolving plasma at the plasma rear and the probe pulses that are
reflected at the critical density surface which moves across the plasma on a
femtosecond scale, as the plasma evolves, providing us an unprecedented
opportunity to capture the beginning and linear stages of the instability. As
earlier studies [11] show, probing the plasma at the target front captures
only the nonlinear stage. The examination of the plasma from the rear of the
target (as shown in the schematic experimental set-up of FIG. 2) provides a
view of the expanding plasma inside the target. This happens as the energetic
electrons from the target front continuously dig into the target ionizing it
to form plasma and the trailing electrons in the beam, in conjunction with
electrons from the newly formed plasma, initiate beam plasma instabilities.
The probe beam reaching the critical density layer from the rear thus always
samples the newly formed plasma where the instability is initiated
(illustrated in FIG. 2). Plasma density profile at two times $t_{1}$ and
$t_{2}$ have been shown. It can be observed that the probe pulse sent from the
rear always samples the plasma activity in a freshly ionized layer. On the
other hand, the probe beam from the front observes a region of the plasma in
which the instability has evolved since its inception and consequently is in
the nonlinear phase. Thus the probes from the rear provide an unprecedented
opportunity to look at the beginning (linear) stages of the instability. By
employing a probe beam from the rear side of the target, we track the
beginning, the linear growth stage, of the instability in a direct,
unambiguous manner. Typically, probing the front surface reveals information
about the nonlinear phase. The information about the nonlinear phase can also
be extracted by probing the evolution of the magnetic field at the target rear
when the electron beam ceases to ionize and form a fresh plasma layer and the
critical density layer does not get pushed further at the rear.

In this paper, Section II explains a new experiment to study the time
evolution of the beam plasma system, Section III displays the results of
numerical simulations showing how the disturbances change over time while
Section IV summarizes the findings of the study.

Figure 1: Pump-probe polarimetric setup to measure spatiotemporal evolution of
ultra-intense laser-generated mega-gauss magnetic fields at the rear of thin
targets. CL: Collimating lens, BS: Beam-splitter, PP: Parallel-polariser, CP:
Cross-polariser, IM1, and IM2: Imaging lens 1 and 2, respectively. Figure 2:
Schematic of the evolution of plasma density at various time delays, both in
front of and deep inside the solid target. In front of the target, the probe
pulse samples the plasma ahead of the critical surface, which has undergone
hydrodynamic expansion. contrarily, at the rear side, the probe pulse
consistently captures the freshly ionized plasma created by the propagation of
the ionization front. Figure 3: Experimentally measured magnetic field profile
at the rear of 50 $\mu$m thick mylar target at a laser intensity of $3\times
10^{19}W/cm^{2}$. Figure 4: Experimentally measured power spectra at the rear
of 50 $\mu$m thick mylar target at a laser intensity of $3\times
10^{19}W/cm^{2}$. The dashed lines indicate the trend.

##  II Experimental setup

The experiment was performed at the Tata Institute of Fundamental Research
(TIFR, Mumbai) with a 100 TW laser system. The experimental setup is shown in
FIG. 1. A p-polarized laser pulse (800 nm, $25$ femtosecond, and 1.0 Joule)
was focused with an f/3 off-axis parabolic mirror to a 10 $\mu$m spot on the
target at 45o angle of incidence, giving a laser intensity of $3\times
10^{19}W/cm^{2}$ on target. We used thin dielectric plastic (50 $\mu$m thick
mylar) as targets. The high contrast of the laser pulse (picosecond pedestal/
peak intensity $10^{-9}$) ensures the creation of minimal pre-plasma during
irradiation of thin targets. The target was translated by a computer-
controlled X-Y stage to ensure a fresh spot during each laser shot. A small
fraction ($5\%$) of the main interaction pulse was extracted with the help of
a beam-splitter and up-converted to its second harmonic ($2\omega$ = $400nm$)
using a $\beta$-barium borate (BBO) crystal, which was then used as the probe
pulse. The $2\omega$ probe was appropriately attenuated and focused to a 75
$\mu$m diameter spot to image the rear of the plasma. The probe pulse was
precisely delayed with respect to the main interaction laser pulse (precision
of $7$ fs) using a computer-controlled delay stage. Note that the transverse
spot of the probe pulse is more than seven times bigger than the transverse
spot size of the pump pulse at the interaction region so that we could capture
the magnetic fields at longer time delays when the plasma can undergo lateral
expansion[14]. Additionally, we measured the space-integrated evolution of the
magnetic fields by replacing the CCDs with two identical photo-diodes (PDs).
We captured the spatio-temporal dynamics of the giant magnetic fields
generated during laser-matter interaction at relativistic laser intensity
($3\times 10^{19}W/cm^{2}$), where the probe laser pulse was sent at the
target rear to image the magnetic field profiles. Probing the target rear with
visible light (400 nm) is very different for transparent dielectrics (mylar)
since the probe pulse can travel through the transparent target material till
the critical density layer in the plasma. We relied on the Cotton-Mouton
effect to measure the magnetic fields, where a normally incident, linearly
polarized probe pulse gets elliptically polarized in the presence of giant
magnetic fields as it transits the plasma and gets reflected at the critical
surface for the 400nm probe. We measured the induced ellipticity for each
pixel in the captured images of the probe pulse and retrieved the strength of
the magnetic fields [13]. Figure 3 illustrates the transverse magnetic field
profile captured from the rear plasma surface at various times.

To identify the prominence of the spatial scales that appear, we carry out the
spatial Fourier transform and evaluate the power spectra of the magnetic field
at various times. This has been shown in FIG.4. The power spectra are observed
to be high at long scales right from the very beginning. At earlier times
there is a clear dip of spectral power at intermediate scales. This particular
feature is in contrast with the data that has been observed when the front
surface of the target is probed [10, 11, 12] where the reflected probe beam
samples the entire depth of the target, up to the critical layer inside the
plasma. Thus, the spatial domain contains the integrated behavior of the
nonlinearly perturbed medium. On the other hand, when the medium is probed
from the rear, it only samples the freshly disturbed target region and shows
the spectral behavior at a very early stage. Thus the difference between the
spectral behavior when probed from the rear and front surfaces essentially
mirrors the differences of excited magnetic scales at the early and later
phases of evolution. We, therefore, infer that at the very early phase of
evolution, the excitations at the skin depth scales are essentially absent and
they appear only in later phases. It should be noted here that the typical
well-known instabilities like Weibel [3] and filamentation[5, 15] have a
dominant growth rate at the skin depth scales. The absence of these scales at
the earliest phase of evolution, then suggests that at the earliest time,
scales are dominant for much longer than the skin depth scales. It is only
later that the magnetic field power spectra as probed from the rear surface of
the target show the appearance of the skin depth scales. In the next section,
we try to understand this feature with the help of Particle-in-cell (PIC)
simulations.

##  III Simulation

We have used the OSIRIS framework to carry out particle - in - cell (PIC)
simulation [16], [17], [18] for studying the propagation of an energetic
electron beam in the plasma medium. We have chosen a 2-D slab geometry with a
$X-Y$ plane as a simulation domain. The simulation box has been chosen as
$L_{x}=30L$ and $L_{y}$ is $30L$ in the simulations. Here $L(=c/\omega_{pe})$
is the skin depth of the plasma. Here the beam propagation direction is along
$y$. The plasma density $n_{0}=10^{21}cm^{-3}$ is chosen to be homogeneous in
the entire simulation domain. As an initial condition $10\%$ of electrons
within a transverse extent of $14L$ and longitudinal extent of $6L$ are given
a directed velocity of $0.9c$ along $\hat{y}$. The simulation box is divided
into $1500\times 1500$ cells with $8\times 8$ particles per cell and the box
has absorbing boundary conditions in both directions for particles and fields.

Figure 5: Simulated power spectra using $1/\omega_{pe}$ as the unit of time.
The power spectra from t=9.6 to t=12.0 represent the linear phase of the
instability, while t=12.8 to t=14.4 display the nonlinear growth of the beam-
plasma instability. Figure 6: Simulated surface plot of the magnetic field.
The linear region is on the right side of the red line.

As the beam propagates in the plasma it disturbs the plasma medium. The
generation and evolution of the magnetic field around the tip of the
propagating beam are of prime interest to us here, as this is the region that
gets captured by the probe beam sent from the rear of the target in the
experiment. The simulation provides for the spatio-temporal profiles of the
magnetic field spectra in the 2-D plane. The surface plot of the magnetic
field at $\omega_{pe}t=9.84$ has been shown in Fig.6. The red line in the
figure demarcates the region in the left and right domains. This line has been
placed to demarcate the position towards its right being chosen around this
time to capture the linear phase which subsequently provides the data also for
the nonlinear phase as the beam moves with time. The location of this
particular line has been chosen by ensuring that the strength of the magnetic
field in this region is typically two orders of magnitude smaller than the
maximum magnetic field strength at this particular time. The Fourier transform
of the magnetic disturbances as a function of $k_{x}$ over this limited region
in $y$ (corresponding to the right of the demarcating red line) has been
evaluated at various times. The spectra for this region have been obtained
after integrating it over $dk_{y}$ to obtain power spectra as a function of
$k_{x}$ in equation 1.

$Q_{x}=|B_{k}|^{2}(k_{x})=\int_{box}{|B_{k}|^{2}(k_{x},k_{y})dk_{y}}$ (1)

The plot of $Q_{x}$ as a function of $k_{x}$ has been shown in FIG.5. At
earlier times (till $t=12.0$) the behavior of the spectra compares well with
the experimental spectra shown in Fig.4 provided by the rear probe beam. The
experimental spectra shown in FIG.4 show a dip at central scales. This dip is
also evident in the simulation observations. It can be observed that
subsequently, the simulation spectra show a characteristic change, and for
time beyond $t=12.8$, the dip disappears, and power law behavior indicates the
onset of a nonlinear stage in this particular spatial domain is apparent. This
characteristic is similar to the experimental data observed when the probe
pulse is sent from the front and has been reported in earlier publications.

It is also interesting to note that in some other studies when energy is
pumped at long scales a dip in spectra is observed. This happens as the
nonlinear cascade rate is unable to push the energy pumped at long scales to
smaller dissipative scales rapidly enough. As an example, the recent active
fluid turbulence studies presented in FIG. 3 of Mukherjee et al., [19] also
show such a dip at intermediate scales for large negative values of $\alpha$
(which implies energy injection at long scales). Magnetic field turbulence in
a beam-driven plasma medium can also be viewed in a sense as active matter.
The energy injection occurs due to inherent instabilities. In our case, the
energy injection occurs due to a boundary-driven mechanism at long scales [4].
Thus, until the nonlinear cascade rates remain small and have not disbursed
energy to smaller scales and/or other instabilities (such as Weibel and Kelvin
-Helmholtz etc.,) which dump energy at skin depth scale do not appear, one
observes a dip at the skin depth scale.

##  IV Results and Discussion

We have observed both in experiments and simulations that the spectral power
of the magnetic field associated with the initial excitations in the medium
through the propagating electron beam has identical features. For both of
them, (a) the power spectra are maximum at the most extended scales of the
beam width dimension, (b) the spectra during the initial phase do not show a
power law form, and (c) during the initial phase there is a dip in the spectra
around electron skin depth scales. We now try to understand the physical
reason for these characteristics of the spectra. The conventional beam plasma
instabilities like filamentation, Weibel, and two-streams have the highest
growth rate at the skin depth scale. Thus the earliest phase that has been
captured in experiments with the probe beam entering from the rear side does
not include the excitations of these modes. Yet the power level at much longer
scales is present during this phase. This can only be understood provided
there is a mechanism of long-scale magnetic field excitations present in the
system. The possibility of such an excitation was discussed theoretically and
demonstrated through simulations by Das et al. [4]. It was shown that the
emission of electromagnetic waves from the boundaries results in the forming
of long-scale (of the size of the beam width) magnetic field structures. These
are essentially the negative energy modes excited by dissipation occurring
from the emission of electromagnetic waves from the boundaries. The
experiments here thus provide clear evidence of such a mechanism being
operative in the context of laser-plasma experiments.

We would also like to point out that in the context of gravitational
interactions, such a mechanism of gravitational wave emission leading to the
excitation of negative energy modes from rotating black holes has already been
described theoretically in detail by Chandrashekar [20, 21, 22]. We are
demonstrating the presence of a similar mechanism for electromagnetic
interaction. For the gravitational case, it is not possible to identify this
mechanism with observational data as yet. We on the other hand are showing the
presence of such a mechanism for electromagnetic interactions clearly through
experiments and PIC simulations. Our experiments are thus a demonstration of
yet another possibility wherein laser-plasma experiments transcend their
limited scope to become relevant to other branches of science.

##  V Concluding remarks

The significant accomplishments of this study are as follows. It identifies
novel diagnostics to fetch the spatial spectra of the magnetic field
disturbances created by an electron beam plasma system during the earliest
temporal phase of evolution. The characteristic features of the magnetic field
spectra are quite distinct from the power law behavior. This represents the
early excitation of magnetic field disturbances which have not yet entered any
nonlinear interaction regime. It is noteworthy that even though the magnetic
field has significant power at long beam width scales, there is a clear
absence of spectral disturbances around skin depth scales. This demonstrates
that the conventional beam plasma instabilities are operative only at later
times. At the beginning of time scales, some other mechanisms that generate
scales much greater than the skin depth are present. Earlier [4], a mechanism
was identified for the generation of these long-scale magnetic field
structures. Its presence in experimental data was shown from the probe
diagnostics involving the front surface. However, its appearance before the
conventional instabilities was demonstrated [4] only through PIC simulations.
Here, on the other hand, the onset of such a mechanism at the very beginning
has been identified clearly through experiments. This reinforces our inference
about the existence of such a mechanism and its appearance much ahead of the
conventional modes often discussed in the context of beam plasma interactions.

It was identified by Das et al [4], that the radiative emission of
electromagnetic waves acts as a dissipation to drive the negative energy modes
forming long-scale magnetic structure from the beam energy. In the context of
rotating black holes it has been shown by Chandrashekar and others [20, 21,
22] that the loss of energy through gravitational radiation acts as a
dissipation and leads to the excitation of negative energy modes. Our work,
therefore, shows for the first time the electromagnetic counterpart of
radiative dissipation leading to the excitation of negative energy modes.

###### Acknowledgements.

G.R.K. acknowledges partial support from J. C. Bose Fellowship Grant
No.JCB-37/2010 of the Department of Science and Technology, Ministry of
Science and Technology, Government of India. A.D. acknowledges support from
the Science and Engineering Research Board (SERB), Department of Science and
Technology (DST) Government of India, and the J. C. Bose Fellowship grant
JCB/2017/000055 and CRG proposals CRG/2018/000624 and CRG/2022/002782. A.D.L.
acknowledges partial support from the Infosys-TIFR Leading Edge Research Grant
(Cycle 2). A.D and D.M would like to acknowledge the OSIRIS Consortium,
consisting of UCLA and IST (Lisbon, Portugal) for providing access to the
OSIRIS framework which is the work supported by NSF ACI-1339893. AD also
acknowledges help from Laxman Prasad Goswami for the re-affirmation of
simulation using various platforms. We thank Prof. K. Subramanian (IUCAA,
Pune) and Prof. Dipankar Bhattacharya (RRI, Bangalore) for insightful
discussions on the negative energy modes from black holes. We thank Ankit
Dulat for his help with the manuscript.

## References

  * Kaw [2017] P. K. Kaw, Nonlinear laser–plasma interaction, Rev. Mod. Plasma Phys. 1 (2017).
  * Das [2020] A. Das, Laser–plasma interaction:aapps-dpp conference, 12- 17 nov 2018, kanazawa, Rev. Mod. Plasma Phys. 4 (2020).
  * Weibel [1959] E. S. Weibel, Spontaneously growing transverse waves in a plasma due to an anisotropic velocity distribution., Phys. Rev. Lett. 2, 83 (1959).
  * Das  _et al._ [2020] A. Das, A. Kumar, C. Shukla, R. K. Bera, D. Verma, D. Mandal, A. Vashishta, B. Patel, Y. Hayashi, K. A. Tanaka, G. Chatterjee, A. D. Lad, G. R. Kumar, and P. Kaw, Boundary driven unconventional mechanism of macroscopic magnetic field generation in beam-plasma interaction, Phys. Rev. Res. 2, 033405 (2020).
  * Bret  _et al._ [2005] A. Bret, M. C. Firpo, and C. Deutsch, Characterization of the initial filamentation of a relativistic electron beam passing through a plasma., Phys. Rev. Lett. 94, 115002 (2005).
  * Remington  _et al._ [2014] B. Remington, D. Arnett, R. Paul, and H. Takabe, Modeling astrophysical phenomena in the laboratory with intense lasers., Science 284, 1488 (2014), https://doi.org/10.1126/science.284.5419.1488 .
  * Califano  _et al._ [1997] F. Califano, F. Pegoraro, and S. V. Bulanov, Spatial structure and time evolution of the weibel instability in collisionless inhomogeneous plasmas., Phys. Rev. E 56, 963 (1997).
  * Califano  _et al._ [1998] F. Califano, R. Prandi, F. Pegoraro, and S. V. Bulanov, Nonlinear filamentation instability driven by an inhomogeneous current in a collisionless plasma., Phys. Rev. E 58, 7837 (1998).
  * Sandhu  _et al._ [2002] A. S. Sandhu, A. K. Dharmadhikari, P. P. Rajeev, G. R. Kumar, S. Sengupta, A. Das, and P. K. Kaw, Laser-generated ultrashort multimegagauss magnetic pulses in plasmas., Phys. Rev. Lett. 89, 225002 (2002).
  * Mondal  _et al._ [2012] S. Mondal, V. Narayanan, A. D. Ding, W. J.and Lad, B. Hao, S. Ahmad, W. M. Wang, Z. M. Sheng, S. Sengupta, P. Kaw, A. Das, and G. R. Kumar, Direct observation of turbulent magnetic fields in hot, dense laser produced plasmas., Proc. Nat. Ac. Sci. 109, 8011 (2012), https://www.pnas.org/doi/pdf/10.1073/pnas.1200753109 .
  * Chatterjee  _et al._ [2017] G. Chatterjee, K. M. Schoeffler, P. K. Singh, A. Adak, A. D. Lad, S. Sengupta, P. Kaw, L. O. Silva, A. Das, and G. R. Kumar, Magnetic turbulence in a table-top laser-plasma relevant to astrophysical scenarios., Nature Communications 8, 15970 (2017).
  * Shaikh  _et al._ [2016] M. Shaikh, A. D. Lad, K. Jana, D. Sarkar, I. Dey, and G. R. Kumar, Megagauss magnetic fields in ultra-intense laser generated dense plasmas., Plasma. Phys. Control. Fusion. 59, 014007 (2016).
  * Chatterjee  _et al._ [2014] G. Chatterjee, P. K. Singh, A. Adak, A. D. Lad, and G. R. Kumar, High-resolution measurements of the spatial and temporal evolution of megagauss magnetic fields created in intense short-pulse laser-plasma interactions., Rev. Sci. Instrum. 85, 013505 (2014).
  * Singh  _et al._ [2013] P. K. Singh, Y. Q. Cui, G. Chatterjee, A. Adak, W. M. Wang, S. Ahmed, A. D. Lad, Z. M. Sheng, and G. R. Kumar, Direct observation of ultrafast surface transport of laser-driven fast electrons in a solid target., Phys. Plasmas. 20, 110701 (2013).
  * Wharton  _et al._ [1998] K. B. Wharton, S. P. Hatchett, S. C. Wilks, M. H. Key, J. D. Moody, V. Yanovsky, A. A. Offenberger, B. A. Hammel, M. D. Perry, and C. Joshi, Experimental measurements of hot electrons generated by ultraintense laser-plasma interactions on solid-density targets., Phys. Rev. Lett. 81, 822 (1998).
  * Hemker [2000] R. G. Hemker, Particle-in-cell modeling of plasma-based accelerators in two and three dimensions., Thesis, University of California, Los Angeles https://doi.org/10.48550/arXiv.1503.00276 (2000).
  * Fonseca  _et al._ [2008] R. A. Fonseca, S. F. Martins, L. O. Silva, J. W. Tonge, F. S. Tsung, and W. Mori, One-to-one direct modeling of experiments and astrophysical scenarios: pushing the envelope on kinetic plasma simulations., Plasma Phys. Control. Fusion. 50, 124034 (2008).
  * Fonseca  _et al._ [2002] R. A. Fonseca, L. O. Silva, F. S. Tsung, V. K. Decyk, W. Lu, C. Ren, W. B. Mori, S. Deng, S. Lee, T. Katsouleas, and J. C. Adam,  _OSIRIS: A Three-Dimensional, Fully Relativistic Particle in Cell Code for Modeling Plasma Based Accelerators._ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2002) pp. 342–351.
  * Mukherjee  _et al._ [2023] S. Mukherjee, R. Singh, and M. James, Fluctuations and maximal chaos in an emergent universal state of active turbulence, Nat. Phys. 19, 891–897 (2023).
  * Chandrasekhar [1970] S. Chandrasekhar, Solutions of two problems in the theory of gravitational radiation, Phys. Rev. Lett. 24, 611 (1970).
  * Friedman and Schutz [1978] J. L. Friedman and B. F. Schutz, Secular instability of rotating newtonian stars., Astrophysical Journal 222, 281 (1978).
  * Ballai  _et al._ [2020] I. Ballai, B. Pintér, R. Oliver, and M. Alexandrou, Dissipative instability in a partially ionised prominence plasma slab., A and A 2, 033405 (2020).