Equal contributions.

Equal contributions.

[1,4,]Chao Wang

[1]Department of Statistics and Data Science, Southern University of Science
and Technology, Shenzhen, P.R. China 2]School of Data Science, City University
of Hong Kong, Hong Kong SAR, P.R. China 3]School of Mathematical Science,
University of Electronic Science and Technology of China, Chengdu, P.R. China
4]National Centre for Applied Mathematics Shenzhen, Shenzhen, P.R. China
5]Department of Mathematics, Hong Kong Baptist University, Hong Kong SAR, P.R.
China

#  Hyperspectral and multispectral image fusion with arbitrary resolution
through self-supervised representations

Ting Wang  Zipei Yan  Jizhou Li  Xile Zhao<EMAIL_ADDRESS>Michael Ng
*  [  [  [  [

###### Abstract

The fusion of a low-resolution hyperspectral image (LR-HSI) with a high-
resolution multispectral image (HR-MSI) has emerged as an effective technique
for achieving HSI super-resolution (SR). Previous studies have mainly
concentrated on estimating the posterior distribution of the latent high-
resolution hyperspectral image (HR-HSI), leveraging an appropriate image prior
and likelihood computed from the discrepancy between the latent HSI and
observed images. Low rankness stands out for preserving latent HSI
characteristics through matrix factorization among the various priors.
However, this method only enhances resolution within the dimensions of the two
modalities. To overcome this limitation, we propose a novel continuous low-
rank factorization (CLoRF) by integrating two neural representations into the
matrix factorization, capturing spatial and spectral information,
respectively. This approach enables us to harness both the low rankness from
the matrix factorization and the continuity from neural representation in a
self-supervised manner. Theoretically, we prove the low-rank property and
Lipschitz continuity in the proposed continuous low-rank factorization.
Experimentally, our method significantly surpasses existing techniques and
achieves user-desired resolutions without the need for neural network
retraining.

###### keywords:

Low-rank factorization, arbitrary resolution, image fusion, continuous
representation.

##  1 Introduction

Figure 1: The pipeline of CLoRF for arbitrary resolution. (A) Train the CLoRF.
(B) Use the trained CLoRF to infer arbitrary resolutions of HSIs with given
spatial and spectral coordinates. (C) An example of CLoRF for super-resolution
on the Pavia University ($336\times 336\times 93$). The CLoRF is trained given
its LR-MSI ($168\times 168\times 4$) and HR-HSI ($42\times 42\times 50$), then
infers the original resolution. Bicubic interpolation directly upsamples from
LR. (D) Visualize the spectrum of a random pixel from the results on (C).

Hyperspectral images (HSIs) have widespread applications across various fields
due to their rich spectral information. The abundant spectral details offered
by HSIs facilitate accurate scene interpretation and enhance the efficacy of
numerous applications, including object classification cite[13 and anomaly
detection cite[15. However, the inherent trade-off between spectral and
spatial resolution in HSI systems, constrained by hardware limitations, often
results in HSIs with lower spatial resolution than RGB, panchromatic (PAN),
and multispectral images (MSI). To enhance the spatial resolution of HSIs, a
natural approach is to fuse LR-HSI and HR-MSI, known as hyperspectral and
multispectral image fusion (HSI-MSI fusion). HSI-MSI fusion resembles MSI
pansharpening, where low spatial resolution MSI is merged with high-resolution
PAN imagery. However, directly applying these pansharpening methods to fuse
HSI and MSI images suffers from challenges, as PAN images have limited
spectral information, leading to spectral distortion cite[21. Consequently,
numerous approaches tailored for HSI-MSI fusion are introduced, which can be
generally categorized into model-based methods and deep learning-based models.

Figure 2: Illustration of the proposed CLoRF for MSI-HSI fusion. The spatial
coordinates and spectral coordinates of HR-MSI and LR-HSI are fed into the
Spatial-INR $\Phi_{\theta}(\cdot)$ and Spectral-INR $\Psi_{\alpha}(\cdot)$ to
generate coefficients and the bases, respectively. Thereafter, the generated
coefficients and the bases are multiplied to recover the HR-HSI.

Model-based methods leverage the low-rank structure of HSIs by characterizing
their low-rankness through matrix factorization, decomposing the HSI matrix
into the basis and coefficients or endmembers and abundances. Therefore, the
principles of matrix factorization-based methods rely on appropriate prior
information and the likelihood determined by the relationships between the
latent HSI and the observed LR-HSI and HR-MSI. Given the known degradation
model, most existing methods focus on modeling prior information for HSIs,
including explicit and implicit methods. The explicit methods employ
handcrafted explicit prior information, such as low-rankness cite[45, 34,
smoothness cite[28, sparsity cite[11, and non-local similarity cite[9, to
depict the prior distribution of the latent HSI for super-resolution (SR).
Implicit smooth regularizations introduce basis functions to parameterize the
prior information, extending the modeling to functional representations. For
example, ? utilize non-negative matrix factorization parameterized by basis
functions to reveal implicit smoothness. Furthermore, with the success of deep
learning, ? implicitly introduce a denoising operator into the optimization
framework to learn the deep image prior, which is shared across all HSIs using
deep neural networks with different structures. Although these methods have
achieved significant success in HSI-MSI fusion, manually designed explicit
prior information may make the regularization terms problematic, resulting in
significant computational complexity. Meanwhile, implicitly designed prior
information sometimes becomes inappropriate for capturing the complex and
detailed structures of HSIs. Moreover, these methods are restricted to fixed
resolution and cannot be instantly applied to recovering HSIs with arbitrary
resolution.

Most deep learning-based methods typically involve supervised training on
large datasets consisting of observed and ground truth data to learn a complex
function mapping. Some representative deep learning methods are based on
convolution neural network (CNN) for HSI-MSI fusion, such as the two-stream
fusion network designed for HR-MSI and LR-HSI cite ?, [38, 18, 17, 36. With
the remarkable performance of implicit neural representations (INRs) in
continuous multi-dimensional data representation, some prior research cite ?,
[37, 16 combined INR and CNN for HSI-MSI fusion. Although these methods have
demonstrated promising results in fusion, they heavily rely on high-quality
training data for supervised learning. Collecting high-quality ground truth
data is extremely time-consuming and costly. In addition, unsupervised fusion
methods cite ?, [46, 24, 40 are introduced, eliminating the need for expensive
training data. However, these methods typically rely on complex network
structures. Besides, they are also confined to fixed resolution, and inferring
results at the user-desired resolutions requires retraining.

To address the aforementioned issues, we propose an unsupervised HSI-MSI
fusion framework for continuous low-rank factorization (CLoRF). Specifically,
CLoRF integrates two implicit neural representations into the low-rank
factorization, capturing continuous spatial and spectral information of HR-MSI
and LR-HSI, respectively. Unlike classic discrete matrix factorization, we
define the continuous matrix function factorization following tensor function
cite ?, [22, which characterizes low-rankness in continuous representation.
Each continuous function is a realization of INR parameterized by multi-layer
perceptrons (MLPs). As MLPs are Lipschitz smooth, further characterizing
smoothness in continuous representation. Compared to previous fusion methods,
as shown in Fig. 1, CLoRF can achieve arbitrary resolution in both spatial and
spectral domains without additional image information and retraining.

Our contributions are summarized as follows:

  1. 1.

We propose a novel unsupervised HSI-MSI fusion method: CLoRF, which represents
HSIs in a continuous representation using low-rank function factorization.
CLoRF can infer with arbitrary resolution without the need for retraining.

  2. 2.

We theoretically prove that the implicit regularization terms of low-rankness
and smoothness are unified in continuous representation, which justifies the
potential effectiveness for HSIs.

  3. 3.

We experimentally demonstrate that CLoRF significantly surpasses existing
techniques, which confirms its wide applicability and superiority.

  4. 4.

We further expand the application of CLoRF to HSI and PAN image fusion (HSI-
PAN fusion).

The structure of the remainder of this paper is outlined as follows. Sec. 2
provides a review of related work. Sec. 3 introduces the proposed CLoRF in
detail. Sec. 4 presents the experimental results and subsequent discussion,
and Sec. 5 concludes this work.

##  2 Related Work

###  2.1 Implicit Neural Representation

INR offers a novel approach to representing implicitly defined, continuous,
differentiable signals parameterized by neural networks cite ?, [29. INR has
demonstrated remarkable performance in representing complex data structures,
such as 3D reconstruction cite ?, [23, 29, 31 and 2D image super-resolution
cite ?, [12, 3 and generation cite ?, [41, 5, etc. Recently, INR-based
approaches have been explored for HSIs, such as HSI SR cite ?, [44, unmixing
cite ?, [35, and fusion cite ?, [37, 8. Despite these commendable efforts, INR
still encounters challenges in HSIs. For instance, HSIs consist of numerous
spectral bands obtained through continuous imaging within a specific spectral
range. However, INR itself may not possess sufficient stability to directly
learn a valid continuous representation from the spectral domain of HSIs. In
the fusion task, learning in both spatial and spectral domains is necessary;
therefore, utilizing a single INR for learning representation may confront
limited representation capabilities.

###  2.2 HSI-MSI Fusion via Continuous Representation

Several recent studies explore INRs for HSI-MSI fusion. For instance, ?
proposes spatial-INR and spectral-INR for spatial and spectral resolution
reconstruction, respectively. Besides, ? proposes an innovative fusion method
that integrates CNN and INR. Based on them, ? develops two spectral-spatial
INRs for arbitrary-resolution hyperspectral pansharpening. Although these
fusion methods utilize spatial-spectral-based INRs, they significantly diverge
from CLoRF. First, these methods rely on local implicit image functions cite
?, [5 for HSI SR. They employ a complex CNN network to encode spatial and
spectral features, then feed these features and the three-dimensional
coordinates of HSIs into an MLP to recover HSIs. Second, these HSI-MSI fusion
methods are supervised learning, which heavily relies on pairs of images for
training. Third, as they do not leverage the low-rank and smooth physical
information inherent in HSIs, the computational cost of these INRs is
intensive, primarily due to the large size of HSIs. Conversely, our method is
model-based and unsupervised, finely encoding the low-rankness and smoothness
into the continuous spatial-spectral factorization function. Consequently,
this improves the stability of the continuous representation in HSIs and
significantly reduces the complexity of the network structure and the
computational cost. Therefore, our method is more feasible and general in
various HSIs.

##  3 Proposed Method

In this section, we first present the problem formulation for HSI-MSI fusion,
then introduce the details of our proposed framework in the following section.

###  3.1 Problem Formulation

Given the HR-MSI and LR-HSI data, we aim to approximate their corresponding
HR-HSI data. Specifically, the HR-HSI, LR-HSI, and HR-MSI data are transformed
into the matrix format along the spectral dimension. The HR-HSI is denoted as
$\mathbf{Z}\in\mathbb{R}^{L\times N}$, where $L$ is the number of spectral
bands, and $N=H*W$ is the total number of pixels, in which $H$ and $W$
indicate spatial resolution. The LR-HSI is denoted as
$\mathbf{X}\in\mathbb{R}^{L\times n}$, where $n$ represents the number of LR
spatial pixels, i.e., $n\ll N$. Finally, the HR-MSI is denoted as
$\mathbf{Y}\in\mathbb{R}^{l\times N}$, where $l\ll L$ signifies that
$\mathbf{Y}$ has fewer spectral bands than $\mathbf{X}$.

The LR-HSI $\mathbf{X}$ can be interpreted as a diminished-quality
representation of HR-HSI $\mathbf{Z}$ in the spatial dimension, which is
formulated as follows:

$\mathbf{X}=\mathbf{Z}\mathbf{B}\mathbf{S}+\mathbf{N}_{h},$ (1)

where $\mathbf{N}_{h}\sim\mathcal{N}(\boldsymbol{0},\sigma_{h}\mathbf{I})$
represents the additive Gaussian noise. Besides,
$\mathbf{B}\in\mathbb{R}^{N\times N}$ is a spatial blurring operator of
$\mathbf{Z}$, representing the point spread function (PSF) of the
hyperspectral sensor. Additionally, $\mathbf{S}\in\mathbb{R}^{N\times n}$ is
the spatial downsampling matrix.

Similarly, the HR-MSI $\mathbf{Y}$ can be considered as a downsampled
realization of HR-HSI $\mathbf{Z}$ in the spectral dimension, which is
formulated as:

$\mathbf{Y}=\mathbf{H}\mathbf{Z}+\mathbf{N}_{m},$ (2)

where $\mathbf{H}\in\mathbb{R}^{l\times L}$ is the spectral response function
(SRF) and $\mathbf{N}_{m}\sim\mathcal{N}(\boldsymbol{0},\sigma_{m}\mathbf{I})$
denotes the additive Gaussian noise.

As HSIs generally have a low-rank structure, thus they lie in a low-
dimensional subspace cite ?, [28, 47. The low-rank factorization aims to
approximate a target matrix $\mathbf{Z}$ as a product of two matrices:

$\mathbf{Z}\approx\mathbf{E}\mathbf{A},$ (3)

where $\mathbf{E}\in\mathbb{R}^{L\times K}$ is a spectral dictionary and
$\mathbf{A}\in\mathbb{R}^{K\times N}$ is a coefficient matrix, respectively.
And $K\ll L$ represents a hyperparameter controlling the number of spectral
bases. The low-rank factorization representation offers three main advantages.
First, it maximizes the utilization of strong correlations among the spectral
bands. Second, by keeping $K$ small (where $K\ll L$), the size of the spectral
mode is reduced, thereby enhancing computational efficiency. Third, each row
of matrix $\mathbf{A}$ can be linearly represented by the rows (bands) of
matrix $\mathbf{Z}$, and the rows of matrix $\mathbf{A}$ maintain the spatial
structures of matrix $\mathbf{Z}$.

By integrating Eq.(3) into Eq.(1) and Eq.(2), $\mathbf{X}$ and $\mathbf{Y}$
are formulated as:

$\mathbf{X}=\mathbf{E}\mathbf{A}\mathbf{B}\mathbf{S}+\mathbf{N}_{h},\quad\mathbf{Y}=\mathbf{H}\mathbf{E}\mathbf{A}+\mathbf{N}_{m}.$
(4)

Thereafter, the fusion problem is transformed into the task of estimating the
spectral dictionary $\mathbf{E}$ and its corresponding coefficient matrix
$\mathbf{A}$ from matrices $\mathbf{X}$ and $\mathbf{Y}$, which follows the
following optimization problem:

$\min_{\mathbf{E},\mathbf{A}}\|\mathbf{X}-\mathbf{E}\mathbf{A}\mathbf{B}\mathbf{S}\|_{\text{F}}^{2}+\lambda\|\mathbf{Y}-\mathbf{H}\mathbf{E}\mathbf{A}\|_{\text{F}}^{2},$
(5)

where $\|\cdot\|_{\text{F}}$ denotes the Frobenius norm, and $\lambda$ denotes
the balancing factor. As there is a lack of specific prior information, the
problem Eq.(5) is undetermined; therefore, existing works focus on exploring
the appropriate prior information. Nonetheless, these methods process HR-HSI
within the dimensions of two modalities and fail to fuse the arbitrary
resolutions of HSIs effectively.

###  3.2 Deep Continuous Low-rank Factorization Model

INR is widely adopted for learning continuous data representation, such as
HSIs cite ?, [37, 8. However, simply utilizing a single INR to represent the
HSI volume results in low efficiency and expensive computation, as it neglects
the distinct low-rank structure of HSIs. Conversely, we propose a continuous
low-rank factorization (CLoRF) model for learning HSI representation
continuously and effectively. Our method fully explores the low-rank structure
of HSIs by simultaneously learning its low-rank continuous representation by
spatial and spectral INR. As a result, our approach effectively captures the
low-rankness and smoothness of HSIs while overcoming the computational burden
associated with existing INR-based methods cite ?, [37, 8 in HSIs.

As illustrated in Fig. 2, we present an overview of the proposed CLoRF.
Specifically, CLoRF consists of two INRs for learning the low-rank continuous
representation of HR-HSI. Inspired by the advance of Sinusoidal Representation
Networks (SIRENs) cite ?, [29, we employ two SIRENs to estimate $\mathbf{E}$
and $\mathbf{A}$ in Eq.(3), respectively. Specifically, we denote one SIREN
$\Psi_{\theta}(\cdot)$ parameterize by $\theta$ for approximating
$\mathbf{E}$, and another SIREN $\Phi_{\alpha}(\cdot)$ parameterized by
$\alpha$ for approximating $\mathbf{A}$, which is defined as follows:

$\begin{split}\mathbf{\hat{E}}(\mathbf{b};\theta)&=[\Psi_{\theta}(b_{1}),\Psi_{\theta}(b_{2}),\dots,\Psi_{\theta}(b_{L})]^{T},\\\
\mathbf{\hat{A}}(\mathbf{O};\alpha)&=[\Phi_{\alpha}(\mathbf{o}_{11}),\Phi_{\alpha}(\mathbf{o}_{12}),\dots,\Phi_{\alpha}(\mathbf{o}_{HW})],\end{split}$

where $\Psi_{\theta}(b_{i}):\mathbb{R}\rightarrow\mathbb{R}^{K}$ is a spectral
basis, with $b_{i}\in\mathbb{R}$ being the 1D coordinate for the $i$-th band
index of the LR-HSI. Besides,
$\Phi_{\alpha}(\mathbf{o}_{ij}):\mathbb{R}^{2}\rightarrow\mathbb{R}^{K}$ is a
spatial basis with the 2D coordinate $\mathbf{o}_{ij}\in\mathbb{R}^{2}$ of the
HR-MSI. And we denote the spectral bases as
$\mathbf{b}=[b_{1},b_{2},\dots,b_{L}]^{T}$ and the spatial bases as
$\mathbf{O}=[\mathbf{o}_{11};\mathbf{o}_{12};\dots;\mathbf{o}_{HW};].$ Both
networks aim to learn how to map from a fixed coordinate to the target
representation. Here, we formalize these networks as follows:

$\displaystyle\Psi_{\theta}(b_{i})=\mathbf{W}^{1}_{d_{1}}(\cdots(\sigma(\mathbf{W}^{1}_{1}b_{i}+\mathbf{c}^{1}_{1}))\cdots)+\mathbf{c}^{1}_{d_{1}},$
$\displaystyle\Phi_{\alpha}(\mathbf{o}_{ij})=\mathbf{W}^{2}_{d_{2}}(\cdots(\sigma(\mathbf{W}^{2}_{1}\mathbf{o}_{ij}^{T}+\mathbf{c}^{2}_{1})\cdots)+\mathbf{c}^{2}_{d_{2}},$

where $\sigma$ denotes the activation function,
$\theta=(\\{\mathbf{W}^{1}_{i}\\}_{i=1}^{d_{1}},\\{\mathbf{c}^{1}_{i}\\}_{i=1}^{d_{1}})$
and
$\alpha=(\\{\mathbf{W}^{2}_{i}\\}_{i=1}^{d_{2}},\\{\mathbf{c}^{2}_{i}\\}_{i=1}^{d_{2}})$
contains weight matrices and bias vectors for spectral- and spatial-INR,
respectively. Our method is a natural progression of low-rank factorization
from discrete mesh grids to the continuous domain. And the target HR-HSI is
approximated as
$\mathbf{\hat{Z}}=\mathbf{\hat{E}}(\mathbf{b};\theta)\mathbf{\hat{A}}(\mathbf{O};\alpha).$

As matrix $\mathbf{\hat{A}}(\mathbf{O};\alpha)$ preserves the spatial
structures of HSIs. Without loss of generality, we consider the spatial
smoothness of HSIs. Moreover, a total variation (TV) loss on the predicted
coefficient matrix $\mathbf{\hat{A}}(\mathbf{O};\alpha)$ is further
incorporated for noise-disruption scenarios. Mathematically, the TV
regularization of $\mathbf{\hat{A}}(\mathbf{O};\alpha)$ is formulated as:

$\sum_{k=1}^{K}\text{TV}(\mathbf{\hat{a}}_{k})=\sum_{k=1}^{K}(\|\mathbf{D}_{h}{\mathbf{\hat{a}}}_{k}\|_{1}+\|\mathbf{D}_{w}{\mathbf{\hat{a}}}_{k}\|_{1}),$
(6)

where $\hat{\mathbf{a}}_{k}$ is the $k$-th row of
$\mathbf{\hat{A}}(\mathbf{O};\alpha)$. $\mathbf{D}_{h}$ and $\mathbf{D}_{w}$
denote the differential operation along the height and width direction in the
matrix form of $\hat{\mathbf{a}}_{k}$, respectively. Here, $\|\cdot\|_{1}$
indicates the $\ell_{1}$ norm. By incorporating the TV loss, we promote
spatial smoothness and improve the overall quality of the fusion.

Therefore, the optimization problem with TV prior can be summarized as
follows:

$\min_{\theta,\alpha}\mathcal{L}_{\text{MSI}}+\lambda\mathcal{L}_{\text{HSI}}+\eta\sum\limits_{k=1}^{K}\text{TV}(\mathbf{\hat{a}}_{k}),$
(7)

where
$\mathcal{L}_{\text{MSI}}=\|\mathbf{X}-\mathbf{\hat{E}}(\mathbf{b};\theta)\mathbf{\hat{A}}(\mathbf{O};\alpha)\mathbf{BS}\|_{\text{F}}^{2},$
$\mathcal{L}_{\text{HSI}}=\|\mathbf{Y}-\mathbf{H}\mathbf{\hat{E}}(\mathbf{b};\theta)\mathbf{\hat{A}}(\mathbf{O};\alpha)\|_{\text{F}}^{2}$,
and $\eta$ is the regularization parameter.

Approximating $\mathbf{\hat{E}}(\mathbf{b};\theta^{\ast})$ and
$\mathbf{\hat{A}}(\mathbf{O};\alpha^{\ast})$ to maintain the low-rank
representation of $\mathbf{\hat{Z}}$ can be achieved after training networks,
with $\theta^{\ast},\alpha^{\ast}$ corresponding to the parameters of well-
trained networks. Recall the networks take the coordinates as the input, and
the optimization in Eq.(7) does not involve the ground-truth HSIs as the
supervision label; therefore, our method is self-supervised. We employ the
Adam optimizer for optimization, which is a stochastic gradient descent
algorithm. Moreover, we can infer an arbitrary-resolution HSI by inputting any
scale coordinates $\\{\tilde{\mathbf{b}},\tilde{\mathbf{O}}\\}$ into the well-
trained network, i.e.,
$\mathbf{\hat{E}}(\tilde{\mathbf{b}};\theta^{\ast})\mathbf{\hat{A}}(\tilde{\mathbf{O}};\alpha^{\ast})$.

Compared to existing INRs-based fusion methods, our method stands out for
several advantages. First, it fully exploits the low-rank and smooth prior of
HSIs through low-rank continuous learning representations. Second, its
computational complexity is significantly reduced through continuous low-rank
factorization. Third, it can achieve the user-desired resolution at arbitrary
locations in HSIs by inputting any scale spatial and spectral coordinates.

###  3.3 Theoretical Analysis

In this section, we theoretically demonstrate that the low-rank and smooth
regularizations are implicitly unified in continuous low-rank matrix
factorization. Our analysis is inspired by the concept of tensor function
factorization in cite ?, [22. Here, we start with rank factorization in the
matrix computation field.

######  Theorem 1 (rank factorization cite ?, [25).

Let $\mathbf{X}\in\mathbb{R}^{n_{1}\times n_{2}}$, where
$\mathrm{rank}(\mathbf{X})=K$, then there exists two matrices
$\mathbf{U}\in\mathbb{R}^{n_{1}\times K}$,
$\mathbf{V}\in\mathbb{R}^{n_{2}\times K}$ such that
$\mathbf{X}=\mathbf{U}\mathbf{V}^{T}$.

Subsequently, we provide a detailed introduction to the proposed continuous
representation of HSIs. Let
$f(\cdot):\mathcal{A}_{f}\times\mathcal{Z}_{f}\rightarrow\mathbb{R}$ be a
bounded real function, where $\mathcal{A}_{f}\subset\mathbb{R}^{2}$,
$\mathcal{Z}_{f}\subset\mathbb{R}$ are definition domains in spatial and
spectral domains, respectively. The function $f$ gives the value of data at
any coordinate in $\mathcal{D}_{f}:=\mathcal{A}_{f}\times\mathcal{Z}_{f}$. We
interpret $f$ as a matrix function since it maps a spatial and spectral
coordinate to the corresponding value, implicitly representing matrix data.

######  Definition 1 (sampled matrix set).

For a matrix function $f(\cdot):\mathcal{D}_{f}\to\mathbb{R}$, we define the
sampled matrix set $\mathcal{S}[f]$ as

$\begin{split}\mathcal{S}[f]:=\\{&\mathbf{M}|\mathbf{M}_{(i,j)}=f(\mathbf{s}_{i},b_{j}),\mathbf{s}_{i}\in\mathcal{A}_{f},\\\
&{b}_{j}\in\mathcal{Z}_{f},\mathbf{M}\in\mathbb{R}^{n_{1}\times
n_{2}},n_{1},n_{2}\in\mathbb{N}_{+}\\},\end{split}$

where $\mathbf{s}_{i},b_{j}$ denote the spatial and spectral coordinates,
respectively.

Regarding the definition of rank, we expect any matrix sampled on
$\mathcal{S}[f]$ to be low-rank. Naturally, we can then define the rank of the
matrix function as follows.

######  Definition 2 (matrix function rank).

Given a matrix function
$f:\mathcal{D}_{f}=\mathcal{A}_{f}\times\mathcal{Z}_{f}\to\mathbb{R}$, we
define a measure of its complexity, denoted by
$\mathrm{MF}\text{-}\mathrm{rank}[f]$ (function rank of $f(\cdot)$), as the
supremum of the matrix rank in the sampled matrix set $\mathcal{S}[f]{:}$

$\mathrm{MF}\text{-}\mathrm{rank}[f]:=\sup_{\mathbf{M}\in\mathcal{S}[f]}\mathrm{rank}(\mathbf{M}).$

We call a matrix function $f(\cdot)$ as a low-rank matrix function if
$K\ll\min\\{n_{1},n_{2}\\}$. When $f(\cdot)$ is defined on a given matrix, we
will show that the MF-rank degenerates into the discrete case, i.e., the
classical matrix rank.

######  Proposition 2.

Consider $\mathbf{X}\in\mathbb{R}^{n_{1}\times n_{2}}$ as an arbitrary matrix.
Let $\mathcal{A}_{f}=\mathcal{N}^{(l_{1})}\times\mathcal{N}^{(l_{2})}\
(l_{1}l_{2}=n_{1})$ represent a two-dimensional discrete set, and
$\mathcal{Z}_{f}=\mathcal{N}^{(n_{2})}$ is an one-dimensional discrete set,
where $\mathcal{N}^{(k)}$ is the set $\\{1,2,\dots,k\\}$. We denote
$\mathcal{D}_{f}=\mathcal{A}_{f}\times\mathcal{Z}_{f}$ and define the matrix
function $f(\cdot):\mathcal{D}_{f}\rightarrow\mathbb{R}$ as
$f(\mathbf{s},b)=\mathbf{X}_{(\mathbf{s},b)}$ for any
$(\mathbf{s},b)\in\mathcal{D}_{f}$. Consequently,
$\mathrm{MF}$-$\mathrm{rank}[f]=\mathrm{rank}(\mathbf{X})$.

The proof of Proposition 2 is shown in Appendix A. Proposition 2 expands the
concept of rank from discrete matrices to matrix functions for continuous
representations. Analogous to classical matrix representations, it is
pertinent to consider whether a low-rank matrix function $f$ can employ
certain matrix factorization strategies to encode its low-rank. We provide an
affirmative response, as stated in the theorem below.

######  Theorem 3 (continuous low-rank factorization).

Let
$f(\cdot):\mathcal{D}_{f}=\mathcal{A}_{f}\times\mathcal{Z}_{f}\to\mathbb{R}$
be a bounded matrix function, where $\mathcal{A}_{f}\subset\mathbb{R}^{2}$,
$\mathcal{Z}_{f}\subset\mathbb{R}$. If $\mathrm{MF}$-$\mathrm{rank}[f]=K$,
then there exist two functions
$f_{\mathrm{spatial}}(\cdot):\mathcal{A}_{f}\rightarrow\mathbb{R}^{K}$,
$f_{\mathrm{spectral}}(\cdot):\mathcal{Z}_{f}\rightarrow\mathbb{R}^{K}$ such
that $f(\mathbf{s},b)=f_{\mathrm{spatial}}(\mathbf{s})\cdot
f^{T}_{\mathrm{spectral}}(b)$ for any pair of inputs
$(\mathbf{s},b)\in\mathcal{D}_{f}$.

The proof of Theorem 3 is illustrated in Appendix B. Theorem 3 is a natural
extension of rank factorization (Theorem 1) from discrete meshgrid to the
continuous domain. Specifically, we employ two MLPs $\Phi_{\alpha}(\cdot)$ and
$\Psi_{\theta}(\cdot)$ with parameters $\theta$ and $\alpha$ to parameterize
the factor functions $f_{\mathrm{spatial}}(\cdot)$ and
$f_{\mathrm{spectral}}(\cdot)$.

Smoothness is another prevalent attribute in HSIs, such as the spatial and
spectral smoothness of HSIs cite ?, [30. Here, we theoretically validate that
our method incorporates implicit smooth regularization derived from the
specific structures of MLPs.

######  Theorem 4 (Lipschitz continuity).

Let $\mathbf{X}\in\mathbb{R}^{n_{1}\times n_{2}}$, and
$\Phi_{\alpha}(\cdot):\mathcal{A}_{f}\rightarrow\mathbb{R}^{K}$,
$\Psi_{\theta}(\cdot):\mathcal{D}_{f}\rightarrow\mathbb{R}^{K}$ be two MLPs
structured with parameters $\alpha,\theta$ where
$\mathcal{A}_{f}\subset\mathbb{R}^{2}$, $\mathcal{Z}_{f}\subset\mathbb{R}$ are
bounded, i.e., $\|\mathbf{s}\|_{1}\leq\zeta,$ $b\leq\zeta$ for any
$\mathbf{s}\in\mathcal{A}_{f}$, $b\in\mathcal{Z}_{f}$. Suppose the MLPs share
the same activation function $\sigma(\cdot)$ and depth $d$ with
$\mathbf{c}^{1}_{i}=\mathbf{c}^{2}_{i}=\mathbf{0},\forall\ i$. Besides, we
assume that

  * •

$\sigma$ is Lipschitz continuous with the Lipschitz constant $\kappa$, and
$\sigma(0)=0$;

  * •

$\|\mathbf{X}\|_{1}$ is bounded by a positive constant $\eta_{1}$;

  * •

$\|\mathbf{W}_{i}^{1}\|_{1},\|\mathbf{W}_{i}^{2}\|_{1}$ are bounded by a
positive constant $\eta_{2}$ for all $i$.

Define a matrix function
$f(\cdot):\mathcal{D}_{f}=\mathcal{A}_{f}\times\mathcal{Z}_{f}\rightarrow\mathbb{R}$
as $f(\mathbf{s},b)=\Phi_{\alpha}(\mathbf{s})\cdot\Psi_{\theta}(b)^{T}$. Then,
the following inequalities hold for any
$(\mathbf{s}_{1},b_{1}),(\mathbf{s}_{2},b_{2})\in\mathcal{D}_{f}$:

$\displaystyle|f(\mathbf{s}_{1},b_{1})-f(\mathbf{s}_{2},b_{2})|$
$\displaystyle\leq\delta\|\mathbf{s}_{1}-\mathbf{s}_{2}\|_{1}+\delta|b_{1}-b_{2}|,$

where $\delta=\eta^{2d+1}\kappa^{2d-2}\zeta$, and
$\eta=\max\\{\eta_{1},\eta_{2}\\}$.

The proof is shown in Appendix C. This smoothness is implicitly encoded with
mild assumptions regarding nonlinear activation functions and weight matrices,
which are readily attainable in a real-world implementation. For instance,
commonly used activation functions such as ReLU, LeakyReLU, Sine, and Tanh are
Lipschitz continuous. Additionally, cite ?, [22 effectively governs the level
of smoothness by managing the structures of MLPs.

##  4 Experiments and Analysis

In this section, we evaluate the performance of our method on five datasets
separately. Additionally, we compare several SOTA methods and evaluate the
fusion results qualitatively and quantitatively. Finally, we expand the
application of CLoRF to HSI-PAN fusion and verify its efficacy on a dataset.

###  4.1 Experimental Details

1) Datasets: Five simulated datasets, including Pavia University cite ?, [7,
Pavia Center cite ?, [7, Indian Pines cite ?, [4, Washington DC cite ?, [48,
47 and University of Houston cite ?, [20 were used for simulations in our
experiments. The five synthetic datasets, each with a simulated PSF and SRF,
are used to generate two observed images: LR-HSI and HR-MSI. More
specifically, we utilize a Gaussian blur of $5\times 5$ pixels with a 0 mean
and 1 standard deviation to simulate the PSF for generating the LR-HSI. The
downsampling ratio was set to 4 for all datasets. The spectral response of the
IKONOS satellite (a Nikon D700 camera) was used to simulate the SRF for
generating the HR-MSI. The i.i.d Gaussian noise was added to HR-MSI and LR-HSI
with signal-to-noise ratios (SNRs) of 30dB, respectively. The details for the
five datasets are summarized as follows.

  * •

Pavia University: The image measures $610\times 340\times 115$ pixels with a
spatial resolution of 1.3 meters and spectral coverage spanning from 0.43 $\mu
m$ to 0.86 $\mu m$. Due to the effects of noise and water vapor absorption, 12
bands were removed. An area covering $336\times 336$ pixels in the lower-left
corner of the image and containing 103 bands was selected for the experiment.

  * •

Pavia Center: The size of the Pavia Center is $1096\times 1096\times 115$ and
spectral ranging from 0.38 to 1.05 $\mu m$. After removing bands caused by
water vapor absorption and low SNRs, the subregion consisting of $336\times
336\times 93$ pixels were chosen from the full dataset as a reference HR-HSI.

  * •

Indian Pines: The Indian Pines has 224 spectral bands with a size of
$145\times 145$ pixels. the spectral wavelength range is from 0.4 to 2.5 $\mu
m$. After removing bands caused by water vapor absorption and low SNRs, the
subregion consisting of $144\times 144\times 191$ pixels was chosen from the
full dataset as a reference HR-HSI.

  * •

Washington DC Mall: The Washington DC Mall has a region of $1280\times 307$
pixels, and the image consists of 210 bands with spectral wavelength ranging
from 0.4 to 2.5 $\mu m$. After removing the low SNRs and atmospheric
absorption bands, 191 bands were kept. The subregion consisting of $304\times
304\times 191$ pixels was clipped from the full dataset as a reference HR-HSI.

  * •

University of Houston: The University of Houston contains $601\times 2384$
pixels and 48 bands ranging from 0.38 to 1.05 $\mu m$. The sub-image
consisting of $320\times 320\times 46$ pixels was chosen from the whole
dataset as a reference HR-HSI for our experiments.

2) Comparison methods: Since our method is an unsupervised algorithm, we
mainly compare it to unsupervised fusion methods; such a choice also aligns
with the practical demands of real-world scenarios. Specifically, we compare
the proposed method with several SOTA unsupervised fusion methods, including
CNMF cite ?, [43, HySure cite ?, [28, HyCoNet cite ?, [46, CNN-FUS cite ?,
[10, MSE-SURE cite ?, [24, and E2E-FUS cite ?, [39. For baselines: CNMF,
HySure, CNN-FUS, and E2E-FUS methods are implemented using MATLAB R2023a on
Windows 11 with 32GB RAM, while HyCoNet, MSE-SURE, and CLoRF are evaluated on
an RTX 4090 GPU with 32GB RAM. Besides, parameters in baselines are manually
fine-tuned to achieve optimal results in all experiments.

3) Evaluation metrics: Four distinct quantitative metrics are employed to
assess the efficacy of fusion results, where ground truth (GT) is given. These
metrics include the mean structured similarity (MSSIM), an extension of SSIM
to assess HSI quality by averaging across all spectral bands; mean peak
signal-to-noise ratio (MPSNR), which is calculated as the average PSNR across
all bands extended for HSI; the spectral angle mapper (SAM) index; and the
relative global dimension error (ERGAS) index.

4) Hyperparameters of CLoRF: For all tasks, both the spatial-INR and spectral-
INR adopt the SIREN network with an initial network parameter $\omega_{0}=30$.
The spatial-INR has 5 hidden layers with a size of 512, while the spectral-INR
has 2 hidden layers with a size of 128 (In Indian Pines, the spatial-MLP has 3
hidden layers with a size of 512, while the spectral-MLP also has 3 hidden
layers with a size of 256). The learning rate is set to 3e-5, and the training
epochs are fixed to 30000. An early stopping strategy is implemented to
prevent overfitting for all datasets. Furthermore, the hyperparameters are
summarized in Table 1.

Table 1: Hyperparameters used for training the proposed model. Hyperparameter | Pavia University | Pavia Center | Indian Pines | Washington DC | University of Houston
---|---|---|---|---|---
$K$ | 9 | 9 | 22 | 11 | 10
$\lambda$ | 1.25 | 1.25 | 0.55 | 1.25 | 1
$\eta$ | 0.0025 | 0.0050 | 0.0060 | 0.0025 | 0.0025
Table 2: Quantitative performance comparison with different algorithms on the different datasets. The best results are bold-faced, and runner-ups are underlined. (MPSNR ↑, MSSIM ↑, SAM ↓, ERGAS ↓) | Metric | CNMF | HySure | HyCoNet | CNN-FUS | MSE-SURE | E2E-FUS | CLoRF
---|---|---|---|---|---|---|---|---
Pavia University | MPSNR | 32.70 | 32.86 | 40.10 | 41.57 | 41.31 | 40.92 | 42.23
MSSIM | 0.92 | 0.93 | 0.97 | 0.98 | 0.98 | 0.98 | 0.99
SAM | 3.55 | 5.81 | 3.10 | 2.38 | 2.20 | 2.31 | 2.05
ERGAS | 3.63 | 3.71 | 1.41 | 1.38 | 1.44 | 1.50 | 1.31
Pavia Center | MPSNR | 38.97 | 33.96 | 41.88 | 42.16 | 42.87 | 43.12 | 43.47
MSSIM | 0.97 | 0.93 | 0.98 | 0.98 | 0.99 | 0.99 | 0.99
SAM | 7.70 | 8.42 | 3.77 | 4.10 | 3.51 | 3.67 | 3.64
ERGAS | 2.35 | 3.56 | 1.45 | 1.47 | 1.37 | 1.33 | 1.26
Indian Pines | MPSNR | 27.48 | 26.54 | 29.44 | 31.07 | 30.85 | 31.68 | 31.96
MSSIM | 0.94 | 0.93 | 0.93 | 0.94 | 0.95 | 0.95 | 0.95
SAM | 3.11 | 3.97 | 2.87 | 2.70 | 2.54 | 2.49 | 2.56
ERGAS | 2.52 | 2.84 | 2.10 | 1.78 | 1.93 | 1.73 | 1.70
Washington DC | MPSNR | 28.47 | 26.37 | 33.35 | 32.24 | 33.52 | 36.15 | 36.75
MSSIM | 0.94 | 0.94 | 0.98 | 0.96 | 0.98 | 0.98 | 0.99
SAM | 4.94 | 7.67 | 3.51 | 4.07 | 2.56 | 2.56 | 2.64
ERGAS | 3.22 | 4.79 | 2.67 | 2.37 | 1.98 | 1.53 | 1.44
University of Houston | MPSNR | 31.71 | 29.41 | 37.02 | 37.38 | 39.75 | 38.84 | 40.55
MSSIM | 0.95 | 0.92 | 0.97 | 0.97 | 0.99 | 0.98 | 0.99
SAM | 2.77 | 4.77 | 2.50 | 2.46 | 1.61 | 1.76 | 1.47
ERGAS | 2.10 | 3.13 | 1.19 | 1.31 | 0.90 | 1.06 | 0.83

Figure 3: The first row shows the Pavia University image (1st band) of the
estimated HR-HSI, and the second row shows the error map between the estimated
image and GT.

Figure 4: The first row shows the Pavia Center image (8th band) of the
estimated HR-HSI, and the second row shows the error map between the estimated
image and GT.

Figure 5: The first row shows the Indian Pines image (160th band) of the
estimated HR-HSI, and the second row shows the error map between the estimated
image and GT.

Figure 6: The first row shows the Washington DC image (5th band) of the
estimated HR-HSI, and the second row shows the error map between the estimated
image and GT.

Figure 7: The first row shows the University of Houston image (46th band) of
the estimated HR-HSI, and the second row shows the error map between the
estimated image and GT.

Figure 8: Reconstructed spectral signatures of two randomly selected locations
at different datasets. From left to right, they are Pavia University, Indian
Pines, Washington DC, and University of Houston. To provide a clearer
comparison, we have compared CLoRF alongside E2E-FUS and GT.

###  4.2 Experimental Results

The fusion performance of all methods on five datasets, evaluated across
MPSNR, MSSIM, SAM, and ERGAS metrics, is shown in Table 2. CLoRF consistently
outperforms others in terms of quantitative results across most scenarios,
affirming the efficacy of the HSI-MSI fusion task. Moreover, CLoRF exhibits
commendable performance in spatial structures. our method outperforms E2E-FUS
by 2.19% in terms of MPSNR. Nonetheless, in certain datasets, particularly
those with a high number of bands, CLoRF slightly lags behind MSE-SURE in SAM
index performance. It is notable that MSE-SURE, throughout its training phase,
is equipped with knowledge of noise levels in LR-HSI and HR-MSI, and leverages
back-projection techniques to enhance spectral detail capture, while our
method does not require this information.

To assess the effectiveness of various methods in preserving spatial
structures, Fig. 3-7 illustrates estimated HR-HSI on five simulated datasets.
The first rows exhibit fusion results virtually, showcasing the 1st band of
Pavia University, the 8th band of Pavia Center, the 160th band of Indian
Pines, the 5th band of Washington DC Mall, and the 46th band of the University
of Houston. Conversely, the second rows depict corresponding error maps,
illustrating the mean squared error (MSE) between GT and the estimated HR-HSI.
While the results for the images in the first rows appear similar across most
methods, with differences almost imperceptible to the naked eye, notable
disparities emerge in the error maps depicted in the second column. It is
evident that CLoRF exhibits minimal noise and preserves superior spatial
structures. As shown in Fig. 8, we plot the two randomly selected spectral
vectors reconstructed by CLoRF and E2E-FUS with four datasets. It is apparent
that CLoRF effectively preserves the high-frequency information compared to
E2E-FUS in some pixels. However, in some specific pixels, the performance of
the E2E-FUS method is superior to CLoRF.

###  4.3 Arbitrary Resolution

Here, we evaluate the performance of CLoRF to fuse HSIs at arbitrary spatial
and spectral resolutions. As shown in Fig. 1, our model is well-trained based
in an unsupervised manner, and we can input any scale spatial and spectral
position coordinates during the inference stage, obtaining arbitrary
resolutions in both spatial and spectral domains. We use PSF and downsampling
ratio of the same size as in Sec. 4.1, but different sizes of SRF (sampled
from SRF in Sec. 4.1) to synthesize the reduced-resolution LR-HSI ($42\times
42\times 50$) and HR-MSI ($168\times 168\times 4$) with Pavia University for
training spatial-INR and spectral-INR. Then, we predict HR-HSIs of arbitrary
resolutions using different scales of spatial and spectral coordinates. For a
fair comparison, our method is directly compared with bicubic interpolation
from the original LR. As shown in Table 3, we can fix spectral resolution to
obtain arbitrary spatial resolution, fix spatial resolution to obtain
arbitrary spectral resolution, or simultaneously achieve resolutions in both
spatial and spectral domains. CLoRF performs nearly as well as bicubic
interpolation in spatial resolution when spectral resolution is fixed, yet it
notably surpasses bicubic when spatial resolution is fixed while achieving
spectral resolution. Moreover, when simultaneously enhancing resolutions in
both spatial and spectral domains, CLoRF consistently outperforms bicubic
interpolation. From Fig. 9-10, it can be seen that when simultaneously
upsampling in both spatial and spectral domains, CLoRF obtains finer details
compared to the bicubic interpolation.

Figure 9: An example of CLoRF for super-resolution on the Pavia Unversity with
a spatial resolution of $336\times 336$ (60th band) and its corresponding
error map while simultaneously achieving arbitrary resolutions in both spatial
and spectral domains. Figure 10: Visualize the spectral signatures of
different pixels in different spectral resolutions while simultaneously
obtaining arbitrary resolutions in both spatial and spectral domains. Table 3:
Experimental results in spatial (fixed spectral), spectral (fixed spatial),
and (spatial, spectral) with arbitrary resolution, metric: MPSNR/MSSIM. The
original dimension of HR-HSI is (168, 168, 50). The best are highlighted.
Arbitrary resolution in the spectral domain
---
Dimension | (168,168,61) | (168,168,72) | (168,168, 83) | (168,168,93)
Bicubic | 37.64 /0.97 | 36.60/0.95 | 38.76/0.96 | 37.23/0.95
CLoRF | 42.42/0.99 | 42.85/0.98 | 42.64/0.98 | 42.56/0.98
Arbitrary resolution in the spatial domain
Dimension | (210,210,50) | (252,252,50) | (294,294,50) | (336,336,50)
Bicubic | 30.07/0.90 | 29.68/0.89 | 30.24/0.89 | 29.29/0.88
CLoRF | 29.86/0.88 | 29.52/0.87 | 30.04/0.86 | 29.22/0.86
Arbitrary resolution in the spectral and spatial domain
Dimension | (210,210,61) | (252,252,72) | (294,294,83) | (336,336,93)
Bicubic | 28.87/0.87 | 28.00/0.85 | 28.56/0.85 | 27.48/0.83
CLoRF | 29.53/0.88 | 29.21/0.87 | 29.60/0.86 | 28.61/0.84

###  4.4 Ablation Study

We conduct a series of ablation studies to verify the effectiveness of CLoRF.

1) Loss function: As shown in Table 4, incorporating TV loss into the loss
function can lead to better recovery quality. On each dataset, the TV loss is
capable of enhancing PSNR by 1-2 dB.

Table 4: Ablation study of the TV loss. | Metric | w/o | w/
---|---|---|---
Pavia University | MPSNR | 41.48 | 42.23
MSSIM | 0.98 | 0.99
University of Houston | MPSNR | 39.04 | 40.55
MSSIM | 0.98 | 0.99
Table 5: Ablation study of the different activation functions. (MPSNR / MSSIM) Act.Func | Pavia University | Houston
---|---|---
ReLU | 25.92/0.65 | 24.58/0.54
ReLU+PE | 36.26/0.93 | 37.56/0.98
Gauss | 35.96/0.90 | 36.10/0.96
Spder | 41.63/0.98 | 39.16/0.98
Sine | 42.23/0.99 | 40.55/0.99

2) Different activation functions: Table 5 displays the results obtained using
different activation functions, such as ReLU, ReLU+Position Encoding (PE),
Gauss cite ?, [26, Spder cite ?, [27, and Sine cite ?, [29. Notably, the Sine
activation function demonstrates outstanding performance across all metrics.
Therefore, we have made it the default activation function for our model.

3) Impact of hyper-parameters: We examine six hyperparameters: $K$, $\lambda$,
$\eta$, the hidden depth of two MLPs, and the learning rate (LR). Each
hyperparameter is explored within a specified range while the others are held
constant. Specifically, for Pavia University, we explore $K$ within the range
of 5 to 17. For $\lambda$, we consider values from the set
$\\{0.5,0.75,1,1.25,1.5,1.75\\}$, For $\eta$, the search space consists of
$\\{10^{-3},2.5\times 10^{-3},5\times 10^{-3},7.5\times 10^{-3},10^{-2}\\}$.
Regarding the hidden depth of the two MLPs, we initially set the hidden depth
of the spectral MLP at 2 and explored the depth of the spatial MLP from 1 to
7. Similarly, we fix the hidden depth of the spatial MLP at 5 and then search
for the hidden depth of the spectral MLP from 1 to 6. For LR, we consider
values from the set $\\{10^{3},5\times 10^{3},10^{4},5\times
10^{4},10^{5},3\times 10^{5},5\times 10^{5}\\}$. The results are displayed in
Fig. 11.

Figure 11: Impact of hyperparameters on MPSNR and MSSIM. The red dot
represents the best result obtained by traversing all values.

###  4.5 PAN-HSI Fusion

In this section, we extend our proposed fusion method to the PAN-HSI fusion.
PAN image, compared to MSI, has fewer bands, making fusion more challenging.
Below are the results and details of some experiments we conducted.

1) Dataset: The dataset acquired for the PRISMA contest1 cite ?, [33, namely
RR1. The image is of size $900\times 900\times 59$. The synthetic datasets,
plus a simulated PSF and SRF, were used to generate two observed images: LR-
HSI and PAN. PSF is set to be the same as the 4.1. The wavelength cite ?, [33
was used to simulate the SRF to generate the PAN image. The i.i.d Gaussian
noise is added to the LR-HSI (30dB) and PAN (30dB) image.

2) Compared Methods: We primarily compare CLoRF with two categories of
baselines: model-based methods include GS cite ?, [19, GSA cite ?, [2, AWLP
cite ?, [32, MTF-GLP cite ?, [1, and HySure cite ?, [28, while the
unsupervised deep learning method is R-PNN cite ?, [14. All training
parameters for CLoRF remain the same for Pavia University.

As shown in Table 6, CLoRF outperforms the other methods in terms of MPSNR and
ERGAS. However, due to spectral distortion in the PAN image, CLoRF fails to
learn a continuous representation in the spectral domain effectively. For
clarity, we depict error maps for two scenarios: noise-free and noisy, as
shown in Fig. 12. The results indicate that CLoRF outperforms other methods in
spatial domains.

Table 6: Quantitative performance comparison with different algorithms on the RR1 dataset. The best results are bold-faced, and runner-ups are underlined. (MPSNR ↑, MSSIM ↑, SAM ↓, ERGAS ↓) | Metric | GS | GSA | AWLP | MTF-GLP | HySure | R-PNN | CLoRF
---|---|---|---|---|---|---|---|---
Noise-free | MPSNR | 27.40 | 32.25 | 29.64 | 30.13 | 29.56 | 30.05 | 32.87
MSSIM | 0.84 | 0.92 | 0.90 | 0.90 | 0.88 | 0.90 | 0.91
SAM | 9.70 | 4.46 | 4.90 | 4.41 | 6.49 | 4.78 | 5.81
ERGAS | 5.27 | 3.06 | 3.92 | 3.75 | 4.27 | 4.28 | 2.95
| Metric | GS | GSA | AWLP | MTF-GLP | HySure | R-PNN | CLoRF
Noisy | MPSNR | 27.25 | 31.90 | 29.23 | 29.74 | 29.44 | 29.33 | 32.61
MSSIM | 0.82 | 0.90 | 0.87 | 0.88 | 0.88 | 0.88 | 0.89
SAM | 10.05 | 5.51 | 5.96 | 5.76 | 7.07 | 6.05 | 6.50
ERGAS | 5.35 | 3.16 | 4.08 | 3.90 | 4.12 | 4.60 | 2.99

Figure 12: The first and second rows show the error images (Noise-free and
Noisy) between the estimated HR-HSI and GT.

##  5 Conclusion

In this work, we introduced an innovative continuous low-rank factorization
representation for HSI-MSI fusion, which incorporates two INRs into the
decomposition to capture spatial and spectral information, respectively.
Theoretical analysis reveals that this continuous function representation
adeptly portrays the low-rank and smoothness priors of HSIs. Extensive
numerical experiments conducted for HSI-MSI fusion confirm its effectiveness
and wide applicability. Nonetheless, our method still faces certain
limitations. Given its unsupervised manner, CLoRF demands an extensive number
of training epochs. Another limitation is that we solely utilize SIREN, which
is widely used in INR. There’s a need to explore diverse continuous
representations for both spatial and spectral domains. The adoption of
continuous low-rank factorization representation for processing and analyzing
HSIs shows potential for future applications across various tasks, e.g., HSI
unmixing and single RGB-HSI super-resolution. Our future efforts will be
dedicated to expanding the versatility of CLoRF to address diverse HSI tasks.

Acknowledgements

This work was partly supported by the Natural Science Foundation of China
(Nos. 12201286, 52303301), the Shenzhen Science and Technology Program
(20231115165836001), HKRGC Grant No.CityU11301120, CityU Grant No. 9229120,
HKRGC GRF 17201020 and 17300021, HKRGC C7004-21GF, and Joint NSFC and RGC
N-HKU769/21, National Key R&D Program of China (2023YFA1011400), and the
Shenzhen Fundamental Research Program (JCYJ20220818100602005).

Data Availability Data will be made available on reasonable request.

##  Appendix A Proof of Proposition 1

Since $\mathbf{X}\in S[f]$, we directly get MF-
rank$[f]\geq\mathrm{rank}({\mathbf{X}})$. Now, we aim to prove the other side:
MF-rank$[f]\leq\mathrm{rank}({\mathbf{X}})$.

Let $\mathbf{M}$ be any matrix within the set $S[f]$. Each column vector of
$\mathbf{M}$ is denoted by $\mathbf{M}_{(:,p)}$ for
$p\in\\{1,2,\ldots,n_{2}\\}$. According to the definition of $S[f]$, there
exists an index $l_{p}\in\\{1,2,\ldots,n_{1}\\}$ dependent on $p$ such that
$\mathbf{M}_{(:,p)}$ is a permutation of the elements in
$\mathbf{X}_{(:,l_{p})}$, allowing for repeated sampling. In other words, for
each $\mathbf{M}_{(:,p)}$, there exists a permutation matrix
$\mathbf{P}\in\\{0,1\\}^{n_{1}\times n_{1}}$ and a corresponding column of
$\mathbf{X}$ depending on $p$ (specifically $\mathbf{X}_{(:,l_{p})}$), such
that $\mathbf{M}_{(:,p)}=\mathbf{P}\mathbf{X}_{(:,l_{p})}$. Additionally, the
permutation matrix $\mathbf{P}$ is consistent across all columns
$\mathbf{M}_{(:,p)}$ for $p=1,2,\ldots,n_{2}$, i.e.,
$\mathbf{M}_{(:,p)}=\mathbf{P}\mathbf{X}_{(:,l_{p})}$ for each $p$.

Define
$\tilde{\mathbf{X}}:=[\mathbf{X}_{(:,l_{1})},\mathbf{X}_{(:,l_{2})},\ldots,\mathbf{X}_{(:,l_{n_{2}})}]\in\mathbb{R}^{m_{1}\times
n_{2}}$, we have that the rank of $\tilde{\mathbf{X}}$ is less than or equal
to the rank of $\mathbf{X}$:
$\mathrm{rank}(\tilde{\mathbf{X}})\leq\mathrm{rank}(\mathbf{X})$. Finally,
since $\mathbf{M}=\mathbf{P}\tilde{\mathbf{X}}$, it follows that
$\mathrm{rank}(\mathbf{M})\leq\mathrm{rank}(\mathbf{X})$, which leads to MF-
rank$[f]\leq\mathrm{rank}({\mathbf{X}})$.

##  Appendix B Proof of Theorem 3

First, we establish a linear representation for each factor function
$f(\mathbf{s},b)$ with fixed inputs $\mathbf{s}$ and $b$, employing a set of
basis functions. Then, we focus on presenting the continuous low-rank
factorization and demonstrating that this factorization preserves the matrix
factorization rank (MF-rank).

Suppose that MF-rank$[f]=K$ with $K<\infty$, thus there exist a matrix
$\mathbf{M}\in\mathbb{R}^{n_{1}\times K}$ with rank$(M)=K$. Denote:
$\mathcal{S}=\\{\mathbf{s}_{i}\ |\
\mathbf{M}_{i,j}=f(\mathbf{s}_{i},b_{j}),i=1,\dots,n_{1}\\},$ and
$\mathcal{T}=\\{b_{j}\ |\
\mathbf{M}_{i,j}=f(\mathbf{s}_{i},b_{j}),j=1,\dots,K\\}.$ It is easy to see
that $\\{\mathbf{M}_{(:,i)}\\}_{i=1}^{K}$ are the column basis of
$S[f]\bigcup\mathbb{R}^{n_{1}}$.

Given any matrix $\mathbf{U}\in\mathbb{R}^{n_{1}\times n_{2}}(n_{2}\geq K)$ as
$\mathbf{U}_{(i,j)}=f(\mathbf{s}_{i},b_{j})$ with
$\mathbf{s}_{i}\in\mathcal{S},b_{j}\in\mathcal{Z}_{f}$, we have $\mathbf{U}\in
S[f]$ and $\mathrm{rank}(\mathbf{U})\leq K$. Furthermore, the column vector in
$\mathbf{U}$ is a linear combination of the column basis
$\\{\mathbf{M}_{(:,i)}\\}_{i=1}^{K}$:

$\mathbf{U}_{(:,j)}=\sum_{k=1}^{K}c_{k}^{(b_{j})}\mathbf{M}_{(:,k)},\text{ for
}j=1,2,\dots,n_{2}.$ (8)

Here, we utilize Eq.(8) and rewrite $f(\mathbf{s},b)$ by
$\mathbf{c}^{(b)}=[c^{(b)}_{1},c^{(b)}_{2},\dots,c^{(b)}_{K}]$:

$f(\mathbf{s},b)=\sum_{k=1}^{K}c_{k}^{(b)}\mathbf{M}_{(i,k)}=\sum_{k=1}^{K}c_{k}^{(b)}f(\mathbf{s},b_{k}),$
(9)

for any $\mathbf{s}\in\mathcal{S},b\in\mathcal{Z}_{f}.$

Next, we will generalize this conclusion from
$\mathbf{\mathbf{missing}}\mathbf{s}\in\mathcal{S}$ to any
$\mathbf{s}\in\mathcal{A}_{f}$. Given
$\mathbf{\tilde{s}}\in\mathcal{A}_{f}/\mathcal{S}$ and we construct a matrix:
$\mathbf{T}\in\mathbb{R}^{(n_{1}+1)\times n_{2}}$, where
$\mathbf{T}_{(i,j)}=f(\mathbf{s}_{i},b_{j})$ and
$\mathbf{s}_{i}\in\mathcal{S}$ for $i=1,2,\dots,n_{1}$ and
$\mathbf{s}_{n_{1}+1}=\mathbf{\tilde{s}}$. Assume there exist $K$ column
vectors $\\{\mathbf{T}_{(:,j_{k})}\\}_{k=1}^{K}$ such that
$\mathbf{T}_{(1:n_{1},j_{k})}=\mathbf{M}_{(:,k)},$ for $k=1,2,\dots,K$. Hence,
we get that $\mathrm{rank}(\mathbf{T})=K$ and for $j=1,2,\dots,n_{2},$

$\mathbf{T}_{(:,j)}=\sum_{k=1}^{K}d_{k}^{(b_{j})}\mathbf{T}_{(:,j_{k})},$

$\mathbf{T}_{(1:n_{1},j)}=\sum_{k=1}^{K}c_{k}^{(b_{j})}\mathbf{M}_{(:,k)}.$

Owing to the uniqueness of the coefficient vector, we get that
$\mathbf{d}^{(b_{j})}=\mathbf{c}^{(b_{j})}$. Hence, we have

$\mathbf{T}_{(n_{1}+1,j)}=\sum_{k=1}^{K}c_{k}^{(b_{j})}\mathbf{T}_{(n_{1}+1,k)},$
(10)

which leads
$f(\mathbf{\tilde{s}},b)=\sum_{k=1}^{K}c_{k}^{(b)}f(\mathbf{\tilde{s}},b_{k})$
for any $\mathbf{\tilde{s}}\in\mathcal{A}_{f}/\mathcal{S}$. This gives the
linear representation form of the factor function $f(\mathbf{s},b)$ (with
fixed $\mathbf{s}$ and $b$) using some basis functions $f(\mathbf{s},b_{k})$
with $b_{k}\in\mathcal{T}$.

We define the factor function
$f_{\mathrm{spatial}}(\cdot):\mathcal{A}_{f}\rightarrow\mathbb{R}^{K}$ as

$f_{\mathrm{spatial}}(\mathbf{\tilde{s}}):=[f(\mathbf{\tilde{s}},b_{1}),f(\mathbf{\tilde{s}},b_{2}),\dots,f(\mathbf{\tilde{s}},b_{K})]^{T}.$

Also, define the matrix function $h(\cdot):\
\mathcal{N}^{(K)}\times\mathcal{Z}_{f}\rightarrow\mathbb{R}$ as

$h(i,b):=c_{i}^{(b)},$

where $\mathcal{N}^{(K)}=\\{1,2,\dots,K\\}.$ From the above analysis, we see
that for any
$(\mathbf{s},b)\in\mathcal{D}_{f}=\mathcal{A}_{f}\times\mathcal{Z}_{f}$, it
holds that

$f(\mathbf{s},b)=\sum_{k=1}^{K}h(k,b)(f_{\mathrm{spatial}}(\mathbf{s}))_{(k)}.$
(11)

Denote $f_{\mathrm{spectral}}(\cdot):\mathcal{Z}_{f}\rightarrow\mathbb{R}^{K}$
as

$f_{\mathrm{spectral}}(b):=[h(1,b),h(2,b),\cdots,h(K,b)]^{T}\in\mathbb{R}^{K\times
1},$

then Eq.(11) is rewritten as

$f(\mathbf{s},b)=f_{\mathrm{spatial}}(\mathbf{s})\cdot
f^{T}_{\mathrm{spectral}}(b).$ (12)

##  Appendix C Proof of Theorem 4

For any $(\mathbf{s}_{1},b_{1}),(\mathbf{s}_{2},b_{2})\in\mathcal{D}_{f}$, we
have

$\displaystyle|f(\mathbf{s}_{1},b_{1})-f(\mathbf{s}_{2},b_{1})|$ (13)
$\displaystyle=$
$\displaystyle|\Phi_{\alpha}(\mathbf{s}_{1})\cdot\Psi_{\theta}^{T}(b_{1})-\Phi_{\alpha}(\mathbf{s}_{2})\cdot\Psi_{\theta}^{T}(b_{1})|$
$\displaystyle\leq$
$\displaystyle|(\Phi_{\alpha}(\mathbf{s}_{1})-\Phi_{\alpha}(\mathbf{s}_{2}))\cdot\Psi_{\theta}^{T}(b_{1})|$
$\displaystyle\leq$
$\displaystyle\|\Phi_{\alpha}(\mathbf{s}_{1})-\Phi_{\alpha}(\mathbf{s}_{2})\|_{1}\|\Psi_{\theta}^{T}(b_{1})\|_{1}.$

Note that $\sigma(\cdot)$ is Lipschitz continuous, i.e.,
$|\sigma(x)-\sigma(y)|\leq\kappa|x-y|$ holds for any $x,y,$ and letting $y=0$
derives $\sigma(x)\leq\kappa|x|$ since $\sigma(0)=0$. On the other hand,
denote $\psi^{(1)}(b)=\mathbf{W}_{1}^{1}t$ and
$\psi^{(k)}(b)=\mathbf{W}_{k}^{1}\sigma(\psi^{(k-1)}(b))$, we get

$\displaystyle\|\Psi_{\theta}(b)\|_{1}=\|\psi^{(d)}(b)\|_{1}$ (14)
$\displaystyle\leq$
$\displaystyle\|\mathbf{W}_{d}^{1}\|_{1}\|\sigma(\psi^{(d-1)}(b))\|_{1}$
$\displaystyle\leq$
$\displaystyle\eta\kappa\|\psi^{(d-1)}(b)\|_{1}\leq\eta^{d}\kappa^{d-1}|b|.$

Meanwhile, we denote $\phi^{(1)}(\mathbf{s})=\mathbf{W}_{1}^{2}\mathbf{s}$ and
$\phi^{(k)}(\mathbf{s})=\mathbf{W}_{k}^{2}\sigma(\phi^{(k-1)}(\mathbf{s}))$.
Then it holds that

$\displaystyle\|\Phi_{\alpha}(\mathbf{s}_{1})-\Phi_{\alpha}(\mathbf{s}_{2})\|_{1}$
(15) $\displaystyle=$
$\displaystyle\|\phi^{(d)}(\mathbf{s}_{1})-\phi^{(k)}(\mathbf{s}_{2})\|_{1}$
$\displaystyle=$
$\displaystyle\|\mathbf{W}_{d}^{2}(\sigma(\phi^{(d-1)}(\mathbf{s}_{1}))-\sigma(\phi^{(d)}(\mathbf{s}_{2}))\|_{1}$
$\displaystyle\leq$
$\displaystyle\eta\kappa\|\phi^{(d-1)}(\mathbf{s}_{1})-\phi^{(d-1)}(\mathbf{s}_{2})\|_{1}$
$\displaystyle\leq$
$\displaystyle\eta^{d}\kappa^{d-1}\|\mathbf{s}_{1}-\mathbf{s}_{2}\|_{1},$

Combining the inequalities Eq.(14) and Eq.(15), we have

$\displaystyle|f(\mathbf{s}_{1},b_{1})-f(\mathbf{s}_{2},b_{1})|$ (16)
$\displaystyle\leq$
$\displaystyle\eta^{2d+1}\kappa^{2d-2}|b_{1}|\|\mathbf{s}_{1}-\mathbf{s}_{2}\|_{1}$
$\displaystyle\leq$
$\displaystyle\eta^{2d+1}\kappa^{2d-2}\zeta\|\mathbf{s}_{1}-\mathbf{s}_{2}\|_{1}.$

Similarly, we prove that

$|f(\mathbf{s}_{2},b_{1})-f(\mathbf{s}_{2},b_{2})|\leq\eta^{2d+1}\kappa^{2d-2}\zeta|b_{1}-b_{2}|.$

Combining the above two inequalities, the triangle inequality completes the
proof.

## References

  * Aiazzi et al. [2006] Aiazzi, B., Alparone, L., Baronti, S., Garzelli, A., Selva, M.  (2006).  MTF-tailored multiscale fusion of high-resolution MS and Pan imagery.  _Photogrammetric Engineering & Remote Sensing_, _72_(5), 591–596.  doi: 10.14358/PERS.72.5.591
  * Aiazzi et al. [2007] Aiazzi, B., Baronti, S., Selva, M.  (2007).  Improving component substitution pansharpening through multivariate regression of MS $+$ pan data.  _IEEE Transactions on Geoscience and Remote Sensing_ , _45_(10), 3230–3239.  doi: 10.1109/TGRS.2007.901007
  * Anokhin et al. [2021] Anokhin, I., Demochkin, K., Khakhulin, T., Sterkin, G., Lempitsky, V., Korzhenkov, D.  (2021).  Image generators with conditionally-independent pixel synthesis.  In _IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ (pp. 14278–14287).  doi: 10.1109/CVPR46437.2021.01405
  * Baumgardner et al. [2015] Baumgardner, M. F., Biehl, L. L., Landgrebe, D. A.  (2015).  220 band aviris hyperspectral image data set: June 12, 1992 indian pine test site 3.  _Purdue University Research Repository_ , _10_(7), 991.
  * Chen et al. [2021] Chen, Y., Liu, S., Wang, X.  (2021).  Learning continuous image representation with local implicit image function.  In _IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ (pp. 8628–8638).  doi: 10.1109/CVPR46437.2021.00852
  * Debals et al. [2017] Debals, O., Van Barel, M., De Lathauwer, L.  (2017).  Nonnegative matrix factorization using nonnegative polynomial approximations.  _IEEE Signal Processing Letters_ , _24_(7), 948–952.  doi: 10.1109/LSP.2017.2697680
  * Dell’Acqua et al. [2004] Dell’Acqua, F., Gamba, P., Ferrari, A., Palmason, J. A., Benediktsson, J. A., Árnason, K.  (2004).  Exploiting spectral and spatial information in hyperspectral urban data with high resolution.  _IEEE Geoscience and Remote Sensing Letters_ , _1_(4), 322–326.  doi: 10.1109/LGRS.2004.837009
  * Deng et al. [2023] Deng, S., Wu, R., Deng, L.-J., Ran, R., Jiang, T.-X.  (2023).  Implicit neural feature fusion function for multispectral and hyperspectral image fusion.  _arXiv preprint arXiv:2307.07288_.
  * Dian et al. [2017] Dian, R., Fang, L., Li, S.  (2017).  Hyperspectral image super-resolution via non-local sparse tensor factorization.  In _IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ (pp. 5344–5353).  doi: 10.1109/CVPR.2017.411
  * Dian et al. [2021] Dian, R., Li, S., Kang, X.  (2021).  Regularizing hyperspectral and multispectral image fusion by CNN denoiser.  _IEEE Transactions on Neural Networks and Learning Systems_ , _32_(3), 1124–1135.  doi: 10.1109/TNNLS.2020.2980398
  * Dong et al. [2016] Dong, W., Fu, F., Shi, G., Cao, X., Wu, J., Li, G., Li, X.  (2016).  Hyperspectral image super-resolution via non-negative structured sparse representation.  _IEEE Transactions on Image Processing_ , _25_(5), 2337–2352.  doi: 10.1109/TIP.2016.2542360
  * Dupont et al. [2021] Dupont, E., Teh, Y. W., Doucet, A.  (2021).  Generative models as distributions of functions.  In _International Conference on Artificial Intelligence and Statistics (AISTATS)_ (pp. 2989–3015).
  * Gao et al. [2014] Gao, L., Li, J., Khodadadzadeh, M., Plaza, A., Zhang, B., He, Z., Yan, H.  (2014).  Subspace-based support vector machines for hyperspectral image classification.  _IEEE Geoscience and Remote Sensing Letters_ , _12_(2), 349–353.  doi: 10.1109/LGRS.2014.2341044
  * Guarino et al. [2024] Guarino, G., Ciotola, M., Vivone, G., Scarpa, G.  (2024).  Band-wise hyperspectral image pansharpening using CNN model propagation.  _IEEE Transactions on Geoscience and Remote Sensing_ , _62_ , 1-18.  doi: 10.1109/TGRS.2023.3339337
  * Guo et al. [2014] Guo, Q., Zhang, B., Ran, Q., Gao, L., Li, J., Plaza, A.  (2014).  Weighted-rxd and linear filter-based rxd: Improving background statistics estimation for anomaly detection in hyperspectral imagery.  _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , _7_(6), 2351–2366.  doi: 10.1109/JSTARS.2014.2302446
  * He et al. [2024] He, L., Fang, Z., Li, J., Chanussot, J., Plaza, A.  (2024).  Two spectral–spatial implicit neural representations for arbitrary-resolution hyperspectral pansharpening.  _IEEE Transactions on Geoscience and Remote Sensing_ , _62_ , 1-21.  doi: 10.1109/TGRS.2024.3380067
  * Jia et al. [2023] Jia, S., Min, Z., Fu, X.  (2023).  Multiscale spatial-spectral transformer network for hyperspectral and multispectral image fusion.  _Information Fusion_ , _96_ , 117–129.  doi: 10.1016/j.inffus.2023.03.011
  * Khader et al. [2023] Khader, A., Yang, J., Xiao, L.  (2023).  Model-guided deep unfolded fusion network with nonlocal spatial-spectral priors for hyperspectral image super-resolution.  _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , _16_ , 4607-4625.  doi: 10.1109/JSTARS.2023.3272370
  * Laben Brower [2000] Laben, C. A., Brower, B. V.  (2000, January 4).  _Process for enhancing the spatial resolution of multispectral imagery using pansharpening._ Google Patents.  (US Patent 6,011,875)
  * Le Saux et al. [2018] Le Saux, B., Yokoya, N., Hansch, R., Prasad, S.  (2018).  2018 ieee GRSS data fusion contest: Multimodal land use classification [technical committees].  _IEEE Geoscience and Remote Sensing Magazine_ , _6_(1), 52–54.  doi: 10.1109/MGRS.2018.2798161
  * Loncan et al. [2015] Loncan, L., De Almeida, L. B., Bioucas-Dias, J. M., Briottet, X., Chanussot, J., Dobigeon, N., … others  (2015).  Hyperspectral pansharpening: A review.  _IEEE Geoscience and Remote Sensing Magazine_ , _3_(3), 27–46.  doi: 10.1109/MGRS.2015.2440094
  * Luo et al. [2024] Luo, Y., Zhao, X., Li, Z., Ng, M. K., Meng, D.  (2024).  Low-rank tensor function representation for multi-dimensional data recovery.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , _46_(5), 3351-3369.  doi: 10.1109/TPAMI.2023.3341688
  * Mildenhall et al. [2021] Mildenhall, B., Srinivasan, P. P., Tancik, M., Barron, J. T., Ramamoorthi, R., Ng, R.  (2021).  NeRF: Representing scenes as neural radiance fields for view synthesis.  _Communications of the ACM_ , _65_(1), 99–106.  doi: 10.1145/3503250
  * Nguyen et al. [2022] Nguyen, H. V., Ulfarsson, M. O., Sveinsson, J. R., Dalla Mura, M.  (2022).  Deep SURE for unsupervised remote sensing image fusion.  _IEEE Transactions on Geoscience and Remote Sensing_ , _60_ , 1–13.  doi: 10.1109/TGRS.2022.3215902
  * Piziak Odell [1999] Piziak, R., Odell, P. L.  (1999).  Full rank factorization of matrices.  _Mathematics Magazine_ , _72_(3), 193–201.  doi: 10.2307/2690882
  * Ramasinghe Lucey [2022] Ramasinghe, S., Lucey, S.  (2022).  Beyond periodicity: Towards a unifying framework for activations in coordinate-mlps.  In _European Conference on Computer Vision (ECCV)_ (Vol. 13693, pp. 142–158).  doi: 10.1007/978-3-031-19827-4_9
  * Shah Sitawarin [2024] Shah, K., Sitawarin, C.  (2024).  SPDER: Semiperiodic damping-enabled object representation.  In _International Conference on Learning Representations (ICLR)._
  * Simoes et al. [2014] Simoes, M., Bioucas-Dias, J., Almeida, L. B., Chanussot, J.  (2014).  A convex formulation for hyperspectral image super-resolution via subspace-based regularization.  _IEEE Transactions on Geoscience and Remote Sensing_ , _53_(6), 3373–3388.  doi: 10.1109/TGRS.2014.2375320
  * Sitzmann et al. [2020] Sitzmann, V., Martel, J., Bergman, A., Lindell, D., Wetzstein, G.  (2020).  Implicit neural representations with periodic activation functions.  In _Advances in Neural Information Processing Systems (NeurIPS)_ (Vol. 33, pp. 7462–7473).
  * Sun et al. [2021] Sun, L., Dong, W., Li, X., Wu, J., Li, L., Shi, G.  (2021).  Deep maximum a posterior estimator for video denoising.  _International Journal of Computer Vision_ , _129_(10), 2827–2845.  doi: 10.1007/S11263-021-01510-7
  * Takikawa et al. [2021] Takikawa, T., Litalien, J., Yin, K., Kreis, K., Loop, C., Nowrouzezahrai, D., … Fidler, S.  (2021).  Neural geometric level of detail: Real-time rendering with implicit 3D shapes.  In _IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ (pp. 11358–11367).  doi: 10.1109/CVPR46437.2021.01120
  * Vivone et al. [2014] Vivone, G., Alparone, L., Chanussot, J., Dalla Mura, M., Garzelli, A., Licciardi, G. A., … Wald, L.  (2014).  A critical comparison among pansharpening algorithms.  _IEEE Transactions on Geoscience and Remote Sensing_ , _53_(5), 2565–2586.  doi: 10.1109/TGRS.2014.2361734
  * Vivone et al. [2023] Vivone, G., Garzelli, A., Xu, Y., Liao, W., Chanussot, J.  (2023).  Panchromatic and hyperspectral image fusion: Outcome of the 2022 WHISPERS hyperspectral pansharpening challenge.  _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , _16_ , 166–179.  doi: 10.1109/JSTARS.2022.3220974
  * K. Wang et al. [2020] Wang, K., Wang, Y., Zhao, X.-L., Chan, J. C.-W., Xu, Z., Meng, D.  (2020).  Hyperspectral and multispectral image fusion via nonlocal low-rank tensor decomposition and spectral unmixing.  _IEEE Transactions on Geoscience and Remote Sensing_ , _58_(11), 7654–7671.  doi: 10.1109/TGRS.2020.2983063
  * T. Wang et al. [2024] Wang, T., Li, J., Ng, M. K., Wang, C.  (2024).  Nonnegative matrix functional factorization for hyperspectral unmixing with nonuniform spectral sampling.  _IEEE Transactions on Geoscience and Remote Sensing_ , _62_ , 1-13.  doi: 10.1109/TGRS.2023.3347414
  * W. Wang et al. [2024] Wang, W., Deng, L.-J., Ran, R., Vivone, G.  (2024).  A general paradigm with detail-preserving conditional invertible network for image fusion.  _International Journal of Computer Vision_ , _132_(4), 1029–1054.  doi: 10.1007/S11263-023-01924-5
  * X. Wang, Cheng, et al. [2023] Wang, X., Cheng, C., Liu, S., Song, R., Wang, X., Feng, L.  (2023).  SS-INR: Spatial-spectral implicit neural representation network for hyperspectral and multispectral image fusion.  _IEEE Transactions on Geoscience and Remote Sensing_ , _61_ , 1-14.  doi: 10.1109/TGRS.2023.3317413
  * X. Wang, Wang, et al. [2023] Wang, X., Wang, X., Song, R., Zhao, X., Zhao, K.  (2023).  MCT-Net: Multi-hierarchical cross transformer for hyperspectral and multispectral image fusion.  _Knowledge-Based Systems_ , _264_ , 110362.  doi: 10.1016/j.knosys.2023.110362
  * Z. Wang et al. [2023] Wang, Z., Ng, M. K., Michalski, J., Zhuang, L.  (2023).  A self-supervised deep denoiser for hyperspectral and multispectral image fusion.  _IEEE Transactions on Geoscience and Remote Sensing_ , _61_ , 1-14.  doi: 10.1109/TGRS.2023.3303921
  * Wu et al. [2024] Wu, H., Wu, S., Zhang, K., Liu, X., Shi, S., Bian, C.  (2024).  Unsupervised blind spectral–spatial cross-super-resolution network for HSI and MSI fusion.  _IEEE Transactions on Geoscience and Remote Sensing_ , _62_ , 1-14.  doi: 10.1109/TGRS.2024.3362862
  * Xu Jiao [2023] Xu, W., Jiao, J.  (2023).  Revisiting implicit neural representations in low-level vision.  In _International Conference on Learning Representations Workshop (ICLR Workshop)._
  * Yokota et al. [2015] Yokota, T., Zdunek, R., Cichocki, A., Yamashita, Y.  (2015).  Smooth nonnegative matrix and tensor factorizations for robust multi-way data analysis.  _Signal Processing_ , _113_ , 234–249.  doi: 10.1016/j.sigpro.2015.02.003
  * Yokoya et al. [2012] Yokoya, N., Yairi, T., Iwasaki, A.  (2012).  Coupled nonnegative matrix factorization unmixing for hyperspectral and multispectral data fusion.  _IEEE Transactions on Geoscience and Remote Sensing_ , _50_(2), 528–537.  doi: 10.1109/TGRS.2011.2161320
  * K. Zhang et al. [2022] Zhang, K., Zhu, D., Min, X., Zhai, G.  (2022).  Implicit neural representation learning for hyperspectral image super-resolution.  _IEEE Transactions on Geoscience and Remote Sensing_ , _61_ , 1–12.  doi: 10.1109/TGRS.2022.3230204
  * L. Zhang et al. [2018] Zhang, L., Wei, W., Bai, C., Gao, Y., Zhang, Y.  (2018).  Exploiting clustering manifold structure for hyperspectral imagery super-resolution.  _IEEE Transactions on Image Processing_ , _27_(12), 5969–5982.  doi: 10.1109/TIP.2018.2862629
  * Zheng et al. [2020] Zheng, K., Gao, L., Liao, W., Hong, D., Zhang, B., Cui, X., Chanussot, J.  (2020).  Coupled convolutional neural network with adaptive response function learning for unsupervised hyperspectral super resolution.  _IEEE Transactions on Geoscience and Remote Sensing_ , _59_(3), 2487–2502.  doi: 10.1109/TGRS.2020.3006534
  * Zhuang Bioucas-Dias [2018] Zhuang, L., Bioucas-Dias, J. M.  (2018).  Fast hyperspectral image denoising and inpainting based on low-rank and sparse representations.  _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , _11_(3), 730–742.  doi: 10.1109/JSTARS.2018.2796570
  * Zhuang Ng [2021] Zhuang, L., Ng, M. K.  (2021).  FastHyMix: Fast and parameter-free hyperspectral image mixed noise removal.  _IEEE Transactions on Neural Networks and Learning Systems_ , _34_(8), 4702–4716.  doi: 10.1109/TNNLS.2021.3112577