# Best of both worlds: Fusing hyperspectral data from two generations of spectro-imagers for X-ray astrophysics

Julia Lascar, Jérôme Bobin, Fabio Acero

► To cite this version:

Julia Lascar, Jérôme Bobin, Fabio Acero. Best of both worlds: Fusing hyperspectral data from two generations of spectro-imagers for X-ray astrophysics. Astronomy & Astrophysics - A&A, 2025, 694, pp.A34. 10.1051/0004-6361/202451608. hal-04667869

HAL Id: hal-04667869

https://hal.science/hal-04667869v1

Submitted on 14 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Best of both worlds: Fusing hyperspectral data from two generations of spectro-imagers for X-ray astrophysics

J. Lascar<fnref n="1" />,*, J. Bobin<fnref n="1" />, and F. Acero<fnref n="1" />,<fnref n="2" /> ID

Received 22 July 2024 / Accepted 25 November 2024

## ABSTRACT

*Context.* With the recent launch of the X-Ray Imaging and Spectroscopy Mission (XRISM) and the advent of microcalorimeter detectors, X-ray astrophysics is entering in a new era of spatially resolved high-resolution spectroscopy. But while this new generation of X-ray telescopes offers much finer spectral resolution than the previous one (e.g. XMM-Newton, Chandra), these instruments also have coarser spatial resolutions, leading to problematic cross-pixel contamination. This issue is currently a critical limitation for the study of extended sources, such as galaxy clusters or supernova remnants (SNRs).

*Aims.* To increase the scientific output of XRISM's hyperspectral data, we propose that it be fused with XMM-Newton data, and seek to obtain a cube with the best spatial and spectral resolution of both generations. This is the aim of hyperspectral fusion. In this article, we explore the potential and limitations of hyperspectral fusion for X-ray astrophysics.

*Methods.* We implemented an algorithm that jointly deconvolves the spatial response of XRISM and the spectral response of XMM-Newton. To do so, we construct a forward model adapted for instrumental systematic degradations and Poisson noise, and then tackle hyperspectral fusion as a regularized inverse problem. We test three methods of regularization: spectral low-rank approximation with a spatial Sobolev regularization; spectral low-rank approximation with a 2D wavelet sparsity constraint; and a 2D-1D wavelet sparsity constraint.

*Results.* We test our method on toy models constructed from hydrodynamic simulations of SNRs. We find that our method reconstructs the ground truth well even when the toy model is complex. For the regularization term, we find that while the low-rank approximation works well as a spectral denoiser in models with less spectral variability, it introduces a bias in models with more spectral variability, in which case the 2D-1D wavelet sparsity regularization works best. Following the present proof of concept, we aim to apply this method to real X-ray astrophysical data in the near future.

**Key words.** methods: data analysis – methods: statistical – ISM: supernova remnants – X-rays: general

## 1. Introduction

With the advent of microcalorimeter detectors, X-ray astrophysics has entered a new era of spatially resolved high-resolution spectroscopy. Indeed, the recent launch of the X-Ray Imaging and Spectroscopy Mission (XRISM) and its microcalorimeter array “Resolve”<fnref n="1" />, and the future Athena telescope and its X-ray Integral Field Unit (X-IFU)<fnref n="2" /> camera, has granted access to spatially resolved spectroscopy at spectral resolutions improved by an order of magnitude. This is set to allow astrophysicists to investigate previously indistinguishable spectral features, such as the emission lines of rare elements, line broadening, or velocity measurements. These characteristics are essential in the study of extended sources such as supernova remnants (SNRs), galaxy clusters, and many of the other exotic sources of the X-ray Universe.

However, while the XRISM/Resolve instrument has a much finer spectral resolution than the Chandra X-ray Observatory's Advanced CCD Imaging Spectrometer (Chandra/ACIS)<fnref n="3" /> or the X-ray Multi-Mirror-Newton's European Photon Imaging Camera (XMM-Newton/EPIC)<fnref n="4" /> cameras, its spatial resolution is

much lower than that of its predecessors. As we demonstrate here, this spatial blur is an important obstacle when analyzing XRISM/Resolve data. Given this trade-off, we set out to establish whether or not it is possible to combine the data from these two generations of X-ray observatories and obtain a result with the best spectral and spatial resolution of each instrument. This is the idea at the heart of hyperspectral fusion.

In this paper, we explore the potential and limitations of hyperspectral fusion for X-ray data taken by pairs of telescopes such as XMM-Newton/EPIC with XRISM/Resolve, or Chandra with Athena-X-IFU.

We implemented a modular optimization algorithm that can use one of three potential regularizations: first, inspired by Guilloteau et al. (2020), and given the widespread usage of low-rank approximations in the literature, we implement a low-rank approximation with Sobolev spatial regularization. Since Sobolev regularization does not take into account multi-scale features, we then considered a low-rank approximation with a 2D wavelet sparsity regularization. Finally, to investigate whether low-rankedness is indeed appropriate, we implemented a heuristic approach using a 2D-1D Wavelet sparsity regularization. This hyperspectral fusion code is openly available on Github<fnref n="5" />.

In Section 2, we present the astrophysical and methodological context behind our work. In Section 3, we present our

<footnote n="1">¹ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France</footnote>

<footnote n="2">² FSLAC IRL 2009, CNRS/IAC, La Laguna, Tenerife, Spain</footnote>

<footnote>* Corresponding author; julia.lascar@cea.fr</footnote>

<footnote n="1">¹ https://xrism.isas.jaxa.jp</footnote>

<footnote n="2">² https://x-ifu.irap.omp.eu/</footnote>

<footnote n="3">³ https://cxc.harvard.edu/</footnote>

<footnote n="4">⁴ https://www.cosmos.esa.int/web/xmm-newton</footnote>

<footnote n="5">⁵ https://github.com/JMLascar/HIFReD_fusion</footnote>
methodology, and describe the forward model we have devel-
oped for X-ray astrophysics, and the proposed regularization
schemes. In Section 4, we test our method on toy models
of varying complexity using synthetic X-ray data cubes con-
structed from hydrodynamic simulation of SNRs. We discuss
the achievements, limitations, and perspective for future work
in Section 5.

## 2. Context

### 2.1. Astrophysical context: X-ray spectro-imaging and microcalorimeters

X-ray spectro-imagers can collect hyperspectral images. These
are 3D data arrays: they have two spatial dimensions, like a
classical image, and a third spectral dimension (energy or wave-
length). In other words, an energy spectrum is associated to each
spatial pixel.

In recent decades, charge-coupled device (CCD) cameras
have been the main spectro-imaging instruments on board X-ray
telescopes such as XMM-Newton, Chandra, Swift, Suzaku, and
eROSITA. When the X-ray photon interacts with the semicon-
ductor material, an electric charge is generated proportional to
the energy of the incident photon. The charge is then transferred
row by row to a readout register. For example, the XMM-Newton
MOS camera can accommodate a large number (600×600) of
small pixels (40 µm, Turner et al. 2001). The energy resolution
is of the order 150 eV<fnref n="6" /> at 6 keV (E/ΔE ~ 40), which is close
to the Fano limit, the physical lower bound of energy resolution
due to the statistical fluctuations in the number of electron-hole
pairs produced in the photon-material interaction.

To overcome this limit, a change of technology is required. Microcalorimeter devices are providing this revolutionary leap forward for X-ray spectro-imaging instruments. The fundamental principle of microcalorimeters is to measure the photon energy via the small increase in temperature caused by the conversion of this energy into heat in the sensor. To be sensitive to small heat increments and reach a spectral resolution on the order of a few electron volts, the device has to be cooled down to 50 mK and remain at very stable temperature. It is important to note that, unlike gratings on board XMM-Newton or Chandra, microcalorimeters are non-dispersive spectrometers and can be used to study diffuse structure.

However, with the current available technology, only a lim-
ited number of pixels can be accommodated on the focal plane.
The XRISM/Resolve instrument is a microcalorimeter array of
6 × 6 sensors (pixel size of 800 µm, Takahashi et al. 2014)
with a constant spectral resolution over the entire bandpass of
5 eV FWHM (XRISM Science Team 2020). While future X-ray
telescopes combining microcalorimeter and high spatial reso-
lution are planned (Athena/XIFU, see Barret et al. 2018, for
more details), they are expected to launch in the late 2030s, and
will still not reach the spatial resolution of the X-ray telescope
Chandra.

In the meantime, mitigation plans to overcome this spatial-
resolution limitation should be investigated to try to partially
recover some of the spatial information. To that end, we explore
in this study the performance of hyperspectral fusion to combine
data sets with high spectral resolution and high spatial resolu-
tion using the example of fusion between XMM-Newton and
XRISM/Resolve data. The final goal of this reconstruction is

to enhance the scientific return of XRISM/Resolve by retriev-
ing maps of physical parameters (plasma temperature, veloc-
ity, abundances, etc.) at higher resolution with reduced cross-
contamination between pixels (spectrospatial mixing).

Indeed, this mixing is an important problem when analyz-
ing hyperspectral data. If an instrument's spatial resolution is
larger than the coherence scale of a physical parameter, the
reconstructed information will be an average of the small-scale
features. To take the example of a velocity map from the ejecta of
a SNR, the small-scale redshifted and blueshifted ejecta features
can mistakenly be measured at small or close-to-zero redshift
values, as was demonstrated by the detailed analysis of Tycho's
SNR in Godinaud et al. (2023). By jointly deconvolving the
data sets from two X-ray instruments, hyperspectral data fusion
could solve this problem and recover higher-fidelity scientific
information at small scale.

2.2. Hyperspectral fusion

Given two hyperspectral (HS) images, the general task of
hyperspectral fusion consists in obtaining a cube with the
best spatial and spectral resolution of each data set. For
instance, Fig. 1 shows HS images of the Perseus galaxy cluster
taken by the X-ray telescopes XMM-Newton and Hitomi/SXS<fnref n="7" />.
XMM-Newton has better spatial resolution, and Hitomi has bet-
ter spectral resolution. The goal of fusing the two datasets would
be to obtain a HS image with XMM-Newton's spatial resolution
and Hitomi's spectral resolution.

This problem has been explored in the literature for applica-
tions to satellite pictures of the earth, and has recently received
some attention in astrophysics with the launch of the James
Webb Space Telescope (JWST). However, X-ray astrophysics
has its own specific challenges. X-ray data are dominated by
Poisson noise, and both instruments produce HS images with
many spectral channels, as opposed to one of the data sets being a
multi-spectral image with a small number of spectral bands typ-
ically obtained with different filters. To use hyperspectral fusion
for X-ray astrophysics, it is thus important to find the fusion
method most appropriate for this case study.

In this section, we provide an overview of the various aspects
to take into account when approaching a problem of hyperspec-
tral fusion in order to find the best way to solve the problem. For
a detailed review of the literature, we refer the interested reader
to Yokoya et al. (2017).

The first aspect to consider is the forward model of the two
instruments. This defines how a source image will be degraded
in the process of image formation and recording, which depends
on the systematic characteristics of the instrument, and the noise
we expect to measure.

When considering systematic instrumental degradation, two
dimensions may be taken into account: the spatial degradation,
and the spectral degradation. Spatially, the degradation will be
due to optical diffraction or other systematic blurs, and this infor-
mation is encapsulated in the point spread function (PSF). Given
a point source, the PSF returns its radiance distribution on the
recorded image. Similarly, if a source emits a Dirac spectrum,
the measured spectrum will be the instrument's spectral response
function, and this will be affected by, for instance, the chemical
and physical properties of the spectro-imager. Depending on the
instrument, the PSF and spectral response may depend on the
energy, on the spatial location of the pixel on the camera, or both.

<footnote n="6">6 https://xmm-tools.cosmos.esa.int/external/xmm_user_support/documentation/uhb/epic_specres.html</footnote>

<footnote n="7">⁷ XRISM/Resolve camera is identical to the SXS on board the short-lived Hitomi telescope launched in 2016.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Comparison of hyperspectral data of the Perseus galaxy cluster as observed by two generations of X-ray spectro-imagers. The XMM-Newton/MOS (top) and Hitomi/SXS (bottom) data have pixels of 8"/30" and spectral channels of 15/1 eV respectively. The left panels show the counts maps summed over all spectral channels and the right panels the spectra summed over all pixels.</figcaption></figure>

Then, for the stochastic noise, while most of the literature on fusion assumes additive Gaussian noise, it is important to realize this is not always appropriate. In X-ray astrophysics, photons have a low flux and are measured individually, and thus follow a Poisson distribution.

Once the instrumental systematic degradation and the stochastic noise have been taken into account, if we take two HS images X and Y of sizes ($M \times N \times l$) and ($m \times n \times L$), the general problem of fusing two HS images can be written as:

$$ \hat{Z} = \arg\min_{Z} D(M_{X}(Z)|X) + D(M_{Y}(Z)|Y), \qquad (1) $$

where $\hat{Z}$ is the fusion result, $M_X$ and $M_Y$ are the model operators that take into account the systematic instrumental degradation (PSF, spectral response, etc.) and the rebinning needed to take account different voxel sizes, and D is a well-chosen function that evaluates the difference between the model $M_X(Z)$ and data X, and between model $M_Y(Z)$ and data Y, given the expected noise distribution. Thus, this problem ensures that Z is a plausible source from which the first instrument can measure X and the second can measure Y. In the case of X-ray astrophysics, the forward model is the one described in Section 3.1 and Fig. 2.

However, this problem is often ill-posed (in the sense of Hadamard 1923) and ill-conditioned, meaning that while there might be a unique solution, the noise will be amplified when inverting the model, corrupting the solution, and slight changes in the data would lead to drastic changes in the solution (Bertero et al. 2022). One well-studied approach to solving ill-posed inverse problems is to regularize the problem, which restricts the space of acceptable solutions. This can be approached by considering the spatial features of Z, its spectral features, or both. To regularize in the spectral dimension,

one classic method is to project the spectra of Z onto a lower-dimensional subspace, as is done with source separation (e.g., Yokoya et al. 2012; Prévost et al. 2022; Pineau et al. 2023) or low-rank approximation (e.g., Simões et al. 2015; Xu et al. 2022; Guilloteau et al. 2020). We can also add a regularization term $\varphi(Z)$ to the function we seek to minimize, which leads us to:

$$ \hat{Z} = \arg\min_{Z} D(M_{X}(Z)|X) + D(M_{Y}(Z)|Y) + \varphi(Z). \qquad (2) $$

This $\varphi(Z)$ will promote certain spatial or spectral structures in the solution. Its choice must thus be appropriate for our data. For instance, in remote sensing (i.e., the study of satellite pictures of the earth), it may be useful to promote spatial piece-wise smoothness using a total variation term, as done by Simões et al. (2015), but this would be ill-advised for astrophysical data, where such features are rarely observed. For our case of X-ray astrophysics, we consider three potential regularization approaches, which are described in Section 3.3.

Finally, to solve the problem in Equation (2), we need to choose an appropriate optimization scheme. In some cases, there might be an analytical solution (as in the work of Pineau et al. 2023 for JWST, when the spectral content is assumed to be a linear combination of known spectra), but most often an iterative algorithm must be constructed. The choice of optimizer will depend on the data size, the specific challenges presented by the forward model, and the chosen regularization terms. It may also be the case that constraints placed on the optimizer (such as insisting that the algorithm reaches convergence as quickly as possible) will lead to constraints on the choice of regularization and cost function.

In astrophysics, work has been done to tackle multispectral/HS fusion for the JWST (see Guilloteau et al. 2022; Pineau et al. 2023) in order to fuse images taken by two of
<figure><img src="image_3.png" /><figcaption>Fig. 2. Simplified diagram of the forward model of two X-ray telescopes (here, XMM-Newton/EPIC and XRISM/Resolve). The problem of fusion consists in inverting this model and using the observed data X and Y to retrieve the deconvoluted cube Z. The mathematical details behind the operators and kernels are detailed in Section 3. Not shown is the exposure time of each data set.</figcaption></figure>

its instruments (NIRCam and NIRSpec), but to our knowledge, no algorithm has been implemented for HS/HS fusion with Poisson statistics, considering two instruments with their PSF and spectral responses, as is needed in X-ray astrophysics.

### 3. Methodology

In this section, we present the methodology of our fusion algorithm for X-ray astrophysics. We present the general problem to be solved in Section 3.1. Section 3.2 presents the rebinning operator needed when the data sets have different voxel sizes. Section 3.3 describes our choices of potential regularization terms. Section 3.4 presents our algorithmic strategy to solve the inverse problem.

#### 3.1. Forward model

Let X be the hyperspectral data from the instrument with higher spatial resolution but lower spectral resolution (e.g., XMM-Newton/EPIC). Let Y be the data from the other instrument, with lower spatial resolution and higher spectral resolution (e.g., XRISM/Resolve). When available, X̄ and Ȳ are the ground truths at the resolutions of the respective instruments.

In this section, we consider the problem where both HS images are of the same size, that is, X and Y ∈ ℝ<sup>M×N×L</sup>, and therefore so does Z (see Section 3.2 for the case where X and Y have different sizes). Our forward model (seen in Fig. 2) must include the spatial degradation of each instrument, which is modeled by their PSFs, as well as the spectral degradation, which is modeled by their spectral responses. We must also account for the effective area of the instrument, and the exposure time. For the former, we assume an effective area with no vignetting; that is, the effective area does not vary spatially<sup>8</sup>. We use A<sub>X</sub> and A<sub>Y</sub> to denote the effective area cubes, which have the same shape as

X and Y. t<sub>X</sub> and t<sub>Y</sub> are scalars, the exposure times of each data set. The spectral responses of the instruments are s<sub>X</sub> and s<sub>Y</sub>, and their PSFs are B<sub>X</sub> and B<sub>Y</sub>.

We define B<sub>X</sub> and B<sub>Y</sub> as the operators which, given a cube Z, apply a 2D convolution on its spatial slices with the PSF B<sub>X</sub> or B<sub>Y</sub>. Similarly, we use S<sub>X</sub> and S<sub>Y</sub> to denote the operators that convolve each spectra of Z with s<sub>X</sub> or s<sub>Y</sub>, respectively. Formally, if Z<sub>k</sub> is the 2D slice of Z at the kth spectral channel, and Z<sub>ij</sub> is the spectra of Z at pixel (i, j):

$$ (\mathcal{B}_Y Z)_k = B_Y \otimes Z_k, \qquad (3) $$

$$ (Z S_X)_{ij} = Z_{ij} \otimes S_X, \qquad (4) $$

where ⊗ represents the convolution operator (for simplicity of notation, this is written in the same way for 2D and 1D convolution).

The aim is to reach, as closely as possible given the noisy data, a fused cube Z that has the spatial resolution of X and the spectral resolution of Y. The general problem is

$$ \hat{Z} = \arg \min_{Z \ge 0} \mathcal{L}_P(X | Z_X) + \mathcal{L}_P(Y | Z_Y) + \varphi(Z), \qquad (5) $$

where

$$ \begin{aligned} Z_X &= (t_X Z \odot A_X) S_X \\ Z_Y &= \mathcal{B}_Y(t_Y Z \odot A_Y), \end{aligned} \qquad (6) $$

where ⊙ is the term-by-term product, the data-fidelity metric L<sub>P</sub> is the negative Poisson log-likelihood, and φ(Z) is the regularization term. We impose Z ≥ 0, as there is no physical interpretation for a negative spectra. The likelihood terms promote faithfulness between our model and the data at comparable resolutions.

It is worth noting that in order to fully invert the forward model, we would need to compare X with B<sub>X</sub>(t<sub>X</sub>Z ⊙ A<sub>X</sub>)S<sub>X</sub>, and likewise for Y. However, our aim is not to retrieve a cube Z with infinite resolution, but rather our intention is for Z to have the spatial resolution of X, and the spectral resolution Y. The

<footnote n="8">8 However, accounting for vignetting would be easy, since the effective area is simply multiplied to the model term by term. Thus whether we consider the effective area as a spatially uniform cube or a vignetted cube has no influence on the algorithmic complexity.</footnote>
problem expressed in Equations (5) and (6) will promote such a solution, since no spatial degradation is applied to go from Z to Z_X, and likewise spectrally for Z_Y with Y.

More explicitly, we obtain

$$ \hat{Z} = \arg \min_{Z \ge 0} \sum_{i=1}^{M} \sum_{j=1}^{N} \sum_{k=1}^{L} \left[ \begin{aligned} & (t_X Z \odot A_X) S_X \Big|_{ijk} \\ & - X_{ijk} \ln((t_X Z \odot A_X) S_X) \Big|_{ijk} \\ & + (\mathcal{B}_Y(t_Y Z \odot A_Y)) \Big|_{ijk} - Y_{ijk} \ln(\mathcal{B}_Y(t_Y Z \odot A_Y)) \Big|_{ijk} \end{aligned} \right] + \varphi(Z), \quad (7) $$

where the shape of Z, X, Y, are all (M, N, L).

## 3.2. Rebinning operator

In most cases, two different hyperspectral instruments will not produce data of the same voxel size. Rather, we will have $X \in \mathbb{R}_{+}^{M \times N \times I}$ and $Y \in \mathbb{R}_{+}^{m \times n \times L}$, where $L > l$, $M > m$, and $N > n$. In other words, X has more pixels than Y, and Y has more spectral channels than X. Our result Z should be of size ($M \times N \times L$). In order to compare Z to the data X and Y in the data fidelity terms of Equation (5), we thus need to change the size of Z to that of X and Y, respectively. This is done by applying a rebinning operator.

For the sake of simplicity, and to keep this operator differentiable, we assume that M is a multiple of m, N is a multiple of n, and L is a multiple of l. Let us denote the multiplying factors ($v_M$, $v_N$, $v_L$), such that $M = v_M m$, and likewise for the other dimensions. We need three rebinning operators: $\mathcal{R}_M$, $\mathcal{R}_N$, and $\mathcal{R}_L$.

In this case, the rebinning operator is that which sums together pixels in nearby groups of v.

For instance, for $\mathcal{R}_M$, a single voxel (s, j, k) of the rebinned tensor is

$$ (\mathcal{R}_M Z)(s, j, k) = \sum_{p=S^*v_M}^{p=(s+1)^*v_M-1} Z(p, j, k). \quad (8) $$

Such an operator can be neatly written as a matrix product between each slice of the tensor Z and a block diagonal matrix filled only with zeroes and ones. For simplicity, we use $\mathcal{R}_{MN}$ to denote the spatial rebinning operator (which applies $\mathcal{R}_M$ on the lines, and $\mathcal{R}_N$ on the columns).

Now, to include rebinning in our optimization problem posed in Equation (5), $Z_X$ and $Z_Y$ become

$$ \begin{aligned} Z_X &= (t_X Z \odot A_X) S_X \mathcal{R}_L \\ Z_Y &= \mathcal{R}_{MN} \mathcal{B}_Y(t_Y Z \odot A_Y), \end{aligned} \quad (9) $$

and plugging this into the Poisson likelihood, we obtain

$$ \hat{Z} = \arg \min_{Z \ge 0} \sum_{i=1}^{M} \sum_{j=1}^{N} \sum_{k=1}^{L} \left[ \begin{aligned} & (t_X Z \odot A_X) S_X \mathcal{R}_L \Big|_{ijk} \\ & - X_{ijk} \ln((t_X Z \odot A_X) S_X \mathcal{R}_L) \Big|_{ijk} + (\mathcal{R}_{MN} \mathcal{B}_Y(t_Y Z \odot A_Y)) \Big|_{ijk} \\ & - Y_{ijk} \ln(\mathcal{R}_{MN} \mathcal{B}_Y(t_Y Z \odot A_Y)) \Big|_{ijk} \end{aligned} \right] + \varphi(Z). \quad (10) $$

## 3.3. Regularization

In X-ray astrophysics, low signal-to-noise ratio presents a major problem, as HS images are dominated by Poisson noise. Thus, we are looking for regularization terms that will be efficient denoisers for Poisson statistics, and that will be appropriate for astrophysical images (often displaying diffuse and/or isotropic features) and X-ray spectra (which can present strong spectral variability).

Our first consideration was sparsity under a 2D-1D starlet transform. Then, inspired by the work of Guilloteau et al. (2022) regarding JWST, we implemented a low-rank approximation with Sobolev spatial regularization. Finally, we wished to test a compromise between the two, with a low-rank approximation and a 2D spatial wavelet sparsity regularization. The following sections describe these three options.

### 3.3.1. 2D-1D wavelet sparsity

A classic regularization method is to choose a term that promotes the sparsity of Z under a well-chosen orthogonal transform – that is to say, that in some domain, Z can be represented by only a small number of non-zero coefficients. Wavelets are useful for capturing multi-scale features. In particular, the isotropic undecimated wavelet transform, dubbed the starlet transform, is appropriate for astrophysical data due to its isotropic structure (Starck & Murtagh 1994). The data are convolved by $B_3$-spline filters, which produce a set of increasingly coarse resolution versions of the data, and the S wavelet scales are defined as the difference between two subsequent resolutions. The coarsest resolution is dubbed the coarse scale, and no assumption of sparsity is placed upon it. This allows the regularization to represent data with a smooth continuum (the coarse scale) and sharp features (the non-zero elements of the fine scales).

More specifically, we chose the 2D-1D starlet transform, which we define as applying the 2D starlet transform channel-wise, then the 1D starlet transform pixel-wise. Thus, in this case, $\varphi(Z)$ is the $l_1$ norm (the minimization of which promotes sparsity; see Bach et al. 2011) of the 2D-1D starlet transform coefficients of Z. For S wavelet scales, we define the 2D-1D starlet transform as $\mathcal{W}_S^{2D1D}$, and the $l_1$ norm regularization term is

$$ \begin{aligned} \varphi(Z) &= \|\mathcal{W}_S(Z)\|_1 \\ &= \mu \sum_{i=1}^{M} \sum_{j=1}^{N} \sum_{k=1}^{L} \sum_{s=1}^{S} |\mathcal{W}_S(Z)|_{i,j,k,s}^{2D1D}, \end{aligned} \quad (11) $$

where $\mu$ is a hyperparameter to control the intensity of regularization.

### 3.3.2. Low rank and Sobolev

This regularization option was inspired by work regarding JWST carried out by Guilloteau et al. (2020), though we designed an implementation adapted for HS-HS fusion and the specifics of X-ray astronomy. The spectral regularization is a low-rank approximation, which corresponds to decomposing Z = WV, where $V \in \mathbb{R}^{L \times r}$ spans the subspace of rank r onto which Z is projected, and $W \in \mathbb{R}^{r \times M \times N}$ is the tensor of representation coefficients in that subspace. V is found by performing principal component analysis on Y (the hyperspectral image with the most spectral information) and keeping only the highest singular values. Therefore, now $Z_X$ and $Z_Y$ from Equation (9) become
variables of W, where

$$ 
\begin{aligned}
Z_X(W) &= (t_X W V \odot A_X) S_X R_L \\
Z_Y(W) &= R_{MN} B_Y (t_Y W V \odot A_Y),
\end{aligned}
\qquad (12) $$

and these can be plugged into Equation (10), now seeking to minimize W. For spatial regularization, the Sobolev regularization (Karl 2005) implies adding the following term to the cost function:

$$ \varphi(W) = \mu \|\mathcal{D}(W)\|_{F}^{2}, \qquad (13) $$

where $\mathcal{D}$ is a first-order 2D finite difference operator. This term will act to minimize the differences between neighboring pixels on the representation coefficients of the subspace spanned by $V$.

Both the Sobolev regularization and the low-rankedness will result in a faster algorithm, due to the dimension reduction and simple expression of the $\mathcal{D}$ operator. However, Sobolev does not capture multi-scale spatial features. Further, as $V$ is extracted from $Y$, it may reproduce features found in $Y$ too strongly. This is acceptable (and even desirable) in MS/HS fusion, where the MS image does not bring a significant amount of spectral information, but in HS/HS fusion, $X$ contains spectral information, in particular because it displays much less spatial mixing than $Y$. The assumption that all spectra in $Z$ can be decomposed as linear combinations of $Y$ is not necessarily true, and forcing a low-rank approximation in this case will bring on spectral distortions. However, if this assumption is valid, low-rankedness will be an undeniable asset.

### 3.3.3. Low rank and 2D wavelet sparsity

As Sobolev regularization does not capture multi-scale information in the same way wavelets do, we considered the option to keep a low-rank approximation as described in the previous section, and then to use a starlet sparsity term for the spatial regularization:

$$ \varphi(W) = \mu \|W_S^{2D}(W)\|_1, \qquad (14) $$

where the operator $W_S^{2D}$ applies the 2D starlet transform to every slice of $W$.

## 3.4. Algorithm description

### 3.4.1. Proximal gradient descent

If we consider the problem posed in Equations (7) or (10), and replace the regularization term by one of the three options we consider above, we note that the cost function is convex, and can be separated into two parts, one differentiable and one non-differentiable. In such cases, proximal gradient methods have been shown to be very effective (Parikh & Boyd 2014).

The proximal gradient descent is an iterative method where, for each iteration, we first undergo a gradient descent step on the differentiable part of the cost function before applying a proximal operator (prox), the purpose of which is to find a solution close to the result obtained by gradient descent while minimizing the function's non-differentiable part. Formally, if the result obtained by the gradient descent step is $x$, and the non-differentiable part of the function is $\varphi$, the operator is:

$$ \text{prox}_{\varphi}(x) = \arg\min_{y}(\varphi(y) + \frac{1}{2}\|y - x\|_{2}^{2}). \qquad (15) $$

For the $l_1$ norm, the proximal operator is a soft thresholding function:

$$ \text{prox}_{l_1, \mu}(x) = \begin{cases} 0 & \text{if } |x| < \mu \\ x - \mu \cdot \text{sign}(x) & \text{if } |x| \ge \mu \end{cases}, \qquad (16) $$

where $\mu$ is the parameter that controls the strength of regularization. For the non-negativity constraint, the prox is simply the operator that forces all negative values to 0.

The gradient of the differentiable part in Equation (5) is

$$ 
\begin{aligned}
\nabla_Z \mathcal{L}_P &= \nabla_Z (\mathcal{L}_P(X | Z_X)) + \nabla_Z (\mathcal{L}_P(Y | Z_Y)) \\
&= \frac{\partial(Z_X)}{\partial Z} \nabla_{Z_X} (\mathcal{L}_P(X | Z_X)) + \frac{\partial(Z_Y)}{\partial Z} \nabla_{Z_Y} (\mathcal{L}_P(Y | Z_Y)). \quad (17)
\end{aligned}
 $$

In the case with rebinning, we have

$$ 
\begin{aligned}
\frac{\partial(Z_X)}{\partial Z} &= t_X A_X R_L^\dagger S_X^\dagger, \\
\frac{\partial(Z_Y)}{\partial Z} &= t_Y A_Y B_Y^\dagger R_{MN}^\dagger,
\end{aligned}
\qquad (18) $$

where $\dagger$ denotes the conjugate transpose of the operators, in the sense that, if we were to vectorize the spatial dimensions of our data and write the operator $S$ as a product with the 2D matrix $S$, $S^\dagger$ would apply the product with matrix $S^\dagger$. In the case with a low-rank approximation, where the gradient descent is done on $W = ZV^T$:

$$ 
\begin{aligned}
\frac{\partial(Z_X)}{\partial W} &= V t_X A_X R_L^\dagger S_X^\dagger, \\
\frac{\partial(Z_Y)}{\partial W} &= V t_Y A_Y B_Y^\dagger R_{MN}^\dagger.
\end{aligned}
\qquad (19) $$

In the case of Sobolev regularization, $\varphi_{Sobolev}(W) = \mu \|\mathcal{D}(W)\|_{F}^{2}$ is differentiable, and so a term is added to the gradient:

$$ \varphi'_{Sobolev}(W) = 2\mu \mathcal{D}^T \mathcal{D}(W). \qquad (20) $$

As is the case for the other operators, $\mathcal{D}$ can be interpreted as a product with the 2D matrix $D$, since the operator $\mathcal{D}$ is a convolution with a kernel of $[1, -1]$ on the lines followed by a convolution with a kernel of $[1, -1]^T$ on the columns. $\mathcal{D}^T$ applies the product with $D^T$.

In the code, the operators $\mathcal{R}$ and $\mathcal{D}$ are written as matrix products. The operators $B_Y$ and $S_X$, however, can be dealt with more efficiently by taking advantage of the diagonalization of the convolution operation in the Fourier domain. To do so, we pad the kernels and $Z$ with zeros so that they are the same size (dimension of $Z$ + dimension of the kernel in question). We then take their Fourier transform, and convolution becomes a Hadamar product in the Fourier space. The Fourier transform of $Z$ will need to be calculated at each iteration, but the transform of the kernels need only be calculated once.

Putting all this together, we obtain Algorithm 22, which we have dubbed HIFReD: Hyperspectral Image Fusion with Regularized Deconvolution. For simplicity, rebinning is not explicit in the pseudo-code.
**Algorithm 1: HIFReD: Hyperspectral Image Fusion with Regularized Deconvolution**

1 **Input:** data $(X, Y)$, instrument response kernels $(S_X, B_Y)$, regularization function $(\varphi)$, Low Rank (True/False), rank $r$, regularization parameter $\mu$, number of wavelet scales $S_{2D}, S_{1D}$, maximum number of iterations $t_{max}$.

2 **PWS:** Proximal Wavelet Step, defined in Alg. 2.

3 **Output:** $Z$

4 **Initialization**
Calculate $\mathcal{F}^{1D}(S_X), \mathcal{F}^{2D}(B_Y)$.
$Z \leftarrow X$

if Low Rank then

5 $\hat{V} \leftarrow \text{eigvectors of } (Y - \bar{Y})(Y - \bar{Y})^T$

6 $V \leftarrow \hat{V}[0: r,]$

7 **else**

8 $\lfloor V \leftarrow 1$

9 $t \leftarrow 0$

10 **while** $t < t_{max}$ **do**

11 $W \leftarrow V^T Z$
$\tilde{W}_t = \mathcal{F}^{1D}(\mathcal{F}^{2D}(W_t))$.
**—Gradient Step—**
**if** $\varphi(W) = \|\mathcal{D}(W)\|_F^2$ **then**
12 $\left| \tilde{W}_{t+1} \leftarrow \tilde{W}_t - \alpha[\nabla_{\tilde{W}}\mathcal{L}_P(\tilde{W}_t) - \mu\varphi'_{Sobolev}(W_t)] \right|$ 
Eqs. (17, 20).

13 **else**

14 $\lfloor \tilde{W}_{t+1} \leftarrow \tilde{W}_t - \alpha[\nabla_{\tilde{W}}\mathcal{L}_P(\tilde{W}_t)] \rfloor$ Eq. (17).

15 $W_{t+1} \leftarrow \mathcal{I}\mathcal{F}^{1D}(\mathcal{I}\mathcal{F}^{2D}(\tilde{W}_{t+1}))$

**—Proximal Step—**

**if** $\varphi(W) = \|\mathcal{W}_S^{2D}(W)\|_1$ **then**

16 $\lfloor W_{t+1} \leftarrow \text{PWS}[W_{t+1}, \mathcal{W}_{\mu,S2D}^{2D}] \rfloor$

17 **if** $\varphi(W) = \|\mathcal{W}_S^{2D-1D}(W)\|_1$ **then**

18 $\lfloor W_{t+1} \leftarrow \text{PWS}[W_{t+1}, \mathcal{W}_{\mu,(S1D,S2D)}^{2D-1D}] \rfloor$

19 $Z_{t+1} \leftarrow W_{t+1}V$

20 $Z_{t+1} \leftarrow Z_{t+1} \odot \delta(Z_{t+1} \ge 0)$

21 $t \leftarrow t + 1$

22 **return** $Z_t$

**Algorithm 2: PWS: Proximal Wavelet Step**

1 **input** cube $Z$, wavelet transform $\mathcal{W}$, regularization parameter $\mu$, number of wavelet scales $S$
$\{w_c, w_1, ..., w_S\} \leftarrow \mathcal{W}(Z)$

$\rho \leftarrow \mu \textbf{mad}(w_{1,1})$ (median absolute deviation)

**for** $i$ **in** $\{1, ..., S\}$ **do**

2 $\lfloor w_i \leftarrow \text{prox}_{l1}(w_i, \rho) \rfloor$ Eq. (16)

3 $\hat{Z} \leftarrow w_{c,c} + \sum_{i=1}^{S} w_i$
Return($\hat{Z}$)

### 3.4.2. Hyperparameters

For the choice of hyperparameters (namely the regularizing factor $\mu$, the gradient step $\alpha$, and the rank $r$), we need to take a heuristic approach. Indeed, the Poisson likelihood is not globally

Lipschitz-continuous, and so we cannot analytically calculate an optimal gradient step, and there is no a priori best choice for the regularization parameters either. Further, we find that one set of hyperparameters that works best for one toy model is not necessarily ideal for another model. The photon count, level of spectral variability, shape of the effective area, and more are all factors that will impact the choice of best parameters. It would be very costly to compute every single combination of hyperparameters to find the best suited for every single model, but below, we describe some guidelines to help pick suitable parameters, if not the exact optimal value.

The gradient step should be taken as high as possible for a faster rate of convergence, but if it is taken too high, the likelihood will diverge dramatically. This is easily diagnosed, and is seen as huge jumps in the cost function, when it should rather decrease monotonously. This usually happens in early iterations, and so we checked the value of the likelihood periodically at the start (first 1000 iterations, which takes less than an hour), and then stopped if it became pathological. If it this was the case, we then lowered the gradient step and tried again with that parameter. As the Poisson likelihood is convex, there is no risk of falling into a local minimum, provided this step is well-chosen as described.

For the $\mu$ parameter of the Sobolev regularization, for the Gaussian toy model, we tested values of between 0.01 and 0.1, and found the best value to be 0.03. This value turned out to be well suited for the other toy models as well. An overly high value of this parameter is detectable as an over-smoothing of spatial features, while an excessively low value will be recognized in the presence of visible noise on spatial slices.

For the rank $r$, we tested values of between 5 and 50, and selected 10. The rank has a similar impact to $\mu$, but on the spectral features. A lower rank might over-smooth and lead to a biased result, while a higher rank might not be enough of a denoiser. On top of that, the rank will impact the algorithm's computational complexity, with a lower rank implying fewer dimensions and thus a faster algorithm.

For the wavelet threshold parameter $\mu$, this is usually taken to be between 0.5 and 3, which in the case of Gaussian noise is comparable to a $0.5 \sigma - 3 \sigma$ interval. There is no such easy interpretation for Poisson noise, but even so, like with the Sobolev, lower $\mu$ will allow more variability and less bias, while higher $\mu$ will lead to more smoothing, and thus potentially more bias. In particular, we find that in dim voxels (close to zero, and dominated by noise), higher values of $\mu$ (close to 3) are effective in smoothing out noise, while in brighter voxels, those high $\mu$ will over-smooth spectral features. We find the opposite to be true with lower values of $\mu$. As a compromise, and to account for the high variabilities of flux between voxels, it is advantageous to apply a varying $\mu$, using a sigmoid function going from 0.5 to 3. Let the coarsest scale of $Z$ in the wavelet domain be $C_Z$. Further, let $\tilde{C} = |C_Z| - |\tilde{C}_Z|$. In this case,

$$ \mu = \frac{(3 - 0.5)}{1 + e^{\tilde{C}}} + 0.5, \qquad (21) $$

which will make $\mu(i, j, k)$ close to 1.25 when $|C_Z|$ is close to $|\tilde{C}_Z|$, and go towards 0.5 for small values of $|C_Z|$, and towards 3 for high values of $|C_Z|$. This simple scheme is however ineffective in cases where the toy model is not very sparse in the wavelet domain, as is the case for the realistic model. In this case, it is better to pick a low, constant value of $\mu$, such as 0.5. Further, one needs to consider the number of wavelet scales defined in the starlet transform, be it the 2D or the 2D-1D transform. As
<table><caption>**Table 1.** Characteristics of the simulated data sets.</caption><thead><tr><th>Model</th><th colspan="2">Gaussian</th><th colspan="2">Gaussian (rebinning)</th><th colspan="2">Realistic (rebinning)</th></tr><tr><th>Data set</th><th>X</th><th>Y</th><th>X</th><th>Y</th><th>X</th><th>Y</th></tr></thead><tbody><tr><td>Number of pixels</td><td colspan="2">45×45</td><td>45×45</td><td>15×15</td><td>45×45</td><td>15×15</td></tr><tr><td>Number of spectral channels</td><td colspan="2">(a) 2500 – (b) 2000</td><td>60</td><td>2400</td><td>60</td><td>2400</td></tr><tr><td>Energy range</td><td colspan="2">(a) 0.5–1.4 keV – (b) 6.2–6.9 keV</td><td colspan="2">0.5–1.4 keV</td><td colspan="2">0.5–1.4 keV</td></tr><tr><td>PSF std</td><td>10"</td><td>30"</td><td>10"</td><td>30"</td><td>11"</td><td>33"</td></tr><tr><td>Spectral response std</td><td>21 eV</td><td>2.5 eV</td><td>21 eV</td><td>2.5 eV</td><td>51 eV</td><td>1.4 eV</td></tr><tr><td>Photon count</td><td>(a) 2.4 × 10<sup>7</sup><br>(b) 3.2 × 10<sup>7</sup></td><td>(a) 2.3 × 10<sup>7</sup><br>(b) 3.1 × 10<sup>7</sup></td><td>2.7 × 10<sup>6</sup></td><td>2.5 × 10<sup>6</sup></td><td>4.5 × 10<sup>7</sup></td><td>3.0 × 10<sup>7</sup></td></tr></tbody></table>

**Notes.** The spectral and spatial responses are given in standard deviation (std).

the coarse scale is left untouched, the number of fine scales will determine the maximum size of the features that can be thresh- olded, as for each scale the size of the filter increases by a factor of 2 (Starck et al. 2010). The number of 2D wavelet scales was chosen to be 2 for all models, and the number of 1D scales was selected to be 3 for the Gaussian model, and 1 for the realistic model, for which the spectral details are much finer due to the finer spectral response.

The last hyperparameter to consider is $t_{max}$, the maximum number of iterations. The user should check the cost function curve at the end to see whether optimization has converged. If not, they may begin another batch of iterations, using the last result as first guess. An automatic stopping criterion could eas- ily be implemented, with iterations stopping at iteration *t* once $\left|\frac{\mathcal{L}_p(Z_t) - \mathcal{L}_p(Z_{t+1})}{\mathcal{L}_p(Z_t)}\right| > \epsilon$ for a small $\epsilon$, but the results in this article were obtained with a manually set $t_{max}$.

# 4. Results on simulated data

In this section, we present the results of our method on toy models derived from hydrodynamics simulations of a SNR. Section 4.1 presents our three considered toy models: the Gaussian model, the Gaussian model with rebinning, and the realistic model with rebinning. Section 4.2 then presents our results for each using the three regularization schemes described in Section 3.3.

## 4.1. Toy models

The toy models we constructed are intended to resemble the case study of SNRs, such as Cassiopeia A. To this end, we used the 3D hydrodynamics simulations of Orlando et al. (2016). For each voxel of the 3D simulation, a synthetic X-ray spectrum was gen- erated using the non-equilibrium ionization Xspec model NEI<fnref n="9" />. In the numerical simulation dataset, the temperature kT ranges from 0.2 keV to 3 keV, the ionization timescale log($\tau$) from 8 to 11.5 cm⁻³ s, and the velocity along the line of sight $V_z$ from

-7000 to 7000 km s⁻¹. The absorption along the line of sight was fixed to $n_H = 0.5 \times 10^{22}$ cm⁻². Then the hypercube (X, Y, Z, E) is projected along the Z axis to emulate a synthetic X-ray observation of a SNR. Due to the projection effect, the spectrum for each pixel can be quite complex, and can contain, for example, double-peaked spectral features associated with the red- and blueshifted half spheres of the SNR or different temperatures and ionization timescales.

We then wished to consider different levels of complexity for our algorithm benchmarking. In particular, we considered three characteristics that can make the problem more complex: the intensity of spectral variability, the inclusion of a rebinning operator (X and Y have different voxel size), and the ratio between the spectral and spatial responses of the two data sets.

The first toy model we considered is dubbed the Gaussian model. The PSFs and spectral responses of each instrument are assumed to be Gaussian, whose parameters are summarized in Table 1, and shown in Figs. 3 and 4. Further, we assumed a flat effective area, and that X has the same size as Y. With these parameters, we generated one set of X and Y between 0.5 and 1.4 keV, and one between 6.2 and 6.9 keV (around the iron emission lines). The dataset around the iron lines displays more spectral variability, as is demonstrated more formally in Fig. 5, which shows the average angular distance $\theta$ between neighboring pixels for each model. For a given pixel's spectra $Z_k$, it is calculated as

$$ \theta(Z_k) = \sum_{Z_{k'} \in N} \arccos \left( \frac{\langle Z_k, Z_{k'} \rangle}{\sqrt{\|Z_k\|^2 \|Z_{k'}\|^2}} \right) / \text{card}(N), \quad (22) $$

where *N* is the set of pixels neighboring $Z_k$, that is, the eight closest pixels to $Z_k$, and card($N$) is the number of elements in that set. This angular distance provides information about the spectral variations across the remnant. As we see, the Gaussian model around the iron line has the strongest variations.

Then, we considered the same Gaussian model (same PSF/spectral response, same flat effective area) between 0.5 and 1.4 keV, but this time with X and Y having different voxel sizes. X has the dimensions (45 × 45 × 60), with square pixels of 10" in width and energy channels of 14 eV in width, and Y has the dimensions (15 × 15 × 2400), with square pixels of 30" in width

<footnote n="9">⁹ https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/node200.html</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Comparison of the four PSFs used to generate the toy models. On the left, the Gaussian approximation of a XRISM-like and XMM-Newton-like instrument. On the right, the actual PSFs of XRISM/Resolve and XMM-Newton/EPIC. All PSFs are shown on images of size 4'50"×4'50". Colors are displayed in square-root scale.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Comparison of the four spectral responses used to generate the toy models. The responses for the XMM-Newton-like instrument are shown on the right, while those for the XRISM-like are shown on the left. The full orange line is for the Gaussian approximations, and the dashed blue line is for the realistic instrument responses.</figcaption></figure>

and energy channels of 0.36 eV in width. This allows us to see the impact of the rebinning operator on our algorithm.

Finally, we set out to obtain a more realistic model, which would mimic the case where we try to fuse data from XMM-Newton/EPIC and XRISM/Resolve. We chose PSFs and spectral responses meant to be as close as possible to those of the instruments in question (see Figs. 3 and 4), and likewise for the effective areas. The energy range was (0.5–1.4 keV), and we chose the same dimensions and voxel size as the Gaussian model with rebinning.

In all of these case studies, we generated toy models with a photon count of around $10^7$. We then applied Poisson noise. Some sample spectra and 2D slices are shown in Fig. 6.

## 4.2. Results

This section presents the results we obtained with the four toy models. We compare the wavelet 2D–1D regularization (W2D1D), the low-rank approximation with Sobolev regularization (LRS), and the low-rank approximation with wavelet 2D regularization (LRW2D). We present our results in

the form of sample pixel fits, amplitude maps, error histograms, error maps, and scalar metrics.

First, the reconstructed spectra obtained on two sample pixels are shown in Figure 7. Once again, for visual clarity, we only show the W2D1D and LRS methods. These results allow us to draw several conclusions. First of all, the Gaussian model around the iron line (c, d) is difficult to reconstruct for both methods, but the LRS method appears more biased. Second, we see that the W2D1D method results in smoother spectra, but this can be accompanied by more bias, as seen in the Gaussian model with rebinning (e, f). This occurs when the underlying assumption of sparsity in the wavelet domain is not appropriate.

The amplitude maps are displayed in Fig. B.1. Overall, for all regularizations, we see that the method is able to reconstruct the general spatial features of the remnant. For the Gaussian model with rebinning, we do notice some amount of leakage on the edges of the remnant for all regularizations, particularly on the western side. For the Gaussian model around 6 keV (b), we also see that LRS has over-smoothed some of the features – this effect is visible in the other toy models but is more subtle. Finally, in models with rebinning (c, d), the LRW2D seems to smooth out finer features and give a biased result.

In Fig. B.2, we present the map of the spectral angle mapper (SAM), an error metric defined as (for a spectra $Z_i$)

$$ SAM(Z_i) = \arccos \left( \frac{\langle Z_i, \hat{Z}_i \rangle}{\|Z_i\|_2 \|\hat{Z}_i\|_2} \right). \qquad (23) $$

The SAM is best if low (close to zero), and is designed to capture spectral distortions. In Fig. B.2, we observe some leakage on the low-count pixels at the edges of the remnant, which is present for all methods, but more so for W2D1D. The results shown in panel a suggest that the wavelet regularization presents slightly more spatial distortion for the simplest model (Gaussian model with low spectral variability), while it appears more successful in panel b for the Gaussian model with high spectral variability. In panel b, the higher SAM pixels in LRS and LRW2D correspond to some local features, in areas of high spectral variability (see the middle of Fig. 5); we highlight, for instance, the cluster of bright spots above the center of the remnant. In panel c, the low rank with Sobolev appears to be doing better in the center. All models display a high error close to the bright spot south of the remnant. In panel d, for the realistic model, LRS obtains better result in the center of the remnant. While all three methods struggle to reach a result close to the ground truth in the bright area to the south of a remnant, LRW2D and W2D1D show the poorest performance, which could suggest that the Sobolev regularization is better at handling outliers. As this bright spot is also seen as a problem in the model shown in panel c, this area of high flux and high spatial variability seems especially difficult to reconstruct in models with rebinning.

In Fig. B.3, we display the relative error maps averaged per pixel. We draw similar conclusions about the Gaussian model without rebinning. The difficulty that low-rank methods have in reconstructing high spectral variabilities areas is once again seen in panel b. For panel c, we again notice the leakage on the edges, with the amplitude being judged too high in red areas, and too low in blue areas. There is a general bias on all models that seems to overestimate the northwest part of the remnant, while underestimating the southeast part. This is likely due to the rebinning, as it is also present in panel d. Further, in panel d, we note that W2D1D and LRW2D underestimated the amplitude of the bright spot in the south of the remnant, while LRS was not
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Average angular distance between neighboring pixels on a window of neighboring pixels, given by Equation (22). Higher values indicate a location of high spectral variability between pixels.</figcaption></figure>

<table><caption>**Table 2.** Error metric results comparing regularizations for the four toy models.</caption><thead><tr><th colspan="5">Gaussian (0.5–1.4 keV)</th></tr><tr><th>Regularization</th><th>NMSE</th><th>aSAM</th><th>acSSIM</th><th>ERGAS</th></tr></thead><tbody><tr><td>W2D1D</td><td><b>21.36</b></td><td>0.571</td><td>0.017</td><td>0.220</td></tr><tr><td>LRS</td><td>20.97</td><td><b>0.557</b></td><td><b>0.011</b></td><td>0.116</td></tr><tr><td>LRW2D</td><td>21.01</td><td>0.560</td><td>0.012</td><td><b>0.114</b></td></tr></tbody><thead><tr><th colspan="5">Gaussian (6.2–6.9 keV)</th></tr><tr><th>Regularization</th><th>NMSE</th><th>aSAM</th><th>acSSIM</th><th>ERGAS</th></tr></thead><tbody><tr><td>W2D1D</td><td><b>18.67</b></td><td><b>0.599</b></td><td><b>0.005</b></td><td><b>0.229</b></td></tr><tr><td>LRS</td><td>14.16</td><td>0.600</td><td>0.011</td><td>0.343</td></tr><tr><td>LRW2D</td><td>16.57</td><td>0.609</td><td>0.010</td><td>0.248</td></tr></tbody><thead><tr><th colspan="5">Gaussian with Rebinning (0.5–1.4 keV)</th></tr><tr><th>Regularization</th><th>NMSE</th><th>aSAM</th><th>acSSIM</th><th>ERGAS</th></tr></thead><tbody><tr><td>W2D1D</td><td>10.36</td><td>0.691</td><td>0.048</td><td>0.422</td></tr><tr><td>LRS</td><td><b>11.01</b></td><td><b>0.670</b></td><td><b>0.045</b></td><td>0.341</td></tr><tr><td>LRW2D</td><td>10.34</td><td>0.690</td><td>0.054</td><td><b>0.317</b></td></tr></tbody><thead><tr><th colspan="5">Realistic with rebinning (0.5–1.4 keV)</th></tr><tr><th>Regularization</th><th>NMSE</th><th>aSAM</th><th>acSSIM</th><th>ERGAS</th></tr></thead><tbody><tr><td>W2D1D</td><td>9.956</td><td>0.700</td><td>0.032</td><td>0.608</td></tr><tr><td>LRS</td><td><b>10.16</b></td><td><b>0.691</b></td><td><b>0.031</b></td><td>0.358</td></tr><tr><td>LRW2D</td><td>9.523</td><td>0.712</td><td>0.038</td><td><b>0.351</b></td></tr></tbody></table>

**Notes.** The NMSE (Equation (A.1)) is better if high, while the aSAM (Equation (A.2)) and acSSIM (Equation (A.3)) are better if low. The best values are marked in bold.

as biased. Overall, W2D1D displays more bias in the southern bright outlier than the other methods.

Finally, we measured scalar error metrics. All results can be found in Table 2. The best results are highlighted in bold. As Zhang (2008) showed for the case of fusion between panchromatic and MS images, scalar metrics do not necessarily agree when comparing the efficacy of fusion methods, and are not always easily interpretable. Thus, we present four metrics in order to provide a more complete comparison. The definitions of the metrics can be found in Appendix A. Namely, we calculate the normalised mean squared error (NMSE), the average spectral angle mapper (aSAM), average complementary structural similarity (acSSIM) index, and the relative dimensionless global error (ERGAS).

Looking at the obtained values in Table 2, we see that the four metrics rarely draw the same conclusion, except in the case of the Gaussian model around 6 keV, which is the model with the highest spectral variability between pixels. This allows us to draw a similar conclusion as for the error maps: the methods using a low-rank approximation have more trouble accounting for high spectral variability. This may be due to the fact that the low-rank approximation is learnt by performing PCA on $Y$, the data set that has the best spectral resolution, but its poor spatial resolution will smooth out spectral variability between pixels. If the low rank could be learnt on $Z$, this would not be a problem, but of course that is not possible.

The results on the Gaussian model around 1 keV show that all three methods are the best in at least one metric, and so it is hard to choose an outright leader. In this case, low rankedness does seem to lower spectral distortion (looking at the aSAM and the ERGAS). For the Gaussian model around 1 keV with rebinning, the low-rank methods are more successful, but the difference is small, and the bright southern spot we note on the maps likely plays a part in the W2D1D being a little less successful. We also observe this for the realistic model with rebinning. Overall, it appears that the low-rank methods are best at dealing with outliers, especially the low rank with Sobolev, but they are not as appropriate for models with high spectral variability.

Convergence took between 30k and 100k iterations, depending on the toy model. On average, one iteration of the algorithm took around 3 s with the method with 2D-1D wavelet regularization, resulting in a computation time of between 25 h and 83 h. With the low-rank approximation with Sobolev, each iteration took around 2 s, resulting in a computation time of between 15 h and 55 h; the low rank with wavelet 2D was around the same.

## 5. Discussion

Overall, the wavelet 2D-1D regularization is more successful than the low-rank approximations at dealing with strong spectral variability across the image. However, the assumption of wavelet sparsity is not necessarily ideal to account for certain structures, such as in the presence of singular bright spots like the one in the southern part of the remnant, and this is especially difficult to treat with the inclusion of a rebinning operator. Nonetheless, we find that for all three regularization methods, the fusion results are satisfactory given the difficulty of the problem; however, they still contain some bias.

The low-rank approximation, while ill-suited for SNRs that display strong spectral variability due to the variations of their physical parameters, was found to accelerate the algorithm by
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Example of spectra and 2D slices from the four toy-model data sets. The ground truth $\hat{Z}$ is shown in black, and the two data sets in purple and orange: X with the higher spatial resolution, and Y with the higher spectral resolution. (a) Gaussian model between 0.5 and 1.4 keV, (b) Gaussian model between 6.2 and 6.9 keV, (c) Gaussian model with rebinning between 0.5 and 1.4 keV, and (d) the Realistic model between 0.5 and 1.4 keV.</figcaption></figure>

half thanks to the dimension reduction. As the data sets are large and the algorithm takes a long time to converge, this is a useful attribute.

Thus, it would be ideal to find a spectral regularization less biased than the wavelet sparsity or the low-rank approximation, but that nonetheless includes a reduction in dimensions. This may be done by exploring a low-rank approximation

more refined than one learnt using PCA on the high-spectral-resolution data set. We could also explore acceleration methods, such as algorithm unrolling, as in Fahes et al. (2022).

Alternatively, a promising avenue would be to include physical information. Indeed, while having a blind fusion result is useful (as it requires no assumptions on the model), in order to obtain physical parameter maps, a model will need to be
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Example of how the fusion output compares to the ground truth for two pixels. From top to bottom, the models are: (a, b) Gaussian model between 0.5 and 1.4 keV, (c, d) Gaussian model between 6.2 and 6.9 keV, (e, f) Gaussian model with rebinning, and (g, h) the Realistic model. The left column shows results for pixel (22, 22), which is lower in brightness than pixel (11, 31) on the right.</figcaption></figure>

fitted to the result. In the case of SNRs, more than one model would be needed, as the emission comes from a mixture composed of thermal radiation from the shocked interstellar media and synchrotron radiation due to particle acceleration.

An algorithm for source separation that includes spectral variations was developed by Lascar et al. (2024). Our future work will focus on combining this approach to HIFReD, we the aim being to develop a method capable of simultaneously fusing data

sets from two generations of telescopes and perform unmixing with a physical spectral model.

On a final note, one limitation of the proposed method is that as it is implemented, it cannot account for spatially varying PSF, or a spectrally varying spectral response. In this case, the convolution is no longer a term-by-term product in Fourier. This is an issue with XMM/Resolve, where a constant spectral response over a wide energy range is not realistic. Therefore, we are
currently working on developing a version of HIFReD than can
handle varying convolution kernels, with promising results so
far. The main obstacle is the computational complexity, which
is greatly increased, but by working on an accelerated version of
the code using the JAX<fnref n="10" /> Python package, we find that we should
be able to achieve even faster performances than those described
in this article, even when using a non-stationary convolution.

6. Conclusions

We have developed an algorithm for fusing hyperspectral images
taken by two generations of X-ray telescopes with complemen-
tary spatial and spectral resolutions. Our method simultaneously
deconvolves the instruments' spatial and spectral responses in
order to return to the best resolutions of the two telescopes.

We implemented three regularization schemes, the first being a term to promote sparsity under a 2D-1D wavelet transform, the second being a low-rank approximation with Sobolev regularization (inspired by Guilloteau et al. 2020), and the third being a low-rank approximation with a 2D wavelet sparsity term.

We tested our method on five toy models constructed from
hydrodynamic simulations of SNRs by Orlando et al. (2016)
and convolved with an increasing complexity of instrumental
response. Spatial features are well reconstructed, apart from a
slight leakage in the bright shell in the toy model with rebin-
ning. The reconstruction of spectral features depends on the
regularizations and model.

We find that while low-rank methods accelerated the algo-
rithm by half, they struggled to account for the strong spectral
variability present around the iron line. The 2D-1D wavelet
method better accounts for this variability, but it is somewhat
biased in reconstructing the realistic model.

In conclusion, we believe that hyperspectral fusion
shows promise and would greatly benefit the analysis of
XRISM/Resolve data; though it would be useful to explore other
regularization terms more suited to the variabilities of the data.
In particular, we intend to combine fusion with a spectral phys-
ical model and source-separation algorithm in future works.
While the purpose of this article is to present a proof of con-
cept, we are now actively working on applying this method to
real data.

**Acknowledgements.** The research leading to these results has received funding from the European Union's Horizon 2020 Programme under the AHEAD2020 project (grant agreement no. 871158). This work was supported by CNES, focused on methodology for X-ray analysis. This work was supported by the Action Thématique Hautes Energies of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. We thank S. Orlando for kindly providing the simulation from Orlando et al. (2016) to generate our toy model.

References

Bach, F., Jenatton, R., Mairal, J., & Obozinski, G. 2011, Found. Trends Mach. Learn., 4, 1

Barret, D., Lam Trong, T., den Herder, J.-W., et al. 2018, SPIE Conf. Ser., 10699, 106991G

Bertero, M., Boccacci P., & De Mol, C. 2022, Introduction to Inverse Problems in Imaging (CRC Press)

Fahes, M., Kervazo, C., Bobin, J., & Tupin, F. 2022, in International Conference on Learning Representations

Godinaud, L., Acero, F., Decourchelle, A., & Ballet, J. 2023, A&A, 680, A80

Guilloteau, C., Oberlin, T., Berné, O., Émilie Habart, & Dobigeon, N. 2020, AJ, 160, 28

Guilloteau, C., Oberlin, T., Berné, O., & Dobigeon, N. 2022, in IEEE International Conference on Image Processing (ICIP 2022), Bordeaux, France, 1

Hadamard, J. 1923, Lectures on the Cauchy Problem in Linear Partial Differential Equations (New Haven: Yale University Press)

Karl, W. C. 2005, in L Handbook of Image and Video Processing (Elsevier)

Lascar, J., Bobin, J., & Acero, F. 2024, A&A, 686, A259

Orlando, S., Miceli, M., Pumo, M. L., & Bocchino, F. 2016, ApJ, 822, 22

Parikh, N., & Boyd, S. 2014, Found. Trends Optim., 1, 127

Pineau, D., Orieux, F., & Abergel, A. 2023, Fast Multispectral and Hyperspectral Image Fusion via Hessian Inversion. hal-04337970

Prévost, C., Borsoi, R. A., Usevich, K., et al. 2022, SIAM J. Imaging Sci., 15, 110

Simões, M., Bioucas-Dias, J., Almeida, L. B., & Chanussot, J. 2015, IEEE Trans. Geosci. Rem. Sens., 53, 3373

Starck, J.-L., & Murtagh, F. 1994, A&A, 342

Starck, J.-L., Murtagh, F., & Fadili, Jalal, M. 2010, Sparse Image and Signal Processing: Wavelets, Curvelets, Morphological Diversity (Cambridge University Press)

Takahashi, T., Mitsuda, K., Kelley, R., et al. 2014, SPIE Conf. Ser., 9144, 914425

Turner, M. J. L., Abbey, A., Arnaud, M., et al. 2001, A&A, 365, L27

Wang, Z., Bovik, A., Sheikh, H., & Simoncelli, E. 2004, IEEE Trans. Image Process., 13, 600

XRISM Science Team 2020, arXiv e-prints [arXiv:2003.04962]

Xu, T., Huang, T.-Z., Deng, L.-J., & Yokoya, N. 2022, IEEE Trans. Geosci. Remote Sens., 60, 1

Yokoya, N., Yairi, T., & Iwasaki, A. 2012, IEEE Trans. Geosci. Rem. Sens., 50, 528

Yokoya, N., Grohnfeldt, C., & Chanussot, J. 2017, IEEE Geosci. Rem. Sens. Mag., 5, 29

Zhang, Y. 2008, The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 37

<footnote n="10">¹⁰ https://jax.readthedocs.io/</footnote>
## Appendix A: Metrics definition

In this section, we define the scalar metrics calculated in Table 2. First, the Normalised Mean Squared Error (NMSE) is defined as such (where $Z_i$ is the value of Z at the voxel i):

$$NMS E(\hat{Z}, Z) = 10 \log_{10} \frac{\sum_{i=1}^{LMN} (\hat{Z}_i^2)}{\sum_{i=1}^{LMN} (\hat{Z}_i - Z_i)^2}, \quad (A.1)$$

By this definition, a higher NMSE implies a result of better quality. Then, the average Spectral Angle Mapper (aSAM) is defined as:

$$aSAM(\hat{Z}, Z) = \frac{1}{L} \sum_{i}^{MN} \arccos \left( \frac{\langle Z_i, \hat{Z}_i \rangle}{\|Z_i\|_2 \|\hat{Z}_i\|_2} \right), \quad (A.2)$$

where $Z_i$ is the spectra for the pixel i. A lower aSAM implies a result of better quality, and an ideal value is 0.

The average complementary Structural Similarity (SSIM) index, is:

$$acSSIM(\hat{Z}, Z) = 1 - \frac{1}{L} \sum_{k=1}^{L} SSIM(\hat{Z}_k, Z_k), \quad (A.3)$$

where SSIM is the Structural Similarity index, as defined by Wang et al. (2004). The aim of the acSSIM is to measure the capacity to reconstruct local image structure, luminance, and contrast. The lower the acSSIM, the better the result. An ideal value is 0.

Finally, the Relative Dimensionless Global Error (ERGAS), is defined as:

$$ERGAS(\hat{Z}, Z, d) = 100d \sqrt{\frac{1}{L} \sum_{k=1}^{L} \left( \frac{\|\hat{Z}_k - Z_k\|_F}{\mu_k \sqrt{MN}} \right)^2}, \quad (A.4)$$

where d is the ratio between the spatial resolution of X and the spatial resolution of Y, and $\mu_k$ is the mean of $\hat{Z}_k$. The lower the ERGAS, the better, and an ideal value is 0.

## Appendix B: Additional Figures

This section presents additional figures for the results of section 4. Fig. B.1 shows the obtained amplitude for the three potential regularization, for each of the four toy models. For these same options, Fig. B.2 shows the spectral angle maps. Fig. B.3 displays the relative error maps averaged per pixel.
<figure><img src="image_9.png" /><figcaption>**Fig. B.1.** Amplitude maps (i.e. sum of counts per pixel over every spectral channel) obtained by the three regularizers and compared to the ground truth.</figcaption></figure>

(a) Gaussian model between 0.5-1.4 keV, (b) Gaussian model between 6.2-6.9 keV, (c) Gaussian model with rebinning between 0.5-1.4 keV, and (d) Realistic model between 0.5-1.4 keV.
<figure><img src="image_10.png" /><figcaption>**Fig. B.2.** Maps of the SAM (Equation 23) for the four toy models, and the three regularizations. A higher SAM means a worse spectral reconstruction.</figcaption></figure>

(a) Gaussian model between 0.5-1.4 keV, (b) Gaussian model between 6.2-6.9 keV, (c) Gaussian model with rebinning between 0.5-1.4 keV, and (d) Realistic model between 0.5-1.4 keV.
<figure><img src="image_11.png" /><figcaption>**Fig. B.3.** Maps of the relative error: $Err(Z) = 100(Z - \hat{Z})/\hat{Z}$, for the four toy models, and the three regularizations. A higher absolute error means a worse spectral reconstruction.</figcaption></figure>

(a) Gaussian model between 0.5-1.4 keV, (b) Gaussian model between 6.2-6.9 keV, (c) Gaussian model with rebinning between 0.5-1.4 keV, and (d), Realistic model between 0.5-1.4 keV.