\section{Introduction}

Hyperspectral (HS) cameras are able to acquire images with very high spectral resolution. However, the fundamental compromise between signal-to-noise ratio, spatial resolution, and spectral resolution means that their spatial resolution is usually low~\cite{shaw2003reviewSpecImaging}. Multispectral (MS) devices, on the other hand, are able to achieve a much higher spatial resolution since they contain only a small number of spectral bands.
An approach that attempts to circumvent the physical limitations of imaging sensors consists in combining HS and MS images (MSI) of the same scene to obtain images with high spatial and spectral resolution, in a multimodal image fusion problem~\cite{yokoya2017HS_MS_fusinoRev}, commonly referred to as \emph{hyperspectral super resolution}.

Different algorithms have been proposed to solve this problem. Early approaches were based on component substitution or on multiresolution analysis, which attempt to extract high-frequency spatial details from the MSI and combine them with the HS image (HSI)~\cite{carper1990componentSubstitutionPansharpening,liu2000MRA_HS_MS_fusion}.
Subspace-based formulations have later received a significant amount of interest as they explore the natural representation of the pixels in an HSI as the linear combination of a small number of spectral signatures~\cite{yokoya2017HS_MS_fusinoRev,Keshava:2002p5667,borsoi2018superpixels1_sparseU}.
Different algorithms have been proposed following this approach using, e.g., Bayesian formulations~\cite{hardie2004hs_ms_fusionMAP} or sparse representations on learned dictionaries~\cite{wei2015hs_ms_fusionBayesianSparse}, and different kinds of matrix factorization formulations employing sparse and spatial regularizations~\cite{kawakami2011hs_ms_fusionMFsparse,simoes2015HySure}, or estimating both the basis vectors and their coefficients blindly/unsupervisedly from the images~\cite{yokoya2012coupledNMF}.

The natural representation of HSIs and MSIs as 3-dimensional tensors has been successfully exploited for hyperspectral unmixing, denoising~\cite{imbiriba2018ULTRA_V,veganzones2015tensorDecompMultisngularHSIs,imbiriba2018_ULTRA,liu2012denoisingHSIsTensorCPD} and super-resolution. Superior super-resolution performance and exact recovery guarantees have been obtained using this formulation~\cite{kanatsoulis2018hyperspectralSRR_coupledCPD,prevost2020coupledTucker_hyperspectralSRR_TSP}.

The image fusion problem was formulated in~\cite{kanatsoulis2018hyperspectralSRR_coupledCPD} as a coupled tensor approximation problem. Assuming that the high resolution (HR) image admits a low-rank canonical polyadic decomposition (CPD), the problem was solved using an alternating optimization strategy. The recovery of the correct HR image (HRI) was shown to be guaranteed provided that the CPD of the MSI is identifiable, and state of the art performance was achieved.
A recent work extended this approach by assuming the high resolution images to follow a block term decomposition (BTD), which shows a closer connection to the physical mixing model when compared to the CPD~\cite{zhang2019hyperspectralSRR_coupledBTD}.
A simpler approach was later proposed in~\cite{kanatsoulis2018hyperspSRR_lowRankMatrixTensor} by requiring only the computation of one CPD of the MSI and a singular value decomposition (SVD) of the HSI.

A Tucker decomposition-based approach was later considered in~\cite{prevost2019coupledTucker_hyperspectralSRR,prevost2020coupledTucker_hyperspectralSRR_TSP}. Closed form SVD-based algorithms were proposed for the image fusion problem, achieving results comparable to~\cite{kanatsoulis2018hyperspectralSRR_coupledCPD} at a very small computational complexity; exact recovery guarantees were also provided. A coupled Tucker approximation was also considered in~\cite{li2018hs_ms_fusionTensorFactorization} using an alternating optimization approach and employing a sparsity regularization on the elements of the core tensor. {This work was later extended by incorporating the piecewise smoothness of the reconstructed image tensor by using a Total Variation regularization along each of its modes~\cite{xu2020HSMSfusion_1DTV}.}
Another approach considered the CPD of non-local similar patch tensors to explore the non-local redundancy of the image~\cite{xu2019nonlocalTensorDecompImageFusion}.
{A different non-local approach was also proposed in~\cite{dian2019hyperspectralSRR_tensorMultiRank} by using a recent definition of the nuclear norm of order-3 tensors.}

Most existing algorithms, however, share a common limitation: they assume that the HSI and the MSI are acquired under the same conditions.
However, despite the short revisit cycles provided by the increasing number of optical satellites orbiting the Earth (e.g. Sentinel, Orbview, Landsat and Quickbird missions), the number of platforms carrying both HS and MS sensors is still considerably limited~\cite{eckardt2015DESIS_satelliteSpecs,kaufmann2006EnMAP_satelliteSpecs}. This makes combining HS and MS observations acquired on board of different satellites of great interest to obtain HRIs~\cite{hilker2009fusion_MODIS_denseTimeSeries_exp2,emelyanova2013fusion_MODIS_denseTimeSeries_exp}.
Images acquired at different time instants can be impacted by, e.g., illumination, atmospheric or seasonal changes. This may result in significant variations between the HSI and the MSI~\cite{borsoi2020variabilityReview}, negatively impacting traditional image fusion algorithms.

Recently, a method was proposed to combine HSIs and MSIs accounting for seasonal {(inter-image)} spectral variability~\cite{Borsoi_2018_Fusion}. Using a low-rank matrix formulation, the set of spectral basis vectors of the HRIs underlying the HS and the MS observations are allowed to be different from each other, with variations introduced by a set of multiplicative scaling factors~\cite{imbiriba2018glmm}. This algorithm led to significant performance improvements when the HSI and MSI are subject to spatially uniform seasonal or acquisition variations. However, it does not account for spatially localized changes commonly seen in practical scenes~\cite{borsoi2020variabilityReview}. Moreover, the algorithm in~\cite{Borsoi_2018_Fusion} presented high computation times and does not offer any theoretical guarantees.

In this paper, we propose a tensor-based image fusion formulation that accounts for localized spatial and spectral changes between the HSI and MSI. A general observation model is considered, in which the HRI underlying the MSI admits an additive variability term to account for changes between the scenes. Studying the general identifiability of this model, we show that this variability term can only be identified in general up to its smooth structure (which is defined according to the degradation operators). To introduce additional a priori information and mitigate the ambiguity associated with the proposed model, both the HRI and the additive perturbations are assumed to have low multilinear rank (i.e., to admit a Tucker decomposition). Two algorithms are then proposed, one totally algebraic and another based on an optimization procedure. Theoretical guarantees for the exact recovery of the HRI are provided for both. Simulation results show that the proposed optimization-based algorithm yields superior performance at a considerably lower computational cost when compared to~\cite{Borsoi_2018_Fusion}, especially when spatially localized variability is considered.

\section{Tensors -- background} \label{sec:background}
\subsection{Notation and definitions}
\label{sec:back_defs}

An order-3 tensor $\mathcal{T}\in\amsmathbb{R}^{N_1\times N_2 \times N_3}$ is an $N_1\times N_2 \times N_3$ array whose elements are indexed by~$[\mathcal{T}]_{n_1,n_2,n_3}$. Each dimension of \textcolor{darkgreen}{an order-3 tensor} is called a \emph{mode}. A mode-$k$ \emph{fiber} of tensor $\mathcal{T}$ is the one-dimensional subset of $\mathcal{T}$ which is obtained by fixing \textcolor{darkgreen}{all but one of the three modes -- the $k$-th dimension}. Similarly, a \emph{slab} or \emph{slice} of a tensor $\mathcal{T}$ is a matrix whose elements are the two-dimensional subset of $\mathcal{T}$ obtained by fixing all but two of its modes. Operator $\operatorname{vec}(\cdot)$ represents the standard matrix column-major vectorization, or tensor vectorization. The (left) pseudo-inverse of matrix $\boldsymbol{X}$ is denoted by $\boldsymbol{X}^\dagger$. We denote scalars by lowercase ($x$) or uppercase ($X$) plain font, vectors and matrices by lowercase ($\boldsymbol{x}$) and uppercase ($\boldsymbol{X}$) bold font, respectively, and order-3 tensors by calligraphic plain font ($\mathcal{T}$) or using the blackboard Greek alphabet ($\mathcal{\Theta}$). In the following, we review some useful operations of multilinear (tensor) algebra that will be used in the rest of the manuscript (see, e.g.,~\cite{cichocki2015tensor,kolda2009tensor} for more details).

\begin{definition}
\noinden
{\color{red}The \textbf{mode-$k$ product} between a tensor $\mathcal{T}$ and a matrix $\cb{B}$ produces a tensor $\mathcal{U}$ that is evaluated such that each mode-$k$ fiber of $\mathcal{T}$ is multiplied by $\cb{B}$.
For instance, the mode-2 product between $\mathcal{T} \in \amsmathbb{R}^{N_1\times N_2\times N_3}$ and $\cb{B} \in \amsmathbb{R}^{M_2\times N_2}$ produces a tensor
$\mathcal{U}\in\amsmathbb{R}^{N_{1}\times M_2\times N_{3}}$, denoted by
$\mathcal{U}=\mathcal{T}\times_2\cb{B}$.}
{\color{magenta}
Its elements are accessed as $[\mathcal{U}]_{n_1,m_2,n_3}=\sum_{i=1}^{N_2}[\mathcal{T}]_{n_1,i,n_3} [\cb{B}]_{m_2,i}\ , \ m_2=1,\ldots,M_2$.}
\end{definition}
{\color{magenta}Note that the mode-$k$ product has the following properties:
\begin{align}
    \mathcal{T} \times_i\boldsymbol{A} \times_j\boldsymbol{B} &= \mathcal{T} \times_j\boldsymbol{B}  \times_i \boldsymbol{A}\,, \ i \neq j \,,
    \label{eq:mode_k_prop1}
    \\
    \mathcal{T} \times_k \boldsymbol{A} \times_k\boldsymbol{B} &= \mathcal{T}\times_k\big(\boldsymbol{A}\boldsymbol{B} \big) \,.
    \label{eq:mode_k_prop2}
\end{align}}
\begin{definition}

\textcolor{darkgreen}{The full \textbf{multilinear product} is denoted by $\big\ldbrack\mathcal{T};\boldsymbol{B}_{1},\boldsymbol{B}_{2},\boldsymbol{B}_{3}\big\rdbrack$, and consists of a series of successive mode-$k$ products, for {$k\in\{1,2,3\}$, between a tensor $\mathcal{T}$ and matrices $\boldsymbol{B}_1$, $\boldsymbol{B}_2$ and $\boldsymbol{B}_3$,} respectively, and {is} expressed as $\mathcal{T}\times_1\boldsymbol{B}_{1}\times_2\boldsymbol{B}_{2}\times_3\boldsymbol{B}_{3}$.}
\end{definition}
\begin{definition}
The \textbf{mode-$k$ matricization} of an order-3 tensor $\mathcal{T}\in\amsmathbb{R}^{N_1\times N_2\times N_3}$, denoted by $\unfold{\boldsymbol{T}}{k}$, arranges its mode-$k$ fibers to be the columns of the resulting matrix $\unfold{\boldsymbol{T}}{k}\in\amsmathbb{R}^{N_k\times N_{\ell}N_{m}}$, $k,\ell,m\in\{1,2,3\}$, $k \neq \ell \neq m$, where the $n_k$-th row of $\unfold{\boldsymbol{T}}{k}$ consists of the vectorization of the slice of $\mathcal{T}$ obtained by fixing the index of the $k$-th mode of $\mathcal{T}$ as $n_k$.
\end{definition}

\begin{definition}
We define by $\operatorname{tSVD}_R(\boldsymbol{X})$ the operator which returns a matrix containing {\color{red} the $R$ left singular vectors associated with the largest singular values of the matrix $\boldsymbol{X}$}
\end{definition}

\subsection{Tensor decompositions}

{The Tucker decomposition is able to represent an order-3 tensor $\mathcal{T}\in\amsmathbb{R}^{N_1\times N_2\times N_3}$ compactly, using a set of factor matrices given by $\boldsymbol{B}_{i}\in\amsmathbb{R}^{N_i\times K_i}$, $i\in\{1,2,3\}$ and a small core tensor $\mathcal{G}\in\amsmathbb{R}^{K_1\times K_2\times K_3}$, as~\cite{kolda2009tensor}
\begin{align}
    \mathcal{T} &= \big\ldbrack \mathcal{G};\boldsymbol{B}_{1},\boldsymbol{B}_{2},\boldsymbol{B}_{3} \big\rdbrack \,.
    \label{eq:tensorDef_Tucker}
\end{align}
The tuple $(K_1,K_2,K_3)$ is called the multilinear rank of $\mathcal{T}$. Each value $K_i$ is also equal to the rank of the mode-$i$ unfolding of~$\mathcal{T}$~\cite{cichocki2015tensor}.}

{\color{magenta}
The Tucker decomposition
allows the rank along each mode of the tensor to be different (i.e., $K_i\neq K_j$)~\cite{sidiropoulos2017tensor}. This property can be very useful since it allows one to set a higher rank to specific modes of the decomposition in order to adequately represent the data diversity while still keeping the \mbox{model low rank.}
}

The matricizations and vectorization of a tensor $\mathcal{T}$ following the Tucker decomposition~\eqref{eq:tensorDef_Tucker} are given by~\cite{kolda2009tensor}:
\begin{align}

    \operatorname{vec}(\mathcal{T}) &= (\boldsymbol{B}_{3} \otimes \boldsymbol{B}_{2} \otimes \boldsymbol{B}_{1} ) \operatorname{vec}(\textcolor{magenta}{\mathcal{G}}) \,,
    \label{eq:gen_tensor_vec}
    \\
    \unfold{\boldsymbol{T}}{1} &= \boldsymbol{B}_{1} \unfold{\textcolor{magenta}{\boldsymbol{G}}}{1} (\boldsymbol{B}_{3} \otimes \boldsymbol{B}_{2})^\top \,,
    \label{eq:gen_tensor_mat_1}
    \\
    \unfold{\boldsymbol{T}}{2} &= \boldsymbol{B}_{2} \unfold{\textcolor{magenta}{\boldsymbol{G}}}{2} (\boldsymbol{B}_{3} \otimes \boldsymbol{B}_{1})^\top \,,
    \label{eq:gen_tensor_mat_2}
    \\
    \unfold{\boldsymbol{T}}{3} &= \boldsymbol{B}_{3} \unfold{\textcolor{magenta}{\boldsymbol{G}}}{3} (\boldsymbol{B}_{2} \otimes \boldsymbol{B}_{1})^\top \,.
    \label{eq:gen_tensor_mat_3}
\end{align}
{The Tucker decomposition can be computed using fast algorithms such as the high-order SVD~\cite{lathauwer2000multilinearHOSVD}.}

{Another classic tensor decomposition is the Canonical Polyadic Decomposition (CPD)~\cite{kolda2009tensor}. The CPD enjoys uniqueness properties under very mild conditions, and its rank can exceed the dimensions of the tensor. However, it is also more difficult to compute, and the same rank value is used to represent all modes of the tensor.}

{\color{magenta}
The Block Term Decomposition (BTD) generalizes CPD and Tucker, and allows us to combine benefits from both approaches~\cite{lathauwer2008tensor_BTD2_uniqueness}.}
Specifically, the BTD of an order-3 tensor $\mathcal{T}$ is defined as a sum of rank-$(K_{1,r},K_{2,r},K_{3,r})$ terms as~\cite{lathauwer2008tensor_BTD2_uniqueness}:
\begin{align} \label{eq:tensor_BTD_i}
    \mathcal{T} & = \sum_{r=1}^R \big\ldbrack\mathcal{G}_r;\boldsymbol{B}_{r,1},\boldsymbol{B}_{r,2},\boldsymbol{B}_{r,3} \big\rdbrack \,,
\end{align}
where $\mathcal{G}_r$, $\boldsymbol{B}_{r,1}$, $\boldsymbol{B}_{r,2}$, and $\boldsymbol{B}_{r,3}$, for $r=1,\ldots,R$, are the core tensors and the factors corresponding to each mode of $\mathcal{T}$. {Differently from the Tucker decomposition, the BTD additionally requires the selection of parameter $R$ (number of blocks), and can be more costly to compute.}
{\color{magenta}
However, the BTD benefits from uniqueness results which, although not as strong than those of the CPD, are still interesting for many applications~\cite[Section~5]{lathauwer2008tensor_BTD2_uniqueness}.}
{These results will prove very important to derive the recoverability guarantees for the algorithm in Section~\ref{sec:alg2_optimizationBased}.
Note that in the CPD can be viewed as a special case of the BTD with $R=1$.}

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{figures/diagram_tmp2b2.pdf}

    \vspace{-0.7cm}
    \caption{{Proposed imaging model. The observed HSI and MSI are acquired under different conditions or at different time instants. Their underlying HRIs, represented by tensors $\mathcal{Z}_h$ and $\mathcal{Z}_m$, respectively, can be different from each other due to the effect of spectral and acquisition variations, as well as scenery changes. The changes occurring between the HRIs $\mathcal{Z}_h$ and $\mathcal{Z}_m$ are represented using an additive tensor $\mathcal{\Psi}$ which captures variability and changes, and both the HRI and the variability tensors are assumed to have low Tucker rank.}}
    \label{fig:gen_diagram}
\end{figure}

\section{Proposed model and its undeterminacies}
\label{sec:imaging_mdl_identif}

\subsection{The imaging model}

Let an HSI with high spectral resolution and low spatial resolution be represented as an order-3 tensor $\mathcal{Y}_h\in\amsmathbb{R}^{N_1\times N_2\times L_h}$, where $N_1$ and $N_2$ are the spatial and $L_h$ the spectral dimensions.
Similarly, an MSI with high spatial and low spectral resolution is denoted by an order-3 tensor $\mathcal{Y}_m\in\amsmathbb{R}^{M_1\times M_2\times L_m}$, where $M_1>N_1$ and $M_2>N_2$ are the spatial and $L_m<L_h$ the spectral dimensions.
Both the HSI and the MSI are assumed to be degraded versions of a tensor $\mathcal{Z}\in\amsmathbb{R}_+^{M_1\times M_2\times L_h}$, with high spectral and spatial resolutions.
This degradation process is commonly described as~\cite{kanatsoulis2018hyperspectralSRR_coupledCPD,prevost2019coupledTucker_hyperspectralSRR,prevost2020coupledTucker_hyperspectralSRR_TSP,li2018hs_ms_fusionTensorFactorization}:
\begin{align}
    \mathcal{Y}_h &= \mathcal{Z} \times_1 \boldsymbol{P}_1 \times_2 \boldsymbol{P}_2 + \mathcal{E}_h \,,
    \label{eq:tensor_obs_mdl_hsi_novariability}
    \\
    \mathcal{Y}_m &= \mathcal{Z} \times_3 \boldsymbol{P}_3 + \mathcal{E}_m  \,,
    \label{eq:tensor_obs_mdl_msi_novariability}
\end{align}
where tensors $\mathcal{E}_m\in\amsmathbb{R}^{M_1\times M_2\times L_m}$ and $\mathcal{E}_h\in\amsmathbb{R}^{N_1\times N_2\times L_h}$ represent additive noise.
The matrix $\boldsymbol{P}_3\in\amsmathbb{R}^{L_m\times L_h}$ contains the spectral response functions (SRF) of each band of the multispectral sensor, and matrices $\boldsymbol{P}_1\in\amsmathbb{R}^{N_1\times M_1}$ and $\boldsymbol{P}_2\in\amsmathbb{R}^{N_2\times M_2}$ represent the spatial blurring and downsampling in the hyperspectral sensor, which we assume to be separable for each spatial dimension as previously done in, e.g.,~\cite{kanatsoulis2018hyperspectralSRR_coupledCPD,prevost2019coupledTucker_hyperspectralSRR,prevost2020coupledTucker_hyperspectralSRR_TSP,li2018hs_ms_fusionTensorFactorization}.
{Note that since the mode-$k$ product obeys~\eqref{eq:mode_k_prop1}, the choice of the ordering of $\boldsymbol{P}_1$ and $\boldsymbol{P}_2$ in~\eqref{eq:tensor_obs_mdl_hsi_novariability} does not affect the result.}
To make notation more convenient, we also denote the (linear) spatial and spectral degradation operators more compactly as
\begin{align}
    \mathscr{P}_{1,2}(\mathcal{T}) &= \mathcal{T} \times_1 \boldsymbol{P}_1 \times_2 \boldsymbol{P}_2 \,,
    \label{eq:operatorNotation_a}
    \\
    \mathscr{P}_{3}(\mathcal{T}) &= \mathcal{T} \times_3 \boldsymbol{P}_3 \,.
    \label{eq:operatorNotation_b}
\end{align}

Most previous works consider that $\mathcal{Y}_h$ and $\mathcal{Y}_m$ are acquired under the same conditions, implicitly assuming that no variability {or changes} occur between the images.
{However, when the HSI and MSI are not acquired from the same mission/instrument and at the same time, the scene which underlies the (degraded) observations $\mathcal{Y}_h$ and $\mathcal{Y}_m$ can be subject to significant changes, referred to as {\color{magenta} inter-image} variability\footnote{\color{magenta} Here, variability should not be confused with the spectral variability considered in \cite{borsoi2020variabilityReview}, which focuses on inter-pixel endmembers variations.}, which include spatial and spectral variations as illustrated in Fig.~\ref{fig:gen_diagram}. Spectral variations originate from different \textcolor{magenta}{chemical}, atmospheric, illumination or seasonal conditions between the scenes~\cite{somers2011variabilityReview,Zare-2014-ID324-variabilityReview,borsoi2020variabilityReview}, and are typical even under short acquisition time differences. Spatial variations, on the other hand, occur due to, e.g., some regions of the scene being affected unequally by seasonal effects (which are strongly material-dependent~\cite{borsoi2020variabilityReview}) or due to the sudden insertion/removal of an object~\cite{liu2019reviewCD_GRSM}. Spatial variations  can be very prominent when large acquisition time differences are considered.
These effects are not accounted for in the majority of the existing algorithms, what motivates the development of more flexible models.
}

Recently, spatially uniform spectral variability has been considered in~\cite{Borsoi_2018_Fusion}. The image fusion problem was formulated as a matrix factorization problem, and the (multiplicative) spectral variability as well as the spatial coefficients were estimated from the observed images.
However, this work still did not address two fundamental problems: 1) How to account for both spatial and spectral variability and 2) what theoretical guarantees can be offered for the recovery of the HRI and (possibly) of the variability factors under these more challenging conditions.

To address these issues, we adopt a more general approach by considering two different HRIs $\mathcal{Z}_h\in\amsmathbb{R}_+^{M_1\times M_2\times L_h}$ and $\mathcal{Z}_m\in\amsmathbb{R}_+^{M_1\times M_2\times L_h}$, both with high spectral and spatial resolutions, underlying the observed HSI and the MSI, respectively. This leads to the following extension of model~\eqref{eq:tensor_obs_mdl_hsi_novariability}--\eqref{eq:tensor_obs_mdl_msi_novariability}:
\begin{align}
    \mathcal{Y}_h &= \mathscr{P}_{1,2}(\mathcal{Z}_h) + \mathcal{E}_h \,,
    \label{eq:tensor_obs_mdl_hsi_0}
    \\
    \mathcal{Y}_m &= \mathscr{P}_{3}(\mathcal{Z}_m) + \mathcal{E}_m  \,.
    \label{eq:tensor_obs_mdl_msi_0}
\end{align}

{\color{magenta} The tensors $\mathcal{Z}_h$ and $\mathcal{Z}_m$ \textcolor{darkgreen}{represent} the underlying HRIs of the observed scene under different acquisition conditions and at possibly different times. To account for inter-image variability, the HRIs are related to each other as follows:}
\begin{align}
    \mathcal{Z}_m = {\mathcal{Z}_h} + \mathcal{\Psi},
    \label{eq:hr_hsi_msi_add_mdl_0}
\end{align}
where $\mathcal{\Psi}\in\amsmathbb{R}^{M_1\times M_2\times L_h}$ is an additive variability tensor representing changes between the scenes.
{\color{magenta}By introducing $\mathcal{\Psi}$, model~\eqref{eq:hr_hsi_msi_add_mdl_0} makes the inter-image variability between the HSI and the MSI explicit.}

Considering the variability model~\eqref{eq:hr_hsi_msi_add_mdl_0} along with~\eqref{eq:tensor_obs_mdl_hsi_0}--\eqref{eq:tensor_obs_mdl_msi_0}, we obtain the following observation model for the acquired HSI and MSI:
\begin{align}
    \mathcal{Y}_h &= \mathscr{P}_{1,2}({\mathcal{Z}_h}) + \mathcal{E}_h \,,
    \label{eq:tensor_obs_mdl_hsi}
    \\
    \mathcal{Y}_m &= \mathscr{P}_{3}({\mathcal{Z}_h} + \mathcal{\Psi}) + \mathcal{E}_m.
    \label{eq:tensor_obs_mdl_msi_add}
\end{align}
{\color{magenta}In the following, $\mathcal{Z}_h$ and $\mathcal{\Psi}$ will be referred to as the HRI and the variability tensor, respectively.}

\subsection{The image fusion problem and its undeterminacies}

The image fusion problem in this case consists in recovering~$\mathcal{\Psi}$ and ${\mathcal{Z}_h}$ from the observed images $\mathcal{Y}_h$ and $\mathcal{Y}_m$. More precisely,
\begin{align} \label{eq:HS_MS_fus_var_prob}
\Bigg\{
\begin{split}
    & \text{find } {\mathcal{Z}_h}\in\Omega_Z \text{ and }  \mathcal{\Psi}\in\Omega_{\Psi},
    \\
    & \text{such that equations~\eqref{eq:tensor_obs_mdl_hsi}--\eqref{eq:tensor_obs_mdl_msi_add} are satisfied.}
\end{split}
\end{align}
The sets $\Omega_Z\subseteq\amsmathbb{R}^{M_1\times M_2\times L_h}$  and $\Omega_{\Psi}\subseteq\amsmathbb{R}^{M_1\times M_2\times L_h}$ denote prior information about the HRI and the variability factor, respectively.

Since the number of unknowns is significantly greater than the number of observations, problem~\eqref{eq:HS_MS_fus_var_prob} is severely ill-posed and additional \emph{a priori} information about the structure of ${\mathcal{Z}_h}$ and $\mathcal{\Psi}$ must be introduced through the sets $\Omega_Z$ and $\Omega_{\Psi}$ in order to obtain a stable recovery.
Common information that has been used to construct $\Omega_Z$ includes spatial \mbox{(piecewise-)} smoothness~\cite{simoes2015HySure}, low matrix (spectral) rank~\cite{yokoya2012coupledNMF,Borsoi_2018_Fusion}, low tensor rank~\cite{kanatsoulis2018hyperspectralSRR_coupledCPD,prevost2020coupledTucker_hyperspectralSRR_TSP,xu2019nonlocalTensorDecompImageFusion}, and non-local spatial information~\cite{li2018hs_ms_fusionTensorFactorization}.

The choice of prior information in $\Omega_Z$ and $\Omega_{\Psi}$ turns to the question of whether assuming additional structure over the pair $({\mathcal{Z}_h},\mathcal{\Psi})$ makes these variables \emph{identifiable} from the observations $(\mathcal{Y}_{h},\mathcal{Y}_{m})$.
Recent works in HS-MS image fusion advocates for a low-rank tensor model~\cite{kanatsoulis2018hyperspectralSRR_coupledCPD,prevost2020coupledTucker_hyperspectralSRR_TSP,prevost2019coupledTucker_hyperspectralSRR}. However, the case at hand is more challenging because of the additional variability $\mathcal{\Psi}$, which makes the model more ambiguous.

In many inverse problems such as matrix or tensor factorization, dictionary learning and blind deconvolution, identifiability of the underlying variables often can only be defined up to some fundamental ambiguities.
Transformation groups and equivalence classes~\cite{li2016identifiabilityBilinearInverseProbs} can be used to precisely define which sets of solutions can generate each possible observations.
These ideas can be leveraged to characterize some of the fundamental ambiguities associated with the model~\eqref{eq:tensor_obs_mdl_hsi}--\eqref{eq:tensor_obs_mdl_msi_add}, and to provide insights into the development of efficient algorithms.
First, we will show that the presence of $\mathcal{\Psi}$ makes the model fundamentally ambiguous, as the content in ${\mathcal{Z}_h}$ cannot be easily distinguished from that of $\mathcal{\Psi}$. Moreover, we will define an equivalence class that characterizes the sets of images ${\mathcal{Z}_h}$ and factors $\mathcal{\Psi}$ which are certain to result in different observed HSI and MSI.
This gives us insight into what kind of structure from these variables can be recovered from $(\mathcal{Y}_{h},\mathcal{Y}_{m})$.

{To proceed, let us first denote by $\mathscr{A}:(\mathcal{Z}_h,\mathcal{\Psi})\mapsto(\mathcal{Y}_h,\mathcal{Y}_m)$ the operator which describes the degradation process in equations~\eqref{eq:tensor_obs_mdl_hsi}--\eqref{eq:tensor_obs_mdl_msi_add}.}
{By representing operators $\mathscr{P}_{1,2}$, $\mathscr{P}_{3}$ and $\mathscr{A}$ in matrix form as $\widetilde{\boldsymbol{P}}_{1,2}\in\amsmathbb{R}^{N_1N_2L_h \times M_1M_2L_h}$, $\widetilde{\boldsymbol{P}}_{3}\in\amsmathbb{R}^{M_1M_2L_m \times M_1M_2L_h}$ and $\boldsymbol{A}\in\amsmathbb{R}^{(N_1N_2L_h+M_1M_2L_m)\times 2M_1M_2L_h}$, respectively, we can write the model~\eqref{eq:tensor_obs_mdl_hsi}--\eqref{eq:tensor_obs_mdl_msi_add} in the noiseless case ($\mathcal{E}_h = \mathcal{E}_m = \mathcal{0}$) equivalently as: }
\begin{align}
    \begin{bmatrix}
    \operatorname{vec}(\mathcal{Y}_h) \\ \operatorname{vec}(\mathcal{Y}_m)
    \end{bmatrix}
    = \underbrace{
    \begin{bmatrix}
    \widetilde{\boldsymbol{P}}_{1,2} & \cb{0} \\
    \widetilde{\boldsymbol{P}}_{3} & \widetilde{\boldsymbol{P}}_{3}
    \end{bmatrix}
    }_{\boldsymbol{A}}
    \begin{bmatrix}
    \operatorname{vec}(\mathcal{Z}_h) \\ \operatorname{vec}(\mathcal{\Psi})
    \end{bmatrix} \,.
    \label{eq:thm1_mtxForm_1}
\end{align}

{Define also the equivalence relation $\sim_{{Z,\Psi}}$ based on operator $\mathscr{A}$ as follows:
\begin{align}
    (\mathcal{Z}_h,\mathcal{\Psi}) \sim_{{Z,\Psi}} (\mathcal{Z}_h',\mathcal{\Psi}')  \iff  \mathscr{A}(\mathcal{Z}_h,\mathcal{\Psi}) - \mathscr{A}(\mathcal{Z}_h',\mathcal{\Psi}') =\mathcal{0}
    \label{eq:thm_prel2_eq_rel}
\end{align}
and its associated equivalence class (EC) $[(\mathcal{Z}_h,\mathcal{\Psi})]_{\sim_{{Z,\Psi}}}$ as
\begin{align}
    [(\mathcal{Z}_h,\mathcal{\Psi})]_{\sim_{{Z,\Psi}}} & = \big\{\mathcal{X}\in\Omega_Z\times\Omega_{\Psi}\,:\,\mathcal{X} \sim_{{Z,\Psi}} (\mathcal{Z}_h,\mathcal{\Psi})\big\}.

    \label{eq:thm_prel2_eq_class}
\end{align}
Now, we are ready to present the following result.}

\begin{theorem} \label{thm:identifiability_ECc_psi}
Suppose that the observation noise is zero (i.e., $\mathcal{E}_h=\mathcal{E}_m=\mathcal{0}$) and that $\Omega_Z=\amsmathbb{R}^{M_1\times M_2\times L_h}_+$ and $\Omega_{\Psi}=\amsmathbb{R}^{M_1\times M_2\times L_h}$.
Then, given a set of HSI and MSI observations $(\mathcal{Y}_h,\mathcal{Y}_m)${, the following is verified:}
\begin{itemize}
    \item[a)] {If operator $\mathscr{A}$ has nontrivial nullspace (e.g., if $L_m<L_h$ or if $N_1N_2<M_1M_2$), then the pair $(\mathcal{Z}_h,\mathcal{\Psi})$ cannot be uniquely identified from the observations $(\mathcal{Y}_h,\mathcal{Y}_m)$.}

    \item[b)] There is only one (unique) equivalence class $[(\mathcal{Z}_0,\mathcal{\Psi}_0)]_{\sim_{{Z,\Psi}}}$ containing HR images and scaling factors $({\mathcal{Z}_h},\mathcal{\Psi})$ that can generate~$(\mathcal{Y}_h,\mathcal{Y}_m)$ according to model~\eqref{eq:tensor_obs_mdl_hsi}--\eqref{eq:tensor_obs_mdl_msi_add}. In other words, $({\mathcal{Z}_h},\mathcal{\Psi})$ can be identified uniquely up to $[(\mathcal{Z}_0,\mathcal{\Psi}_0)]_{\sim_{{Z,\Psi}}}$.
\end{itemize}

\end{theorem}

\begin{proof}
{\textbf{Proof of a):}} {Due to the special structure of matrix $\boldsymbol{A}$ in~\eqref{eq:thm1_mtxForm_1}, is is clear that $\operatorname{rank}(\boldsymbol{A})=\operatorname{rank}(\widetilde{\boldsymbol{P}}_{1,2})+\operatorname{rank}(\widetilde{\boldsymbol{P}}_{3})\leq 2M_1M_2L_h$. Thus, if $\mathscr{A}$ has nontrivial nullspace then either $\mathscr{P}_{1,2}$ or $\mathscr{P}_{3}$ have nontrivial nullspace.}
If the operator $\mathscr{P}_{1,2}$ has nontrivial nullspace, then we can find ${\mathcal{Z}_h}$, ${\mathcal{Z}_h'}$, different from one another, such that
\begin{align}
    \mathscr{P}_{1,2}({\mathcal{Z}_h}) & = \mathscr{P}_{1,2}({\mathcal{Z}_h'}) \,,
    \label{eq:thm1_i}
\end{align}
implying that $\mathcal{Y}_h=\mathcal{Y}_h'$. Now, we can always find $\mathcal{\Psi}$ and $\mathcal{\Psi}'$ satisfying
\begin{align}
    {\mathcal{Z}_h} + \mathcal{\Psi} & = {\mathcal{Z}_h'} + \mathcal{\Psi}',
    \label{eq:thm1_ii}
\end{align}
which implies that $\mathcal{Y}_m=\mathcal{Y}_m'$ (i.e., the model is not identifiable).

Similarly, if operator $\mathscr{P}_{3}$ has nontrivial nullspace, then suppose we select ${\mathcal{Z}_h}={\mathcal{Z}_h'}$. This makes $\mathcal{Y}_h=\mathcal{Y}_h'$. Then, we can select $\mathcal{\Psi}$ and $\mathcal{\Psi}'$, distinct from one another, satisfying
\begin{align}
    \mathscr{P}_{3}(\mathcal{\Psi}) - \mathscr{P}_{3}(\mathcal{\Psi}') = \mathcal{0} \,,
    \label{eq:thm1_iii}
\end{align}
where $\mathcal{0}$ is the tensor of zeros. Since ${\mathcal{Z}_h}={\mathcal{Z}_h'}$, this leads to:
\begin{align}
    \mathscr{P}_{3}({\mathcal{Z}_h}) & = \mathscr{P}_{3}({\mathcal{Z}_h'})
    \nonumber\\
    & = \mathscr{P}_{3}({\mathcal{Z}_h'}) + \mathscr{P}_{3}(\mathcal{\Psi}) - \mathscr{P}_{3}(\mathcal{\Psi}') \,,
    \label{eq:thm1_iv}
\end{align}
which also implies that $\mathcal{Y}_m=\mathcal{Y}_m'$ (i.e., the model is not identifiable).

{\textbf{Proof of b):}} Note that since there is no additive noise, $\mathscr{A}$ is a linear operator and, thus, satisfies the following relation:
\begin{align}
    \mathscr{A}({\mathcal{Z}_h},\mathcal{\Psi}) - \mathscr{A}({\mathcal{Z}_h'},\mathcal{\Psi}') =
    \mathscr{A}({\mathcal{Z}_h}-{\mathcal{Z}_h'},\mathcal{\Psi}-\mathcal{\Psi}') \,.
\end{align}
By inspecting the definition of the equivalence relation in~\eqref{eq:thm_prel2_eq_rel}, it can be seen that the equivalence class in~\eqref{eq:thm_prel2_eq_class} is characterized by the
kernel of the operator $\mathscr{A}$, and can also be written as:
\begin{align}
    [({\mathcal{Z}_h},\mathcal{\Psi})]_{\sim_{{Z,\Psi}}}
    & = \big\{({\mathcal{Z}_h},\mathcal{\Psi}) + \mathcal{X} \,:\, \mathcal{X} \in \ker(\mathscr{A})\big\}.
    \label{eq:thm_prel2_eq_class2}
\end{align}

Now, suppose that we have two sets of HR images and variability factors belonging to different ECs, i.e., $({\mathcal{Z}_h},\mathcal{\Psi})\in[(\mathcal{Z}_0,\mathcal{\Psi}_0)]_{\sim_{{Z,\Psi}}}$, $({\mathcal{Z}_h'},\mathcal{\Psi}')\in[(\mathcal{Z}_0',\mathcal{\Psi}_0')]_{\sim_{{Z,\Psi}}}$, with $[(\mathcal{Z}_0,\mathcal{\Psi}_0)]_{\sim_{{Z,\Psi}}}\neq[(\mathcal{Z}_0',\mathcal{\Psi}_0')]_{\sim_{{Z,\Psi}}}$. Comparing the observations $(\mathcal{Y}_h,\mathcal{Y}_m)$ and $(\mathcal{Y}_h',\mathcal{Y}_m')$ generated by elements of each EC, we have:
\begin{align}
    (\mathcal{Y}_h,\mathcal{Y}_m)
    &= \mathscr{A}(({\mathcal{Z}_h},\mathcal{\Psi}) + \mathcal{X})
    \nonumber \\
    & \neq \mathscr{A}(({\mathcal{Z}_h'},\mathcal{\Psi}')  + \mathcal{X}')
    =  (\mathcal{Y}_h',\mathcal{Y}_m') \,,
\end{align}
for all $\mathcal{X},\mathcal{X}'\in\ker(\mathscr{A})$. Thus, elements selected from different ECs will always lead to different observations, ensuring that the EC is identifiable.
\end{proof}

Intuitively, {item a)} shows that variations in $\mathcal{\Psi}$ that occur in the nullspace of operators $\mathscr{P}_{1,2}$ and $\mathscr{P}_{3}$ are not reflected in the corresponding observations $(\mathcal{Y}_h,\mathcal{Y}_m)$ in such a way that they can be differentiated from possible changes in ${\mathcal{Z}_h}$.
More generally, changes of $\mathcal{\Psi}$ and ${\mathcal{Z}_h}$ that occur in the nullspace of matrix $\boldsymbol{A}$ do not affect $(\mathcal{Y}_h,\mathcal{Y}_m)${, which is clear from~\eqref{eq:thm1_mtxForm_1}}. This notion can be further extended by noting that changes occurring in the column space of $\boldsymbol{A}$ will certainly lead to different observations, which is made precise in {item b)}.

{Item b) in} Theorem~\ref{thm:identifiability_ECc_psi} allows us to characterize the ambiguities in the model {in more detail}. However, it is important to consider the characteristics of $\ker(\mathscr{A})$ in our problem to better understand the recoverability of the variability factor $\mathcal{\Psi}$. Let us consider the model in~\eqref{eq:thm1_mtxForm_1}, and two sets of variables $({\mathcal{Z}_h},\mathcal{\Psi})\sim_{{Z,\Psi}}({\mathcal{Z}_h'},\mathcal{\Psi}')$, belonging to the same EC.
It can be seen that to generate the same observations, the HR images need to satisfy ${\mathcal{Z}_h}-{\mathcal{Z}_h'}\in\ker(\mathscr{P}_{1,2})$, while the variability factors have to satisfy $\mathscr{P}_{3}({\mathcal{Z}_h}-{\mathcal{Z}_h'})=-\mathscr{P}_{3}(\mathcal{\Psi}-\mathcal{\Psi}')$. Therefore, the general form of the difference between the variability factors inside each equivalence class is of the form:
\begin{align}
    \mathcal{\Psi}-\mathcal{\Psi}' = -\underbrace{({\mathcal{Z}_h}-{\mathcal{Z}_h'})}_{\in\ker(\mathscr{P}_{1,2})} + \underbrace{\mathcal{X}}_{\in\ker(\mathscr{P}_{3})} \,.
    \label{eq:psi_EC_general_form}
\end{align}
The set of all possible $\mathcal{\Psi}-\mathcal{\Psi}'$ satisfying~\eqref{eq:psi_EC_general_form} is the sum of $\ker(\mathscr{P}_{3})$ and $\ker(\mathscr{P}_{1,2})$, which is given by $\ker(\mathscr{P}_{3}\circ\mathscr{P}_{1,2})$. We can readily see that $\mathcal{\Psi}$ cannot be recovered from the observations.
Only the spectrally degraded variability factors $\mathscr{P}_{3}(\mathcal{\Psi})$ can be uniquely recovered (which comes ``for free'' with the recovery of ${\mathcal{Z}_h}$ since it can be computed as $\mathscr{P}_{3}(\mathcal{\Psi})=\mathcal{Y}_m-\mathscr{P}_{3}({\mathcal{Z}_h})$). This makes it sufficient to study the capability of an algorithm to recover ${\mathcal{Z}_h}$ in our model. Since the matrices $\boldsymbol{P}_i$, $i\in\{1,2,3\}$ are essentially low-pass filtering and downsampling operators, their nullspaces intuitively encode high-frequency information along each tensor mode. Thus, only the smooth structure of $\mathcal{\Psi}$ can be identified uniquely from observations~$(\mathcal{Y}_h,\mathcal{Y}_m)$, since otherwise we cannot separate the effects of $\mathcal{\Psi}$ from ${\mathcal{Z}_h}$.

We also note that each EC in~\eqref{eq:thm_prel2_eq_class}, which contains all factors $\mathcal{\Psi}$ whose difference lies in the nullspace of the combined operator $\mathscr{P}_{1,2}\circ\mathscr{P}_{3}$, is strictly larger than if we considered changes that occur in the nullspace of each of these operators individually (i.e., $\mathscr{P}_{1,2}$ and $\mathscr{P}_{3}$).

Theorem~\ref{thm:identifiability_ECc_psi} guarantees that tensors belonging to different ECs will result in different observations, which is the minimal requirement for having identifiable ${\mathcal{Z}_h}$ and $\mathscr{P}_{3}(\mathcal{\Psi})$. However, the coresponding inverse problem still remains ill-posed as the number of {\color{red}unknowns} is greater than the number of observations.
Thus, stronger identifiability conditions cannot be obtained unless we provide stricter a priori characterizations of the sets $\Omega_Z$ and $\Omega_{\Psi}$.

\subsection{A Low-Multilinear-Rank Model}
\label{sec:low_rank_model}

One possible condition that can be imposed on the structures of both $\Omega_{\Psi}$ and $\Omega_Z$ is the low-rank tensor model. This kind of structure makes it possible to obtain identifiability and exact recovery guarantees for problem~\eqref{eq:HS_MS_fus_var_prob}, where spatial and spectral variabilities are present. Moreover, it also makes the problem well-posed and easier to solve since the number of unknowns becomes smaller than the amount of available data.

Suppose that ${\mathcal{Z}_h}$ and $\mathcal{\Psi}$ have multilinear ranks $(K_{Z,1},K_{Z,2},K_{Z,3})$ and $(K_{\Psi,1},K_{\Psi,2},K_{\Psi,3})$, respectively. This means that they can be represented as
\begin{align}
    {\mathcal{Z}_h} &= \big\ldbrack\mathcal{G}_{Z}; \boldsymbol{B}_{Z,1}, \boldsymbol{B}_{Z,2}, \boldsymbol{B}_{Z,3} \big\rdbrack \,,
    \label{eq:tucker_mdl_Z}
    \\
    \mathcal{\Psi} &= \big\ldbrack\mathcal{G}_{\Psi}; \boldsymbol{B}_{\Psi,1}, \boldsymbol{B}_{\Psi,2}, \boldsymbol{B}_{\Psi,3} \big\rdbrack \,,
    \label{eq:tucker_mdl_Psi}
\end{align}
where $\boldsymbol{B}_{Z,i}\in\amsmathbb{R}^{M_i\times K_{Z,i}}$, $\boldsymbol{B}_{\Psi,i}\in\amsmathbb{R}^{M_i\times K_{\Psi,i}}$, $i\in\{1,2\}$, $\boldsymbol{B}_{Z,3}\in\amsmathbb{R}^{L_h\times K_{Z,3}}$, $\boldsymbol{B}_{\Psi,3}\in\amsmathbb{R}^{L_h\times K_{\Psi,3}}$ are the factor matrices and $\mathcal{G}_{Z}\in\amsmathbb{R}^{K_{Z,1}\times K_{Z,2}\times K_{Z,3}}$, $\mathcal{G}_{\Psi}\in\amsmathbb{R}^{K_{\Psi,1}\times K_{\Psi,2}\times K_{\Psi,3}}$ are the core tensors.

Our objective is to study the identifiability and exact recovery of these variables given the observation model in~\eqref{eq:tensor_obs_mdl_hsi}--\eqref{eq:tensor_obs_mdl_msi_add}.
Using this model{, and applying the definition of the multilinear product and the properties of the mode-$k$ product defined in Section~\ref{sec:back_defs}}, the noiseless case of the degradation model~\eqref{eq:tensor_obs_mdl_hsi}--\eqref{eq:tensor_obs_mdl_msi_add} can be written as
\begin{align}
    \mathcal{Y}_h ={}&  \big\ldbrack\mathcal{G}_{Z}; \boldsymbol{P}_1 \boldsymbol{B}_{Z,1}, \boldsymbol{P}_2 \boldsymbol{B}_{Z,2}, \boldsymbol{B}_{Z,3} \big\rdbrack \,
    \label{eq:obs_mdl_hsi_lowrank}
    \\
    \mathcal{Y}_m ={}& \big\ldbrack\mathcal{G}_{Z}; \boldsymbol{B}_{Z,1}, \boldsymbol{B}_{Z,2}, \boldsymbol{P}_3 \boldsymbol{B}_{Z,3} \big\rdbrack
    \nonumber \\
    & + \big\ldbrack\mathcal{G}_{\Psi}; \boldsymbol{B}_{\Psi,1}, \boldsymbol{B}_{\Psi,2}, \boldsymbol{P}_3 \boldsymbol{B}_{\Psi,3} \big\rdbrack \,.
    \label{eq:obs_mdl_msi_lowrank_btd}
\end{align}

Note that we can represent the multispectral image model in~\eqref{eq:obs_mdl_msi_lowrank_btd} equivalently using a standard Tucker model as:
\begin{align}
    \mathcal{Y}_m ={}& \big\ldbrack \mathcal{C}_{m}; \boldsymbol{C}_{m,1}, \boldsymbol{C}_{m,2}, \boldsymbol{P}_3 \boldsymbol{C}_{m,3} \big\rdbrack \,,
    \label{eq:obs_mdl_msi_lowrank_hightucker}
\end{align}
where $\boldsymbol{C}_{m,i}$, $i\in\{1,2,3\}$ and $\mathcal{C}_m$ are the factor matrices and the core tensor of the MSI, which satisfy:
\begin{align}
    \mathcal{C}_{m} &= \mathcal{G}_{Z}\oplus\mathcal{G}_{\Psi} \,,   \label{eq:tucker_mdl_msi_var_ext_core}
    \\
    \boldsymbol{C}_{m,i} &= \big[\boldsymbol{B}_{Z,i} \, \boldsymbol{B}_{\Psi,i}\big],\,\,\, i\in\{1,2,3\} \,, \label{eq:tucker_mdl_msi_var_ext_factors}
\end{align}
{\color{red}where for two tensors $\mathcal{A}$ and $\mathcal{B}$, the binary operation $\mathcal{A}\oplus\mathcal{B}$ returns a block-diagonal tensor whose diagonal blocks are $\mathcal{A}$ and $\mathcal{B}$.}

{The model in equations~\eqref{eq:tucker_mdl_Z}--\eqref{eq:tucker_mdl_msi_var_ext_factors} will be subsequently used in Sections~\ref{sec:alg1_algebraic} and~\ref{sec:alg2_optimizationBased} to develop two image fusion algorithms, one algebraic (faster, but with stringent rank constraints) and another based on an optimization procedure (which allows for higher rank values). In each case, a new algorithm will be presented followed by its recoverability guarantees. It should be noted that in practice, the HRI $\mathcal{Z}_h$ and the variability tensor $\mathcal{\Psi}$ can have high rank. Nevertheless, we will perform only a coupled tensor approximation, with which we are able to capture most of the energy of the data even with insufficient ranks. In practice, however, higher rank models will be preferred to ensure the data is well represented and to avoid the presence artifacts in the reconstructed HRI (which will be achieved with the method of Section~\ref{sec:alg2_optimizationBased}).}

\section{An algebraic algorithm}
\label{sec:alg1_algebraic}

Considering the model in Section~\ref{sec:low_rank_model}, the image fusion problem consists in estimating the factors and core tensor $\mathcal{G}_Z$, $\boldsymbol{B}_{Z,i}$, $i\in\{1,2,3\}$. However, if the values composing the multilinear rank of ${\mathcal{Z}_h}$ are sufficiently low, those variables can be computed by solving the following coupled system of equations:
\begin{align}
\begin{cases}
    \mathcal{Y}_h &= \big\ldbrack\mathcal{G}_{Z}; \boldsymbol{C}_{h,1}, \boldsymbol{C}_{h,2}, \boldsymbol{B}_{Z,3} \big\rdbrack
    \\

    \mathcal{Y}_m &= \big\ldbrack\mathcal{C}_{m}; \boldsymbol{C}_{m,1}, \boldsymbol{C}_{m,2}, \boldsymbol{P}_3 \boldsymbol{C}_{m,3} \big\rdbrack
    \\
    \boldsymbol{C}_{h,i} &= \boldsymbol{P}_i \boldsymbol{B}_{Z,i}, \,\,\, i\in\{1,2\}
    \\
    \boldsymbol{C}_{m,i} &= \big[\boldsymbol{B}_{Z,i},\boldsymbol{B}_{\Psi,i}\big], \,\,\, i\in\{1,2,3\}
\end{cases},
\label{eq:alg1_coupled_eqs_formulation}
\end{align}
where $\boldsymbol{C}_{h,i}$, $i\in\{1,2\}$ denote the spatial factor matrices of the HSI, and the HRI is obtained from the solution of~\eqref{eq:alg1_coupled_eqs_formulation} as ${\mathcal{Z}_h}=\ldbrack\mathcal{G}_{Z};\boldsymbol{B}_{Z,1},\boldsymbol{B}_{Z,2},\boldsymbol{B}_{Z,3}\rdbrack$.

If we suppose that $K_{Z,i}+K_{\Psi,i}\leq N_i$, $i\in\{1,2\}$,~\eqref{eq:alg1_coupled_eqs_formulation} can be solved using an efficient, algebraic approach detailed in Alg.~\ref{alg:alg1}, which we call CT-STAR (Coupled Tucker decompositions for hyperspectral Super-resoluTion with vARiability). {The basic intuition behind this algorithm is to use the correspondence between the mode-1 and mode-2 matricizations of the HSI and MSI in order to separate the HRI from the variability tensor when computing its factor matrices. More details will be provided in the following.}

\begin{algorithm} [thb]
\footnotesize
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\caption{\mbox{Algebraic image fusion method (CT-STAR)}\label{alg:alg1}}
\Input{\mbox{Images $\mathcal{Y}_h$, $\mathcal{Y}_m$ ranks $K_{Z,i}$, $K_{\Psi,i}$, $i\in\{1,2,3\}$}}
\Output{HRI ${\widehat{\mathcal{Z}}_h}$, spectrally degraded variability factors $\mathscr{P}_{3}(\widehat{\mathcal{\Psi}})$}

Check if $K_{Z,i}+K_{\Psi,i}\leq N_i$, $i\in\{1,2\}$ \;

Compute $\widehat{\boldsymbol{C}}_{h,3} = \operatorname{tSVD}_{K_{Z,3}}(\unfold{\boldsymbol{Y}_{\!h}}{3})$ \;

Compute $\widehat{\boldsymbol{C}}_{m,i} = \operatorname{tSVD}_{K_{Z,i}+K_{\Psi,i}}(\unfold{\boldsymbol{Y}_{\!m}}{i})$ for $i\in\{1,2\}$ \;

Compute $\widetilde{\boldsymbol{Q}}_i$, for $i\in\{1,2\}$, as $\widetilde{\boldsymbol{Q}}_i = \big(\boldsymbol{P}_i \widehat{\boldsymbol{C}}_{m,i}\big)^{\dagger} \operatorname{tSVD}_{K_{Z,i}}(\unfold{\boldsymbol{Y}_{\!h}}{i})$\;

Compute $\widetilde{\boldsymbol{C}}_{m,i}=\widehat{\boldsymbol{C}}_{m,i} \widetilde{\boldsymbol{Q}}_i$, for $i\in\{1,2\}$ \;

Compute $\widehat{\mathcal{G}}_Z$ by solving
$(\widehat{\boldsymbol{C}}_{h,3} \otimes \boldsymbol{P}_2\widetilde{\boldsymbol{C}}_{m,2} \otimes \boldsymbol{P}_1\widetilde{\boldsymbol{C}}_{m,1} )\operatorname{vec}(\mathcal{G}_Z)=\operatorname{vec}(\mathcal{Y}_h)$\;

Compute ${\widehat{\mathcal{Z}}_h}=\ldbrack\widehat{\mathcal{G}}_{Z}; \widetilde{\boldsymbol{C}}_{m,1}, \widetilde{\boldsymbol{C}}_{m,2}, \widehat{\boldsymbol{C}}_{h,3}\rdbrack$ \;

Compute $\mathscr{P}_{3}(\widehat{\mathcal{\Psi}})=\mathcal{Y}_m-{\widehat{\mathcal{Z}}_h}\times_3\boldsymbol{P}_3$ \;
\end{algorithm}

It is important to note that CT-STAR does not \textcolor{magenta}{enforce} the block diagonal structure of the core tensor of the MSI (described in~\eqref{eq:tucker_mdl_msi_var_ext_core}). The following theorem gives a constructive proof of exact recovery conditions from which Alg.~\ref{alg:alg1} is derived.

\begin{theorem} \label{thm:exacr_rec_alg1}
Suppose that the HRI ${\mathcal{Z}_h}$ and the variability tensor $\mathcal{\Psi}$ have multilinear ranks $(K_{Z,1},K_{Z,2},K_{Z,3})$ and $(K_{\Psi,1},K_{\Psi,2},K_{\Psi,3})$, respectively, that $\mathcal{Y}_h$ and {\color{red}$\mathcal{Y}_m$} admit Tucker decompositions as denoted in~\eqref{eq:alg1_coupled_eqs_formulation}, that the observation noise is zero (i.e. $\mathcal{E}_h=\mathcal{0}$, $\mathcal{E}_m=\mathcal{0}$), and that
\begin{align}
    & \operatorname{rank}(\boldsymbol{P}_i \boldsymbol{B}_{Z,i}) = K_{Z,i}\,, \,\,\, i\in\{1,2\}
    \label{eq:thm3_rank_cond1}
    \\
    & \operatorname{rank}(\boldsymbol{P}_i \boldsymbol{B}_{\Psi,i}) \leq K_{\Psi,i}\,, \,\,\, i\in\{1,2\}
    \label{eq:thm3_rank_cond2}
    \\
    & \operatorname{rank}(\unfold{\boldsymbol{Y}_{\!h}}{i}) = K_{Z,i}\,, \,\,\, i\in\{1,2,3\}
    \label{eq:thm3_rank_cond3}
    \\
    & \operatorname{rank}(\unfold{\boldsymbol{Y}_{\!m}}{i}) = K_{Z,i}+K_{\Psi,i} \leq N_i\,, \,\,\, i\in\{1,2\}
    \label{eq:thm3_rank_cond4}
\end{align}
Then, if all columns in $\boldsymbol{P}_i \boldsymbol{B}_{Z,i}$ are linearly independent from those in $\boldsymbol{P}_i \boldsymbol{B}_{\Psi,i}$, for $i\in\{1,2\}$, Algorithm~\ref{alg:alg1} exactly recovers ${\mathcal{Z}_h}$ from the observations.
\end{theorem}

\begin{proof}
Let us compute matrices $\widehat{\boldsymbol{C}}_{m,i}$, $i\in\{1,2\}$, and $\widehat{\boldsymbol{C}}_{h,i}$, $i\in\{1,2,3\}$ as the left-singular vectors associated with the non-zero singular values of $\unfold{\boldsymbol{Y}_{\!m}}{i}$, $i\in\{1,2\}$ and $\unfold{\boldsymbol{Y}_{\!h}}{i}$, $i\in\{1,2,3\}$, respectively.
Then, due to~\eqref{eq:thm3_rank_cond3}--\eqref{eq:thm3_rank_cond4} and to the non-uniqueness of matrix decomposition, these matrices satisfy:
\begin{align}
    \widehat{\boldsymbol{C}}_{h,i}  &= \boldsymbol{P}_i \boldsymbol{B}_{Z,i} \boldsymbol{Q}_{h,i}, \,\,\, i\in\{1,2\}
    \label{eq:thm3_4a}
    \\
    \widehat{\boldsymbol{C}}_{h,3}  &= \boldsymbol{B}_{Z,3} \boldsymbol{Q}_{h,3}
    \label{eq:thm3_4b}
    \\
    \widehat{\boldsymbol{C}}_{m,i}  &= \boldsymbol{C}_{m,i} \boldsymbol{Q}_{m,i}, \,\,\, i\in\{1,2\}
    \label{eq:thm3_4c}
\end{align}
for invertible matrices $\boldsymbol{Q}_{h,i}\in\amsmathbb{R}^{K_{Z,i}\times K_{Z,i}}$, $i\in\{1,2,3\}$ and $\boldsymbol{Q}_{m,i}\in\amsmathbb{R}^{(K_{Z,i}+K_{\Psi,i})\times (K_{Z,i}+K_{\Psi,i})}$, $i\in\{1,2\}$.

Now, the main problem caused by the presence of variability is that matrices $\boldsymbol{Q}_{m,1}$ and $\boldsymbol{Q}_{m,2}$ preclude us from distinguishing the factors $\boldsymbol{B}_{Z,1}$ and $\boldsymbol{B}_{Z,2}$, associated with ${\mathcal{Z}_h}$, from $\boldsymbol{B}_{\Psi,1}$ and $\boldsymbol{B}_{\Psi,2}$, associated with $\mathcal{\Psi}$, using only information available in the MSI. These two become mixed in the spatial factors $\widehat{\boldsymbol{C}}_{m,i}$.
Nonetheless, consider the relationship between spatial degradation of the factors estimated from the MSI and the spatial factors of the HSI:
\begin{align}
    \boldsymbol{P}_i \widehat{\boldsymbol{C}}_{m,i} &= \boldsymbol{P}_i \boldsymbol{C}_{m,i} \boldsymbol{Q}_{m,i}
    \label{eq:thm3_5a}
    \\
    &= \big[\boldsymbol{P}_i \boldsymbol{B}_{Z,i}, \, \boldsymbol{P}_i \boldsymbol{B}_{\Psi,i}\big] \boldsymbol{Q}_{m,i} \,.
    \label{eq:thm3_5b}
\end{align}
Now, let us compute matrices $\widetilde{\boldsymbol{Q}}_i\in\amsmathbb{R}^{(K_{Z,i}+K_{\Psi,i})\times K_{Z,i}}$, $i\in\{1,2\}$ such that
\begin{align}
    & \widehat{\boldsymbol{C}}_{h,i} = \boldsymbol{P}_i \widehat{\boldsymbol{C}}_{m,i} \widetilde{\boldsymbol{Q}}_i \,.
    \label{eq:thm3_6a}
\end{align}
By partitioning the following matrix product as $\boldsymbol{Q}_{m,i}\widetilde{\boldsymbol{Q}}_i=[\overline{\boldsymbol{Q}}_{Z,i}^\top,\overline{\boldsymbol{Q}}_{\Psi,i}^\top]^\top$,~\eqref{eq:thm3_6a} can be represented as
\begin{align}
    \widehat{\boldsymbol{C}}_{h,i} &= \boldsymbol{P}_i \boldsymbol{B}_{Z,i} \boldsymbol{Q}_{h,i}
    \nonumber \\
    &= \boldsymbol{P}_i \boldsymbol{B}_{Z,i} \overline{\boldsymbol{Q}}_{Z,i} + \boldsymbol{P}_i \boldsymbol{B}_{\Psi,i} \overline{\boldsymbol{Q}}_{\Psi,i} \,.
    \label{eq:thm3_7a}
\end{align}

Since all columns in $\boldsymbol{P}_i \boldsymbol{B}_{Z,i}$ are linearly independent from those in $\boldsymbol{P}_i \boldsymbol{B}_{\Psi,i}$, equality~\eqref{eq:thm3_6a} will be satisfied if and only if the result of the product $\widehat{\boldsymbol{C}}_{m,i} \widetilde{\boldsymbol{Q}}_i$ (i.e., the right hand side of~\eqref{eq:thm3_7a}) does not contain any nontrivial linear combination of the columns of $\boldsymbol{P}_i \boldsymbol{B}_{\Psi,i}$. Thus, $\overline{\boldsymbol{Q}}_{\Psi,i}=\cb{0}$ and $\overline{\boldsymbol{Q}}_{Z,i}=\boldsymbol{Q}_{h,i}$ due to~\eqref{eq:thm3_rank_cond1} and~\eqref{eq:thm3_rank_cond3}.
This allows us to ``separate'' the variability and image factors as
\begin{align}
    \widetilde{\boldsymbol{C}}_{m,i} &= \widehat{\boldsymbol{C}}_{m,i} \widetilde{\boldsymbol{Q}}_i
    \nonumber\\
    & = \boldsymbol{B}_{Z,i} \boldsymbol{Q}_{h,i} \,,
\end{align}
for $i\in\{1,2\}$.
Now, consider the vectorization of the HSI as:
\begin{align}

    (\widehat{\boldsymbol{C}}_{h,3} \otimes \boldsymbol{P}_2\widetilde{\boldsymbol{C}}_{m,2} \otimes \boldsymbol{P}_1\widetilde{\boldsymbol{C}}_{m,1} )\operatorname{vec}(\mathcal{G}_Z) = \operatorname{vec}(\mathcal{Y}_h) \,.
    \label{eq:thm3_7b}
\end{align}
Since {\color{magenta} the matrix in the  left hand side of~\eqref{eq:thm3_7b}  has full column rank}, $\widehat{\mathcal{G}}_Z$ can be uniquely recovered from this equation, and will satisfy $\widehat{\mathcal{G}}_Z=\mathcal{G}_Z\times_1\boldsymbol{Q}_{h,1}^{-1}\times_2\boldsymbol{Q}_{h,2}^{-1}\times_3\boldsymbol{Q}_{h,3}^{-1}$.
The HRI and the spectrally degraded scaling factors are then finally recovered as:
\begin{align}
\begin{split}
    {\widehat{\mathcal{Z}}_h} &= \ldbrack\widehat{\mathcal{G}}_{Z}; \widetilde{\boldsymbol{C}}_{m,1}, \widetilde{\boldsymbol{C}}_{m,2}, \widehat{\boldsymbol{C}}_{h,3}\rdbrack
    \\
    &= \ldbrack \mathcal{G}_{Z}; \boldsymbol{B}_{Z,1}, \boldsymbol{B}_{Z,2}, \boldsymbol{B}_{Z,3}\rdbrack = \mathcal{Z}
\end{split}
\end{align}
and
\begin{align} \label{eq:recovering_low_psi1}
    \widehat{\mathcal{\Psi}\times_3\boldsymbol{P}_3}=\mathcal{Y}_m-{\widehat{\mathcal{Z}}_h}\times_3\boldsymbol{P}_3 \,,
\end{align}
which completes the proof.
\end{proof}

{Note that CT-STAR does not use the spectral degradation operation $\boldsymbol{P}_3$ to recover the HRI $\mathcal{Z}_h$, what makes it a spectrally ``blind'' algorithm.}
The CT-STAR algorithm is {also} fast (see Section~\ref{sec:comp}), but only works for the cases where the ranks of the spatial modes of ${\mathcal{Z}_h}$ are smaller than the dimensions of the HSI, which is quite restrictive. {This can make CT-STAR unsuited to process real images which can have high rank values, and motivates the search for a method with more flexibility in the selection of the ranks.} Moreover, both Alg.~\ref{alg:alg1} and Theorem~\ref{thm:exacr_rec_alg1}, in considering model~\eqref{eq:alg1_coupled_eqs_formulation}, made no assumptions about the (block diagonal) structure of the core tensor of the MSI. Although this led to more freedom from a modeling perspective, the recoverability conditions turned out to be restrictive.
{In the following section, we will explore the block diagonal structure of $\mathcal{C}_m$ using an optimization-based algorithm to address these limitations.}

\section{An optimization-based algorithm}
\label{sec:alg2_optimizationBased}

In this section, we pursue a different approach. Assume that model~\eqref{eq:obs_mdl_hsi_lowrank}--\eqref{eq:obs_mdl_msi_lowrank_btd} holds and that the values forming the multilinear ranks of both ${\mathcal{Z}_h}$ and $\mathcal{\Psi}$ are sufficiently low so that $\mathcal{Y}_m$ admits a block term decomposition (BTD) in the noiseless case~\cite{lathauwer2008tensor_BTD2_uniqueness}. We can then use uniqueness results thereof to guarantee the identifiability of ${\mathcal{Z}_h}$ under less restrictive conditions.

Let us consider the image fusion problem as the solution to the following optimization problem:
\begin{align}
     \mathop{\min}_{\boldsymbol{\Theta}} & \,\, J(\boldsymbol{\Theta}) \triangleq \Big\|\mathcal{Y}_h - \big\ldbrack {\mathcal{G}}_{Z}; \boldsymbol{P}_1 {\boldsymbol{B}}_{Z,1}, \boldsymbol{P}_2 {\boldsymbol{B}}_{Z,2}, {\boldsymbol{B}}_{Z,3} \big\rdbrack \Big\|_F^2
    \nonumber\\
    & \hspace{-1ex} + {\lambda}\,\bigg\|\mathcal{Y}_m - \sum_{\iota\in\{Z,\Psi\}} \big\ldbrack {\mathcal{G}}_{\iota}; {\boldsymbol{B}}_{\iota,1}, {\boldsymbol{B}}_{\iota,2}, \boldsymbol{P}_3 {\boldsymbol{B}}_{\iota,3} \big\rdbrack \bigg\|_F^2
    \label{eq:opt_prob_fus_btd}
\end{align}
where $\boldsymbol{\Theta}=\{{\mathcal{G}}_{\iota},{\boldsymbol{B}}_{\iota,i}:\iota\in\{Z,\Psi\}, \, i\in\{1,2,3\}\}$ {and $\lambda\in\amsmathbb{R}_+$ is a fixed parameter which balances the contribution of each term in the cost function. In the following, we first describe a procedure to solve problem~\eqref{eq:opt_prob_fus_btd} in Section~\ref{sec:CBSTAR_opt}, and later provide exact recovery guarantees in Section~\ref{sec:CBSTAR_recovery}.}

\begin{algorithm} [thb]
\footnotesize
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\caption{\mbox{Optimization-based image fusion (CB-STAR)}\label{alg:alg2}}
\Input{\mbox{Images $\mathcal{Y}_h$,$\mathcal{Y}_m$ ranks $K_{Z,i}$,$K_{\Psi,i}$, $i\in\{1,2,3\}${, iterations~$F$}}}
\Output{HRI ${\widehat{\mathcal{Z}}_h}$, spectrally degraded variability factors $\mathscr{P}_{3}(\widehat{\mathcal{\Psi}})$}

Initialize $\boldsymbol{\Theta}^{(0)}$ according to Section~\ref{sec:alg2_initialization}\;

\While{Stopping criteria is not satisfied}{

{Compute $\mathcal{G}_Z$ and $\boldsymbol{B}_{Z,i}$, $i\in\{1,2,3\}$ by solving~\eqref{eq:opt_subp_Z} with Algorithm~\ref{alg:alg3} in Appendix~\ref{sec:appendix1}, using $F$ iterations\;}

Compute $\boldsymbol{B}_{\Psi,1}$, $\boldsymbol{B}_{\Psi,2}$, $\boldsymbol{P}_3\boldsymbol{B}_{\Psi,3}$ and $\mathcal{G}_{\Psi}$ by solving \eqref{eq:HOSVD_optPsi} using the high-order SVD with rank $(K_{\Psi,1},K_{\Psi,2},K_{\Psi,3})$\;
}
Compute ${\widehat{\mathcal{Z}}_h}=\ldbrack\widehat{\mathcal{G}}_{Z}; \widehat{\boldsymbol{B}}_{Z,1}, \widehat{\boldsymbol{B}}_{Z,2}, \widehat{\boldsymbol{B}}_{Z,3}\rdbrack$ \;

Compute $\mathscr{P}_{3}(\widehat{\mathcal{\Psi}})=\mathcal{Y}_m-{\widehat{\mathcal{Z}}_h}\times_3\boldsymbol{P}_3$ \;
\end{algorithm}

\subsection{Optimization}
\label{sec:CBSTAR_opt}

{In order to minimize the cost function in~\eqref{eq:opt_prob_fus_btd}, we consider a block coordinate descent strategy, which successively minimizes $J$ w.r.t. $\mathcal{Z}_h$ (i.e., $\mathcal{G}_Z$ and $\boldsymbol{B}_{Z,i}$, $i\in\{1,2,3\}$) and w.r.t. $\mathcal{\Psi}$ (i.e., $\mathcal{G}_{\Psi}$ and $\boldsymbol{B}_{\Psi,i}$, $i\in\{1,2,3\}$), while keeping the remaining variables fixed.}
The optimization procedure is detailed in Alg.~\ref{alg:alg2}, which we call CB-STAR (Coupled Block term decompositions for hyperspectral Super-resoluTion with vARiability).

\begin{color}{orange}

\subsubsection{Optimizing w.r.t. $\mathcal{Z}_h$}
The optimization problem w.r.t. $\mathcal{Z}_h$ can be written as:
\begin{align} \label{eq:opt_subp_Z}
    \min_{\mathcal{G}_Z,\boldsymbol{B}_{Z,i}} \,\,
    & \Big\|\mathcal{Y}_h - \big\ldbrack {\mathcal{G}}_{Z}; \boldsymbol{P}_1{\boldsymbol{B}}_{Z,1}, \boldsymbol{P}_2{\boldsymbol{B}}_{Z,2}, {\boldsymbol{B}}_{Z,3}\big\rdbrack \Big\|_F^2
    \nonumber \\
    +\lambda & \Big\|\textcolor{magenta}{\mathcal{Y}}_0 - \big\ldbrack {\mathcal{G}}_{Z}; {\boldsymbol{B}}_{Z,1}, {\boldsymbol{B}}_{Z,2}, \boldsymbol{P}_3{\boldsymbol{B}}_{Z,3} \big\rdbrack \Big\|_F^2 \,,
\end{align}
where $\textcolor{magenta}{\mathcal{Y}}_0=\mathcal{Y}_m-\big\ldbrack {\mathcal{G}}_{\Psi}; {\boldsymbol{B}}_{\Psi,1}, {\boldsymbol{B}}_{\Psi,2}, \boldsymbol{P}_3{\boldsymbol{B}}_{\Psi,3}\big\rdbrack$. This is a variability-free, Tucker-based image fusion problem. We propose to solve~\eqref{eq:opt_subp_Z} using a block coordinate descent strategy w.r.t.
$\mathcal{G}_Z$ and $\boldsymbol{B}_{Z,i}$, $i\in\{1,2,3\}$, with a small number of iterations~$F$. This procedure is detailed in Alg.~\ref{alg:alg3} and in Appendix~\ref{sec:appendix1}.
An approximate closed form solution to~\eqref{eq:opt_subp_Z} can also be computed efficiently using the SCOTT algorithm~\cite{prevost2020coupledTucker_hyperspectralSRR_TSP}.
\end{color}

\subsubsection{Optimizing w.r.t. $\mathcal{\Psi}$}
This optimization problem can be written equivalently as
\begin{align} \label{eq:HOSVD_optPsi}
    \min_{\mathcal{G}_{\Psi},\boldsymbol{B}_{\Psi,i},\boldsymbol{X}_2} \,\, \Big\|\textcolor{magenta}{\mathcal{Y}}_1 - \big\ldbrack {\mathcal{G}}_{\Psi}; {\boldsymbol{B}}_{\Psi,1}, {\boldsymbol{B}}_{\Psi,2}, \boldsymbol{X}_2 \big\rdbrack \Big\|_F^2 \,,
\end{align}
where $\textcolor{magenta}{\mathcal{Y}}_1=\mathcal{Y}_{m}-\ldbrack {\mathcal{G}}_{Z}; {\boldsymbol{B}}_{Z,1}, {\boldsymbol{B}}_{Z,2}, \boldsymbol{P}_3 {\boldsymbol{B}}_{Z,3}\rdbrack$ and $\boldsymbol{X}_2=\boldsymbol{P}_3 {\boldsymbol{B}}_{\Psi,3}$. This problem can be solved by computing the high-order SVD of $\textcolor{magenta}{\mathcal{Y}}_1$ with rank $(K_{\Psi,1},K_{\Psi,2},K_{\Psi,3})$~\cite{lathauwer2000multilinearHOSVD}. Note that problem~\eqref{eq:HOSVD_optPsi} only returns $\boldsymbol{X}_2=\boldsymbol{P}_3 {\boldsymbol{B}}_{\Psi,3}$ instead of $\boldsymbol{B}_{\Psi,3}$. This is not a problem since the variations of $\boldsymbol{B}_{\Psi,3}$ in the nullspace of $\boldsymbol{P}_3$ are not identifiable.

\subsubsection{Initialization}
\label{sec:alg2_initialization}

Since this optimization problem is non-convex, the choice of initialization can have a significant impact on the performance of the algorithm. This can be particularly prominent in this algorithm since the model considered in~\eqref{eq:tensor_obs_mdl_hsi}--\eqref{eq:tensor_obs_mdl_msi_add} allows for a significant amount of ambiguity.
Fortunately, for practical scenes, we can consider a simple strategy to provide a reasonably accurate initial guess.

In the noiseless case, the following relation is satisfied:
\begin{align} \label{eq:initialization_psi_1}
    \mathscr{P}_{1,2}(\mathcal{Y}_m)
    & = \mathscr{P}_{1,2}\big(\mathscr{P}_{3}({\mathcal{Z}_h} + \mathcal{\Psi})\big)
    \\ \nonumber
    & = \mathscr{P}_{3}(\mathcal{Y}_h) + \mathscr{P}_{1,2}\big(\mathscr{P}_{3}(\mathcal{\Psi})\big) \,.
\end{align}
Thus, we can obtain a spatially and spectrally degraded version of $\mathcal{\Psi}$ directly from the HSI and MSI simply as:
\begin{align}
    \widetilde{\mathcal{\Psi}}
    & = \mathscr{P}_{1,2}(\mathcal{Y}_m) - \mathscr{P}_{3}(\mathcal{Y}_h)
    \nonumber \\
    & = \mathscr{P}_{1,2}\big(\mathscr{P}_{3}(\mathcal{\Psi})\big) \,.
    \label{eq:initialization_psi_2}
\end{align}
Then, if $\mathcal{\Psi}$ is smooth, we can spatially upscale $\widetilde{\mathcal{\Psi}}$ using some form of interpolation (e.g., bicubic), leading to $\mathscr{P}_{3}(\mathcal{\Psi}^{(0)})$.
Finally, we can use the Tucker decomposition of $\mathscr{P}_{3}(\mathcal{\Psi}^{(0)})$ to initialize $\mathcal{G}_{\Psi}$, $\boldsymbol{B}_{\Psi,i}$, $i\in\{1,2,3\}$. We call this the \emph{interpolation} initialization.

Another option is to try to invert~\eqref{eq:initialization_psi_2} using the pseudoinverse of operator $\mathscr{P}_{1,2}$, which can be computed using properties of the tensor vectorization and Kronecker product:
\begin{align}
    \big(\mathscr{P}_{1,2}\big)^{\dagger} =
    (\cdot) \times_1 \boldsymbol{P}_1^\dagger \times_2 \boldsymbol{P}_2^\dagger \,,
\end{align}
where $\boldsymbol{X}^{\dagger}$ denotes the pseudoinverse of $\boldsymbol{X}$. the initialization can be then computed as $\mathscr{P}_{3}(\mathcal{\Psi}^{(0)})=\widetilde{\mathcal{\Psi}}\times_1 \boldsymbol{P}_1^\dagger \times_2 \boldsymbol{P}_2^\dagger$. We call this the \emph{pseudoinverse} initialization.

The initialization of $\boldsymbol{B}_{Z,i}$, $i\in\{1,2,3\}$ can then be performed as ${\boldsymbol{B}}_{Z,3} = \operatorname{tSVD}_{K_{Z,3}}(\unfold{\boldsymbol{Y}_{\!h}}{3})$ and ${\boldsymbol{B}}_{Z,i} = \operatorname{tSVD}_{K_{Z,i}}(\unfold{\boldsymbol{X}}{i})$ for $i\in\{1,2\}$, where $\mathcal{X}=\mathcal{Y}_m-\mathscr{P}_{3}(\mathcal{\Psi}^{(0)})$.

\subsection{Exact Recovery}
\label{sec:CBSTAR_recovery}

Suppose that $K_{Z,i}=K_{\Psi,i}$, $i\in\{1,2,3\}$, without loss of generality, so that the MSI follows a standard BTD as considered in~\cite{lathauwer2008tensor_BTD2_uniqueness}. Then we have the following result regarding the identifiability of the proposed algorithm.

\begin{theorem} \label{thm:exacr_rec_alg2}
Suppose that $K_i\equiv K_{Z,i}=K_{\Psi,i}$, $i\in\{1,2,3\}$, that the observations are noise free (i.e., $\mathcal{E}_h=\mathcal{0}$, $\mathcal{E}_m=\mathcal{0}$), that $\{\mathcal{G}_{\iota},\boldsymbol{B}_{\iota,i}:\iota\in\{Z,\Psi\},\,i\in\{1,2,3\}\}$ are drawn from some joint absolutely continuous distribution, and that the following conditions on the dimensions hold:
\begin{align}
    M_1 \ge 2 K_1  & \mbox{ and } M_2 \ge 2 K_2
    \label{eq:dimension_conditions_thm_alg2_c}

    \\
    K_{Z,3} \leq \min & \big\{ N_1 N_2, \, K_{Z,1}K_{Z,2} \big\}
    \label{eq:dimension_conditions_thm_alg2_ef}
\end{align}
and either one of the following:
\begin{align}
\begin{cases}
    \begin{cases}
    K_3 > K_1+K_2-2 \,, &  \text{or} \\

    |K_1-K_2| > K_3-2 \,, & \\
    \end{cases}
    \\
    \mbox{ and } L_m \ge 2 K_3
\end{cases}
\label{eq:dimension_conditions_thm_alg2_a}
\end{align}
or
\begin{align}
\begin{split}
    K_1=K_2 \,,\,\,\, K_3 \geq 3 \,\,\,\mbox{and}\,\,\, K_3 < L_m \,,
\end{split}
\label{eq:dimension_conditions_thm_alg2_b}
\end{align}
is satisfied. Then, the solution to optimization problem~\eqref{eq:opt_prob_fus_btd} satisfies {$\widehat{\mathcal{Z}}_h=\mathcal{Z}_h$} almost surely.
\end{theorem}

\begin{proof}
Since there is no additive noise, the optimal solution to optimization problem~\eqref{eq:opt_prob_fus_btd} will necessarily make both terms in the cost function equal to zero. This implies that
\begin{align}
    \mathcal{Y}_h & = \ldbrack \widehat{\mathcal{G}}_{Z}; \boldsymbol{P}_1\widehat{\boldsymbol{B}}_{Z,1}, \boldsymbol{P}_2\widehat{\boldsymbol{B}}_{Z,2}, \widehat{\boldsymbol{B}}_{Z,3} \rdbrack \,,
    \\
    \mathcal{Y}_m & = \sum_{\iota\in\{Z,\Psi\}} \ldbrack \widehat{\mathcal{G}}_{\iota}; \widehat{\boldsymbol{B}}_{\iota,1}, \widehat{\boldsymbol{B}}_{\iota,2}, \boldsymbol{P}_3 \widehat{\boldsymbol{B}}_{\iota,3} \rdbrack \,,
    \label{eq:recovered_MSI_BTD_opt}
\end{align}
where $\{\widehat{\mathcal{G}}_{\iota},\widehat{\boldsymbol{B}}_{\iota,i}:\iota\in\{Z,\Psi\},\,i\in\{1,2,3\}\}$ denotes a solution to~\eqref{eq:opt_prob_fus_btd}.

Since the dimension conditions~\eqref{eq:dimension_conditions_thm_alg2_c} and either~\eqref{eq:dimension_conditions_thm_alg2_a} or~\eqref{eq:dimension_conditions_thm_alg2_b} are satisfied and the core tensor and factor matrices are drawn from joint absolutely continuous distributions, the BTD decomposition of the MSI in~\eqref{eq:recovered_MSI_BTD_opt} is essentially unique according to Theorems~5.1 and~5.5 in~\cite{lathauwer2008tensor_BTD2_uniqueness}. This means that the following conditions are satisfied:
\begin{align}
    \widehat{\boldsymbol{B}}_{Z,i} & = \boldsymbol{B}_{\iota_1,i} \boldsymbol{Q}_{\iota_1,i}, \,\, i\in\{1,2\}
    \label{eq:ambig_thm2_1}
    \\
    \widehat{\boldsymbol{B}}_{\Psi,i} & = \boldsymbol{B}_{\iota_2,i} \boldsymbol{Q}_{\iota_2,i}, \,\, i\in\{1,2\}
    \label{eq:ambig_thm2_2}
    \\
    \boldsymbol{P}_3\widehat{\boldsymbol{B}}_{Z,3} & = \boldsymbol{P}_3 \boldsymbol{B}_{\iota_1,3} \boldsymbol{Q}_{\iota_1,3}
    \label{eq:ambig_thm2_3}
    \\
    \boldsymbol{P}_3\widehat{\boldsymbol{B}}_{\Psi,3} & = \boldsymbol{P}_3 \boldsymbol{B}_{\iota_2,3} \boldsymbol{Q}_{\iota_2,3}
    \label{eq:ambig_thm2_4}
    \\
    \widehat{\mathcal{G}}_{Z} & = \mathcal{G}_{\iota_1} \times_1 \boldsymbol{Q}_{\iota_1,1}^{-1} \times_2 \boldsymbol{Q}_{\iota_1,2}^{-1} \times_3 \boldsymbol{Q}_{\iota_1,3}^{-1}
    \label{eq:ambig_thm2_5}
    \\
    \widehat{\mathcal{G}}_{\Psi} & = \mathcal{G}_{\iota_2} \times_1 \boldsymbol{Q}_{\iota_2,1}^{-1} \times_2 \boldsymbol{Q}_{\iota_2,2}^{-1} \times_3 \boldsymbol{Q}_{\iota_2,3}^{-1}
    \label{eq:ambig_thm2_6}
\end{align}
{\color{magenta} To account for the possible permutation of the two BTD terms,  indexes $\iota_1$ and $\iota_2$ can be either $(\iota_1,\iota_2)=(Z,\Psi)$ or $(\iota_1,\iota_2)=(\Psi,Z)$,} and $\boldsymbol{Q}_{Z,i}$, $\boldsymbol{Q}_{\Psi,i}$, $i\in\{1,2,3\}$ are invertible matrices of appropriate size, {\color{magenta} which account for rotational {and scaling} ambiguities of the model.}

{\color{magenta} Let us consider the mode-3 unfolding of the spectrally degraded HSI:}
\begin{align*}
    \boldsymbol{P}_3\unfold{{\boldsymbol{Y}_{\!h}}}{3} & = \boldsymbol{P}_3\widehat{\boldsymbol{B}}_{Z,3} \underbrace{\unfold{\mathop{\widehat{\boldsymbol{G}}_{Z}}}{3} \big(\boldsymbol{P}_2\widehat{\boldsymbol{B}}_{Z,2} \otimes \boldsymbol{P}_1\widehat{\boldsymbol{B}}_{Z,1}\big)^\top}_{\,\widehat{\!\boldsymbol{X}}}
    \\
    & = \boldsymbol{P}_3{\boldsymbol{B}}_{Z,3} \underbrace{\unfold{\mathop{{\boldsymbol{G}}_{Z}}}{3} \big(\boldsymbol{P}_2{\boldsymbol{B}}_{Z,2} \otimes \boldsymbol{P}_1{\boldsymbol{B}}_{Z,1}\big)^\top}_{\boldsymbol{X}} \,,
\end{align*}
where $\,\widehat{\!\boldsymbol{X}}$ and ${\boldsymbol{X}}$ are generically full row rank by the assumption~\eqref{eq:dimension_conditions_thm_alg2_ef} on the ranks and on the dimensions.
Therefore, we have
\begin{align}
    \operatorname{span}\big( \boldsymbol{P}_3\widehat{\boldsymbol{B}}_{Z,3} \big) = \operatorname{span}\big( \boldsymbol{P}_3{\boldsymbol{B}}_{Z,3} \big) \,.
    \label{eq:thm2_btd_span_eq}
\end{align}

However, since $K_3<L_m$ and due to the distributional assumptions on the factor matrices, we have that, generically, matrices $\boldsymbol{P}_3\boldsymbol{B}_{Z,3}$ and $\boldsymbol{P}_3\boldsymbol{B}_{\Psi,3}$ both have rank $K_3$, and matrix $\begin{bmatrix}
\boldsymbol{P}_3 {\boldsymbol{B}}_{Z,3}  &  \boldsymbol{P}_3\boldsymbol{B}_{\Psi,3}
\end{bmatrix}$ is full column rank (i.e., it has rank greater than $K_3$). Therefore, the subspaces spanned by $\boldsymbol{P}_3\boldsymbol{B}_{Z,3}$ and $\boldsymbol{P}_3\boldsymbol{B}_{\Psi,3}$ are different, and it is not possible to have $\boldsymbol{P}_3 {\boldsymbol{B}}_{Z,3}  =  \boldsymbol{P}_3\boldsymbol{B}_{\Psi,3} \boldsymbol{S}$, for any matrix $\boldsymbol{S}$.
Thus,~\eqref{eq:thm2_btd_span_eq} implies that the equation $\boldsymbol{P}_3\widehat{\boldsymbol{B}}_{Z,3}=\boldsymbol{P}_3\boldsymbol{B}_{\Psi,3}\boldsymbol{Q}_{\Psi,3}$ is not possible (i.e., it cannot be satisfied for any matrix $\boldsymbol{Q}_{\Psi,3}$). This ensures, due to \eqref{eq:ambig_thm2_3} (and to the essential uniqueness of the MSI BTD), that we must have $\iota_1=Z$ in~\eqref{eq:ambig_thm2_1},~\eqref{eq:ambig_thm2_3} and~\eqref{eq:ambig_thm2_5}, (i.e., the factor matrices related to $\mathcal{\Psi}$ can not fit the HSI), which shows that the correct permutation of the BTD terms is selected.

Finally, since $\,\widehat{\!\boldsymbol{X}}$ is generically full row rank, (or since $\unfold{{\boldsymbol{Y}_{\!h}}}{3}$ has rank $K_3$), we have that $\widehat{\boldsymbol{B}}_{Z,3}=\boldsymbol{B}_{Z,3}\boldsymbol{S}$ for some $\boldsymbol{S}$, and~\eqref{eq:ambig_thm2_3} ensures that $\boldsymbol{S}=\boldsymbol{Q}_{Z,3}$. This means that $\widehat{\boldsymbol{B}}_{Z,3}=\boldsymbol{B}_{Z,3}\boldsymbol{Q}_{Z,3}$ (almost surely), and, using~\eqref{eq:ambig_thm2_1} and~\eqref{eq:ambig_thm2_5} and $\iota_1=Z$, the reconstructed image consequently satisfies
\begin{align}
    {\widehat{\mathcal{Z}}_h} & = \big\ldbrack \widehat{\mathcal{G}}_{Z}; \widehat{\boldsymbol{B}}_{Z,1}, \widehat{\boldsymbol{B}}_{Z,2}, \widehat{\boldsymbol{B}}_{Z,3} \big\rdbrack
    \nonumber \\
    & = \big\ldbrack {\mathcal{G}}_{Z}; {\boldsymbol{B}}_{Z,1}, {\boldsymbol{B}}_{Z,2}, {\boldsymbol{B}}_{Z,3} \big\rdbrack
    = {\mathcal{Z}_h},
\end{align}
(almost surely), which concludes the proof.
\end{proof}

By taking the block diagonal structure of $\mathcal{C}_m$ into account, Theorem~\ref{thm:exacr_rec_alg2} obtains generally less restrictive recovery conditions. Comparing Theorems~\ref{thm:exacr_rec_alg1} and~\ref{thm:exacr_rec_alg2}, we can see that: 1) conditions~\eqref{eq:thm3_rank_cond3} and~\eqref{eq:dimension_conditions_thm_alg2_ef} are equivalent; 2) the conditions for the spectral ranks are not directly comparable but are similarly restrictive for both theorems; and, most notably, 3) The constraint~\eqref{eq:thm3_rank_cond4} on spatial ranks is much more restrictive when compared to the one in~\eqref{eq:dimension_conditions_thm_alg2_c}, required by Theorem~\ref{thm:exacr_rec_alg2}. {This shows that, although computationally more demanding, CB-STAR has more flexibility and may be able to deliver \textcolor{darkgreen}{better performance} when compared to CT-STAR for images with complex spatial content.}

\section{{Computational complexity}}\label{sec:comp}

{The computational complexity of the algorithms is given as follows. The total cost involved with the main operations in CT-STAR (Alg.~\ref{alg:alg1}) are the following: 1) computation of the truncated SVDs in steps~2 and~3, which requires $\mathcal{O}(\max\{K_{Z,1},K_{Z,2}\}M_1M_2L_m+K_{Z,3}N_1N_2L_h)$ flops, 2) computation of the equations in steps~4 and~5, which requires~$\mathcal{O}(\max\{K_{Z,1},K_{Z,2}\}N_1N_2L_h+M_1N_1K_{Z,1}+M_2N_2K_{Z,2})$ flops, and solution of the linear equation in step~6, which requires~$\mathcal{O}\big(L_hN_1N_2(K_{Z,1}K_{Z,2}K_{Z,3})^2\big)$ flops.}

{ For CB-STAR (Alg.~\ref{alg:alg2}), the total cost involved in each iteration is 1) solving problem~\eqref{eq:opt_subp_Z} in step~3 using Alg.~\ref{alg:alg3}, whose main costs are based on the solution to Sylvester equations, and are given by $\mathcal{O}\big(F(M_1^3+M_2^3+L_h^3+(K_{Z,2}K_{Z,3})^3)\big)$ (where we assume $K_{Z,1}\geq K_{Z,3}$ for simplicity), and 2) computing the high-order SVD in step~4, which costs~$\mathcal{O}\big(\max_i\{K_{\Psi,i}\}M_1M_2L_m\big)$ flops.
}

\section{Experiments}
\label{sec:results}

In this section, the performance of the proposed approach is illustrated through numerical experiments considering both synthetic and real data containing spatial and spectral variability.
All simulations were coded in MATLAB and run on a desktop with a 4.2 GHz Intel Core i7 and 16GB RAM.

\subsection{Experimental Setup} \label{sec:sim_setup}

We compared CT-STAR and CB-STAR to both matrix and tensor factorization-based algorithms. Among the matrix factorization-based methods, we considered the HySure~\cite{simoes2015HySure} and CNMF~\cite{yokoya2012coupledNMF} methods, the FuVar~\cite{Borsoi_2018_Fusion} method, which accounts for spectral variability, and the multiresolution analysis-based GLP-HS algorithm~\cite{aiazzi2006GLP_HS}. We also considered the {LTMR~\cite{dian2019hyperspectralSRR_tensorMultiRank},} STEREO~\cite{kanatsoulis2018hyperspectralSRR_coupledCPD} and SCOTT~\cite{prevost2020coupledTucker_hyperspectralSRR_TSP} algorithms, {which are tensor-based image fusion methods.}

The real {HRIs and MSIs}, which were acquired at different time instants but at the same spatial resolution, were pre-processed as described in~\cite{simoes2015HySure}. {This consisted} in the manual removal of water absorption and low-SNR bands, followed by the normalization of all bands of the {HRIs and MSIs} such that the 0.999 intensity quantile {corresponded} to a value of~1.
Afterwards, the {HRIs were} denoised (as described in~\cite{roger1996denoisingHSI}) to yield the high-SNR reference image~${\mathcal{Z}_h}$~\cite{yokoya2017HS_MS_fusinoRev}.
The observed HSIs $\mathcal{Y}_{h}$ were then generated from ${\mathcal{Z}_h}$ by applying a separable degradation operator, with $\boldsymbol{P}_1=\boldsymbol{P}_2$ (a Gaussian filter with unity variance followed by a subsampling with a decimation factor of two\footnote{Details on how to construct $\boldsymbol{P}_1$ and $\boldsymbol{P}_2$ can be found in~\cite{prevost2020coupledTucker_hyperspectralSRR_TSP}.}). Gaussian noise was also added to obtain an SNR of 30dB.
The observed MSIs $\mathcal{Y}_{m}$ were generated by adding noise to the reference MSI to obtain an SNR of 40dB. The spectral response function $\boldsymbol{P}_3$ was obtained from calibration measurements and known a priori\footnote{Available for download \href{https://earth.esa.int/web/sentinel/user-guides/sentinel-2-msi/document-library/-/asset_publisher/Wk0TKajiISaR/content/sentinel-2a-spectral-responses}{here}.}.
{Note that the HRI $\mathcal{Z}_h$ is not known by the algorithms, and is just used to assess their performance during the experiments using quantitative metrics.}

The parameters of the algorithms were selected as follows.
We selected the ranks and regularization parameters for HySure, CNMF, FuVar and LTMR according to the original works~\cite{simoes2015HySure,yokoya2012coupledNMF,Borsoi_2018_Fusion,dian2019hyperspectralSRR_tensorMultiRank}. For STEREO, we selected the rank in the interval $[5,80]$ which led to the best reconstruction results. Similarly, for SCOTT and for the proposed algorithms, we selected the spatial ranks in the intervals $[10,80]$ and the spectral ranks in the interval $[2,30]$, which led to the best reconstruction results. {For simplicity, we also set $\lambda=1$ for CB-STAR.} The spatial and spectral degradation operators (or, equivalently, the blurring kernels for HySure, CNMF and FuVar) were assumed to be known a priori for all methods.
The BCD procedure in Alg.~\ref{alg:alg2} was performed until the relative change in the objective function value was smaller than~$10^{-3}$. {At each iteration of Alg.~\ref{alg:alg2}, \textcolor{darkgreen}{we ran Alg.~\ref{alg:alg3} for} one single inner iteration (i.e., $F=1$), which resulted in good experimental performance with moderate execution times.} Both the \emph{interpolation} and the \emph{pseudoinverse} initializations described in Section~\ref{sec:alg2_initialization} were considered, but only the first one (which performed better) is shown in the visual results.

To evaluate the quality of the reconstructed images ${\widehat{\mathcal{Z}}_h}$, we considered four quantitative metrics, which were previously used in~\mbox{\cite{yokoya2017HS_MS_fusinoRev,simoes2015HySure,Borsoi_2018_Fusion}}.
The first metric is the peak signal to noise ratio (PSNR), defined as
\begin{align}
	{\rm PSNR}

	{}={} \frac{10}{L} \sum_{\ell=1}^L
    \log_{10} \Bigg(
    \frac{M_1 M_2 \operatorname{E}\!\big\{\!\max\big({[\mathcal{Z}_h]_{:,:,\ell}}\big)\big\}}{\big\|{[\mathcal{Z}_h]_{:,:,\ell}  -[\widehat{\mathcal{Z}}_h]_{:,:,\ell}\big\|_F^2}} \Bigg)
    \nonumber\,,
\end{align}
where $\operatorname{E}\{\cdot\}$ denotes the expectation operator.

The second metric is the Spectral Angle Mapper (SAM):
\begin{align}
	{\rm SAM}

	{}={} \frac{1}{M_1 M_2} \sum_{n,m} \arccos \Bigg( \frac{{[\mathcal{Z}_h]_{:,n,m}^\top[\widehat{\mathcal{Z}}_h]_{:,n,m}}}
    {\big\|{[\mathcal{Z}_h]_{:,n,m}}\big\|_2\big\|{[\widehat{\mathcal{Z}}_h]_{:,n,m}}\big\|_2} \Bigg)
    \nonumber\,.
\end{align}

The ERGAS~\cite{wald2000qualityERGAS} metric provides a global statistical measure of the quality of the fused data, and is defined as:
\begin{align}
	{\rm ERGAS}

	{}={} \frac{M_1M_2}{N_1N_2} \sqrt{\!\frac{10^4}{L_h} \!
    \sum_{\ell=1}^{L_h}  \frac{\big\|{[\mathcal{Z}_h]_{:,:,\ell}
    -[\widehat{\mathcal{Z}}_h]_{:,:,\ell}}\big\|_F^2}{(\cb{1}^\top{[\mathcal{Z}_h]_{:,:,\ell}}\cb{1}/(M_1M_2))^2}}
    \nonumber \,.
\end{align}

The last metric is the average of the bandwise UIQI~\cite{wang2002qualityUIQI}, which evaluates image distortions caused by loss of correlation and by luminance and contrast distortion, with value approaching one as ${\widehat{\mathcal{Z}}_h}$ approaches ${\mathcal{Z}_h}$.

We also evaluate the reconstructed images visually, by displaying true- and pseudo-color representations of the visual and infrared spectra of~${\widehat{\mathcal{Z}}_h}$ (corresponding to the wavelengths $0.45$, $0.56$ and $0.66~\mu m$, and $0.80$, $1.50$ and~$2.20~\mu m$, respectively). Due to space limitations, we only display the results of FuVar, {LTMR,} STEREO, SCOTT, CT-STAR and CB-STAR, since these are the methods which performed best, and the ones which were conceptually closest to our approach.
The spectrally degraded additive factors $\mathscr{P}_{3}(\widehat{\mathcal{\Psi}})$ estimated by CT-STAR and CB-STAR are also evaluated visually, through pseudo-color representations of its visible and infrared spectra, and by the norm (over all bands) of each of its pixels.

\begin{table}[ht!]
\caption{Results - synthetic example (``rk'' stands for ``ranks'')} \vspace{-5pt}
\centering
\renewcommand{\arraystretch}{1.15}
\setlength\tabcolsep{3.5pt}
\resizebox{\linewidth}{!}{
\pgfplotstabletypeset[header=false,
col sep = tab,
	columns/0/.style={string type},
every head row/.style={before row={\hline}, after row=\hline},
columns/0/.style={string type, column name={Algorithm}, column type/.add={@{}@{\,}}{}},
columns/1/.style={column name={SAM}, column type/.add={@{\,}|@{\,}}{}},
columns/2/.style={column name={ERGAS}, column type/.add={@{\,}|@{\,}}{}},
columns/3/.style={column name={PSNR}, column type/.add={@{\,}|@{\,}}{}},
columns/4/.style={column name={UIQI}, column type/.add={@{\,}|@{\,}}{}},
columns/5/.style={column name={time}, fixed, precision=2, column type/.add={@{\,}|@{\,}}{}},
columns ={0,1,2,3,4,5},
every row 9 column 1/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 13 column 1/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 13 column 2/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 13 column 3/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 7 column 4/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 8 column 4/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 9 column 4/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 10 column 4/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 11 column 4/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 12 column 4/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 13 column 4/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 6 column 5/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
]{tables/tab_synthetic_MC.txt}}
\label{tab:ex0_bis}
\end{table}

\subsection{Examples -- Synthetic data}

To evaluate the {proposed} algorithms in a controlled scenario, we first considered a simulation with synthetic data.
The tensors ${\mathcal{Z}_h}$ and $\mathcal{\Psi}$, of dimensions $100\times100\times200$, were generated following the Tucker model, with uniformly distributed entries on the interval $[0,1]$ and ranks $(10,10,5)$ and $(5,5,3)$, respectively.
The spectral response function $\boldsymbol{P}_3\in\amsmathbb{R}^{10\times200}$ was constructed by uniformly averaging groups of $20$ bands, and the rest of the simulation setup was the same as described in Section~\ref{sec:sim_setup}. For this experiment, we initialized CB-STAR with the results of CT-STAR.
{\textcolor{darkgreen}{We considered} two examples. First, \textcolor{darkgreen}{we compared} the proposed method to other state-of-the-art algorithms, for different choices of rank. Afterwards, \textcolor{darkgreen}{we evaluated} the sensitivity of CT-STAR and CB-STAR to the presence of additive noise. In both cases, we report average results of a Monte Carlo simulation with 100 \textcolor{magenta}{noise} realizations.}

\subsubsection{Comparison to other algorithms}

{For this comparison, we set the} ranks of STEREO and SCOTT were $50$ and $(60,60,5)$, respectively. {\textcolor{darkgreen}{We also ran} the proposed methods with four different rank values (smaller, equal, and larger than the true data ranks)}, indicated in Table~\ref{tab:ex0_bis}.
The results in Table~\ref{tab:ex0_bis} show that the proposed methods yielded significant improvements when compared to the other algorithms, which is expected since this dataset was generated according to the model~\eqref{eq:tensor_obs_mdl_hsi}--\eqref{eq:tensor_obs_mdl_msi_add}. Moreover, the performance of both CT-STAR and CB-STAR as a function of the ranks was similar, with the best results {observed} when the rank was the same as the ground truth, but with similar performance when the rank was underestimated. {This indicates that CT-STAR and CB-STAR can still perform well when the data rank is higher than the one specified for the model.} When the selected rank was overestimated, the performance of the proposed methods degraded more sharply (with a more prominent decrease for CB-STAR), indicating that the ranks should not be much greater than the true values in order to obtain the best performance.

\begin{color}{orange}
\subsubsection{Effect of noise}
The exact recovery results obtained in Theorems~\ref{thm:exacr_rec_alg1} and~\ref{thm:exacr_rec_alg2} assume a noiseless observation model, which is not the case in real applications.
To illustrate how the performance of the proposed methods is affected by the presence of additive noise~$\mathcal{E}_h$ and~$\mathcal{E}_m$, \textcolor{darkgreen}{we evaluated} the quantitative reconstruction metrics for different SNRs, varying from 0dB up to the noiseless case (SNR = $\infty$). For simplicity, we set $K_Z$ and $K_{\Psi}$ as the correct ranks of $\mathcal{Z}_h$ and $\mathcal{\Psi}$, respectively. The results are shown in Table~\ref{tab:exp_synth_SNRs}. It can be seen that the performance of both CT-STAR and CB-STAR decreased for low SNRs. Moreover, CB-STAR achieved better results for high SNRs ($\geq$40dB), but was surpassed by CT-STAR for low SNR scenarios ($\leq$30dB), whose results were considerably better for extremely noisy scenarios (0dB). Moreover, for SNRs found in typical scenes ($\geq$20dB), \textcolor{darkgreen}{these results indicate that} both methods are able to obtain satisfactory results.
Finally, we note that since the ranks of $\mathcal{Z}_h$ and $\mathcal{\Psi}$ \textcolor{darkgreen}{satisfied the} requirements of Theorems~\ref{thm:exacr_rec_alg1} and~\ref{thm:exacr_rec_alg2}, exact reconstruction was obtained for SNR~=~$\infty$ (i.e., the original and reconstructed images were equal up to machine precision).

The theoretical results presented in this work \textcolor{darkgreen}{made use of two hypotheses} to guarantee the recovery of the HRI, namely, that the image and variability factors have low multilinear rank, and that there is no observation noise. However, those hypotheses are often not strictly satisfied in practice. Nonetheless, CT-STAR and CB-STAR only perform a low-rank approximation of the data. Thus, they can still perform well when the true rank of the image is larger than the one specified to the algorithm and under the presence of noise, as illustrated in the previous simulations (even though exact recovery is not guaranteed in those cases). In the following section, we will evaluate the performance of the proposed methods with real HSIs and MSIs, which do not necessarily come from the specified low-rank models.

\end{color}

\begin{table}[th]
\renewcommand{\arraystretch}{1.15}
\setlength\tabcolsep{3.5pt}
\centering
\caption{{Performance of the proposed methods for different SNRs}}  \vspace{-5pt}
\resizebox{\linewidth}{!}{
\begin{tabular}{c|cccccccccc} \hline
\multicolumn{10}{c}{CT-STAR} \\\hline
SNR	&	$0$	&	$10$	&	$20$	&	$30$	&	$40$	&	$60$	&	$80$	&	$100$	&	$\infty$	\\\hline
SAM	&	$9.622$	&	$3.191$	&	$1.197$	&	$0.647$	&	$0.370$	&	$0.033$	&	$0.003$	&	$0$	&	$0$	\\
ERGAS	&	$15.33$	&	$6.948$	&	$2.005$	&	$0.887$	&	$0.436$	&	$0.040$	&	$0.004$	&	$0$	&	$0$	\\
PSNR	&	$17.18$	&	$23.76$	&	$34.59$	&	$41.72$	&	$48.35$	&	$69.23$	&	$89.35$	&	$109.4$	&	$297.4$	\\
UIQI	&	$0.501$	&	$0.856$	&	$0.987$	&	$0.998$	&	$0.999$	&	$1$	&	$1$	&	$1$	&	$1$	\\
\hline
\multicolumn{10}{c}{CB-STAR} \\\hline
SNR	&	$0$	&	$10$	&	$20$	&	$30$	&	$40$	&	$60$	&	$80$	&	$100$	&	$\infty$	\\\hline
SAM	&	$36.89$	&	$14.47$	&	$4.567$	&	$1.272$	&	$0.295$	&	$0.023$	&	$0.002$	&	$0$	&	$0$	\\
ERGAS	&	$52.83$	&	$17.26$	&	$5.130$	&	$1.433$	&	$0.342$	&	$0.028$	&	$0.003$	&	$0$	&	$0$	\\
PSNR	&	$7.74$	&	$17.32$	&	$27.66$	&	$38.17$	&	$50.54$	&	$71.99$	&	$91.95$	&	$111.9$	&	$265.9$	\\
UIQI	&	$0.157$	&	$0.590$	&	$0.927$	&	$0.994$	&	$1$	&	$1$	&	$1$	&	$1$	&	$1$	\\
\hline
\end{tabular}}
\label{tab:exp_synth_SNRs}
\end{table}

\begin{figure}[thb]
	\centering
    \begin{minipage}{0.36\linewidth}
        \centering
        \includegraphics[width=1\linewidth]{figures/HSI_MSI_4b_IsabellaLake_sub1.pdf}
        Lake Isabella
    \end{minipage}
    \,\,
    \begin{minipage}{0.42\linewidth}
        \centering
        \includegraphics[width=1\linewidth]{figures/HSI_MSI_lockwood_s1.pdf}
        Lockwood
    \end{minipage}
    \caption{Hyperspectral and multispectral images with a small acquisition time difference used in the experiments.}
    \label{fig:ex3_HSIs_MSIs_a}
\end{figure}
\begin{figure}[thb]
	\centering
    \begin{minipage}{0.28\linewidth}
        \centering
        \includegraphics[width=1\linewidth]{figures/HSI_MSI_Playa-crop}
        Ivanpah Playa
    \end{minipage}
    \,\,
    \begin{minipage}{0.32\linewidth}
        \centering
        \includegraphics[width=1\linewidth]{figures/HSI_MSI_Tahoe_better_reg.pdf}
        Lake Tahoe A
    \end{minipage}
    \,\,
    \begin{minipage}{0.32\linewidth}
        \centering
        \includegraphics[width=1\linewidth]{figures/HSI_MSI_Tahoe_changed.pdf}
        Lake Tahoe B
    \end{minipage}
    \caption{Hyperspectral and multispectral images with a large acquisition time difference used in the experiments.}
    \label{fig:ex3_HSIs_MSIs_b}
\end{figure}

\subsection{Example -- Real data}

In this example, we evaluated the algorithms using real HS and MS images acquired at different time instants, thus presenting different acquisition and seasonal conditions.
The reference hyperspectral and multispectral images, with a pixel size of 20~m, acquired by the AVIRIS and by the Sentinel-2A instruments, respectively, were originally considered in~\cite{Borsoi_2018_Fusion}.
Four sets of image pairs were available. Two of which contained images acquired less than three months apart (thus containing moderate variability). The other two contained images acquired with a time difference of more than one year (thus containing more significant variability).
The HSI and MSI contained $L_h=173$ and $L_m=10$ bands, respectively.
The selected ranks for the tensor-based methods are shown in Table~\ref{tab:exp_ranks_values}.
{Although the ranks of CB-STAR \textcolor{darkgreen}{satisfied Theorem~\ref{thm:exacr_rec_alg2} only for} the Ivanpah Playa image pair, \textcolor{darkgreen}{this did not have} a negative impact on its performance, as will be shown in the following.}

\subsubsection{\textcolor{magenta}{Rank} {sensitivity analysis}}
{
Before proceeding, we evaluate how sensitive the performance of the proposed methods is to the selection of the ranks by plotting the PSNR as a function of each of the ranks $K_{Z,i}$ and $K_{\Psi,i}$, while keeping the others fixed.
For simplicity, we kept the spatial ranks equal to each other (i.e., $K_{Z,1}=K_{Z,2}$ and $K_{\Psi,1}=K_{\Psi,2}$), and only show the results for the Lockwood image (to be described in Section~\ref{sec:res_realdata_moderate}) due to space limitations\footnote{Additional results are available on the supplemental material.}. The results, shown in Fig.~\ref{fig:sensitivity_ex6}, indicate that CT-STAR performs well when $K_{\Psi,1}=K_{\Psi,2}$ are small, and performs well for values of $K_{Z,3}$ which are not small. The optimal $K_{Z,1}=K_{Z,2}$ were relatively large, but a drop in performance was observed when they approach their upper limit $N_i-K_{\Psi,i}$, $i\in\{1,2,\}$. The performance of CB-STAR increased steadily with $K_{Z,1}=K_{Z,2}$ and for small values of $K_{Z,3}$, but decreased more sharply when all values $K_{Z,i}$, $i\in\{1,2,3\}$ were large. The variability ranks $K_{\Psi,i}$, on the other hand, did not affect the results too much when $K_{\Psi,1}=K_{\Psi,2}$ were sufficiently large.
}

\begin{figure}
    \centering
    \begin{minipage}{0.45\linewidth} \centering
    \includegraphics[width=\linewidth]{figures/sensit/ex6_sensit_CTSTAR_rankZ.pdf}

    \end{minipage}
    \begin{minipage}{0.45\linewidth} \centering
    \includegraphics[width=\linewidth]{figures/sensit/ex6_sensit_CTSTAR_rankPsi.pdf}

    \end{minipage}
    \begin{minipage}{0.45\linewidth} \centering
    \includegraphics[width=\linewidth]{figures/sensit/ex6_sensit_CBSTAR_rankZ.pdf}

    \end{minipage}
    \begin{minipage}{0.45\linewidth} \centering
    \includegraphics[width=\linewidth]{figures/sensit/ex6_sensit_CBSTAR_rankPsi.pdf}

    \end{minipage}
    \caption{{Sensitivity analysis of the proposed methods Lockwood image. Top: PSNR of CT-STAR as a function of the ranks of $\mathcal{Z}_h$ (left) and $\mathcal{\Psi}$ (right). Bottom: PSNR of CB-STAR as a function of the ranks of $\mathcal{Z}_h$ (left) and $\mathcal{\Psi}$ (right)}}
    \label{fig:sensitivity_ex6}
\end{figure}

\begin{table}[ht]
\renewcommand{\arraystretch}{1.15}
\setlength\tabcolsep{3.5pt}
\centering
\caption{Ranks of the tensor-based algorithms used in the experiments}  \vspace{-5pt}
{\small
    \begin{tabular}{c|cccc} \hline
    & CT-STAR: & CB-STAR: & SCOTT: & STEREO: \\
    & $K_Z,K_{\Psi}$ & $K_Z,K_{\Psi}$ & $K$ & $K$ \\ \hline
    Lockwood      & \makecell{$(30,30,8)$,\\$(3,3,2)$} & \makecell{$(70,70,5)$,\\$(40,40,3)$} & $(60,60,5)$ & $50$\\
    Lake Isabella & \makecell{$(30,30,8)$,\\$(3,3,2)$} & \makecell{$(50,50,5)$,\\$(40,40,3)$} & $(60,60,5)$ & $50$\\
    Lake Tahoe    & \makecell{$(30,30,10)$,\\$(3,3,1)$} & \makecell{$(35,35,9)$,\\$(50,50,4)$} & $(40,40,7)$ & $30$\\
    Ivanpah Playa & \makecell{$(16,16,8)$,\\$(3,3,2)$} & \makecell{$(40,40,4)$,\\$(40,40,5)$} & $(30,30,30)$ & $10$\\ \hline
\end{tabular}}
\label{tab:exp_ranks_values}
\end{table}

\subsubsection{Moderate variability}
\label{sec:res_realdata_moderate}

The first pair of images considered in this example contained $80\times80$ pixels and were acquired over the region surrounding Lake Isabella, on 2018-06-27 and on 2018-08-27. The second pair of images contained $80\times100$ pixels and was acquired near Lockwood, on 2018-08-20 and on 2018-10-19. A true color representation of the HSI and MSI for this example can be seen in Fig.~\ref{fig:ex3_HSIs_MSIs_a}.
Due to the relatively small difference between the acquisition dates of both images, the HSI and MSI look similar. However, there are slight differences between them, as seen in the overall color hue of the images and in the upper right part of the Lake Isabella HSI.
The quantitative performance metrics of all algorithms are shown in Tables~\ref{tab:ex_lockwood} and~\ref{tab:ex_isabellalake}, while the reconstructed images are presented in Figs.~\ref{fig:ex_lockwood_Z_Psi_est} and~\ref{fig:ex_isabellalake_Z_Psi_est}.

\begin{table}[ht!]
\caption{Results - Lockwood} \vspace{-5pt}
\centering
\renewcommand{\arraystretch}{1.15}
\setlength\tabcolsep{3.5pt}
{\small
\pgfplotstabletypeset[header=false,
col sep = tab,
	columns/0/.style={string type},
every head row/.style={before row={\hline}, after row=\hline},
columns/0/.style={string type, column name={Algorithm}, column type/.add={@{}@{\,}}{}},
columns/1/.style={column name={SAM}, column type/.add={@{\,}|@{\,}}{}},
columns/2/.style={column name={ERGAS}, column type/.add={@{\,}|@{\,}}{}},
columns/3/.style={column name={PSNR}, column type/.add={@{\,}|@{\,}}{}},
columns/4/.style={column name={UIQI}, column type/.add={@{\,}|@{\,}}{}},
columns/5/.style={column name={time}, fixed, precision=2, column type/.add={@{\,}|@{\,}}{}},
columns ={0,1,2,3,4,5},
every row 8 column 1/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 3 column 2/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 8 column 3/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 8 column 4/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 9 column 4/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 6 column 5/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
]{tables/tab_lockwood.txt}}
\label{tab:ex_lockwood}
\end{table}

\begin{table}[ht!]
\caption{Results - Isabella Lake} \vspace{-5pt}
\centering
\renewcommand{\arraystretch}{1.15}
\setlength\tabcolsep{3.5pt}
{\small
\pgfplotstabletypeset[header=false,
col sep = tab,
	columns/0/.style={string type},
every head row/.style={before row={\hline}, after row=\hline},
columns/0/.style={string type, column name={Algorithm}, column type/.add={@{}@{\,}}{}},
columns/1/.style={column name={SAM}, column type/.add={@{\,}|@{\,}}{}},
columns/2/.style={column name={ERGAS}, column type/.add={@{\,}|@{\,}}{}},
columns/3/.style={column name={PSNR}, column type/.add={@{\,}|@{\,}}{}},
columns/4/.style={column name={UIQI}, column type/.add={@{\,}|@{\,}}{}},
columns/5/.style={column name={time}, fixed, precision=2, column type/.add={@{\,}|@{\,}}{}},
columns ={0,1,2,3,4,5},
every row 8 column 1/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 3 column 2/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 8 column 3/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 8 column 4/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 6 column 5/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
]{tables/tab_isabellalake.txt}}
\label{tab:ex_isabellalake}
\end{table}

The quantitative results show that CB-STAR achieved the overall best results for this example, outperforming the other methods in all metrics except in ERGAS, where it performed very similarly to FuVar (which yielded the best results for this metric). CT-STAR, on the other hand, performed similarly to STEREO and SCOTT, being limited by the more stringent constraints on the image ranks.
The visual inspection of the results indicates that CB-STAR provides reconstructions closest to the ground truth when compared to the remaining methods. Although FuVar also provided good results, it yielded a slightly worse representation of the road in the left part of the Lockwood HSI, as well as more aberrations in the color of the light-brown regions in the middle of the Isabella Lake scene (which are not seen in the results of CB-STAR). {LTMR,} STEREO and SCOTT, not being able to account for variability, yielded slight color aberrations in the reconstructions, which are most clearly seen in the central part of the Isabella Lake image, while CT-STAR produced significant artifacts due to the stringent rank constraints.
The estimated factors $\mathscr{P}_{3}(\widehat{\mathcal{\Psi}})$ are in agreement with the localized changes observed in Fig.~\ref{fig:ex3_HSIs_MSIs_a}, particularly in the upper-central area of the Isabella Lake image pair, which is subject to local illumination changes.
The computation times of the algorithms show a large difference between that of FuVar and those of the other algorithms, which indicates that CB-STAR achieves better results at a significantly smaller computational complexity.

\subsubsection{Significant variability}

The remaining image pairs used in this example were acquired over the Ivanpah Playa and over Lake Tahoe area. The Ivanpah Playa image pair contained $80\times128$ pixels and was acquired on 2015-10-26 and on 2017-12-17. For the Lake Tahoe region, we considered two different image pairs (``A'' and ``B''), both with $100\times80$ pixels, the first one acquired on 2014-10-04 and on 2017-10-24, and the second one acquired on 2014-09-19 and on 2017-10-24. A true color representation of the HSI and MSI for this example can be seen in Fig.~\ref{fig:ex3_HSIs_MSIs_b}.
Due to the considerable difference between the acquisition date/time of the HSI and MSI, significant differences can be found between them. For the Ivanpah Playa images, there are large variations between the sand colors in the central part of the image. For the Lake Tahoe region, significant differences are observed in both image pairs, with differences in the color hues of the ground and of the crop circles for the image pair A, and also a large change in the water level of the lake in the image pair B.
The quantitative performance metrics of all algorithms are shown in Tables~\ref{tab:ex_playa},~\ref{tab:ex_tahoe}, and~\ref{tab:ex_tahoeChanged}, while the reconstructed images are presented in Figs.~\ref{fig:ex_playa_Z_Psi_est},~\ref{fig:ex_tahoe_Z_Psi_est} and~\ref{fig:ex_tahoeChanged_Z_Psi_est}.

\begin{table}[ht!]
\caption{Results - Ivanpah Playa} \vspace{-5pt}
\centering
\renewcommand{\arraystretch}{1.15}
\setlength\tabcolsep{3.5pt}
{\small
\pgfplotstabletypeset[header=false,
col sep = tab,
	columns/0/.style={string type},
every head row/.style={before row={\hline}, after row=\hline},
columns/0/.style={string type, column name={Algorithm}, column type/.add={@{}@{\,}}{}},
columns/1/.style={column name={SAM}, column type/.add={@{\,}|@{\,}}{}},
columns/2/.style={column name={ERGAS}, column type/.add={@{\,}|@{\,}}{}},
columns/3/.style={column name={PSNR}, column type/.add={@{\,}|@{\,}}{}},
columns/4/.style={column name={UIQI}, column type/.add={@{\,}|@{\,}}{}},
columns/5/.style={column name={time}, fixed, precision=2, column type/.add={@{\,}|@{\,}}{}},
columns ={0,1,2,3,4,5},
every row 3 column 1/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 8 column 2/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 8 column 3/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 3 column 4/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 7 column 5/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
]{tables/tab_playa.txt}}
\label{tab:ex_playa}
\end{table}

\begin{table}[ht!]
\caption{Results - Lake Tahoe A} \vspace{-5pt}
\centering
\renewcommand{\arraystretch}{1.15}
\setlength\tabcolsep{3.5pt}
{\small
\pgfplotstabletypeset[header=false,
col sep = tab,
	columns/0/.style={string type},
every head row/.style={before row={\hline}, after row=\hline},
columns/0/.style={string type, column name={Algorithm}, column type/.add={@{}@{\,}}{}},
columns/1/.style={column name={SAM}, column type/.add={@{\,}|@{\,}}{}},
columns/2/.style={column name={ERGAS}, column type/.add={@{\,}|@{\,}}{}},
columns/3/.style={column name={PSNR}, column type/.add={@{\,}|@{\,}}{}},
columns/4/.style={column name={UIQI}, column type/.add={@{\,}|@{\,}}{}},
columns/5/.style={column name={time}, fixed, precision=2, column type/.add={@{\,}|@{\,}}{}},
columns ={0,1,2,3,4,5},
every row 3 column 1/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 8 column 2/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 8 column 3/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 8 column 4/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 9 column 4/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 5 column 5/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
]{tables/tab_tahoe.txt}}
\label{tab:ex_tahoe}
\end{table}

\begin{table}[ht!]
\caption{Results - Lake Tahoe B} \vspace{-5pt}
\centering
\renewcommand{\arraystretch}{1.15}
\setlength\tabcolsep{3.5pt}
{\small
\pgfplotstabletypeset[header=false,
col sep = tab,
	columns/0/.style={string type},
every head row/.style={before row={\hline}, after row=\hline},
columns/0/.style={string type, column name={Algorithm}, column type/.add={@{}@{\,}}{}},
columns/1/.style={column name={SAM}, column type/.add={@{\,}|@{\,}}{}},
columns/2/.style={column name={ERGAS}, column type/.add={@{\,}|@{\,}}{}},
columns/3/.style={column name={PSNR}, column type/.add={@{\,}|@{\,}}{}},
columns/4/.style={column name={UIQI}, column type/.add={@{\,}|@{\,}}{}},
columns/5/.style={column name={time}, fixed, precision=2, column type/.add={@{\,}|@{\,}}{}},
columns ={0,1,2,3,4,5},
every row 8 column 1/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 8 column 2/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 8 column 3/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 8 column 4/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
every row 5 column 5/.style={postproc cell content/.style={ @cell content/.add={$\bf}{$} }},
]{tables/tab_tahoeChanged.txt}}
\label{tab:ex_tahoeChanged}
\end{table}

The quantitative results show that \mbox{CB-STAR} achieved again the overall best results for this example, outperforming the remaining algorithms in most metrics, except in the SAM and UIQI for the Ivanpah Playa HRI and in the SAM of the Lake Tahoe A HRI. Moreover, there was a stronger gap between the performance of the methods that consider variability and the remaining algorithms. \mbox{CT-STAR}, although better than {LTMR,} STEREO and SCOTT, performed significantly worse than CB-STAR due to its stringent constraints on the image ranks.
The visual inspection of the results again indicates that CB-STAR provides reconstructions closest to the ground truth when compared to the remaining methods. Although FuVar also provided good results {(as it accounts for spectral variability)}, the reconstructions by CB-STAR were closer to the ground truth, as can be observed in the color shades of the upper part of the Ivanpah Playa image and of the crop circles of the Lake Tahoe~A image, and especially in the overall colors in the more uniform regions containing soil and water and vegetation in the Lake Tahoe B image. {However, FuVar results showed slightly sharper edges in some regions (e.g., around the crop circles in the Lake Tahoe images), which occurs due to the use of a Total Variation spatial regularization. Nonetheless, a spatial regularization can  also be incorporated to the CB-STAR cost function in~\eqref{eq:opt_prob_fus_btd} to achieve a similar effect.}
{LTMR,} STEREO and SCOTT, on the other hand, produced significant artifacts in all reconstructed images. CT-STAR also produced significant artifacts, which as in the previous example are due to the stringent rank constraints{, with reconstructions visually worse than those of CB-STAR and FuVar}.
The estimated factors $\mathscr{P}_{3}(\widehat{\mathcal{\Psi}})$ were in close agreement with the variability seen in the scenes, notably in the sand region of the central part of the Ivanpah Playa image (which lies at the botton of a hill), and in the regions near the lake in the Lake Tahoe A and B images, which undergo variations in the water level. Moreover, the overall amplitude of the variables was significantly larger than in the previous example, in which the differences between the images were more moderate.
The computation times of all algorithms were similar to those observed in the previous example, except for that of CB-STAR, which was higher since it underwent a larger number of iterations for the data in this example. Nonetheless, the computation times of CB-STAR were still considerably smaller than those of FuVar.

\begin{figure*}
    \centering
    \resizebox{0.85\linewidth}{!}
    {

    \includegraphics[height=10ex]{figures/Zhat/ex6_lockwood_recIm_vis_IR.pdf}
    \includegraphics[height=9.85ex]{figures/estimPsis/ex6_lockwood_both_algs_P3Psi_vis_norms2.pdf}
    }
    \vspace{-0.3cm}
    \caption{Left: Visible (top) and infrared (bottom) representation of the true and estimated versions of the Lockwood HSI. Right: Spectrally degraded additive scaling factors $\mathscr{P}_{3}(\mathcal{\Psi})$ estimated by CT-STAR and CB-STAR.}
    \label{fig:ex_lockwood_Z_Psi_est}
\end{figure*}

\begin{figure*}
    \centering
    \resizebox{0.85\linewidth}{!}
    {
    \includegraphics[height=10ex]{figures/Zhat/ex5_isabellalake_recIm_vis_IR.pdf}
    \includegraphics[height=9.9ex]{figures/estimPsis/ex5_isabellalake_both_algs_P3Psi_vis_norms2.pdf}
    }
    \vspace{-0.3cm}
    \caption{Left: Visible (top) and infrared (bottom) representation of the true and estimated versions of the Isabella Lake HSI. Right: Spectrally degraded additive scaling factors $\mathscr{P}_{3}(\mathcal{\Psi})$ estimated by CT-STAR and CB-STAR.}
    \label{fig:ex_isabellalake_Z_Psi_est}
\end{figure*}

\begin{figure*}
    \centering
    \resizebox{0.85\linewidth}{!}
    {
    \includegraphics[height=10ex]{figures/Zhat/ex3_playa_recIm_vis_IR.pdf}
    \includegraphics[height=9.95ex]{figures/estimPsis/ex3_playa_both_algs_P3Psi_vis_norms2.pdf}
    }
    \vspace{-0.3cm}
    \caption{Left: Visible (top) and infrared (bottom) representation of the true and estimated versions of the Ivanpah Playa HSI. Right: Spectrally degraded additive scaling factors $\mathscr{P}_{3}(\mathcal{\Psi})$ estimated by CT-STAR and CB-STAR.}
    \label{fig:ex_playa_Z_Psi_est}
\end{figure*}

\begin{figure*}
    \centering
    \resizebox{0.85\linewidth}{!}
    {
    \includegraphics[height=10ex]{figures/Zhat/ex4_tahoeA_recIm_vis_IR.pdf}
    \includegraphics[height=9.9ex]{figures/estimPsis/ex4_tahoe_both_algs_P3Psi_vis_norms2.pdf}
    }
    \vspace{-0.3cm}
    \caption{Left: Visible (top) and infrared (bottom) representation of the true and estimated versions of the Lake Tahoe A HSI. Right: Spectrally degraded additive scaling factors $\mathscr{P}_{3}(\mathcal{\Psi})$ estimated by CT-STAR and CB-STAR.}
    \label{fig:ex_tahoe_Z_Psi_est}
\end{figure*}

\begin{figure*}[thb]
    \centering
    \resizebox{0.85\linewidth}{!}
    {
    \includegraphics[height=10ex]{figures/Zhat/ex4_tahoeB_recIm_vis_IR.pdf}
    \includegraphics[height=9.9ex]{figures/estimPsis/ex4_tahoe_changed_both_algs_P3Psi_vis_norms2.pdf}
    }
    \vspace{-0.3cm}
    \caption{Left: Visible (top) and infrared (bottom) representation of the true and estimated versions of the Lake Tahoe B HSI. Right: Spectrally degraded additive scaling factors $\mathscr{P}_{3}(\mathcal{\Psi})$ estimated by CT-STAR and CB-STAR.}
    \label{fig:ex_tahoeChanged_Z_Psi_est}
\end{figure*}

\section{Conclusions}
\label{sec:conclusions}

In this paper, we proposed a novel framework for multimodal (hyperspectral and multispectral) image fusion accounting for spatially and spectrally localized changes. We first studied the general identifiability of the considered model, which becomes more ambiguous due to the presence of changes. Then, assuming that the high resolution image and the variation factors admit a Tucker decomposition, two new algorithms were proposed -- one being purely algebraic (which was computationally more efficient), and another based on an optimization procedure (which allowed for more relaxed specification of the multilinear ranks). Theoretical guarantees for the exact recovery of the high resolution image were provided for both algorithms. The proposed optimization-based algorithm achieved superior experimental performance in the presence of spectral and spatial variations between the images, while also exhibiting a smaller computational cost.

\appendices

\section{Optimization of cost function~\eqref{eq:opt_subp_Z}}
\label{sec:appendix1}

{To solve~\eqref{eq:opt_subp_Z}, we consider a block coordinate descent strategy (as in, e.g.,~\cite{lathauwer2008tensor_BTD3_algorithms}), which successively minimizes the cost function in~\eqref{eq:opt_subp_Z} with respect to each of the variables $\mathcal{G}_Z$, $\boldsymbol{B}_{Z,i}$, $i\in\{1,2,3\}$,} while keeping the remaining ones fixed.
Let us denote by $J\big(\boldsymbol{X}\,|\,\boldsymbol{\Theta}_{\setminus\{\boldsymbol{X}\}}\big)$ the cost function $J$ in which all variables but $\boldsymbol{X}$ are fixed.
Note that we apply the QR factorization after computing each of the factor matrices $\boldsymbol{B}_{\iota,j}$ to constrain them to be unitary at all the iterations,  as performed in~\cite{lathauwer2008tensor_BTD3_algorithms}. This normalization prevents convergence issues by avoiding under/over-flow and keeping these matrices well-conditioned. {Starting from an initialization $\mathcal{G}_Z$, $\boldsymbol{B}_{Z,i}$, $i\in\{1,2,3\}$, this procedure is repeated for $F$ inner iterations.}

\begin{algorithm} [thb]
\footnotesize
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\caption{\mbox{Block coordinate descent to solve~\eqref{eq:opt_subp_Z}}
\label{alg:alg3}}
\Input{\mbox{Tensors $\mathcal{Y}_h$, $\mathcal{X}_0$, $\mathcal{G}_Z$ matrices $\boldsymbol{B}_{Z,i}$,$i\in\{1,2,3\}${, iterations~$F$}}}
\Output{$\mathcal{G}_Z$ matrices $\boldsymbol{B}_{Z,i}$,$i\in\{1,2,3\}$}

\For{$i=1,\ldots,F$}{
Minimize $J(\mathcal{G}_Z|\boldsymbol{\Theta}_{\setminus\{\mathcal{G}_Z\}}^{(i)})$ w.r.t. $\mathcal{G}_Z$\;
Minimize $J(\boldsymbol{B}_{Z,1}|\boldsymbol{\Theta}_{\setminus\{\boldsymbol{B}_{Z,1}\}}^{(i)})$ w.r.t. $\boldsymbol{B}_{Z,1}$\;
Minimize $J(\boldsymbol{B}_{Z,2}|\boldsymbol{\Theta}_{\setminus\{\boldsymbol{B}_{Z,2}\}}^{(i)})$ w.r.t. $\boldsymbol{B}_{Z,2}$\;
Minimize $J(\boldsymbol{B}_{Z,3}|\boldsymbol{\Theta}_{\setminus\{\boldsymbol{B}_{Z,3}\}}^{(i)})$ w.r.t. $\boldsymbol{B}_{Z,3}$\;}
\end{algorithm}

\subsubsection{Optimizing w.r.t. $\boldsymbol{B}_{Z,i}$, $i\in\{1,2,3\}$}

To save space, we present only the case where $i=1$. The extension to $i\in\{2,3\}$ is straightforward.
Note that the cost function $J\big(\boldsymbol{B}_{Z,1}\,|\,\boldsymbol{\Theta}_{\setminus\{\boldsymbol{B}_{Z,1}\}}\big)$ can be equivalently reformulated using the mode-$1$ matricization as
\begin{align} \label{eq:opt_subp_1_i}
    J\big(& \boldsymbol{B}_{Z,1}\,|\,\boldsymbol{\Theta}_{\setminus\{\boldsymbol{B}_{Z,1}\}}\big) = \big\|\unfold{\boldsymbol{Y}_{\!h}}{1}^\top - \boldsymbol{X}_2 \boldsymbol{B}_{Z,1}^\top\boldsymbol{P}_1^\top\big\|_F^2
    \nonumber\\ &
    + {\lambda} \big\|\unfold{\boldsymbol{Y}_{\!m}}{1}^\top - \boldsymbol{X}_3 - \boldsymbol{X}_1 \boldsymbol{B}_{Z,1}^\top\big\|_F^2 \,,
\end{align}
where matrices $\boldsymbol{X}_1$, $\boldsymbol{X}_2$ and $\boldsymbol{X}_3$ are given by
\begin{align}
    \boldsymbol{X}_1 &= \big({\boldsymbol{B}}_{Z,3} \otimes \boldsymbol{P}_3{\boldsymbol{B}}_{Z,2}\big) \unfold{\boldsymbol{G}_Z}{1}^\top \,,
    \\
    \boldsymbol{X}_2 &= \big(\boldsymbol{P}_2{\boldsymbol{B}}_{Z,3} \otimes {\boldsymbol{B}}_{Z,2}\big) \unfold{\boldsymbol{G}_Z}{1}^\top \,,
    \\
    \boldsymbol{X}_3 & = \big( \boldsymbol{P}_3 {\boldsymbol{B}}_{\Psi,3} \otimes {\boldsymbol{B}}_{\Psi,2}\big)\unfold{\boldsymbol{G}_{\Psi}}{1}^\top {\boldsymbol{B}}_{\Psi,1}^\top \,.
\end{align}

Computing the derivative of~\eqref{eq:opt_subp_1_i} and setting it equal to zero results in the following expression:
\begin{align}
\begin{split}
    {\lambda}\boldsymbol{X}_1^\top\boldsymbol{X}_1 {\boldsymbol{B}}_{Z,1}^\top
    & + \boldsymbol{X}_2^\top\boldsymbol{X}_2 {\boldsymbol{B}}_{Z,1}^\top\boldsymbol{P}_1^\top\boldsymbol{P}_1
    \\
    = {\lambda} & \boldsymbol{X}_1^\top(\unfold{\boldsymbol{Y}_{\!m}}{1}^\top - \boldsymbol{X}_3)
    + \boldsymbol{X}_2^\top\unfold{\boldsymbol{Y}_{\!h}}{1}^\top\boldsymbol{P}_1 \,.
\end{split}
\end{align}
This is a Sylvester equation that can be directly solved using existing software with, e.g., the Hessenberg-Schur or the Bartels-Stewart algorithms (see~\cite{simoncini2016solversSylvesterEq} and references therein).

\subsubsection{Optimizing w.r.t. $\mathcal{G}_{Z}$}

Cost function $J\big(\mathcal{G}_{Z}\,|\,\boldsymbol{\Theta}_{\setminus\{\mathcal{G}_{Z}\}}\big)$ can be equivalently reformulated using the tensor vectorization as
\begin{align} \label{eq:opt_subp_4_i}
    J\big(\mathcal{G}_{Z}\,|\,\boldsymbol{\Theta}_{\setminus\{\mathcal{G}_{Z}\}}\big)
    ={} & \big\|\operatorname{vec}(\mathcal{Y}_h) - \boldsymbol{X}_2 \,\boldsymbol{g}_Z \,\big\|_F^2
    \nonumber \\ &
    + {\lambda} \big\|\operatorname{vec}(\mathcal{Y}_m) - \boldsymbol{x}_3 - \boldsymbol{X}_1 \, \boldsymbol{g}_Z \, \big\|_F^2 \,,
\end{align}
where $\boldsymbol{g}_Z=\operatorname{vec}(\mathcal{G}_{Z})$ is the vectorization of the core tensor, and $\boldsymbol{X}_1$, $\boldsymbol{X}_2$ and $\boldsymbol{x}_3$ are given by
\begin{align}
    \boldsymbol{X}_1 & = \big(\boldsymbol{P}_3 \boldsymbol{B}_{Z,3} \otimes \boldsymbol{B}_{Z,2} \otimes \boldsymbol{B}_{Z,1} \big) \,,
    \\
    \boldsymbol{X}_2 & = \big(\boldsymbol{B}_{Z,3} \otimes \boldsymbol{P}_2 \boldsymbol{B}_{Z,2} \otimes \boldsymbol{P}_1 \boldsymbol{B}_{Z,1} \big) \,,
    \\
    \boldsymbol{x}_3 & = \big(\boldsymbol{P}_3\boldsymbol{B}_{\Psi,3} \otimes \boldsymbol{B}_{\Psi,2} \otimes \boldsymbol{B}_{\Psi,1} \big) \operatorname{vec}(\mathcal{G}_{\Psi}) \,.
\end{align}
The solution that minimizes~\eqref{eq:opt_subp_4_i} can be computed through the normal equations, which can be written as
\begin{align}
    ({\lambda}\boldsymbol{X}_1^\top\boldsymbol{X}_1 + \boldsymbol{X}_2^\top\boldsymbol{X}_2) \boldsymbol{g}_Z {}={} &  {\lambda}\boldsymbol{X}_1^\top(\operatorname{vec}(\mathcal{Y}_m)-\boldsymbol{x}_3)
    \nonumber \\
    & + \boldsymbol{X}_2^\top\operatorname{vec}(\mathcal{Y}_h) \,.
\end{align}
As shown in~\cite{prevost2020coupledTucker_hyperspectralSRR_TSP}, this set of equations can be alternatively interpreted as a generalized Sylvester equation, for which efficient solvers can be used.

\bibliographystyle{IEEEtran}