1]School of Mathematics and Statistics, Northwestern Polytechnical University,
Xi’an, 710021, Shaanxi, China

2]School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an,
710049, Shaanxi, China

3]School of Electronic and Information Engineering, Xi’an Jiaotong University,
Xi’an, 710049, Shaanxi, China

# Haar Nuclear Norms with Applications to Remote Sensing Imagery Restoration

Shuang Xu  Chang Yu  Jiangjun Peng  Xiangyong Cao  [  [  [

###### Abstract

Remote sensing image restoration aims to reconstruct missing or corrupted
areas within images. To date, low-rank based models have garnered significant
interest in this field. This paper proposes a novel low-rank regularization
term, named the Haar nuclear norm (HNN), for efficient and effective remote
sensing image restoration. It leverages the low-rank properties of wavelet
coefficients derived from the 2-D frontal slice-wise Haar discrete wavelet
transform, effectively modeling the low-rank prior for separated coarse-
grained structure and fine-grained textures in the image. Experimental
evaluations conducted on hyperspectral image inpainting, multi-temporal image
cloud removal, and hyperspectral image denoising have revealed the HNN’s
potential. Typically, HNN achieves a performance improvement of 1-4 dB and a
speedup of 10-28x compared to some state-of-the-art methods (e.g., tensor
correlated total variation, and fully-connected tensor network) for inpainting
tasks.

###### keywords:

Nuclear norm, tensor nuclear norm, tensor completion, hyperspectral image
denoising, remote sensing image cloud removal.

##  1 Introduction

Remote sensing imagery restoration aims to reconstruct missing or corrupted
areas of images [12], including hyperspectral images (HSIs) and multi-temporal
images (MTIs), typically represented as multi-dimensional arrays. Inpainting
and denoising are key tasks in this restoration process [28, 39]. The
mathematical formulations for these tasks are as follows:

$\min_{\mathcal{X}}R(\mathcal{X})\quad\text{s.
t.}\quad\mathscr{P}_{\Omega}(\mathcal{M})=\mathscr{P}_{\Omega}(\mathcal{X}),$
(1)

and

$\min_{\mathcal{X}}R(\mathcal{X})\quad\text{s.
t.}\quad\mathcal{M}=\mathcal{X}+\mathcal{E},$ (2)

Here $\mathcal{E}$ represents noise, $\mathcal{M}$ is the corrupted image, and
$\mathcal{X}$ is the recovered one. Here, $\Omega$ indicates the set of
indices for observable pixels, $\mathscr{P}_{\Omega}(\cdot)$ is the orthogonal
projection onto the observable pixels, and $R(\cdot)$ is the regularization
term applied to the recovered images to ensure desired properties such as
smoothness or low-rankness.

This problem is typically addressed by low-rank (LR) models, which exploits
the intrinsic structure and correlations in the data, under the assumption
that data lie in a low-dimensional subspace [50, 30]. These images tend to
exhibit spatial and channel redundancies, making LR models particularly
effective. The nuclear norm (NN) and its variants, such as the tensor nuclear
norm (TNN) [19, 20], are commonly used as LR regularization. TNN, for
instance, is extensively applied in various tensor processing problems and is
defined as the sum of singular values from the first frontal slice of the
singular value tensor, derived from tensor singular value decomposition (TSVD)
[14]. Moreover, adjacent pixels in images often have similar values, leading
to gradient sparsity. This property gives rise to a new optimization problem:

$\min_{\mathcal{X}}R(\mathcal{X})+\lambda\|\nabla\mathcal{X}\|_{1},\quad\text{s.
t.}\quad\mathscr{P}_{\Omega}(\mathcal{M})=\mathscr{P}_{\Omega}(\mathcal{X}),$
(3)

where $\|\nabla\mathcal{X}\|_{1}$ represents the anisotropic total variation
(TV) [10], which is the $\ell_{1}$ norm of the gradient map
$\nabla\mathcal{X}$, and $\lambda>0$ is the tuning parameter that balances the
regularization terms. TV improves image restoration but is sensitive to the
$\lambda$ parameter selection, limiting its practicality.

The recently proposed correlated TV (CTV) [24] and tensor CTV (TCTV) [31]
address this issue by imposing NN and TNN on the gradient map. For instance,
TCTV, defined as the sum of TNNs of gradients along each mode, inherits the LR
property from the original tensor and promotes gradient sparsity. In essence,
TCTV simultaneously encourages LR and smoothness priors without parameter
tuning. Most significantly, minimizing TCTV can achieve exact data recovery
with high probability. In contrast to NN and TNN, which models LR prior in the
image domain where low- and high-frequencies are entangled, the efficacy of
CTV and TCTV stems from its ability to model the high-frequency LR prior,
specifically the low-rankness of gradient maps.

However, TCTV and CTV face certain challenges. Firstly, its computation is
time-consuming, particularly for high-order images with a significantly larger
number of channels than RGB images. For instance, processing an HSI of size
$256\times 256\times 193$ can take TCTV approximately 2000 seconds, which is
10 times longer than the processing time for TNN. Secondly, while TCTV and CTV
effectively models the high-frequency LR prior, neglecting the low-frequency
information may result in the loss of crucial data and the color/spectral
distortion. This study aims to address these critical issues, and its main
contributions are outlined as follows:

(1) Drawing inspiration from the Haar wavelet transform (HWT), we propose a
two-dimensional frontal HWT tailored for high-order images to disentangle low-
and high-frequencies. Furthermore, we introduce a novel regularizer, the Haar
nuclear norm (HNN), which models the LR property for both disentangled low-
and high-frequency information, thereby combining LR and smoothness priors at
a reduced scale.

(2) The paper establishes the theory of exact recovery for HNN-regularized
inpainting and denoising models, providing theoretical guarantees for the
effectiveness of the proposed method.

(3) The paper develops an algorithm based on the Alternating Direction Method
of Multipliers (ADMM) for solving the HNN-based restoration models. This
algorithm ensures efficient computation, making HNN practical for real-world
applications. For example, in inpainting experiments, HNN achieves a 1-4dB
improvement in performance over state-of-the-art methods with a speedup of
10-28x.

##  2 Related work

###  2.1 Low-rank models

Low-rank models can be broadly classified into two categories: matrix/tensor
decompositions and low-rank regularization techniques.

####  2.1.1 Low-rank decomposition models

Matrix decompositions, with SVD as a prime example, are fundamental in low-
rank models. SVD represents a matrix as a product of three matrices, enabling
it to be expressed as a sum of rank-one matrices, which is key to many low-
rank methods. For high-order images, tensor decompositions extend this
concept, providing a more suitable framework to capture complex data
relationships.

The Canonical Polyadic (CP) decomposition [11] factorizes a tensor into a sum
of rank-one tensors, similar to how SVD works for matrices. The Tucker
decomposition [29] offers a more flexible representation by expressing a
tensor as a core tensor multiplied by factor matrices along each mode,
capturing varying correlations. Tensor Train (TT) [21, 15] and Tensor Ring
(TR) [40, 18] decompositions represent tensors as chains or rings of third-
order tensors, effectively completing images. The recently proposed Fully
Connected Tensor Network (FCTN) [48, 49, 41] decomposes a tensor into a
sequence of small tensors of the same order, with these tensors interacting to
model the relationships between every pair of modes. FCTN holds promise for
superior reconstruction of missing pixels [42, 43] but may require careful
rank selection. Recently, an architecture searching strategy has enhanced
FCTN’s data adaptability [49]. Besides exploring advanced decomposition
models, the noise modeling techniques, such as non-independent and identically
distributed (non-i.i.d.) [3, 34] or asymmetric [35, 13] noise models, have
shown promising results for hyperspectral image processing.

####  2.1.2 Low-rank regularizations

LR regularization encourages solutions with a low-rank structure. For
matrices, the NN, the sum of singular values, is a convex relaxation of rank
and is used in matrices. Nonetheless, the application of LR regularization to
tensors remains an unresolved issue. The Sum of Nuclear Norm (SNN) [17] is the
sum of the nuclear norms of all possible matrix representations
(matricizations) of a tensor, aiming for low ranks in all two-dimensional
unfoldings, thus implicitly encouraging a low-rank structure for the entire
tensor. Building upon the recently proposed tensor-tensor product, the TNN
[19, 20] is derived and demonstrated to be the convex envelope of the tensor
average rank within the unit ball of the tensor spectral norm. In contrast to
SNN, TNN provides theoretical guarantees of recovery. Recently, some more
advanced norms for tensors are developed for denoising and inpainting [46,
45], such as the non-convex [6, 22, 5] and non-local [4, 37, 36] variants.

###  2.2 Low-rank models combined with smoothness

In addition to the LR prior, integrating local smoothness priors or
regularizations [44, 32] can greatly improve image reconstruction. Models such
as TV-regularized SNN (SNNTV) [16], TV-regularized TNN (TNNTV) [27], and
Smooth PARAFAC Completion with TV (SPCTV) [38] are commonly used. However,
these methods are sensitive to the TV weight, often requiring extensive
parameter tuning [25]. CTV and TCTV address this by modeling LR priors in the
gradient domain, promoting both LR and smooth recovery in a regularization.
However, TCTV is computationally intensive, taking over 50 minutes for a
$512\times 512\times 63$ image. Furthermore, CTV and TCTV focus on high-
frequency LR information, potentially overlooking important LR details in the
image domain.

##  3 Preliminaries

###  3.1 Tensor operations

The Frobenius norm of a tensor
$\mathcal{A}\in\mathbb{R}^{I_{1}\times\cdots\times I_{N}}$ is defined as
$\|\mathcal{A}\|_{F}=\sqrt{\sum_{i_{1},\cdots,i_{N}}|a_{i_{1}\cdots
i_{N}}|^{2}}$. The mode-$n$ unfolding of a tensor results in a matrix
representation, denoted by $\bm{A}_{(n)}\in\mathbb{R}^{I_{n}\times\prod_{i\neq
n}I_{i}}$. The mode-$n$ product of a tensor
$\mathcal{A}\in\mathbb{R}^{I_{1}\times\cdots\times I_{N}}$ and a matrix
$\bm{B}\in\mathbb{R}^{J\times I_{n}}$ results in a new tensor
$\mathcal{C}=\mathcal{A}\times_{n}\bm{B}\in\mathbb{R}^{I_{1}\times\cdots\times
I_{n-1}\times J\times I_{n+1}\times\cdots I_{N}}$, given by

$c_{i_{1},\cdots,i_{n-1},j,i_{n+1},\cdots,i_{N}}=\sum_{i_{k}=1}^{I_{k}}a_{i_{1},\cdots,i_{n-1},i_{n},i_{n+1},\cdots,i_{N}}b_{ji_{n}}.$
(4)

This can be represented in matrix form as $\bm{C}_{(n)}=\bm{B}\bm{A}_{(n)}.$

######  Definition 1 ($n$-Rank [7]).

The $n$-rank of a tensor $\mathcal{A}$ is the rank of the tensor $\mathcal{A}$
unfolded along mode $n$, given by
$\mathrm{rank}_{n}(\mathcal{A})=\mathrm{rank}\left(\bm{A}_{(n)}\right)$.

######  Definition 2 (Tucker-Rank [8]).

The Tucker-rank is a tuple of $n$-ranks,
$\mathrm{rank}_{T}(\mathcal{A})=\left(\mathrm{rank}_{1}(\mathcal{A}),\cdots,\mathrm{rank}_{N}(\mathcal{A})\right)$.

###  3.2 Haar wavelet transform

######  Definition 3 (1-D Haar Discrete Wavelet Transform).

The 1-D Haar discrete wavelet transform (HWT) of a vector
$\bm{a}\in\mathbb{R}^{N}$ is defined by

$\bm{b}=f(\bm{a})=\bm{W}_{N}\bm{a},$ (5)

where the $N$-order orthogonal projection matrix $\bm{W}_{N}$ is composed of
an $\frac{N}{2}\times N$ block $\bm{H}_{N/2}$ and an $\frac{N}{2}\times N$
block $\bm{G}_{N/2}$, and the 1-D inverse HWT (IHWT) is given by

$\bm{a}=f^{-1}(\bm{b})=\bm{W}_{N}^{-1}\bm{b}=\bm{W}_{N}^{T}\bm{b}.$ (6)

$\bm{H}_{N/2}$ and $\bm{G}_{N/2}$ are generally deemed as the average filter
and discrete gradient operator, respectively; and $\bm{W}_{N}$ is defined by

$\displaystyle\bm{W}_{N}=$
$\displaystyle\left[\begin{array}[]{c}\bm{H}_{N/2}\\\ \bm{G}_{N/2}\\\
\end{array}\right]$ (7) $\displaystyle=$
$\displaystyle\underbrace{\left[\begin{matrix}\frac{\sqrt{2}}{2}&\frac{\sqrt{2}}{2}&0&0&\cdots&0&0\\\
0&0&\frac{\sqrt{2}}{2}&\frac{\sqrt{2}}{2}&\cdots&0&0\\\
\vdots&&&&\ddots&&\vdots\\\
0&0&0&0&\cdots&\frac{\sqrt{2}}{2}&\frac{\sqrt{2}}{2}\\\
\hdashline\frac{\sqrt{2}}{2}&-\frac{\sqrt{2}}{2}&0&0&\cdots&0&0\\\
0&0&\frac{\sqrt{2}}{2}&-\frac{\sqrt{2}}{2}&\cdots&0&0\\\
\vdots&&&&\ddots&&\vdots\\\
0&0&0&0&\cdots&\frac{\sqrt{2}}{2}&-\frac{\sqrt{2}}{2}\\\
\end{matrix}\right]}_{N},$

######  Definition 4 (2-D Haar Discrete Wavelet Transform).

The 2-D Haar discrete wavelet transform (2-D HWT) for a matrix
$\bm{A}\in\mathbb{R}^{M\times N}$ is defined as

$\bm{B}=F(\bm{A})=\bm{W}_{M}\bm{AW}_{N}^{T},$ (8)

where $\bm{W}_{M}$ and $\bm{W}_{N}$ are the $M$-order and $N$-order projection
matrices, respectively. The 2-D inverse HWT (2-D IHWT) is given by

$\bm{A}=F^{-1}(\bm{B})=\bm{W}_{M}^{T}\bm{B}\bm{W}_{N}.$ (9)

2-D HWT processes a matrix $\bm{A}\in\mathbb{R}^{M\times N}$ by applying the
1-D HWT to its row vectors and column vectors separately. Consequently, as
depicted in Eq. (7), $\bm{B}$ can be represented as a block matrix:

$\displaystyle\bm{B}=$ $\displaystyle\bm{W}_{M}\bm{A}\bm{W}_{N}^{T}$ (10)
$\displaystyle=$ $\displaystyle\left[\begin{array}[]{c}\bm{H}_{M/2}\\\
\bm{G}_{M/2}\\\
\end{array}\right]\bm{A}\left[\begin{array}[]{c}\bm{H}_{N/2}\\\
\bm{G}_{N/2}\\\ \end{array}\right]^{T}$ $\displaystyle=$
$\displaystyle\left[\begin{matrix}\bm{H}_{M/2}\bm{A}\bm{H}_{N/2}^{T}&\bm{H}_{M/2}\bm{A}\bm{G}_{N/2}^{T}\\\
\bm{G}_{M/2}\bm{A}\bm{H}_{N/2}^{T}&\bm{G}_{M/2}\bm{A}\bm{G}_{N/2}^{T}\\\
\end{matrix}\right]$ $\displaystyle\triangleq$
$\displaystyle\left[\begin{matrix}F_{1}(\bm{A})&F_{2}(\bm{A})\\\
F_{3}(\bm{A})&F_{4}(\bm{A})\\\
\end{matrix}\right]\triangleq\left[\begin{matrix}\bm{B}_{1}&\bm{B}_{2}\\\
\bm{B}_{3}&\bm{B}_{4}\\\ \end{matrix}\right],$

where $\bm{B}_{1},\bm{B}_{2},\bm{B}_{3},\bm{B}_{4}$ represent the blocks known
as approximation, horizontal, vertical, and diagonal wavelet coefficients,
respectively.

The computational complexity of the 2-D HWT, as expressed by
$\bm{B}=\bm{W}_{M}\bm{A}\bm{W}_{N}^{T}$, is $O(MN(M+N))$, which involves two
instances of matrix multiplication. However, $\bm{W}_{M}$ and $\bm{W}_{N}$ are
highly sparse matrices with distinctive patterns. The wavelet coefficients can
be verified through the following expressions:

$\displaystyle\bm{B}_{1}(i,j)$
$\displaystyle=\frac{a_{2i-1,2j-1}+a_{2i-1,2j}+a_{2i,2j-1}+a_{2i,2j}}{2},$
(11) $\displaystyle\bm{B}_{2}(i,j)$
$\displaystyle=\frac{a_{2i-1,2j}+a_{2i,2j}-a_{2i-1,2j-1}-a_{2i,2j-1}}{2},$
$\displaystyle\bm{B}_{3}(i,j)$
$\displaystyle=\frac{a_{2i,2j-1}+a_{2i,2j}-a_{2i-1,2j-1}-a_{2i-1,2j}}{2},$
$\displaystyle\bm{B}_{4}(i,j)$
$\displaystyle=\frac{a_{2i-1,2j-1}+a_{2i,2j}-a_{2i-1,2j}+a_{2i,2j-1}}{2},$

implicating that each wavelet coefficient is the sum of four $M/2\times N/2$
matrices. Consequently, the computational complexity of the 2-D HWT via this
method is diminished to $O(\frac{3}{4}MN)$, with the 2-D IHWT exhibititing the
same complexity of $O(\frac{3}{4}MN)$ upon analogous considerations.

Figure 1: (a) Typical remote sensing imagery
$\mathcal{A}\in\mathbb{R}^{M\times N\times S}$ and its wavelet coefficients,
which exhibit rapid decay in singular values. $\mathscr{F}(\cdot)$ represents
the 2D FHWT. (b) The CE curves for the original image and wavelet
coefficients. (c) Comparison of TNN, TCTV, and HNN, where $\|\cdot\|_{*}$
denotes the matrix nuclear norm and $\|\cdot\|_{\rm TNN}$ represents the TNN.

##  4 Haar Nuclear Norm

###  4.1 Motivation

The 2-D HWT decomposes an image into an approximation coefficient and three
detail coefficients, and it has achieved significant advancements in image
processing. However, it’s unclear if these coefficients maintain low-rank
properties for low-rank data, which this section explores.

The image considered in this paper is structured as a 3-order tensor,
$\mathcal{A}\in\mathbb{R}^{M\times N\times S}$, with $M\times N$ pixels and
$S$ channels. Our focus is on characterizing spatial contents, which can be
achieved through wavelet decomposition along the first and second modes.
Consequently, the 2-D HWT for tensor images is defined as follows:

######  Definition 5 (2-D Frontal Slice-wise Haar Discrete Wavelet Transform).

Given a 3-order tensor $\mathcal{A}\in\mathbb{R}^{M\times N\times S}$, the 2-D
frontal slice-wise HWT (FHWT) is defined as:

$\mathcal{B}=\mathscr{F}(\mathcal{A})=\mathcal{A}\times_{1}\bm{W}_{M}\times_{2}\bm{W}_{N}.$
(12)

Equivalently,
$\mathcal{B}(:,:,k)=F\left(\mathcal{A}(:,:,k)\right)=\bm{W}_{M}\mathcal{A}(:,:,k)\bm{W}_{N}^{T}.$
That is, 2-D FHWT is equivalent to applying 2-D HWT to each frontal slice of
the tensor. The 2-D inverse FHWT (IFHWT) is given by:

$\mathcal{A}=\mathscr{F}^{-1}(\mathcal{B})=\mathcal{B}\times_{1}\bm{W}_{M}^{T}\times_{2}\bm{W}_{N}^{T}.$
(13)

The wavelet coefficient $\mathcal{B}$ can be represented as a block tensor
comprising four wavelet coefficient tensors
$\mathcal{B}_{i}\in\mathbb{R}^{M/2\times N/2\times S}(i=1,2,3,4)$:

$\displaystyle\mathcal{B}=$
$\displaystyle\mathcal{A}\times_{1}\bm{W}_{M}\times_{2}\bm{W}_{N}$ (14)
$\displaystyle=$
$\displaystyle\mathcal{A}\times_{1}\left[\begin{array}[]{c}\bm{H}_{M/2}\\\
\bm{G}_{M/2}\\\
\end{array}\right]\times_{2}\left[\begin{array}[]{c}\bm{H}_{N/2}\\\
\bm{G}_{N/2}\\\ \end{array}\right]$ $\displaystyle\triangleq$
$\displaystyle\left[\begin{matrix}\mathscr{F}_{1}(\mathcal{A})&\mathscr{F}_{2}(\mathcal{A})\\\
\mathscr{F}_{3}(\mathcal{A})&\mathscr{F}_{4}(\mathcal{A})\\\
\end{matrix}\right]\triangleq\left[\begin{matrix}\mathcal{B}_{1}&\mathcal{B}_{2}\\\
\mathcal{B}_{3}&\mathcal{B}_{4}\\\ \end{matrix}\right].$

In this equation, $\mathscr{F}_{i}(\mathcal{A})$ represents a mapping from
tensor $\mathcal{A}$ to the $i$-th wavelet coefficient tensor. It is
noteworthy that the wavelet coefficients can be considered as the sum of four
$M/2\times N/2\times S$ tensors (as detailed in the analysis near Eq. 11),
thus the computational complexity of 2-D FHWT (or IFHWT) is
$O(\frac{3}{4}MNS)$.

If the original image is LR, its wavelet coefficients are intuitively also LR.
This is supported by Fig. 1(a), showing fast-decaying singular values across
all wavelet coefficients. Furthermore, the following theorem delineates the
low-rank property of the wavelet coefficients $\mathcal{B}_{i}$:

######  Theorem 6 (Low-rank Properties for Wavelet Coefficients).

Consider a 3-order tensor $\mathcal{A}\in\mathbb{R}^{M\times N\times S}$. Let
$[\mathcal{B}_{1},\mathcal{B}_{2};\mathcal{B}_{3},\mathcal{B}_{4}]=\mathrm{FHWT}_{2}(\mathcal{A})$,
and then there is

$\mathrm{rank}_{n}(\mathcal{B}_{i})\leq\mathrm{rank}_{n}(\mathcal{A}),$ (15)

for $n=1,2,3$ and $i=1,2,3,4$.

The proof is provided in the supplementary materials. The theorem states that
the Tucker rank of the wavelet coefficients does not exceed that of the
original data, confirming that wavelet coefficients remain LR if the original
data is LR.

Cumulative energy (CE) curves provide more insights into the singular value
distributions of wavelet coefficients. For a singular value sequence of length
$n$, $(\lambda_{1},\ldots,\lambda_{n})$, the $k$-th CE is calculated as the
ratio of the cumulative sum of the first $k$ singular values to the total sum
of singular values,
$\text{CE}_{k}=\frac{\sum_{i=1}^{k}\lambda_{i}}{\sum_{i=1}^{n}\lambda_{i}}$.
Fig. 1(b) shows that the CE curve for the approximation coefficient
$\mathcal{B}_{1}$ rises faster than those for the detail coefficients
$\mathcal{B}_{i}(i=2,3,4)$. This suggests that the singular values of low-
frequencies are less dispersed, allowing for reconstruction with fewer
principal components. In contrast, the singular values of high-frequencies are
more dispersed, needing more components for reconstruction. Consequently,
coarse-grained structures are more low-rank, while fine-grained details are
less so. The original image, a mix of low- and high-frequencies, has a CE
curve that falls between those of the low- and high-frequencies, indicating it
is of intermediate low-rank.

NN/TNN applies LR regularization to the original image without considering the
distinct singular value distributions of low- and high-frequencies,
potentially leading to insufficient detail preservation. CTV/TCTV focuses on
high-frequencies through gradient operators, ensuring their preservation but
possibly inaccurately modeling the LR prior for coarse structures and risking
color or spectral distortions.

###  4.2 Haar nuclear norm

This paper introduces the Haar Nuclear Norm (HNN) to model the disentangled LR
properties of low- and high-frequency wavelet coefficients.

######  Definition 7.

Given a 3-order tensor $\mathcal{A}\in\mathbb{R}^{M\times N\times S}$, let
$[\mathcal{B}_{1},\mathcal{B}_{2};\mathcal{B}_{3},\mathcal{B}_{4}]=\mathrm{FHWT}_{2}(\mathcal{A})$,
and HNN is defined as the sum of the nuclear norms of the mode-3 unfoldings of
all wavelet coefficients,

$\|\mathcal{A}\|_{\mathrm{HNN}}=\sum_{i=1}^{4}\|\left(\bm{B}_{i}\right)_{(3)}\|_{*}.$
(16)

It is crucial to clarify the distinction between HNN and existing methods.
Firstly, as illustrated in Fig. 1, the existing tools primarily concentrate on
LR recovery either in the original image domain (e.g., TNN) or the transformed
domain (e.g., TCTV). The former is good at recovering coarse-grained structure
but the entangled LR modeling may neglect detail preservation, while the
latter maintains fine-grained details but may lack robust structure recovery.
HNN differs from existing methods by explicitly modeling the LR properties of
disentangled coarse-grained structures (approximation coefficients) and fine-
grained textures (detail coefficients). It offers a balance between structure
recovery and detail preservation, potentially outperforming methods like TNN
and TCTV.

Despite involving the NN of four wavelet coefficients, HNN does not
significantly increase computational burden. The most intensive operations are
SVD and 2-D FHWT/IFHWT, with complexities of $O(MNS^{2})$ and
$O(\frac{3}{4}MNS)$, respectively. Since $\frac{3}{4}MNS\ll MNS^{2}$, the 2-D
FHWT/IFHWT is negligible compared to SVD, making HNN only slightly slower than
the NN in terms of computational complexity. However, TNN and TCTV involve
additional complexities due to the fast Fourier transform, $O(MNS\log S)$ and
$O(MNS\log MNS)$, respectively. Consequently, HNN is expected to be more
computationally efficient than TNN and TCTV.

###  4.3 HNN based matrix completion

The HNN based matrix completion (HNN-MC) problem is formulated as

$\min_{\mathcal{X}}\|\mathcal{X}\|_{\textrm{HNN}},\text{ s.t.
}\mathscr{P}_{\Omega}(\mathcal{M})=\mathscr{P}_{\Omega}(\mathcal{X}).$ (17)

Note that remote sensing imagery is often written in a tensor format By
introducing auxiliary variables $\mathcal{E}$, it can be rewritten as

$\left\\{\begin{aligned}
\min_{\mathcal{X}}&\sum_{i=1}^{4}\|\left(\bm{B}_{i}\right)_{(3)}\|_{*},\\\
\text{ s.t.
}&\mathcal{M}=\mathcal{X}+\mathcal{E},\mathscr{P}_{\Omega}(\mathcal{E})=0,\mathcal{B}_{i}=\mathscr{F}_{i}(\mathcal{X}).\end{aligned}\right.$
(18)

Recall that $\mathscr{F}_{i}(\mathcal{X})$ denotes a mapping from a tensor
$\mathcal{A}$ into the $i$-th wavelet coefficient, defined in Eq. (14).

The ADMM algorithm is then exploited to solve the issue, and the corresponding
augmented Lagrangian function is

$\displaystyle\mathscr{L}=$
$\displaystyle\left(\sum_{i=1}^{4}\|\left(\bm{B}_{i}\right)_{(3)}\|_{*}+\frac{\mu_{b}}{2}\left\|\mathcal{B}_{i}-\mathscr{F}_{i}(\mathcal{X})+\frac{\Gamma_{i}}{\mu_{b}}\right\|_{F}^{2}\right)$
(19)
$\displaystyle+\frac{\mu_{a}}{2}\left\|\mathcal{M}-\mathcal{X}-\mathcal{E}+\frac{\Gamma_{5}}{\mu_{a}}\right\|_{F}^{2},$

with a constraint $\mathscr{P}_{\Omega}(\mathcal{E})=0,$ where $\mu_{a}$ and
$\mu_{b}$ are penalty parameters, and $\Gamma_{i}\in\mathbb{R}^{M/2\times
N/2\times S}(i=1,2,3,4)$ and $\Gamma_{5}\in\mathbb{R}^{M\times N\times S}$ are
Lagrangian multipliers. In the ADMM, each unknown variable is updated
iteratively by fixing others, and the updating equations are discussed as
follows:

Update $\mathcal{E}$: Fixing other variables, the optimization problem of
$\mathcal{E}$ is given by

$\min_{\mathcal{E}}\left\|\mathcal{E}-\left(\mathcal{M}-\mathcal{X}+\frac{\Gamma_{5}}{\mu_{a}}\right)\right\|_{2}^{2},\text{
s.t. }\mathscr{P}_{\Omega}(\mathcal{E})=0,$ (20)

This is a constrained least squares problem, and the solution is given by

$\mathcal{E}=\mathscr{P}_{\Omega^{\perp}}\left(\mathcal{M}-\mathcal{X}+\frac{\Gamma_{5}}{\mu_{a}}\right),$
(21)

where $\Omega^{\perp}$ is the complement set of $\Omega$.

Update $\mathcal{B}_{i}$: Fixing other variables, the optimization problem of
$\mathcal{B}_{i}$ is written as

$\min_{\mathcal{B}_{i}}\|\left(\bm{B}_{i}\right)_{(3)}\|_{*}+\frac{\mu_{b}}{2}\left\|\mathcal{B}_{i}-\mathscr{F}_{i}(\mathcal{A})+\frac{\Gamma_{i}}{\mu_{b}}\right\|_{F}^{2}.$
(22)

The solution is given by singular value thresholding,

$\mathcal{B}_{i}=\mathrm{fold}_{3}\left(\bm{U}\mathscr{G}_{1/\mu_{b}}(\bm{\Sigma})\bm{V}^{T}\right),$
(23)

where
$[\bm{U},\bm{\Sigma},\bm{V}^{T}]\triangleq\mathrm{svd}\left(\mathscr{F}_{i}(\mathcal{A})-\Gamma_{i}/\mu_{b}\right)$,
and $\mathscr{G}_{\gamma}(a)=\mathrm{sign}(a)\max(|a|-\gamma,0)$ denotes the
soft-thresholding function.

Update $\mathcal{X}$: Fixing other variables, the optimization problem of
$\mathcal{X}$ is given by

$\min_{\mathcal{X}}\frac{\mu_{b}}{2}\left\|\mathcal{B}-\mathscr{F}(\mathcal{X})+\frac{\Gamma}{\mu_{b}}\right\|_{F}^{2}+\frac{\mu_{a}}{2}\left\|\mathcal{M}-\mathcal{X}-\mathcal{E}+\frac{\Gamma_{5}}{\mu_{a}}\right\|_{F}^{2}.$
(24)

This is a weighted least squares problem, and it is easy to compute the
gradient of this loss function regarding $\mathcal{X}$,

$\displaystyle\frac{\partial\mathscr{L}(\mathcal{X})}{\partial\mathcal{X}}$
$\displaystyle=$
$\displaystyle\mu_{a}\left(\mathcal{X}+\mathcal{E}-\mathcal{M}-\Gamma_{5}/\mu_{a}\right)+\mu_{b}(\mathcal{X}\times_{1}\bm{W}_{M}\times_{2}\bm{W}_{N}$
$\displaystyle-\mathcal{B}-\frac{\Gamma}{\mu_{b}})\times_{1}\bm{W}_{M}^{T}\times_{2}\bm{W}_{N}^{T}$
$\displaystyle=$
$\displaystyle\mu_{a}\left(\mathcal{X}+\mathcal{E}-\mathcal{M}-\Gamma_{5}/\mu_{a}\right)+\mu_{b}\mathcal{X}-\mathscr{F}^{-1}\left(\mu_{b}\mathcal{B}+\Gamma\right).$

Let the gradient be zero, and the solution is

$\mathcal{X}=\frac{\mu_{a}(\mathcal{M}+\Gamma_{5}/\mu_{a}-\mathcal{E})+\mathscr{F}^{-1}\left(\mu_{b}\mathcal{B}+\Gamma\right)}{\mu_{a}+\mu_{b}}.$
(25)

Except for the aforementioned variables, the multipliers and penalty
parameters are updated by

$\displaystyle\Gamma_{i}$
$\displaystyle=\Gamma_{i}+\mu_{b}\left(\mathcal{B}_{i}-\mathscr{F}_{i}\left(\mathcal{X}\right)\right),\,\,\left(i=1,2,...,4\right)$
(26) $\displaystyle\Gamma_{5}$
$\displaystyle=\Gamma_{5}+\mu_{a}\left(\mathcal{M}-\mathcal{X}-\mathcal{E}\right),\,\mu_{a}=\mu_{a}\rho\,\,,\,\,\mu_{b}=\mu_{b}\rho,$

where $\rho$ is constant greater than 1.

###  4.4 HNN based robust principal component analysis

Given a noisy observation $\mathcal{M}$, it is assumed that $\mathcal{M}$ can
be decomposed a low-rank component $\mathcal{X}$ (i.e., clean data) and a
sparse noise component $\mathcal{S}$. The HNN based robust principal component
analysis (HNN-RPCA) problem tries to recover $\mathcal{X}$, formulated as

$\min_{\mathcal{X}}\|\mathcal{X}\|_{\textrm{HNN}}+\lambda\|\mathcal{E}\|_{1},\text{
s.t. }\mathcal{M}=\mathcal{X}+\mathcal{E}.$ (27)

With the definition of HNN, it is rewritten as

$\left\\{\begin{aligned}
\min_{\mathcal{X}}&\sum_{i=1}^{4}\|\left(\bm{B}_{i}\right)_{(3)}\|_{*}+\lambda\|\mathcal{E}\|_{1},\\\
\text{ s.t.
}&\mathcal{M}=\mathcal{X}+\mathcal{E},\mathcal{B}_{i}=\mathscr{F}_{i}(\mathcal{X}).\end{aligned}\right.$
(28)

Its augmented Lagrangian function is

$\displaystyle\mathscr{L}=$
$\displaystyle\left(\sum_{i=1}^{4}\|\left(\bm{B}_{i}\right)_{(3)}\|_{*}+\frac{\mu_{b}}{2}\left\|\mathcal{B}_{i}-\mathscr{F}_{i}(\mathcal{X})+\frac{\Gamma_{i}}{\mu_{b}}\right\|_{F}^{2}\right)$
(29)
$\displaystyle+\lambda\|\mathcal{E}\|_{1}+\frac{\mu_{a}}{2}\left\|\mathcal{M}-\mathcal{X}-\mathcal{E}+\frac{\Gamma_{5}}{\mu_{a}}\right\|_{F}^{2},$

where $\mu_{a}$ and $\mu_{b}$ are penalty parameters, and
$\Gamma_{i}\in\mathbb{R}^{M/2\times N/2\times S}(i=1,2,3,4)$ and
$\Gamma_{5}\in\mathbb{R}^{M\times N\times S}$ are Lagrangian multipliers. The
solution is very similar to that of HNN based image inpainting, as follows:

$\left\\{\begin{aligned}
\mathcal{E}&=\mathscr{G}_{\lambda/\mu_{a}}\left(\mathcal{M}-\mathcal{X}+\frac{\Gamma_{5}}{\mu_{a}}\right),\\\
\mathcal{B}_{i}&=\mathrm{fold}_{3}\left(\bm{U}\mathscr{G}_{1/\mu_{b}}(\bm{\Sigma})\bm{V}^{T}\right),\\\
\mathcal{X}&=\frac{\mu_{a}(\mathcal{M}+\Gamma_{5}/\mu_{a}-\mathcal{E})+\mathscr{F}^{-1}\left(\mu_{b}\mathcal{B}+\Gamma\right)}{\mu_{a}+\mu_{b}},\\\
\Gamma_{i}&=\Gamma_{i}+\mu_{b}\left(\mathcal{B}_{i}-\mathscr{F}_{i}\left(\mathcal{X}\right)\right),\,\,\left(i=1,2,...,4\right)\\\
\Gamma_{5}&=\Gamma_{5}+\mu_{a}\left(\mathcal{M}-\mathcal{X}-\mathcal{E}\right),\\\
\mu_{a}&=\mu_{a}\rho\,\,,\,\,\mu_{b}=\mu_{b}\rho.\end{aligned}\right.$ (30)

##  5 Theory of Recovery Guarantee

The section develops the theory of recovery guarantee for HNN-RPCA and HNN-TC
models. To ease representation, this section would use the matrix notations.
Firstly, incoherence conditions on the transformed maps after Haar
transformation are defined. Specifically, assuming that each transformed map
$(\bm{B}_{i})_{(3)}\in\mathbb{R}^{MN/4\times S},(i=1,\cdots,4)$ of
$\bm{X}_{0}\in\mathbb{R}^{MN\times S}$ with rank of $r$ has the singular value
decomposition $\bm{U}_{i}\bm{\Sigma}_{i}\bm{V}_{i}^{T}$, where
$\bm{U}_{i}\in\mathbb{R}^{MN\times r}$ and $\bm{V}_{i}\in\mathbb{R}^{S\times
r}$, the incoherence condition with a fixed constant $\mu$ is defined as
follows:

$\max_{k}\left\|\bm{U}_{i}^{T}\hat{e}_{k}\right\|^{2}\leq\frac{\mu
r}{n_{1}},i=1,\cdots,4.$ (31)
$\max_{k}\left\|\bm{V}_{i}^{T}\hat{e}_{k}\right\|^{2}\leq\frac{\mu
r}{n_{2}},i=1,\cdots,4.$ (32)
$\|\bm{U}_{i}\bm{V}_{i}^{T}\|_{\infty}\leq\sqrt{\frac{\mu
r}{n_{1}n_{2}}},i=1,\cdots,4.$ (33)

Here, $\|\bm{X}\|_{\infty}=\max\limits_{i,j}|\bm{X}_{i,j}|$, and $\hat{e}_{k}$
is the standard orthonormal basis. The exact decomposition theorem for HNN-
RPCA model is derived as follows, where we now write it in matrix notations,

$\min_{\bm{X}}\|\bm{X}\|_{\textrm{HNN}}+\lambda\|\bm{E}\|_{1},\textrm{ s.t.
}\bm{M}=\bm{X}+\bm{E}.$ (34)

(a) RPCA (b) MC

Figure 2: Phase transitions of (a) RPCA models and (b) MC models. The success
rates are shown in the title of each panel.  Figure 3: The false-color
denoised images (band 24-17-1) of compared methods on the PC dataset for case
6. Table 1: HSI denoising metrics on the BA dataset. The best and the second
best results are highlighted in bold and underline, respectively.

Case | Metrics | Noisy | LLRT | TDL | LTDL | NGMeet | CTV | 3DTNN | 3DLogTNN | WNLRATV | RCILD | HNN [b]
---|---|---|---|---|---|---|---|---|---|---|---|---
1 | PSNR$\uparrow$ | 13.06 | 15.6 | 22.9 | 23.04 | 23.15 | 27.95 | 26.7 | 24.42 | 28.37 | 23.22 |  31.89 [t]
SSIM$\uparrow$ | 0.0308 | 0.062 | 0.6801 | 0.7085 | 0.7408 | 0.4991 | 0.5886 | 0.3614 | 0.7549 | 0.7381 | 0.7864
ERGAS$\downarrow$ | 1114.99 | 838.17 | 416.67 | 414.47 | 411.94 | 196.53 | 256.8 | 303.19 | 210.78 | 411.8 | 126.5
SAM$\downarrow$ | 0.269 | 0.217 | 0.083 | 0.082 | 0.081 | 0.091 | 0.080 | 0.108 | 0.087 | 0.082 |  0.045 [b]
2 | PSNR$\uparrow$ | 14.3 | 16.63 | 24.18 | 25.66 | 24.38 | 29.17 | 28.15 | 26.36 | 28.84 | 26.13 |  32.16 [t]
SSIM$\uparrow$ | 0.0467 | 0.0847 | 0.5333 | 0.748 | 0.7006 | 0.5677 | 0.6416 | 0.4666 | 0.7405 | 0.7892 | 0.8064
ERGAS$\downarrow$ | 1046.25 | 818.79 | 412.75 | 365.74 | 406.52 | 174.26 | 230.81 | 277.38 | 219.36 | 361.53 | 129.73
SAM$\downarrow$ | 0.266 | 0.221 | 0.099 | 0.086 | 0.100 | 0.080 | 0.076 | 0.106 | 0.092 | 0.085 |  0.047 [b]
3 | PSNR$\uparrow$ | 13.33 | 15.77 | 22.82 | 23.77 | 23.12 | 29.23 | 25.05 | 21.35 | 28.48 | 24.91 |  32.02 [t]
SSIM$\uparrow$ | 0.0408 | 0.0772 | 0.4622 | 0.7111 | 0.6552 | 0.5706 | 0.329 | 0.1923 | 0.7041 | 0.7612 | 0.802
ERGAS$\downarrow$ | 1195.65 | 932.48 | 471.94 | 436.95 | 463.06 | 170.52 | 295.23 | 476.74 | 248.07 | 409.96 | 130.87
SAM$\downarrow$ | 0.278 | 0.232 | 0.108 | 0.096 | 0.107 | 0.081 | 0.119 | 0.178 | 0.106 | 0.094 |  0.049 [b]
4 | PSNR$\uparrow$ | 14.35 | 16.68 | 23.96 | 25.36 | 24.6 | 29.1 | 27.45 | 24.89 | 29.58 | 26.52 |  32.26 [t]
SSIM$\uparrow$ | 0.0473 | 0.0862 | 0.5004 | 0.742 | 0.7169 | 0.564 | 0.5438 | 0.3648 | 0.7839 | 0.7972 | 0.8079
ERGAS$\downarrow$ | 1024.37 | 797.55 | 401.08 | 362.24 | 386.97 | 174.07 | 234.16 | 313.58 | 198.36 | 336.48 | 125.28
SAM$\downarrow$ | 0.265 | 0.220 | 0.099 | 0.084 | 0.096 | 0.080 | 0.084 | 0.122 | 0.084 | 0.082 |  0.046 [b]
5 | PSNR$\uparrow$ | 14.86 | 17.22 | 21.92 | 23.48 | 25.74 | 27.35 | 26.74 | 27.42 | 26.34 | 26.72 |  29.16 [t]
SSIM$\uparrow$ | 0.0525 | 0.0971 | 0.2834 | 0.6506 | 0.7111 | 0.5363 | 0.6656 | 0.6634 | 0.567 | 0.7472 | 0.6868
ERGAS$\downarrow$ | 978.93 | 766.22 | 472.17 | 273.05 | 353.64 | 235.23 | 267.21 | 256.9 | 314.17 | 309.07 | 219.12
SAM$\downarrow$ | 0.276 | 0.233 | 0.154 | 0.084 | 0.111 | 0.123 | 0.108 | 0.110 | 0.148 | 0.095 | 0.105 [b]
6 | PSNR$\uparrow$ | 13.47 | 15.97 | 20.75 | 21.46 | 24.1 | 27.25 | 26.45 | 25.35 | 25.78 | 25.11 |  28.99 [t]
SSIM$\uparrow$ | 0.0417 | 0.0817 | 0.2425 | 0.5546 | 0.6633 | 0.5288 | 0.5819 | 0.4671 | 0.5201 | 0.7145 | 0.67
ERGAS$\downarrow$ | 1154.8 | 892.02 | 524.67 | 332.13 | 422.28 | 234.78 | 260.3 | 298.07 | 368.66 | 380.11 | 224.18
SAM$\downarrow$ | 0.288 | 0.242 | 0.153 | 0.096 | 0.115 | 0.122 | 0.110 | 0.131 | 0.171 | 0.111 |  0.108 [b]

######  Theorem 8 (HNN-RPCA Theorem).

Suppose that the transformed map $(\bm{B}_{i})_{(3)}\in\mathbb{R}^{MN/4\times
S},(i=1,\cdots,4)$ of $\bm{X}_{0}\in\mathbb{R}^{MN\times S}$ with rank $r$
satisfies the incoherence conditions (31)-(33), and the support set $\Omega$
of $\bm{E}_{0}$ is uniformly distributed over all sets of cardinality $m$.
Then, with a probability of at least $1-cn_{(1)}^{-10}$ (for the choice of
support set exceeding $\bm{E}_{0}$), the HNN-RPCA model with
$\lambda=4/\sqrt{n_{(1)}}$ guarantees an exact recovery, i.e., the
minimization in Eq. (34) yields $\hat{\bm{X}}=\bm{X}_{0}$ and
$\hat{\bm{E}}=\bm{E}_{0}$, provided that:

$\operatorname{rank}\left(\bm{X}_{0}\right)\leq\rho_{r}n_{2}\mu^{-1}\left(\log
n_{(1)}\right)^{-2},m\leq\rho_{s}n_{1}n_{2}.$ (35)

Here, $\rho_{r}$ and $\rho_{s}$ are numerical constants, where $\rho_{s}$
determines the sparsity of $\bm{E}_{0}$, $\rho_{r}$ is a small constant
related to the rank of $\bm{X}_{0}$, and $n_{(1)}=\max\\{MN/4,S\\}$.

Similarly, for the HNN-MC model which is written in matrix notations as:

$\min_{\bm{X}}\ \|\bm{X}\|_{\textrm{HNN}},\text{ s.t.
}\mathscr{P}_{\Omega}(\bm{M})=\mathscr{P}_{\Omega}(\bm{X}),$ (36)

there is the following theorem.

######  Theorem 9 (HNN-MC Theorem).

Suppose that $\bm{B}_{i}\in\mathbb{R}^{MN/4\times S},(i=1,\cdots,4)$ obey the
incoherence conditions (31)-(32), $\Omega\sim\mbox{Ber}(p)$ and $m$ is the
number of $\Omega$, where $\mbox{Ber}(p)$ represents the Bernoulli
distribution with parameter $p$. Let $n_{1}=MN/4$ and $n_{2}=S$, and without
loss of generality, suppose $n_{1}\geq n_{2}$. Then, there exist universal
constants $c_{0},c_{1}>0$ such that $\bm{X}_{0}$ is the unique solution to
MNN-MC model (36) with probability at least $1-c_{1}n_{1}^{-3}\log n_{1}$,
provided that

$m\geq c_{0}\mu rn_{1}^{5/4}\log(n_{1}).$ (37)

The proofs of the above theorems are placed in the supplementary material.
These theorems show that the HNN-RPCA/-MC model has the exact recovery with a
high probability.

##  6 Experiments

The proposed model’s effectiveness is tested through extensive experiments
using four quality metrics: Peak Signal-to-Noise Ratio (PSNR), Structural
Similarity (SSIM), Erreur Relative Globale Adimensionnelle de Synthese
(ERGAS), and Spectral Angle Mapper (SAM). PSNR, SSIM, and ERGAS evaluate
spatial distortion, while SAM measures spectral distortion. Better image
quality is indicated by higher PSNR and SSIM scores, and lower ERGAS and SAM
scores. Experiments were conducted on a Windows 11 desktop with an 8-core
R7-4800H CPU and 16GB RAM.

Table 2: HSI denoising metrics on the PC dataset. The best and the second best
results are highlighted in bold and underline, respectively.

Case | Metrics | Noisy | LLRT | TDL | LTDL | NGMeet | CTV | 3DTNN | 3DLogTNN | WNLRATV | RCILD | HNN [b]
---|---|---|---|---|---|---|---|---|---|---|---|---
1 | PSNR$\uparrow$ | 12.54 | 16.96 | 23.35 | 23.04 | 23.92 | 25.95 | 25.12 | 27.04 | 24.02 | 24.04 |  27.59 [t]
SSIM$\uparrow$ | 0.117 | 0.255 | 0.7348 | 0.7085 | 0.7947 | 0.7116 | 0.7288 | 0.8074 | 0.7155 | 0.8061 | 0.8099
ERGAS$\downarrow$ | 865.04 | 521.08 | 255.24 | 414.47 | 241.74 | 185.14 | 204.08 | 165.22 | 233.17 | 238.13 | 152.69
SAM$\downarrow$ | 0.243 | 0.155 | 0.047 | 0.082 | 0.043 | 0.081 | 0.042 | 0.039 | 0.047 | 0.043 | 0.044 [b]
2 | PSNR$\uparrow$ | 14.17 | 17.98 | 22.66 | 25.66 | 24.37 | 27.6 | 27.13 | 28.16 | 27.45 | 25.94 |  28.45 [t]
SSIM$\uparrow$ | 0.1789 | 0.3117 | 0.5653 | 0.748 | 0.7478 | 0.783 | 0.7706 | 0.7953 | 0.8035 | 0.8471 | 0.8452
ERGAS$\downarrow$ | 773.26 | 505.14 | 302.96 | 365.74 | 257.31 | 155.13 | 166.15 | 157.8 | 178.52 | 207.37 | 139.07
SAM$\downarrow$ | 0.234 | 0.163 | 0.091 | 0.086 | 0.073 | 0.070 | 0.054 | 0.062 | 0.056 | 0.049 |  0.044 [b]
3 | PSNR$\uparrow$ | 13.5 | 17.39 | 22.05 | 23.77 | 23.44 | 27.86 | 27.31 | 28.36 | 27.06 | 24.86 |  28.24 [t]
SSIM$\uparrow$ | 0.1609 | 0.2901 | 0.535 | 0.7111 | 0.7168 | 0.794 | 0.7805 | 0.8013 | 0.7725 | 0.8227 | 0.8404
ERGAS$\downarrow$ | 842.82 | 554.27 | 326.73 | 436.95 | 316.19 | 150.55 | 162.93 | 161.15 | 254.39 | 252.96 | 145.68
SAM$\downarrow$ | 0.242 | 0.173 | 0.093 | 0.096 | 0.086 | 0.069 | 0.054 | 0.065 | 0.079 | 0.061 |  0.048 [b]
4 | PSNR$\uparrow$ | 13.12 | 17.28 | 22.08 | 25.36 | 23.7 | 26.45 | 25.17 | 27.03 | 25.25 | 25.73 |  27.76 [t]
SSIM$\uparrow$ | 0.1369 | 0.2721 | 0.5492 | 0.742 | 0.7483 | 0.7374 | 0.7312 | 0.8083 | 0.7293 | 0.8391 | 0.8186
ERGAS$\downarrow$ | 835.9 | 521.79 | 304.17 | 362.24 | 259.7 | 175.11 | 203.36 | 165.42 | 214.67 | 214.09 | 149.84
SAM$\downarrow$ | 0.241 | 0.161 | 0.081 | 0.084 | 0.064 | 0.078 | 0.044 | 0.042 | 0.056 | 0.052 | 0.045 [b]
5 | PSNR$\uparrow$ | 13.52 | 17.6 | 19.86 | 23.55 | 24.02 | 24.13 | 23.91 | 25.11 | 25.16 | 25.6 |  25.43 [t]
SSIM$\uparrow$ | 0.1468 | 0.28 | 0.3895 | 0.5351 | 0.7315 | 0.662 | 0.6469 | 0.6957 | 0.7094 | 0.7664 | 0.7447
ERGAS$\downarrow$ | 813.08 | 510.51 | 399.92 | 406.29 | 255.25 | 263.17 | 266.18 | 252.27 | 242.79 | 227.82 | 241.2
SAM$\downarrow$ | 0.252 | 0.173 | 0.136 | 0.124 | 0.080 | 0.115 | 0.092 | 0.096 | 0.087 | 0.080 | 0.094 [b]
6 | PSNR$\uparrow$ | 12.96 | 17.12 | 19.21 | 22.3 | 23.56 | 24.02 | 23.29 | 24.23 | 23.64 | 23.79 |  24.88 [t]
SSIM$\uparrow$ | 0.1353 | 0.2645 | 0.3528 | 0.472 | 0.7248 | 0.6633 | 0.6065 | 0.6518 | 0.5969 | 0.7171 | 0.7288
ERGAS$\downarrow$ | 875.19 | 552.58 | 433.47 | 460.26 | 287.06 | 270.3 | 280.32 | 271.1 | 287.3 | 259.8 | 260.4
SAM$\downarrow$ | 0.260 | 0.183 | 0.145 | 0.133 | 0.092 | 0.119 | 0.098 | 0.105 | 0.092 | 0.085 | 0.108 [b]

###  6.1 Experiments on phase transition

The recovery phenomenon was analyzed across varying Tucker ranks and noise
intensities of $\bm{E}$ for RPCA problems, or varying sampling rates of
$\bm{M}$ for MC problems. A simulated tensor
$\mathcal{X}\in\mathbb{R}^{30\times 30\times 30}$ was generated using the
equation
$\mathcal{X}=\mathcal{C}\times_{1}\bm{U}_{1}\times_{2}\bm{U}_{2}\times_{3}\bm{U}_{3}$,
where $\mathcal{C}\in\mathbb{R}^{R\times R\times R}$ and
$\bm{U}_{i}\in\mathbb{R}^{30\times R}$, for $i=1,2,3$. The rank $R$ varied
from 1 to 30, the noise intensity ranged from 0.01 to 0.5, and the sampling
rate varied from 0.01 to 0.99. For each case, the experiment was repeated 10
times, and a trial was deemed successful if the average relative
reconstruction error was less than 0.1. Figure 2 illustrates the phase
transition maps (yellow = 100%, blue = 0%), indicating that HNN achieved a
higher success rate than CTV and TCTV.

Table 3: Metrics for HSI inpainting on the BA dataset. The best and the second
best results are highlighted in bold and underline, respectively.

SR | Metrics | Observed | SNN | KBR | TNN | SPCTV | TNNTV | FCTNTC | FCTNFR | TCTV | HNN
---|---|---|---|---|---|---|---|---|---|---|---
7% | PSNR$\uparrow$ | 20.09 | 31.82 |  44.88 | 39.47 | 40.07 | 34.05 | 37.06 | 39.89 | 44.43 | 48.46
SSIM$\uparrow$ | 0.0864 | 0.8254 | 0.9742 | 0.9314 | 0.9492 | 0.8542 | 0.8920 | 0.9385 |  0.9749 | 0.9898
ERGAS$\downarrow$ | 535.81 | 138.69 |  31.90 | 61.77 | 52.10 | 108.97 | 73.48 | 52.66 | 33.26 | 22.30
SAM$\downarrow$ | 1.3075 | 0.1374 |  0.0397 | 0.0857 | 0.0632 | 0.0962 | 0.0859 | 0.0643 | 0.0456 | 0.0288
5% | PSNR$\uparrow$ | 20.00 | 31.10 | 41.46 | 36.83 | 37.67 | 33.15 | 36.89 | 39.40 |  42.40 | 45.55
SSIM$\uparrow$ | 0.0708 | 0.8145 | 0.9501 | 0.8908 | 0.9215 | 0.8412 | 0.8883 | 0.9316 |  0.9616 | 0.9824
ERGAS$\downarrow$ | 541.57 | 150.60 | 47.38 | 85.19 | 69.13 | 121.79 | 74.92 | 55.64 |  43.01 | 30.77
SAM$\downarrow$ | 1.3515 | 0.1489 | 0.0534 | 0.1173 | 0.0802 | 0.1059 | 0.0871 | 0.0669 |  0.0576 | 0.0396
4% | PSNR$\uparrow$ | 19.95 | 30.74 | 39.20 | 35.17 | 36.14 | 32.70 | 36.72 | 39.05 |  39.76 | 43.00
SSIM$\uparrow$ | 0.0623 | 0.8097 | 0.9267 | 0.8587 | 0.8993 | 0.8343 | 0.8852 | 0.9258 |  0.9375 | 0.9688
ERGAS$\downarrow$ | 544.37 | 156.67 | 61.55 | 103.90 | 82.73 | 128.49 | 76.36 | 57.91 |  55.64 | 43.63
SAM$\downarrow$ | 1.3764 | 0.1546 |  0.0665 | 0.1425 | 0.0927 | 0.1125 | 0.0890 | 0.0691 | 0.0722 | 0.0576

Table 4: Metrics for HSI inpainting on the PC dataset. The best and the second
best results are highlighted in bold and underline, respectively.

SR | Metrics | Observed | SNN | KBR | TNN | SPCTV | TNNTV | FCTNTC | FCTNFR | TCTV | HNN
---|---|---|---|---|---|---|---|---|---|---|---
7% | PSNR$\uparrow$ | 14.68 | 21.72 |  36.95 | 31.17 | 30.34 | 22.44 | 28.94 | 33.56 | 32.79 | 38.37
SSIM$\uparrow$ | 0.0380 | 0.3756 |  0.9678 | 0.8822 | 0.8749 | 0.3975 | 0.8165 | 0.9288 | 0.9228 | 0.9799
ERGAS$\downarrow$ | 722.79 | 319.96 |  56.13 | 118.83 | 119.25 | 294.82 | 139.44 | 83.19 | 92.52 | 49.24
SAM$\downarrow$ | 1.3109 | 0.1961 |  0.0806 | 0.1845 | 0.1369 | 0.1757 | 0.1742 | 0.1213 | 0.1513 | 0.0559
5% | PSNR$\uparrow$ | 14.59 | 21.25 | 31.38 | 28.69 | 27.72 | 21.89 | 28.60 |  32.53 | 30.86 | 35.04
SSIM$\uparrow$ | 0.0287 | 0.3399 | 0.8912 | 0.8157 | 0.7863 | 0.3529 | 0.8041 |  0.9114 | 0.8861 | 0.9594
ERGAS$\downarrow$ | 730.54 | 337.91 | 105.51 | 152.20 | 160.69 | 313.87 | 144.82 |  92.98 | 113.63 | 72.15
SAM$\downarrow$ | 1.3557 | 0.1931 |  0.1141 | 0.2003 | 0.1556 | 0.1703 | 0.1777 | 0.1272 | 0.1763 | 0.0717
4% | PSNR$\uparrow$ | 14.54 | 21.04 | 28.58 | 26.91 | 25.97 | 21.67 | 28.31 |  31.68 | 29.75 | 32.10
SSIM$\uparrow$ | 0.0240 | 0.3255 | 0.8091 | 0.7446 | 0.7007 | 0.3344 | 0.7929 |  0.8947 | 0.8578 | 0.9225
ERGAS$\downarrow$ | 734.34 | 346.05 | 145.46 | 183.00 | 196.45 | 322.02 | 149.69 |  102.09 | 128.31 | 103.87
SAM$\downarrow$ | 1.3823 | 0.1896 |  0.1295 | 0.2039 | 0.1656 | 0.1681 | 0.1831 | 0.1358 | 0.1911 | 0.0952

###  6.2 Experiments on HSI denoising

This section carries out HSI denoising experiments to assess the performance
of HNN based RPCA. The compared methods include non-local meets global
(NGMeet) [9], hyper-Laplacian regularized unidirectional low-rank tensor
recovery (LLRT) [1], tensor dictionary learning (TDL) [26], CTV [24], three-
directional TNN (3DTNN) [46], 3D log-based TNN (3DLogTNN) [46], and weighted
non-local low-rank model with adaptive TV regularization (WNLRATV) [2].
Besides them, it is also compared with a deep learning based method,
representative coefficient image with a learnable denoiser (RCILD) [23]. The
experiments simulate six scenarios:

  1. 1.

Case 1: I.i.d. Gaussian noise with $\sigma=75$;

  2. 2.

Case 2: Non-i.i.d. Gaussian noise with band-varying standard deviations
$\sigma\in[30,100]$;

  3. 3.

Case 3: Building on Case 2, 1/3 of the bands are corrupted by impulse noise
with varying ratios $p\in[5\%,20\%]$;

  4. 4.

Case 4: Building on Case 2, 1/3 of the bands are corrupted by stripe noise
with varying ratios $p\in[5\%,20\%]$;

  5. 5.

Case 5: Building on Case 2, 1/3 of the bands are corrupted by deadline noise
with varying ratios $p\in[5\%,20\%]$;

  6. 6.

Case 6: Building on Case 2, the data is corrupted by a mixture of impulse,
stripe, and deadline noise as described in the preceding cases.

Tables 1 and 2 list the HSI denoising metrics on the BA and PC datasets,
respectively. The noise intensity in this experiment configuration is severe,
but HNN still exhibits robust performance. Notably, it surpasses CTV in terms
of PSNR, with an improvement of 1.74 dB on the BA datasets for case 6. Visual
inspection of the results in Fig. 3 reveals that methods including TDL, LTDL,
CTV, 3DTNN and 3DLogTNN exhibit noticeable noise. Although NGMeet and WNLRATV
effectively remove noise, they suffer from spectral distortion, due to the
color mismatch. HNN stands out by generating an image that closely aligns with
the ground truth.

###  6.3 Experiments on HSI inpainting

This subsection focuses on HSI inpainting, and we conduct experiments using
the Bay Area (BA) and Pavia Center (PC) datasets. The BA dataset has
dimensions of $256\times 256\times 193$, while the PC dataset has dimensions
of $200\times 200\times 80$. We compare the HNN method with eight methods,
including SNN [17], TNN [19], TCTV [31], Kronecker-Basis-Representation (KBR)
[33], TNNTV [27], FCTN based tensor completition (FCTNTC) [48], FCTN with
factor regularization (FCTNFR) [47], and SPCTV [38], by varying the sampling
rates (SR) at 4%, 5%, and 7%. Tables 3 and 4 present the performance metrics
on the BA and PC datasets, respectively.

Figure 4: False-color images (band: 193-97-19) of all compared methods on the
BA dataset.

Table 3 shows that the HNN method outperforms other advanced methods in all
evaluation metrics on the BA dataset. In terms of computational efficiency,
the running times for the less effective KBR and TCTV methods are 1310.40s and
2020.76s, respectively, which are significantly higher than the HNN method’s
running time of 177.13s. This indicates that the HNN method not only achieves
best performance but also offers improved efficiency. As depicted in Table 4,
the HNN approach consistently outperforms other methods on the PC dataset. For
example, at a SR of 7%, it is evident that HNN shows a superior PSNR value
compared to KBR and FCTBFR by 1.42 dB and 4.81 dB, respectively.

Fig. 4 depicts the false-color images on the BA dataset with a 5% SR. A
notable observation is that the SNN method simply filled in the missing areas,
while the TNN, SPCTV, FCTNTC, and FCTCFR methods introduced grid-like stripes.
Only the KBR, TCTV, and HNN successfully restored the HSIs in a visually
appealing manner. Among these, the HNN method stands out with its good
results, for example, sharper textures and colors that closely resemble the
original clear images. This confirms the effectiveness of the HNN method in
the task of HSI inpainting.

Table 5:  Metrics on the Jizzakh mountain dataset. The best and the second
best results are highlighted in bold and underline, respectively.

Time mode | Metrics | Obs | SNN | KBR | TNN | SPCTV | TNNTV | FCTNTC | FCTNFR | TCTV | HNN
---|---|---|---|---|---|---|---|---|---|---|---
#1 | PSNR$\uparrow$ | 12.39 | 23.09 | 23.46 | 27.20 | 26.79 | 22.60 | 26.84 | 27.38 |  28.16 | 28.24
SSIM$\uparrow$ | 0.5187 | 0.7961 | 0.8726 | 0.8875 | 0.8717 | 0.7746 | 0.8753 | 0.8916 |  0.9089 | 0.9323
ERGAS$\downarrow$ | 764.33 | 182.84 | 173.29 | 113.07 | 117.77 | 191.43 | 117.39 | 110.69 |  101.99 | 100.84
SAM$\downarrow$ | 0.6122 | 0.0751 | 0.0594 | 0.0593 | 0.0634 | 0.0673 | 0.0619 | 0.0603 |  0.0554 | 0.0343
#2 | PSNR$\uparrow$ | 12.28 | 24.96 | 29.68 | 29.23 | 30.44 | 24.47 | 29.54 |  32.09 | 30.73 | 32.48
SSIM$\uparrow$ | 0.5531 | 0.8299 | 0.9472 | 0.9141 | 0.9176 | 0.8009 | 0.9054 |  0.9505 | 0.9386 | 0.9659
ERGAS$\downarrow$ | 801.33 | 152.17 | 89.77 | 92.49 | 80.35 | 160.63 | 88.65 |  66.09 | 78.44 | 64.46
SAM$\downarrow$ | 0.5067 | 0.0420 |  0.0213 | 0.0332 | 0.0406 | 0.0403 | 0.0335 | 0.0245 | 0.0294 | 0.0156
#3 | PSNR$\uparrow$ | 14.01 | 27.60 | 33.81 | 32.49 | 33.40 | 27.26 | 32.22 |  35.21 | 34.13 | 36.23
SSIM$\uparrow$ | 0.6819 | 0.9017 |  0.9749 | 0.9511 | 0.9534 | 0.8892 | 0.9418 | 0.9739 | 0.9675 | 0.9859
ERGAS$\downarrow$ | 594.71 | 110.81 | 55.57 | 63.09 | 56.60 | 114.86 | 64.70 |  46.07 | 52.40 | 41.85
SAM$\downarrow$ | 0.3363 | 0.0265 | 0.0114 | 0.0207 | 0.0215 | 0.0231 | 0.0255 |  0.0167 | 0.0180 | 0.0080
#4 | PSNR$\uparrow$ | 13.24 | 26.34 | 31.83 | 30.02 | 32.01 | 25.64 | 30.83 |  33.28 | 31.81 | 34.24
SSIM$\uparrow$ | 0.7193 | 0.8857 |  0.9726 | 0.9374 | 0.9506 | 0.8613 | 0.9369 | 0.9697 | 0.9571 | 0.9802
ERGAS$\downarrow$ | 663.22 | 128.05 | 68.36 | 82.72 | 65.55 | 138.33 | 75.18 |  56.89 | 67.37 | 51.29
SAM$\downarrow$ | 0.3519 | 0.0245 | 0.0095 | 0.0206 | 0.0158 | 0.0240 | 0.0190 | 0.0137 | 0.0179 | 0.0082
#5 | PSNR$\uparrow$ | 10.54 | 23.02 | 27.59 | 28.58 | 30.26 | 23.85 | 29.28 |  31.72 | 29.42 | 32.59
SSIM$\uparrow$ | 0.4817 | 0.7669 | 0.9340 | 0.8807 | 0.8992 | 0.7331 | 0.8765 |  0.9386 | 0.9146 | 0.9610
ERGAS$\downarrow$ | 1052.98 | 185.63 | 108.97 | 95.38 | 78.63 | 165.94 | 87.76 |  66.56 | 87.25 | 61.15
SAM$\downarrow$ | 0.6842 | 0.0475 |  0.0238 | 0.0386 | 0.0297 | 0.0448 | 0.0346 | 0.0259 | 0.0370 | 0.0159
#6 | PSNR$\uparrow$ | 12.39 | 24.03 | 24.99 | 28.15 | 31.16 | 23.45 | 30.42 |  31.95 | 27.90 | 32.80
SSIM$\uparrow$ | 0.6147 | 0.8421 | 0.9229 | 0.9025 | 0.9247 | 0.8154 | 0.9120 |  0.9415 | 0.9153 | 0.9660
ERGAS$\downarrow$ | 699.76 | 159.04 | 140.87 | 98.41 | 68.07 | 169.59 | 74.07 |  62.35 | 100.54 | 58.96
SAM$\downarrow$ | 0.4981 | 0.0535 | 0.0362 | 0.0438 | 0.0312 | 0.0551 | 0.0343 |  0.0279 | 0.0385 | 0.0187
#7 | PSNR$\uparrow$ | 20.65 | 31.62 | 31.33 | 32.98 | 32.83 | 31.61 | 32.18 | 30.86 | 33.82 |  33.21
SSIM$\uparrow$ | 0.9162 | 0.9647 | 0.9709 | 0.9733 | 0.9655 | 0.9638 | 0.9646 | 0.9632 | 0.9774 |  0.9756
ERGAS$\downarrow$ | 248.40 | 68.02 | 73.98 | 59.48 | 60.15 | 68.41 | 63.95 | 75.15 | 54.14 |  57.84
SAM$\downarrow$ | 0.0956 | 0.0095 | 0.0098 | 0.0099 | 0.0121 |  0.0082 | 0.0122 | 0.0125 | 0.0092 | 0.0081
#8 | PSNR$\uparrow$ | 14.17 | 25.66 | 26.79 | 29.07 |  29.97 | 25.19 | 29.29 | 29.93 | 29.94 | 31.07
SSIM$\uparrow$ | 0.7535 | 0.8883 | 0.9289 | 0.9320 | 0.9268 | 0.8765 | 0.9247 | 0.9411 |  0.9437 | 0.9613
ERGAS$\downarrow$ | 541.15 | 130.14 | 119.73 | 89.92 |  79.98 | 137.46 | 85.99 | 80.78 | 82.59 | 71.08
SAM$\downarrow$ | 0.2950 | 0.0352 |  0.0296 | 0.0304 | 0.0563 | 0.0311 | 0.0370 | 0.0377 | 0.0299 | 0.0199
#9 | PSNR$\uparrow$ | 14.61 | 26.19 | 27.66 | 29.64 | 30.50 | 25.47 | 29.79 | 30.08 |  30.91 | 31.47
SSIM$\uparrow$ | 0.7734 | 0.9105 | 0.9439 | 0.9426 | 0.9369 | 0.9024 | 0.9353 | 0.9488 |  0.9560 | 0.9664
ERGAS$\downarrow$ | 514.77 | 124.67 | 108.15 | 83.85 | 75.66 | 135.52 | 81.86 | 79.63 |  72.40 | 68.25
SAM$\downarrow$ | 0.2575 | 0.0261 |  0.0201 | 0.0227 | 0.0363 | 0.0236 | 0.0287 | 0.0290 | 0.0218 | 0.0140
Mean | PSNR$\uparrow$ | 13.81 | 25.83 | 28.57 | 29.71 | 30.82 | 25.51 | 30.04 |  31.39 | 30.76 | 32.48
SSIM$\uparrow$ | 0.6681 | 0.8651 | 0.9409 | 0.9246 | 0.9274 | 0.8463 | 0.9192 |  0.9465 | 0.9421 | 0.9661
ERGAS$\downarrow$ | 686.20 | 142.25 | 110.09 | 87.94 | 77.73 | 146.45 | 83.56 |  73.64 | 79.27 | 65.81
SAM$\downarrow$ | 0.5160 | 0.1333 | 0.1047 | 0.0853 | 0.0766 | 0.1300 | 0.0787 |  0.0684 | 0.0766 | 0.0550

Table 6: Metrics on the Jizzakh urban dataset. The best and the second best
results are highlighted in bold and underline, respectively.

Time mode | Metrics | Obs | SNN | KBR | TNN | SPCTV | TNNTV | FCTNTC | FCTNFR | TCTV | HNN
---|---|---|---|---|---|---|---|---|---|---|---
#1 | PSNR$\uparrow$ | 10.71 | 24.92 | 26.41 | 24.69 | 25.16 | 25.03 | 24.84 | 24.50 | 25.32 |  26.30
SSIM$\uparrow$ | 0.4982 | 0.7665 |  0.8487 | 0.7846 | 0.7877 | 0.7506 | 0.7880 | 0.8070 | 0.8174 | 0.8515
ERGAS$\downarrow$ | 892.96 | 137.10 | 114.27 | 140.64 | 131.42 | 135.40 | 136.45 | 141.58 | 131.03 |  116.67
SAM$\downarrow$ | 0.6135 | 0.0403 |  0.0365 | 0.0457 | 0.0400 | 0.0401 | 0.0441 | 0.0449 | 0.0417 | 0.0311
#2 | PSNR$\uparrow$ | 12.31 | 24.39 | 26.89 | 25.56 | 27.43 | 24.38 | 26.95 |  27.73 | 27.15 | 28.76
SSIM$\uparrow$ | 0.5487 | 0.8036 |  0.9062 | 0.8408 | 0.8737 | 0.7844 | 0.8650 | 0.8985 | 0.8806 | 0.9279
ERGAS$\downarrow$ | 753.35 | 156.03 | 115.24 | 134.02 | 108.08 | 156.19 | 114.02 |  104.07 | 112.21 | 93.24
SAM$\downarrow$ | 0.5071 | 0.0463 | 0.0356 | 0.0475 | 0.0315 | 0.0463 | 0.0353 |  0.0308 | 0.0410 | 0.0215
#3 | PSNR$\uparrow$ | 14.08 | 26.91 | 30.99 | 28.56 | 30.90 | 26.89 | 30.43 |  31.47 | 30.86 | 32.88
SSIM$\uparrow$ | 0.6720 | 0.8837 |  0.9589 | 0.9109 | 0.9322 | 0.8686 | 0.9234 | 0.9474 | 0.9387 | 0.9680
ERGAS$\downarrow$ | 566.35 | 115.47 | 71.75 | 94.38 | 71.86 | 116.48 | 75.89 |  67.32 | 73.07 | 57.39
SAM$\downarrow$ | 0.3369 | 0.0265 |  0.0162 | 0.0271 | 0.0180 | 0.0264 | 0.0203 | 0.0178 | 0.0224 | 0.0109
#4 | PSNR$\uparrow$ | 14.11 | 26.79 | 29.53 | 27.63 | 30.22 | 26.59 | 29.79 |  30.64 | 29.80 | 32.63
SSIM$\uparrow$ | 0.7192 | 0.8707 |  0.9498 | 0.8906 | 0.9242 | 0.8561 | 0.9192 | 0.9426 | 0.9241 | 0.9652
ERGAS$\downarrow$ | 559.53 | 115.78 | 85.12 | 104.00 | 77.00 | 119.17 | 80.97 |  73.42 | 81.27 | 58.71
SAM$\downarrow$ | 0.3519 | 0.0283 | 0.0217 | 0.0304 | 0.0193 | 0.0282 | 0.0217 |  0.0192 | 0.0256 | 0.0128
#5 | PSNR$\uparrow$ | 10.51 | 23.69 | 28.00 | 25.48 | 28.28 | 23.86 | 27.83 |  28.72 | 26.87 | 29.91
SSIM$\uparrow$ | 0.4747 | 0.7441 |  0.9092 | 0.8088 | 0.8617 | 0.7160 | 0.8485 | 0.8976 | 0.8619 | 0.9339
ERGAS$\downarrow$ | 965.49 | 163.06 | 97.45 | 130.22 | 94.18 | 159.82 | 98.63 |  89.64 | 110.73 | 79.05
SAM$\downarrow$ | 0.6844 | 0.0439 | 0.0303 | 0.0494 |  0.0299 | 0.0458 | 0.0340 | 0.0311 | 0.0456 | 0.0201
#6 | PSNR$\uparrow$ | 11.22 | 25.84 |  29.78 | 24.83 | 29.17 | 25.75 | 28.52 | 28.32 | 25.17 | 30.64
SSIM$\uparrow$ | 0.6041 | 0.8245 |  0.9271 | 0.8323 | 0.8877 | 0.8073 | 0.8747 | 0.8961 | 0.8562 | 0.9378
ERGAS$\downarrow$ | 772.65 | 120.82 |  77.19 | 137.24 | 81.76 | 121.98 | 88.08 | 90.01 | 133.95 | 69.32
SAM$\downarrow$ | 0.4985 | 0.0246 |  0.0148 | 0.0394 | 0.0197 | 0.0267 | 0.0240 | 0.0242 | 0.0378 | 0.0142
#7 | PSNR$\uparrow$ | 18.74 | 30.30 |  32.69 | 29.02 | 31.97 | 30.19 | 30.21 | 29.98 | 29.64 | 33.05
SSIM$\uparrow$ | 0.9108 | 0.9581 | 0.9723 | 0.9502 | 0.9631 | 0.9595 | 0.9573 | 0.9588 | 0.9568 |  0.9721
ERGAS$\downarrow$ | 292.35 | 78.54 |  60.25 | 90.79 | 64.51 | 79.57 | 79.13 | 80.73 | 84.09 | 57.42
SAM$\downarrow$ | 0.0958 | 0.0114 | 0.0083 | 0.0128 |  0.0083 | 0.0111 | 0.0104 | 0.0100 | 0.0117 | 0.0073
#8 | PSNR$\uparrow$ | 13.42 | 26.36 |  28.07 | 25.59 | 27.53 | 26.45 | 26.74 | 26.69 | 26.26 | 28.11
SSIM$\uparrow$ | 0.7359 | 0.8743 |  0.9181 | 0.8768 | 0.8935 | 0.8674 | 0.8883 | 0.9040 | 0.8933 | 0.9203
ERGAS$\downarrow$ | 553.45 | 114.75 | 92.56 | 126.79 | 98.76 | 114.00 | 108.27 | 108.49 | 118.66 |  92.74
SAM$\downarrow$ | 0.2967 | 0.0211 |  0.0162 | 0.0257 | 0.0186 | 0.0213 | 0.0211 | 0.0217 | 0.0245 | 0.0155
#9 | PSNR$\uparrow$ | 14.50 | 26.48 | 28.45 | 26.67 | 27.83 | 26.56 | 26.82 | 26.68 | 27.82 |  28.21
SSIM$\uparrow$ | 0.7625 | 0.8965 | 0.9333 | 0.9028 | 0.9109 | 0.8921 | 0.9055 | 0.9173 | 0.9218 |  0.9329
ERGAS$\downarrow$ | 502.23 | 117.70 | 92.13 | 114.61 | 99.04 | 117.36 | 111.52 | 113.68 | 100.09 |  95.54
SAM$\downarrow$ | 0.2585 | 0.0213 |  0.0163 | 0.0235 | 0.0176 | 0.0213 | 0.0204 | 0.0217 | 0.0211 | 0.0145
Mean | PSNR$\uparrow$ | 13.29 | 26.19 |  28.98 | 26.45 | 28.72 | 26.19 | 28.02 | 28.30 | 27.65 | 30.05
SSIM$\uparrow$ | 0.6584 | 0.8469 |  0.9248 | 0.8664 | 0.8927 | 0.8335 | 0.8855 | 0.9077 | 0.8945 | 0.9344
ERGAS$\downarrow$ | 680.84 | 126.60 |  91.21 | 120.51 | 93.88 | 126.50 | 101.02 | 98.96 | 107.02 | 82.41
SAM$\downarrow$ | 0.5059 | 0.0861 |  0.0641 | 0.0884 | 0.0673 | 0.0871 | 0.0722 | 0.0702 | 0.0791 | 0.0537

Figure 5: The false-color decloud images (the 5th time node, band 1-4-7) of
all compared methods on the Jizzakh mountain dataset.

###  6.4 Experiment on multi-temporal cloud removal

In this section, we concentrate on a specialized inpainting task: multi-
temporal cloud removal. In this scenario, areas covered by clouds are
considered missing pixels, and the objective is to reconstruct these pixels.
However, clouds have distinctive shape patterns, making this task more
challenging than randomly pixel missing cases.

A series of comprehensive synthetic experiments are conducted to evaluate the
effectiveness of the HNN method. The synthetic datasets comprise two remote
sensing images captured over Jizzakh by the Landsat-8 satellite, encompassing
mountain and urban areas, with seven channels and nine time nodes. The size of
these data is $512\times 512\times 7\times 9$. To ensure the fidelity of the
synthetic data to real-world conditions, we selected nine real cloud masks
with varying levels of cloud coverage from the remote sensing cloud detection
dataset.

Table 5 summarizes the performance metrics of all the compared methods on the
Jizzakh mountain dataset. HNN demonstrates the most favorable performance. It
consistently achieves the highest PSNR and SSIM scores while maintaining the
lowest ERGAS and SAM scores across nearly all time nodes.

Table 6 presents the performance metrics on the Jizzakh urban dataset. As with
the results obtained from the Jizzakh mountain dataset, it is evident that the
HNN consistently achieves the best performance metrics across the majority of
time nodes, with the second-best results observed at time nodes 1 and 9. In
comparison to other methods, the results of HNN demonstrate both stability and
superiority.

To further highlight the superiority of the HNN method, Fig. 5 displays false-
color declouded images. To facilitate a detailed comparison, specific regions
within each image are magnified for enhanced visibility. In Fig. 5, the SNN,
SPCTV, and FCTNTC methods are observed to fill the cloud regions rather
haphazardly, while the TNN and TCTV methods do not successfully reconstruct
the original texture details beneath the clouds. The KBR and FCTNFR methods
exhibit color discrepancies compared to the original images. Conversely, the
HNN method effectively restores the intricate texture details, resulting in
the restored images that closely resemble the original cloud-free images.

###  6.5 Disccusion

####  6.5.1 Running time

Table 7 compares the execution times of all methods. For inpainting, HNN is
notably faster than FCTNFR and TCTV. SNN is the only method faster than HNN,
but SNN does not yield satisfactory restoration results. Fig. 6 presents a
scatter plot of the Jizzakh mountain dataset, clearly illustrating the
performance-time relationship. Remarkably, HNN achieves the highest PSNR while
maintaining relatively short execution times. In contrast, although TCTV and
FCTNFR yield impressive results, their lengthy execution times are a concern.
Similarly, while SNN and TNN have shorter execution times, their PSNR does not
meet expectations.

Table 7: The running time of compared methods (in seconds).

Task | Datasets | SNN | KBR | TNN | SPCTV | TNNTV | FCTNTC | FCTNFR | TCTV | HNN
---|---|---|---|---|---|---|---|---|---|---
HSI Inpainting | BA | 71.19 | 1310.4 | 219.42 | 3272.72 | 1108.82 | 244.23 | 285.69 | 2020.76 |  177.13 [t]
PC | 18.45 | 276.63 |  39.38 | 885.58 | 211.36 | 68.22 | 101.13 | 275.12 | 60.13 [b]
Cloud Removal | Mountain | 109.54 | 2099.33 | 417.33 | 4643.82 | 1858.6 | 1320.26 | 5474.12 | 3229.33 |  190.36 [t]
Urban | 71.71 | 1730.17 | 348.06 | 3953.04 | 1828.23 | 968.53 | 4682.67 | 3105.27 |  183.75
Real-world |  89.33 | 727.18 | 139.41 | 1445.49 | 501.31 | 373.3 | 1317.16 | 891.96 |  43.71 [b]
Task | Datasets | LLRT | TDL | LTDL | NGMeet | CTV | 3DTNN | 3DLogTNN | WNLRATV | HNN
HSI Denoising | BA | 1456.47 | 32.73 | 16160.75 |  75.23 | 118.97 | 259.68 | 387.73 | 646.28 | 84.85 [t]
PC | 399.67 | 6.8 | 15060.76 | 35.03 | 21.48 | 46.25 | 70.40 | 200.43 |  19.53 [b]

####  6.5.2 Numerical Analysis of Convergence

The algorithm’s convergence is guaranteed by the fact that the current model
satisfies the sufficient conditions for ADMM convergence. To further validate
its convergence, we conduct a numerical analysis of the HNN’s convergence.
Fig. 7 depicts the curves of PSNR, SSIM, ERGAS, and SAM values versus the
number of iterations on the Jizzakh mountain dataset. Notably, after 40
iterations, their values reach a stable state. This observation provides
further evidence for the convergence of the ADMM algorithm employed in this
study for solving the HNN model.

####  6.5.3 Influence of the Temporal Number

This section explores the influence of temporal number on HNN’s performance,
by varying the temporal number on the Jizzakh mountain dataset. Recalling the
results with 9 time nodes (listed in Table 5), it is shown that, if not
considering HNN, FCTNFC and TCTV are two best performer among the existing
methods, achieving PSNR values of 31.39 dB and 30.76 dB, respectively.
However, Table 8 shows that the HNN achieves PSNR values of 31.80 dB with only
6 time nodes, surpassing the performance of the existing methods by a
considerable improvement. This result further confirms the effectiveness of
the HNN method.

Figure 6: The scatter plot of PSNR versus time for the Jizzakh mountain dataset. Table 8:  Metrics on the Jizzakh mountain dataset with varying temporal numbers. Metrics | 3 | 4 | 5 | 6 | 7 | 8 | 9
---|---|---|---|---|---|---|---
PSNR | 24.33 | 28.86 | 30.58 | 31.80 | 32.21 | 32.09 | 32.48
SSIM | 0.9142 | 0.9488 | 0.9527 | 0.9609 | 0.9645 | 0.9619 | 0.9661
ERGAS | 161.97 | 98.36 | 84.04 | 73.06 | 68.93 | 71.31 | 65.81
SAM | 0.0847 | 0.0606 | 0.0551 | 0.0539 | 0.0526 | 0.0578 | 0.0550
Figure 7: Metric curves versus iteration number.

##  7 Conclusion

This paper introduces the HNN as a regularization term for remote sensing
image restoration. HNN effectively combines low-rank and smoothness priors and
is computationally efficient. Its superior performance in inpainting,
denoising, and cloud removal tasks makes it a valuable tool for enhancing the
quality and utility of remote sensing imagery. Future research can explore the
application of HNN in other remote sensing image processing tasks, such as HSI
unmixing and classification.

## References

  * [1] Yi Chang, Luxin Yan, and Sheng Zhong.  Hyper-laplacian regularized unidirectional low-rank tensor recovery for multispectral image denoising.  In 2017 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2017, Honolulu, HI, USA, July 21-26, 2017, pages 5901–5909. IEEE Computer Society, 2017.
  * [2] Yang Chen, Wenfei Cao, Li Pang, and Xiangyong Cao.  Hyperspectral image denoising with weighted nonlocal low-rank model and adaptive total variation regularization.  IEEE Trans. Geosci. Remote. Sens., 60:1–15, 2022.
  * [3] Yang Chen, Xiangyong Cao, Qian Zhao, Deyu Meng, and Zongben Xu.  Denoising hyperspectral image with non-i.i.d. noise structure.  IEEE Trans. Cybern., 48(3):1054–1066, 2018.
  * [4] Yong Chen, Wei He, Naoto Yokoya, Ting-Zhu Huang, and Xi-Le Zhao.  Nonlocal tensor-ring decomposition for hyperspectral image denoising.  IEEE Trans. Geosci. Remote. Sens., 58(2):1348–1362, 2020.
  * [5] Yong Chen, Ting-Zhu Huang, Wei He, Xi-Le Zhao, Hongyan Zhang, and Jinshan Zeng.  Hyperspectral image denoising using factor group sparsity-regularized nonconvex low-rank approximation.  IEEE Trans. Geosci. Remote. Sens., 60:1–16, 2022.
  * [6] Yongyong Chen, Yanwen Guo, Yongli Wang, Dong Wang, Chong Peng, and Guoping He.  Denoising of hyperspectral images using nonconvex low rank matrix approximation.  IEEE Trans. Geosci. Remote. Sens., 55(9):5366–5380, 2017.
  * [7] Lieven De Lathauwer, Bart De Moor, and Joos Vandewalle.  A multilinear singular value decomposition.  SIAM Journal on Matrix Analysis and Applications, 21(4):1253–1278, 2000.
  * [8] Silvia Gandy, Benjamin Recht, and Isao Yamada.  Tensor completion and low-n-rank tensor recovery via convex optimization.  Inverse Problems, 27(2):025010, jan 2011.
  * [9] Wei He, Quanming Yao, Chao Li, Naoto Yokoya, Qibin Zhao, Hongyan Zhang, and Liangpei Zhang.  Non-local meets global: An iterative paradigm for hyperspectral image restoration.  IEEE Trans. Pattern Anal. Mach. Intell., 44(4):2089–2107, 2022\.
  * [10] Wei He, Hongyan Zhang, Liangpei Zhang, and Huanfeng Shen.  Total-variation-regularized low-rank matrix factorization for hyperspectral image restoration.  IEEE Trans. Geosci. Remote. Sens., 54(1):178–188, 2016.
  * [11] F. L. Hitchcock.  The expression of a tensor or a polyadic as a sum of products.  J. Math. Phys., 6(1-4):164–189, 1927.
  * [12] Danfeng Hong, Wei He, Naoto Yokoya, Jing Yao, Lianru Gao, Liangpei Zhang, Jocelyn Chanussot, and Xiaoxiang Zhu.  Interpretable hyperspectral artificial intelligence: When nonconvex modeling meets hyperspectral remote sensing.  IEEE Geoscience and Remote Sensing Magazine, 9(2):52–87, 2021.
  * [13] Xingfang Huang, Shuang Xu, Chunxia Zhang, and Jiangshe Zhang.  Robust CP tensor factorization with skew noise.  IEEE Signal Process. Lett., 27:785–789, 2020.
  * [14] M. E. Kilmer, K. Braman, N. Hao, and R. C. Hoover.  Third-order tensors as operators on matrices: A theoretical and computational framework with applications in imaging.  SIAM J. Matrix Anal. Appl., 34(1):148–172, 2013.
  * [15] Ching-Yun Ko, Kim Batselier, Luca Daniel, Wenjian Yu, and Ngai Wong.  Fast and accurate tensor completion with total variation regularized tensor trains.  IEEE Trans. Image Process., 29:6918–6931, 2020.
  * [16] Xutao Li, Yunming Ye, and Xiaofei Xu.  Low-rank tensor completion with total variation for visual data inpainting.  In Satinder Singh and Shaul Markovitch, editors, Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence, February 4-9, 2017, San Francisco, California, USA, pages 2210–2216. AAAI Press, 2017\.
  * [17] Ji Liu, Przemyslaw Musialski, Peter Wonka, and Jieping Ye.  Tensor completion for estimating missing values in visual data.  IEEE Trans. Pattern Anal. Mach. Intell., 35(1):208–220, 2013\.
  * [18] Zhen Long, Ce Zhu, Jiani Liu, and Yipeng Liu.  Bayesian low rank tensor ring for image recovery.  IEEE Trans. Image Process., 30:3568–3580, 2021.
  * [19] Canyi Lu, Jiashi Feng, Zhouchen Lin, and Shuicheng Yan.  Exact low tubal rank tensor recovery from gaussian measurements.  In Jérôme Lang, editor, International Joint Conference on Artificial Intelligence, IJCAI, July 13-19, 2018, Stockholm, Sweden, pages 2504–2510, 2018.
  * [20] Canyi Lu, Xi Peng, and Yunchao Wei.  Low-rank tensor completion with a new tensor nuclear norm induced by invertible linear transforms.  In Computer Vision and Pattern Recognition, CVPR, Long Beach, CA, USA, June 16-20, 2019, pages 5996–6004, 2019.
  * [21] I. V. Oseledets.  Tensor-train decomposition.  SIAM J. Sci. Comput., 33(5):2295–2317, 2011.
  * [22] Chong Peng, Yang Liu, Kehan Kang, Yongyong Chen, Xinxing Wu, Andrew Cheng, Zhao Kang, Chenglizhao Chen, and Qiang Cheng.  Hyperspectral image denoising using nonconvex local low-rank and sparse separation with spatial-spectral total variation regularization.  IEEE Trans. Geosci. Remote. Sens., 60:1–17, 2022.
  * [23] Jiangjun Peng, Hailin Wang, Xiangyong Cao, Qian Zhao, Jing Yao, Hongying Zhang, and Deyu Meng.  Learnable representative coefficient image denoiser for hyperspectral image.  IEEE Trans. Geosci. Remote. Sens., 62:1–16, 2024.
  * [24] Jiangjun Peng, Yao Wang, Hong-Ying Zhang, Jianjun Wang, and Deyu Meng.  Exact decomposition of joint low rankness and local smoothness plus sparse matrices.  IEEE Trans. Pattern Anal. Mach. Intell., 45(5):5766–5781, 2023\.
  * [25] Jiangjun Peng, Qi Xie, Qian Zhao, Yao Wang, Yee Leung, and Deyu Meng.  Enhanced 3dtv regularization and its applications on HSI denoising and compressed sensing.  IEEE Trans. Image Process., 29:7889–7903, 2020.
  * [26] Yi Peng, Deyu Meng, Zongben Xu, Chenqiang Gao, Yi Yang, and Biao Zhang.  Decomposable nonlocal tensor dictionary learning for multispectral image denoising.  In 2014 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2014, Columbus, OH, USA, June 23-28, 2014, pages 2949–2956. IEEE Computer Society, 2014.
  * [27] Duo Qiu, Minru Bai, Michael K. Ng, and Xiongjun Zhang.  Robust low-rank tensor completion via transformed tensor nuclear norm with total variation regularization.  Neurocomputing, 435:197–215, 2021.
  * [28] Yuning Qiu, Guoxu Zhou, Qibin Zhao, and Shengli Xie.  Noisy tensor completion via low-rank tensor ring.  IEEE Trans. Neural Networks Learn. Syst., 35(1):1127–1141, 2024\.
  * [29] L. R. Tucker.  Implications of factor analysis of three-way matrices for measurement of change.  Probl. Meas. Change, 15:122–137, May 1963.
  * [30] Hailin Wang, Jiangjun Peng, Xiangyong Cao, Jianjun Wang, Qibin Zhao, and Deyu Meng.  Hyperspectral image denoising via nonlocal spectral sparse subspace representation.  IEEE J. Sel. Top. Appl. Earth Obs. Remote. Sens., 16:5189–5203, 2023.
  * [31] Hailin Wang, Jiangjun Peng, Wenjin Qin, Jianjun Wang, and Deyu Meng.  Guaranteed tensor recovery fused low-rankness and smoothness.  IEEE Trans. Pattern Anal. Mach. Intell., 45(9):10990–11007, 2023\.
  * [32] Yao Wang, Jiangjun Peng, Qian Zhao, Yee Leung, Xi-Le Zhao, and Deyu Meng.  Hyperspectral image restoration via total variation regularized low-rank tensor decomposition.  IEEE J. Sel. Top. Appl. Earth Obs. Remote. Sens., 11(4):1227–1243, 2018.
  * [33] Qi Xie, Qian Zhao, Deyu Meng, and Zongben Xu.  Kronecker-basis-representation based tensor sparsity and its applications to tensor recovery.  IEEE Trans. Pattern Anal. Mach. Intell., 40(8):1888–1902, 2018\.
  * [34] Shuang Xu, Xiangyong Cao, Jiangjun Peng, Qiao Ke, Cong Ma, and Deyu Meng.  Hyperspectral image denoising by asymmetric noise modeling.  IEEE Trans. Geosci. Remote. Sens., 60:1–14, 2022.
  * [35] Shuang Xu, Chunxia Zhang, and Jiangshe Zhang.  Adaptive quantile low-rank matrix factorization.  Pattern Recognit., 103:107310, 2020.
  * [36] Jize Xue, Yongqiang Zhao, Wenzhi Liao, and Jonathan Cheung-Wai Chan.  Nonconvex tensor rank minimization and its applications to tensor recovery.  Inf. Sci., 503:109–128, 2019.
  * [37] Jize Xue, Yongqiang Zhao, Wenzhi Liao, and Jonathan Cheung-Wai Chan.  Nonlocal low-rank regularized tensor decomposition for hyperspectral image denoising.  IEEE Trans. Geosci. Remote. Sens., 57(7):5174–5189, 2019.
  * [38] Tatsuya Yokota, Qibin Zhao, and Andrzej Cichocki.  Smooth PARAFAC decomposition for tensor completion.  IEEE Trans. Signal Process., 64(20):5423–5436, 2016.
  * [39] Qiang Zhang, Yaming Zheng, Qiangqiang Yuan, Meiping Song, Haoyang Yu, and Yi Xiao.  Hyperspectral image denoising: From model-driven, data-driven, to model-data-driven.  IEEE Transactions on Neural Networks and Learning Systems, in press:1–21, 2023.
  * [40] Qibin Zhao, Guoxu Zhou, Shengli Xie, Liqing Zhang, and Andrzej Cichocki.  Tensor ring decomposition.  CoRR, abs/1606.05535, 2016.
  * [41] Wen-Jie Zheng, Xi-Le Zhao, Yu-Bang Zheng, and Ting-Zhu Huang.  Provable stochastic algorithm for large-scale fully-connected tensor network decomposition.  Journal of Scientific Computing, 98(1):1–16, 2024.
  * [42] Wen-Jie Zheng, Xi-Le Zhao, Yu-Bang Zheng, Jie Lin, Lina Zhuang, and Ting-Zhu Huang.  Spatial-spectral-temporal connective tensor network decomposition for thick cloud removal.  ISPRS Journal of Photogrammetry and Remote Sensing, 199:182–194, 2023.
  * [43] Wen-Jie Zheng, Xi-Le Zhao, Yu-Bang Zheng, and Zhi-Feng Pang.  Nonlocal patch-based fully connected tensor network decomposition for multispectral image inpainting.  IEEE Geoscience and Remote Sensing Letters, 19:1–5, art no. 8025105, 2022.
  * [44] Yu-Bang Zheng, Ting-Zhu Huang, Xi-Le Zhao, Yong Chen, and Wei He.  Double-factor-regularized low-rank tensor factorization for mixed noise removal in hyperspectral image.  IEEE Transactions on Geoscience and Remote Sensing, 58(12):8450–8464, 2020.
  * [45] Yu-Bang Zheng, Ting-Zhu Huang, Xi-Le Zhao, Tai-Xiang Jiang, Teng-Yu Ji, and Tian-Hui Ma.  Tensor N-tubal rank and its convex relaxation for low-rank tensor recovery.  Information Sciences, 532:170–189, 2020.
  * [46] Yu-Bang Zheng, Ting-Zhu Huang, Xi-Le Zhao, Tai-Xiang Jiang, Tian-Hui Ma, and Teng-Yu Ji.  Mixed noise removal in hyperspectral image via low-fibered-rank regularization.  IEEE Transactions on Geoscience and Remote Sensing, 58(1):734–749, 2020.
  * [47] Yu-Bang Zheng, Ting-Zhu Huang, Xi-Le Zhao, and Qibin Zhao.  Tensor completion via fully-connected tensor network decomposition with regularized factors.  Journal of Scientific Computing, 92:1–35, 2022.
  * [48] Yu-Bang Zheng, Ting-Zhu Huang, Xi-Le Zhao, Qibin Zhao, and Tai-Xiang Jiang.  Fully-connected tensor network decomposition and its application to higher-order tensor completion.  In Thirty-Fifth AAAI Conference on Artificial Intelligence Virtual Event, February 2-9, pages 11071–11078, 2021.
  * [49] Yu-Bang Zheng, Xi-Le Zhao, Junhua Zeng, Chao Li, Qibin Zhao, Heng-Chao Li, and Ting-Zhu Huang.  Svdinstn: A tensor network paradigm for efficient structure search from regularized modeling perspective.  In 2024 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2024.
  * [50] Lina Zhuang and Michael K. Ng.  Hyperspectral mixed noise removal by $\ell_{1}$-norm-based subspace representation.  IEEE J. Sel. Top. Appl. Earth Obs. Remote. Sens., 13:1143–1157, 2020.


Further author information: (Send correspondence to Nemanja Jovanovic.)
Nemanja Jovanovic: E-mail<EMAIL_ADDRESS>Telephone: +1 626 395 1214

# An all-photonic, dynamic device for flattening the spectrum of a laser
frequency comb for precise calibration of radial velocity measurements

Nemanja Jovanovic  Department of Astronomy, California Institute of
Technology, Pasadena, CA, 91125, USA  Pradip Gatkine  Department of Astronomy,
California Institute of Technology, Pasadena, CA, 91125, USA  Boqiang Shen  T.
J. Watson Laboratory of Applied Physics, California Institute of Technology,
Pasadena, CA, 91125, USA  Maodong Gao  T. J. Watson Laboratory of Applied
Physics, California Institute of Technology, Pasadena, CA, 91125, USA  Nick
Cvetojevic  Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS,
Laboratoire Lagrange, France  Katarzyna Ławniczuk  Bright Photonics BV,
Horsten 1, 5612 AX Eindhoven, The Netherlands  Ronald Broeke  Bright Photonics
BV, Horsten 1, 5612 AX Eindhoven, The Netherlands  Charles Beichman  Jet
Propulsion Laboratory, 4800 Oak Grove Drive, Pasadena, CA 91109, USA
Stephanie Leifer  The Aerospace Corporation, 2310 E. El Segundo Blvd., El
Segundo, CA 90245  Jeffery Jewell  Jet Propulsion Laboratory, 4800 Oak Grove
Drive, Pasadena, CA 91109, USA  Gautam Vasisht  Jet Propulsion Laboratory,
4800 Oak Grove Drive, Pasadena, CA 91109, USA  Dimitri Mawet  Department of
Astronomy, California Institute of Technology, Pasadena, CA, 91125, USA  Jet
Propulsion Laboratory, 4800 Oak Grove Drive, Pasadena, CA 91109, USA

###### Abstract

Laser frequency combs are fast becoming critical to reaching the highest
radial velocity precisions. One shortcoming is the highly variable brightness
of the comb lines across the spectrum (up to 4-5 orders of magnitude). This
can result in some lines saturating while others are at low signal and lost in
the noise. Losing lines to either of these effects reduces the precision and
hence effectiveness of the comb. In addition, the brightness of the comb lines
can vary with time which could drive comb lines with initially reasonable
SNR’s into the two regimes described above. To mitigate these two effects,
laser frequency combs use optical flattener’s.

Flattener’s are typically bulk optic setups that disperse the comb light with
a grating, and then use a spatial light modulator to control the amplitude
across the spectrum before recombining the light into another single mode
fiber and sending it to the spectrograph. These setups can be large (small
bench top), expensive (several hundred thousand dollars) and have limited
stability. To address these issues, we have developed an all-photonic spectrum
flattener on a chip. The device is constructed from optical waveguides on a
SiN chip. The light from the laser frequency comb’s output optical fiber can
be directly connected to the chip, where the light is first dispersed using an
arrayed waveguide grating. To control the brightness of each channel, the
light is passed through a Mach-Zehnder interferometer before being recombined
with a second arrayed waveguide grating. Thermo-optic phase modulators are
used in each channel before recombination to path length match the channels as
needed.

Here we present the results from our first generation prototype. The device
operates from 1400-1800 nm (covering the H band), with 20, 20 nm wide
channels. The device was mounted on a PCB board to enable electrical control
of the active elements and tested in the laboratory. It was demonstrated that
the Mach-Zehnder’s allowed for nearly 40 dBs of dynamic modulation of the
spectrum, which is greater than that offered by most spatial light modulators.
With a smooth spectrum light source (superluminescent light source), we
reduced the spectral variation to 3 dBs, limited by the properties of the
components used. On a laser frequency comb which had strong modulations at
high spatial frequencies, we still managed to reduce the modulation to 5 dBs.
These devices are of the order of a US quarter and could play a significant
role in future PRV and EPRV initiatives.

###### keywords:

laser frequency combs, spectrographs, calibration, precision radial velocity,
spectral shaping, spectral flattening, integrated photonics, exoplanets

##  1 INTRODUCTION

Laser frequency combs (LFC) are becoming essential for the calibration of
spectrographs to enable the detection and characterization of exoplanets via
the precision radial velocity (PRV) technique [1]. In order to best utilize an
LFC, and hence maximize the accuracy of the wavelength calibration, the
largest possible number of comb lines within the observational band pass would
be used to derive the wavelength solution. However, the intensity of the lines
across the spectrum can vary by many orders magnitude making it impossible to
extract all lines with an optimum signal-to-noise ratio (SNR). An example of
the issue can be seen in the left panel of Fig. 1. This shows a trace from the
LFC spectrum used with the PARVI spectrograph at Palomar Observatory [2],
which exhibits amplitude differences of up to 45 dB across the 400 nm window,
from 1400 to 1800 nm used by the instrument currently.

Figure 1: (Left) a spectrum of the PARVI LFC showing the strong intensity
modulations across the spectrum. (Right) PARVI’s bulk optic flattener. The
flattener occupies a volume of $1\times 0.5\times 0.5$ m3.

To combat this problem, spectral flatteners are typically constructed from
bulk optics, where light is dispersed with a grating and a spatial light
modulator is used to control the amplitude across the spectrum, before
recombining the light into a fiber and sending it to the spectrograph [3]. As
seen in the right panel of Fig. 1, these setups can be large (small bench
tops) and expensive (several hundred thousand dollars based on two flatteners
purchased by this team around the time of this article). Given that the final
comb spectrum is typically generated inside a waveguide or optical fiber,
finding an all photonic solution to flattening the spectrum would be the most
elegant way of miniaturizing this necessary component.

In addition, the brightness of the comb lines can vary with time which could
drive comb lines into either the noise floor or non-linear/saturation regions
of the detectors dynamic range. The evolution of the amplitude in the spectrum
of the LFC is typically slow (on timescales of many minutes to hours) and
significantly smaller (1-10 dB) than the amplitude differences across the
spectrum. Nonetheless, these variations necessitate that the flattener be a
dynamic element that can track and control the evolution in the spectrum with
time.

To address both of these issues, we have developed a broadband, all-photonic
spectrum shaper/flattener (BAPSS) on-a-chip. The device was constructed from
optical waveguides on a silicon nitride (SiN) chip and operates from 1400-1800
nm (covering the astronomical H band), with 20, 20 nm wide channels. The
device was mounted on a PCB board to enable electrical control of the active
elements and tested in the laboratory. Here we present the laboratory
characterization of a prototype BAPSS device.

##  2 Circuit layout

The architecture of the circuit we choose for this application is shown in
Fig. 2 and consists of the following elements. Firstly, an arrayed waveguide
grating (AWGs) is used to disperse the light into 20 discrete, narrow band (20
nm) channels. The channel spacing was chosen to match the bandwidth over which
the LFC amplitude changes in unison, in addition to reducing the overall
number of channels. Next, a series of Mach-Zehnder Interferometers (MZIs), one
in each AWG output, are used to adjust the power transmitted through that
channel, by modifying the path length on one arm via Thermo-Optic Phase
Modulators (TOPMs). TOPMs following each MZI are used to adjust the path
length of a given channel with respect to its neighbor to be able to rephase
the spectrum. Finally, another AWG with the same properties as the first is
used to recombine the channels into a single spectrum.

Figure 2: Schematic layout of the circuit architecture used for the photonic
spectral shaper.

This development leveraged significant advances in commercial lithographic
fabrication techniques over the last few decades. Indeed, SiN fabrication via
Multi-Project Wafer (MPW) runs (fabrication runs where customers pay for only
a portion of a wafer), are now producing low loss ($<$0.5 dB/cm at 1550 nm
propagation loss) waveguides and devices that are highly reproducible, which
was ideal for rapid prototyping of concepts like this. The concept was
designed using Nazca Design (https://nazca-design.org/) with the elements
outlined above.

##  3 Experimental procedures

The circuit was fabricated by LioniX International as part of MPW run
$\\#$21\. The SiN MPW offering consists of a dual stripe waveguide geometry
embedded in a silica cladding. For details about the waveguide geometry refer
to the MPW overview manual (https://www.lionix-
international.com/photonics/mpw-services/). The waveguides were optimized to
allow for optimal guiding around 1550 nm with bend radii as tight as 100
$\mu$m without substantial losses. To improve coupling to optical fibers SSCs
were used at the edges of the chip. These consisted of tapers which would
expand the mode to $\sim$10 $\mu$m around 1550 nm.

After the devices were fabricated they were diced, polished and packaged to
simplify testing. This included bonding a 48-fiber vgroove array to one face
of the chip. The fiber used for the vgroove was SMF28. All inputs and outputs
to and from the various devices were accessed through this v-groove. The
circuit was mounted onto a PCB and the DC and ground lines connected so they
could be accessed via ribbon cables from the top and bottom of the device. The
entire assembly was mounted on a sub-mount that included a themo-electric
cooler (TEC). An image of the packaged device is shown in Fig. 3.

Figure 3: Experimental setup for characterizing the broadband all-photonic
spectrum shapers. The light source for the tests was either a SLD or an LFC. A
polarization controller was used to align the polarization of the source with
the TE mode of the waveguides. Light was injected via a v-groove array. The
fully packaged PIC is shown. The electrodes on the circuit (chromium plated)
as well as the AWGs can easily be seen in the photo. The output was routed to
an OSA.

To test the devices the setup shown in Fig. 3 was used. Initial testing was
conducted by injecting light from a super-luminescent diode (SLD, Thorlabs,
S5FC1550P-A2). Since the fibers bonded to the device under test were SMF28, we
used a polarization controller (PCON, Thorlabs, FPC561) to orient the
polarized signal of the source with the TE mode of the waveguides on the chip,
as the waveguides were highly polarization dependent. The output was connected
directly to an optical spectrum analyzer (OSA, Thorlabs, OSA202C).

To test each device, the output of the PCON was connected to a given device,
the output of that device to the OSA, and then the PCON was adjusted to
maximize the signal on the OSA. This ensured that the polarization of the
light source was aligned with the TE mode of the waveguides, which was most
efficient.

To control the TOPMs we used a multi-channel controller. Since the full device
consisted of 20 MZIs and 20 TOPMs, the driver had to support $>40$ active
channels simultaneously. To achieve a maximum of $3\pi$ phase shifts, we
required a driver capable of up to 20 V and 50 mA per line. For this reason we
used a multi-channel (120) driver from Nicslab (XPOW-120AX-CCvCV-U). A linear
power supply was used to power the multi-channel controller. A computer was
used to operate the controller. To connect the relevant pins of the controller
to the PCB and hence the device under test, we used an electrical breadboard
(Nicslab, M6 multiconnector).

Once testing with the SLD was completed, we undertook tests with an LFC.
However, the comb we had access to was pre-broadened, which meant that it had
a triangular shaped spectral profile, not representative of full broadened
combs, which are typically flatter, and was only several hundred nanometers
wide. Data were acquired with 0.2 nm resolution, where the lines of the comb
were not resolved as well as with 0.02 and 0.05 nm resolution where they were.

The SLD had sufficient flux to test the device down to 1400 nm and up to 1650
nm. The OSA limited testing at the long wavelength end to 1700 nm.

##  4 Results

Using the methods described in the previous section, the spectral response of
the BAPSS device was measured. Light was injected into the input of the full
BAPSS device and the output routed to the OSA. Spectra were taken with the
MZI’s turned off (blue) as well as turned on and adjusted for maximum
throughput while maintaining flatness (orange), minimum throughput (red) and
some arbitrary flat level in between (red) and are shown in Fig. 4.

Figure 4: Spectral response of the full device. The blue trace shows the
transmission of the device with no power applied. The orange is when the MZIs
and TOPMS are adjusted for maximum throughput while maintaining flatness. The
red trace is when they are adjusted for minimum throughput, and the green for
an arbitrary level in the middle of the range. The data were normalized to the
reference waveguide/loop back. This process removes the contribution of
coupling losses to/from the chip leaving only the losses of the
component/device.

The minimum losses occur between 1600-1650 nm, and are as low as 9.5 dBs. The
minimum theoretical losses of the AWG and MZI amount to 7.2-7.5 dB, which
don’t quite account for those measured and could be the result of
imperfections in manufacturing over a larger device size (i.e. non-
uniformities for example). The spectral ripple on the various flat states
(orange and green) have an amplitude of 2-3 dB, which is consistent with what
we would expect from the back-to-back AWG. Dashed blue lines indicate the
approximate average level for the maximum and minimum throughput states
indicating that the power can be adjusted over a 38 dB range. The short and
long wavelength ends of the red spectrum were impacted by the noise floor of
the OSA (purple trace).

The next set of tests involved the LFC. Figure 5 shows the results. The top
left panel shows the LFC spectrum through the reference waveguide (blue) as
well as the BAPSS device (orange, green, red and purple). It can be seen that
the LFC spectrum is triangular in shape on a logarithmic scale, due to the
pre-broadened nature of the output of the LFC as outlined in Section 3 above,
which is what we had access to. There is a large dip in the region around 1560
nm because a fiber Bragg grating (FBG) was used to suppress the pump region.
The grating was not optimized for the spectral profile of the pump or the
amplitude and so unfortunately suppressed a lot of light in that region for
these tests.

The orange trace shows the same spectrum through the BAPSS device, with the
device set to maximum throughput. The green and red traces show the spectrum
flattened at the -50 and -60 dBm levels. The flattening was done manually by
eye. The blue shaded boxes are 5 dB high, showing that we managed to flatten
the spectrum to within 5 dBs, across hundreds of nanometers of range as was
the case with the SLD tests. The only exception to this is the region around
the pump, which has a series of sharp spectral features too narrow for the
device to control with the current channel spacing’s used. Finally, the purple
trace shows how much of the LFC spectrum can be suppressed if we put the
device into the minimum transmission state. The brown trace is the noise floor
of the OSA indicating that in large sections of the spectrum we are limited by
the OSA.

Figure 5: Spectral response of the full device with an LFC. (left) LFC
spectrum through the reference waveguide (blue), and through the BAPSS device
with maximum throughput (orange), minimum throughput (purple), and two
arbitrary flat levels (green and red). The OSA noise floor is shown in brown.
(Right) the LFC spectrum taken through the references waveguide with and
without optical attenuation.

Finally, we tested the power handling capabilities of the device. To ensure we
did not damage the device, we connected the LFC to the reference waveguide and
took a spectrum. We then replaced the 20 dB attenuator, used in all LFC
experiemnts till this point, with a 10 and 5 dB version and finally removed
the attenuator altogether. The spectra are shown in the top right panel of
Fig. 5. The spectra are largely identical, and are simply offset vertically by
the amount corresponding to the attenuation factor. We measured the total
integrated power injected into the device with the 20 dB attenuator to be 3-4
mW. As we could not directly measure the power with no attenuator, we can
infer it would be a 20 dB brighter signal that is achromatic (all traces are
parallel to one another) and therefore assume that in the full power mode,
300-400 mW of power were being injected into the full device.

After confirming there were no deleterious effects with these levels of power
on the reference waveguide, we connected the LFC to the BAPSS device without
any attenuation and conducted some basic tests. The BAPSS device performed
identically to how it had with reduced power and is therefore capable of
working with power levels of typical LFCs.

Data were also acquired with resolutions of 0.02 and 0.05 nm where the comb
lines were resolved, and were deemed to be consistent with the lower
resolution data presented here. For simplicity we only show the low resolution
data in this paper.

##  5 Discussion

The results above demonstrate that the concept outlined works extremely well
for spectral shaping/flattening. Specifically, the device was shown to operate
across 250 nm of spectral range, the largest demonstrated from such a layout,
with losses consistent with theory. The MZIs were capable of modulating a
single spectral channel by 40 dBs, which indicates that the amplitude in the
arms of the interferometer is very well balanced, ruling out large
manufacturing imperfections. In addition, this modulation range exceeds nearly
all bulk optic flatterers which rely on SLMs, which are typically limited to
20-30 dBs only.

The circuit demonstrated here is actually designed for two purposes: static
compensation of the native LFC spectrum as well as for tracking the dynamic
aspects. In this work we have mostly focused on compensating for the static
profile of the spectrum. It should be noted that around the pump region of an
LFC there can be sharp features (as seen in Fig. 5), which would require a
photonic device with much narrower channels. This can be done, but more
channels means more electrodes and increases the complexity of the circuit.
Regardless, the ideal thing to do is to use custom FBGs to compensate for the
sharp static aspects of the spectrum, offloading that from the flattener,
leaving that to work on the broader parts of the spectrum. If however one
wishes to do all of this in the flattener itself, it is possible to use a
cascaded AWG approach whereby certain low resolution spectral channels are
sent to higher resolution second stage AWGs that provide the sampling needed,
but only in a limited number of channels. This is currently being explored and
will be focus of future work.

Although we did not demonstrate the ability of the device to track dynamic
changes to the spectrum, the device is more than capable to correct the slow
changes expected in an LFC spectrum: with a modulation amplitude of up to 40
dBs and thermal response time as stated by the vendor (and confirmed in a
recent publication [4]) of 1 ms it will be sufficient to correct for the slow
evolving, low amplitude changes typical of LFC spectra. How to implement the
control law for closed loop operation will require some consideration and if
the focus of future work.

Astronomical spectroscopy typically requires large bandwidths. We demonstrated
successful operation of the BAPSS over a single band (1490-1800 nm) in this
work. It may be possible to re-engineer the device to increase the bandwidth
to capture several bands and/or split the light at the output of the LFC and
send it to several devices each optimized for operation in a different band.
This is another avenue that can be explored in future.

###### Acknowledgements.

This work was produced with support from the Keck Institute for Space Studies.
N. Cvetojevic acknowledges the funding from the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation program
(grant agreement CoG - 683029).

## References

  * [1] Fischer, D. A., Anglada-Escude, G., Arriagada, P., Baluev, R. V., Bean, J. L., Bouchy, F., Buchhave, L. A., Carroll, T., Chakraborty, A., Crepp, J. R., Dawson, R. I., Diddams, S. A., Dumusque, X., Eastman, J. D., Endl, M., Figueira, P., Ford, E. B., Foreman-Mackey, D., Fournier, P., Fűrész, G., Gaudi, B. S., Gregory, P. C., Grundahl, F., Hatzes, A. P., Hébrard, G., Herrero, E., Hogg, D. W., Howard, A. W., Johnson, J. A., Jorden, P., Jurgenson, C. A., Latham, D. W., Laughlin, G., Loredo, T. J., Lovis, C., Mahadevan, S., McCracken, T. M., Pepe, F., Perez, M., Phillips, D. F., Plavchan, P. P., Prato, L., Quirrenbach, A., Reiners, A., Robertson, P., Santos, N. C., Sawyer, D., Segransan, D., Sozzetti, A., Steinmetz, T., Szentgyorgyi, A., Udry, S., Valenti, J. A., Wang, S. X., Wittenmyer, R. A., and Wright, J. T., “State of the field: Extreme precision radial velocities,” Publications of the Astronomical Society of the Pacific 128, 066001 (may 2016).
  * [2] Gibson, R. K., Oppenheimer, R., Matthews, C. T., and Vasisht, G., “Characterization of the C-RED 2: A High Frame Rate Near-Infrared Camera,” arXiv e-prints , arXiv:1911.04567 (Nov 2019).
  * [3] Probst, R. A., Steinmetz, T., Wilken, T., Wong, G. K. L., Hundertmark, H., Stark, S. P., Russell, P. S. J., Hänsch, T. W., Holzwarth, R., and Udem, T., “Spectral flattening of supercontinua with a spatial light modulator,” in [Techniques and Instrumentation for Detection of Exoplanets VI ], Shaklan, S., ed., 8864, 706 – 713, International Society for Optics and Photonics, SPIE (2013).
  * [4] Cvetojevic, N., Martinache, F., Chingaipe, P., Laugier, R., N’Diaye, M., Ligi, R., and Mary, D., “Closed-loop on-chip phase control and fringe tracking in photonic interferometers.,” in [Optical and Infrared Interferometry and Imaging VIII ], 12183, 45, SPIE (2022).