11institutetext: Department of Astronomy, School of Physical Sciences,
University of Science and Technology of China, Hefei, Anhui 230026, China.
e-mail<EMAIL_ADDRESS>22institutetext: CAS Key Laboratory for Research in
Galaxies and Cosmology, University of Science and Technology of China, Hefei,
Anhui 230026, China 33institutetext: School of Astronomy and Space Science,
University of Science and Technology of China, Hefei, Anhui 230026, China

# A hard spectrum diffuse $\gamma$-ray component associated with H ii gas in
the Galactic plane

Bing Liu 112233 Rui-zhi Yang 112233

We analyzed 12-year Fermi Large Area Telescope $\gamma$-ray data in the inner
Galaxy centered at ($l=30^{\circ}$, $b=0^{\circ}$) and ($l=330^{\circ}$,
$b=0^{\circ}$). We found significant hardening of the spectrum of the diffuse
$\gamma$-ray emission in these regions as previously reported. We further
deduced that the diffuse $\gamma$ rays can be divided into two components from
the likelihood analysis. One component is associated with the total gas column
density and reveals a soft spectrum, while the other is associated with the H
ii gas and presents a hard spectrum. Assuming the diffuse $\gamma$-ray
emissions are mainly produced through the interaction between cosmic rays
(CRs) and the ambient gas, these two components are produced by the CR
populations with spectral indices of $2.8$ (”soft”) and $2.3$ (”hard”),
respectively. We argue that the hard CR population may come from the vicinity
of the CR accelerators. The soft CR population has a similar spectral shape
and density as measured in the solar neighborhood, which implies a uniform CR
”sea” with a similar density and spectral shape in the Galaxy.

###### Key Words.:

$\gamma$ rays: diffuse, cosmic rays, interstellar medium

##  1 Introduction

Cosmic rays (CRs) are believed to be well mixed in the magnetic field in our
Galaxy and interact with the cold gas in the interstellar medium (ISM) to
produce diffuse $\gamma$-ray emissions. Thanks to the progress in the
$\gamma$-ray astronomy, such diffuse $\gamma$-ray emissions have been well
studied in recent years (Acero et al. 2016; Yang et al. 2016). Taking into
account the gas content derived from molecular and atomic line spectrometry,
one can obtain information on the CR distribution on the Galactic scale. By
assuming a cylindrically symmetric distribution, Yang et al. (2016) and Acero
et al. (2016) found that CRs are strongly inhomogeneous in our Galaxy. The
radial distribution of CRs peaks at the galactocentric distance of about 4
kpc, and the derived CR spectra in the inner Galaxy are significantly harder.

The spatial distribution of CRs has direct implications on the understanding
of their origin and propagation. Generally speaking, the diffusion process
will make the CR distribution smoother than the distribution of CR sources,
but it cannot erase all hints of the source distributions (see, e.g., Strong &
Moskalenko 1998). Thus, the sharp peak of the derived CR distributions
requires a significant modification to the ”standard” CR propagation model in
our Galaxy. As an example, Recchia et al. (2016) invoked the CR-driven wind in
the Galactic plane to explain the observed CR radial profile.

However, whether such inhomogeneity is a global effect or caused by individual
sources is still unclear. Aharonian et al. (2020) used the $\gamma$-ray
observations of individual molecular clouds to derive the CR density therein,
and unveiled a homogeneous “sea” of CRs with a constant density and spectral
shape, especially for the galactocentric distances exceeding 8 kpc, as well as
the Sagittarius B complex, in the region of the Galactic center (GC). However,
in the region with galactocentric distances between 4 and 6 kpcs, the derived
CR densities reveal significant deviation from the ”sea” and each other. One
possible explanation is that larger amounts of particle accelerators within
this region cause the enhancement of CR density in their vicinity. The
distributions of potential CR accelerators, such as supernova remnants (SNRs),
pulsars, and OB stars, do reveal strong galactocentric radius dependence
(Green 2015; Yusifov & Küçük 2004; Bronfman et al. 2000). Particularly, the
distribution of OB stars presents a similar peak at about 4 kpc (Bronfman et
al. 2000), which makes these massive star clusters an attractive source
population for addressing the higher CR density and harder CR spectra in this
region. Indeed, young massive star clusters (YMCs) were recently identified as
a new population of $\gamma$-ray sources and are believed to be an alternative
source of the Galactic CRs (Aharonian et al. 2019). In several such systems,
including the Cygnus cocoon (Ackermann et al. 2011; Aharonian et al. 2019),
Westerlund 1 (Abramowski et al. 2012), NGC 3603 (Yang & Aharonian 2017),
Westerlund 2 (Yang et al. 2018) and 30 Dor C (H.E.S.S. Collaboration et al.
2015), extended $\gamma$-ray emissions with hard spectra have been detected
from GeV to TeV. Besides this, most of the observed $\gamma$ rays reveal a
good spatial correlation with the H ii region, which is formed due to the
photoionization of these massive stars (Yang & Aharonian 2017; Aharonian et
al. 2019; Sun et al. 2020a, b). Due to strong extinction in the Galactic
plane, there may be numerous unknown YMCs exist in the Galaxy. Thus, if there
are additional contributions from the unknown YMCs, we expect the $\gamma$-ray
emission to present similar spatial distribution as the H ii gas.

To test such a hypothesis, we chose two regions centered at
$(l,b)=(30^{\circ},0^{\circ})$ (region I) and $(l,b)=(330^{\circ},0^{\circ})$
(region II) to investigate the diffuse $\gamma$-ray emissions as well as the
distribution of gases in different phases. We chose these regions because they
are near the tangent point of the 4-kpc ring, and the aforementioned CR excess
should have a strong impact on the total $\gamma$-ray fluxes in these
directions. In the following, we give a detailed analysis of the 12-year
Fermi-LAT data toward these regions and discuss the gas and CR content
therein. The paper is organized as follows. In Sect.2, we investigate the gas
distribution in this region. In Sect.3, we present details of the $\gamma$-ray
data analysis. In Sect.4, we discuss the implications of our analysis results.

##  2 Gas tracers

The 21 cm H i line and 2.6 mm CO line are widely used tracers for atomic
hydrogen and molecular hydrogen, respectively. In this work, we used the data
from the Galactic CO survey of Dame et al. (2001) with the CfA
1.2m-millimeter-wave Telescope, and the HI4PI Survey on HI gas (HI4PI
Collaboration et al. 2016). For the CO data, we used the standard assumption
of a linear relationship between the velocity-integrated CO intensity, $W_{\rm
CO}$, and the column density of molecular hydrogen, N(H2). The conversion
factor $X_{\rm CO}$ may be different in different parts of the Galaxy;
therefore, we used the values derived from Acero et al. (2016) for different
galactocentric distances, such as $X_{\rm CO}=0.5\times
10^{20}\leavevmode\nobreak\ \rm cm^{-2}(K\leavevmode\nobreak\ km/s)^{-1}$
below 1.5 kpc and $X_{\rm CO}=1.5\times 10^{20}\leavevmode\nobreak\ \rm
cm^{-2}(K\leavevmode\nobreak\ km/s)^{-1}$ above.

For the HI data, we use the following equation:

$N_{HI}(v,T_{s})=-log\left(1-\frac{T_{B}}{T_{s}-T_{bg}}\right)T_{s}C_{i}\Delta
v\ ,$ (1)

where $T_{bg}\approx 2.66$ K is the brightness temperature of the cosmic
microwave background radiation at 21cm, and $C_{i}=1.83\times 10^{18}\rm
cm^{2}$. In the case when $T_{B}>T_{s}-5\leavevmode\nobreak\ \rm K$, we
truncated $T_{B}$ to $T_{s}-5\leavevmode\nobreak\ \rm K$; $T_{s}$ was chosen
to be 150 K. The systematic uncertainties due to the different spin
temperatures are discussed in Ackermann et al. (2012) and Acero et al. (2016).
The effect, however, is quite small in most regions of the sky.

In order to investigate the galactocentric radial distribution of CRs, we need
to divide the gas distribution into different rings around the GC. For this
purpose, we used the following relation:

$V_{LSR}=R_{\odot}(\frac{V(R)}{R}-\frac{V_{\odot}}{R_{\odot}})sin(l)cos(b),$
(2)

where $V_{LSR}$ is the radial velocity with respect to the local standard of
rest, $R$ is the galactocentric distance, $V(R)$ is the Galactic rotational
curve, and $l$ and $b$ are the galactic coordinates in the line of sight
(LOS). We adopted the rotational curve parameterized in Clemens (1985);
$V_{\odot}$ and $R_{\odot}$ are fixed to 220 km/s and 8.5 kpc, respectively.
By applying this to both CO and HI data, we can transform the velocity
information into the galactocentric distance of the gas. According to Fig. 8
of Acero et al. (2016), we split the gas into galactocentric distance bins;
for example, below and above 4.5 kpc in the analysis below. It should be noted
that Eq.2 allows emission from the forbidden velocity zones in the CO and HI
data. The reason could be the non-circular motion of the gas. This component
contains only a small fraction of total gas; therefore, for simplicity, we
assigned it to the local rings.

For different reasons, the neural gas cannot always be traced by CO and H i
observations (Grenier et al. 2005). In such cases (e.g. in optically thick
clouds), the infrared emission from cold interstellar dust provides an
alternative and independent measurement of the gas column density. In the
analysis below, we also used the dust column map as gas template for cross-
checking. According to Eq. (4) of Planck Collaboration et al. (2011), the
relation between the dust opacity and the column density can be approximated
as

$\tau_{\rm M}(\lambda)=\left(\frac{\tau_{\rm D}(\lambda)}{N_{\rm
H}}\right)^{\rm dust}[N_{{\rm HI}}+2X_{\rm CO}W_{\rm CO}],$ (3)

where $\tau_{\rm M}$ is the dust opacity as a function of the wavelength
$\lambda$, $(\tau_{\rm D}/N_{\rm H})^{\rm dust}$ is the reference dust
emissivity measured in low-$N_{\rm H}$ regions, $W_{\rm CO}$ is the integrated
brightness temperature of the CO emission, and $X_{\rm CO}=N_{\rm
H_{2}}/W_{\rm CO}$ is the $\rm{H_{2}/CO}$ conversion factor. The substitution
of the latter into Eq. (3) gives

$N_{\rm H}=N_{\rm HI}+2N_{\rm H_{2}}=\tau_{\rm
M}(\lambda)\left[\left(\frac{\tau_{\rm D}(\lambda)}{N_{\rm H}}\right)^{\rm
dust}\right]^{-1}.$ (4)

Here for the dust emissivity at $353\leavevmode\nobreak\ \rm GHz$, we used
$(\tau_{\rm D}/N_{\rm H})^{\rm dust}_{353{\rm\leavevmode\nobreak\
GHz}}=(1.18\pm 0.17)\times 10^{-26}$ cm2 taken from Table 3 of Planck
Collaboration et al. (2011). The derived total gas column density
distributions in region I and region II are shown in the top panel of Fig.1.

We also used the dust data to derive the residual templates for the missing
”dark” gas component. We fit the dust opacity maps as a linear combination of
HI and CO maps, and, in this way, find the residual map to be the dark gas
template. Then, we iterated this fit by including the dark gas template until
convergence is achieved (Yang et al. 2016). This method is similar to the
derivation of the $E(B-V)_{res}$ templates used by the Fermi-LAT collaboration
(Acero et al. 2016), where instead of the extinction maps we used the dust
opacity maps. Indeed, $E(B-V)$ has a nearly perfect linear correlation with
the dust opacity, especially in higher column regions (Planck Collaboration et
al. 2014). Thus, our method and the method used in Acero et al. (2016) should
give similar results.

We used the free-free emission map obtained from the joint analysis of Planck,
WMAP, and 408 MHz observations (Planck Collaboration et al. 2016) to derive
the map of H ii column density. First, we converted the emission measure into
free-free intensity ($I_{\nu}$) by using the conversion factor at 353-GHz in
Table 1 of Finkbeiner (2003). Then, we applied Eq. (5) from Sodroski et al.
(1997),

$\displaystyle N_{{H\textsc{ii}}}=$ $\displaystyle 1.2\times 10^{15}\ {\rm
cm^{-2}}\left(\frac{T_{\rm e}}{1\ \rm K}\right)^{0.35}\left(\frac{\nu}{1\ \rm
GHz}\right)^{0.1}\left(\frac{n_{\rm e}}{1\ \rm cm^{-3}}\right)^{-1}$ (5)
$\displaystyle\times\frac{I_{\nu}}{1\ \rm Jy\ sr^{-1}},$

to convert the free-free intensity into column density in each pixel, with a
frequency at $\nu=\rm 353\ GHz$, and an electron temperature of $T_{e}=\rm
8000\ K$. This equation also shows that the H ii column density is inversely
proportional to the effective density of electrons $n_{\rm e}$. Here, we chose
$n_{\rm e}=2\ \rm cm^{-3}$ (Sodroski et al. 1997) to estimate the upper limit
of the H ii column density. We note that $n_{\rm e}$ varies in different
positions of the Galaxy; however, in this study, most of the H ii column
densities are from the inner Galaxy. As discussed in Sodroski et al. (1997),
the assumption of using $n_{\rm e}=2\ \rm cm^{-3}$ as the effective electron
density in the H ii regions in the inner Galaxy produce the same gas-to-dust
mass ratio for the H ii component as the mean values of the gas-to-dust mass
ratio for the H2 and HI components within the same ranges of radial distance.
Thus, in the following study, we fixed $n_{\rm e}=2\ \rm cm^{-3}$ to derive
the column density of the H ii component. The derived H ii column density maps
of both regions are shown in the bottom panel of Fig. 1.

Figure 1: Dust and H ii templates used in Fermi-LAT data analysis. Top panel:
Derived total gas column density maps (in units of ${\rm cm}^{-2}$) from
Planck dust opacity for region I (left) and region II (right). Bottom panel:
Derived H ii gas density maps (in units of ${\rm cm}^{-2}$) from the joint
analysis of Planck, WMAP, and 408 MHz observations for region I (left) and
region II (right).

Figure 2: Gas column density maps (in units of ${\rm cm}^{-2}$) below 4.5 kpc
(top panel) and above 4.5 kpc (bottom panel) for region I (left) and region II
(right), which are derived from CO and HI observation data (Dame et al. 2001;
HI4PI Collaboration et al. 2016) using method described in Sect.2.

Figure 3:  Dark gas templates for region I (left) and region II (right), which
are derived from dust data using the method described in Sect.2.

##  3 Fermi-LAT data analysis

We collected the Fermi-LAT Pass 8 data from August 4, 2008 (MET 239557417)
until October 28, 2020 (MET 625621368), and used the Fermitools from the Conda
distribution111https://github.com/fermi-lat/Fermitools-conda/ together with
the latest version of the instrument response functions (IRFs) P8R3_SOURCE_V3
for the $\gamma$-ray emission analysis. Two $20\hbox{${}^{\circ}$}\times
20\hbox{${}^{\circ}$}$ square regions centered at the position of
($l=30^{\circ}$, $b=0^{\circ}$) and ($l=330^{\circ}$, $b=0^{\circ}$) were
chosen as two regions of interest (ROI). Here, we selected the ”source” class
events, then applied the recommended data cut expression
$\rm(DATA\\_QUAL>0)\&\&(LAT\\_CONFIG==1)$ to exclude time periods when some
spacecraft events affected the data quality. Moreover, to reduce the
background contamination from the Earth’s albedo, only the events with zenith
angles under $90\hbox{${}^{\circ}$}$ were included in the analysis.

Figure 4:  Observed counts map of $\gamma$ rays in the energy range of 0.2–200
GeV for region I (left) and region II (right).

For each region, the background model includes the sources in the Fermi-LAT
ten-year Source Catalog (4FGL-DR2, Ballet et al. 2020) within the ROI enlarged
by 10∘, and the normalizations and spectral parameters were left free for
sources with significance $>5\sigma$ within the ROI.

###  3.1 Diffuse $\gamma$-ray emission

As for the diffuse background components, we first used the Fermi-LAT Galactic
diffuse background model (gll_iem_v07.fits), hereafter referred to as the
”Fermi” model. Since we are interested in the diffuse $\gamma$-ray emissions,
we also built our own Galactic diffuse background models for cross-checking.
The diffuse $\gamma$-ray emission mainly comes from the pion decay process
induced by the inelastic collision between CR protons and ambient gas, as well
as the inverse Compton (IC) scattering of CR electrons in the interstellar
radiation fields (ISRFs). We calculated the IC component using
GALPROP222http://galprop.stanford.edu/webrun/ (Vladimirov et al. 2011), which
uses information regarding CR electrons and ISRFs. For the GALPROP webrun we
used the parameter set ${}^{S}Y^{Z}4^{R}30^{T}150^{C}2$ in Ackermann et al.
(2012) as a fiducial model. We also used the 128 parameter sets in Ackermann
et al. (2012) to generate the corresponding IC emission templates; and for
each region, we chose the best template by comparing the maximum likelihood
achieved by each likelihood fitting process.

For the pion decay component, as a zero-order approximation, the CRs are
homogeneously distributed in the CR sea in our Galaxy. Thus, we can use the
gas distribution map to model the pion decay component of the diffuse
$\gamma$-ray emission. As mentioned above, to avoid the uncertainties of dark
gas and H i spin temperatures, we used the gas column density maps derived
from the Planck dust opacity. Instead of the template gll_iem_v07.fits in the
Fermi models, models using total gas distribution maps derived from Planck
dust opacity ( as shown in the top panels of Fig. 1), and IC emission
templates are referred to as ”dust” models.

However, since the metallicity varies with galactocentric distance, the dust-
to-gas ratio may also vary in the Galaxy. Furthermore, as derived in Acero et
al. (2016), the spectra and density of CRs vary in galactocentric distance.
According to the results in Acero et al. (2016), both the spectra and density
of CRs and the $X_{CO}$ factor vary significantly below and above 4.5 kpc.
Thus, we split the gas into to rings with galactocentric radius below and
above 4.5 kpc and include both of them in the pion decay component, and the
gas column density distributions within different galactocentric distances are
shown in Fig.2. In this model, as mentioned above, a dark gas component must
be invoked to account for the gas that cannot be traced by CO and
$21\leavevmode\nobreak\ \rm cm$ observations. Thus, we also added a dark gas
component as described in the last section. Models that replace the total gas
distribution map with split gas distribution maps and the dark gas component
maps (as shown in Fig.3) are referred to as ”gas” models. Finally, as we
investigated the possible contribution of an H ii component in the diffuse
$\gamma$-ray emissions, we also added the H ii column density map (as shown in
the bottom panel of Fig.1) as a diffuse template in the following analysis. In
the following analysis, we used a log-parabola function ($dN/dE=N_{0}(E/E_{\rm
b})^{-\alpha-\beta log(E/E_{\rm b})}$) as the energy spectrum of the gas
template and left both the normalization and spectral parameters free in the
likelihood fitting below.

###  3.2 Results

Using photons of energies within 0.2–200 GeV, we performed the standard
likelihood analysis of Fermi-LAT data with the 4FGL-DR2 sources and diffuse
components mentioned above. To check the possible contribution from the H ii
component, we performed the analysis with and without a template derived from
H ii column density, respectively. To compare the goodness of fit in these two
models, we also calculated the Akaike information criterion (AIC) value for
each model. AIC was defined as $AIC=-2\rm log(likelihood)+2k$, where k is the
number of free parameters in the model. The log(likelihood) and corresponding
AIC value for each fit are listed in Table.1. In addition, we compared the
overall maximum likelihood of the model with H ii template (alternative
hypothesis, $L$ ) and that of the model without (null hypothesis, $L_{0}$) to
obtain the test statistics (TS) of the H ii component which is defined as
$-2{\rm log}(L_{0}/L)$ following Lande et al. (2012). Then, the corresponding
significance ($\sigma_{{H\textsc{ii}}}$) of the H ii component can be
estimated as $\sim\sqrt{\rm TS_{{H\textsc{ii}}}}$. The TS value and the
significance of the additional H ii component for each model are also listed
in Table.1.

We found that in both regions the inclusion of the H ii templates for dust and
gas models improves the fit dramatically, and the results are even improved
for the Fermi models. However, the overall results reveal that the gas models
cannot fit the data as well as the Fermi or dust models. Furthermore, by
subtracting the fit model map calculated using gtmodel from the observed
$\gamma$-ray counts map (as shown in Fig.4), we generated the residual maps of
these different models. Then, we divided the residual counts by the square
root of the observed $\gamma$-ray counts in each pixel, which should be equal
to the significance (in $\sigma$) of each pixel approximately. We used these
residual significance (signal-to-noise, S/N) maps to compare the fitting
quality of these models. The improvement can be also seen from the comparisons
of the results with and without the H ii templates, which are shown in Fig.10
and Fig.11. In addition, we used a Gaussian function to fit the distribution
of the S/N ratios of each model to check whether the model is a statistically
good description of the data. As shown in the bottom panel of Fig.10 and
Fig.11, each distribution can be well described by the Gaussian function with
a mean of $\sim$0.0 and standard deviation of $\sim$1.0, which indicates an
unbiased significance distribution of the background.

Table 1: Fitting results for different models of region I and region II Region | Model |  free parameters (k) |  log(likelihood) |  AIC | ${\rm TS}_{\rm{H\textsc{ii}}}$ | $\sigma_{{H\textsc{ii}}}$
---|---|---|---|---|---|---
I | Fermi |  506 |  15876971.5 |  -31752931.0 |  |
Fermi+H ii |  509 |  15877085.4 |  -31753152.8 |  227.8 | 15.1
dust |  509 |  15873005.8 |  -31744993.6 |  |
dust+H ii |  512 |  15876432.2 |  -31751840.4 |  6852.8 | 82.8
gas |  515 |  15870347.9 |  -31739665.8 |  |
gas+H ii |  518 |  15872922.1 |  -31744808.2 |  5148.4 | 71.8
II | Fermi |  592 |  8009736.8 |  -16018289.6 |  |
Fermi+H ii |  595 |  8009738.5 |  -16018287.0 |  3.4 | 1.8
dust |  595 |  8006053.7 |  -16010917.4 |  |
dust+H ii |  598 |  8008499.1 |  -16015802.2 |  4890.8 | 69.9
gas |  601 |  8005227.7 |  -16009253.4 |  |
gas+H ii |  604 |  8006935.4 |  -16012662.8 |  3415.4 | 58.4

Figure 5: SEDs of diffuse $\gamma$-ray emissions associated with different
components for region I (left panel) and region II (right panel). The black
and red points show the $\gamma$ rays associated with total gas (the dust
template) for both dust models and dust+H ii models, respectively. The blue
points represent the $\gamma$ rays related to H ii gas (the H ii template)
obtained for dust+H ii models.

Figure 6: Same as Fig.5, but the SEDs have been normalized to the $\gamma$-ray
emissivity per H atom, which is proportional to the CR density. Also plotted
as black curves are the predicted $\gamma$-ray emissivities for the local CR
spectra (LIS).

###  3.3 Spectral analysis

The TS values of the H ii templates in the Fermi model are much smaller than
those in the other two models. This is not unexpected, since in the standard
Fermi diffuse background model the artificial ”patches” are added to absorb
the residual in the Galactic plane
333https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html, and
the standard diffuse model is not suggested to study the diffuse emission
regarding the interstellar medium. And the IC contribution is also fixed.
Thus, in the spectral analysis, we used only the gas and dust templates
described above. We then divided the energy range 0.2 GeV - 200 GeV into nine
logarithmically spaced energy bins and extracted the spectral energy
distribution (SED) of diffuse $\gamma$ rays in these regions via the maximum
likelihood analysis in each energy bin. The uncertainties include 68%
statistical errors for the energy flux and systematic errors due to the
uncertainties in LAT effective areas. The results for the dust model and
dust+H ii models, as well as gas model and gas+H ii models, are shown in
Fig.5, Fig.6, Fig.7, and Fig.8, respectively.

Figure 7: $\gamma$-ray emissivities per H atom of diffuse $\gamma$-ray
emission associated with the two gas rings in gas models for region I (left
panel) and region II (right panel), respectively.

Figure 8:  $\gamma$-ray emissivities per H atom of diffuse $\gamma$-ray
emission associated with the two gas rings in gas +H ii models for region I
(left panel) and region II (right panel), respectively.

In Fig.6, Fig.7, and Fig.8, we have divided the SEDs with their corresponding
gas column densities derived from Sect.2 to obtain the $\gamma$-ray
emissivities per H atom, which is proportional to the CR density. The derived
emissivities are compared with the one predicted by using the local
interstellar spectrum (LIS) of CRs (black line, Casandjian 2015). We found
that in the dust model without including H ii templates, the derived CR
spectrum from the dust model of each region is significantly harder than the
LIS, which is consistent with the results for the inner Galaxy in Yang et al.
(2016) and Acero et al. (2016). However, for the dust+H ii models, the derived
$\gamma$-ray emissivities associated with the dust templates are in good
agreement with the local emissivities, which imply similar CR spectra in these
regions and the solar neighborhood. Meanwhile, the derived emissivities
associated with the H ii component are much harder, with an index of about
$-2.1$.

For the gas model without the H ii templates, the derived $\gamma$-ray
emissivities for the inner rings (below 4.5 kpc) are significantly higher than
the LIS value and the outer rings (above 4.5 kpc), while the spectra in the
inner rings is significantly harder. These results are also in good agreement
with the results of Yang et al. (2016) and Acero et al. (2016). When including
H ii templates, for region I the $\gamma$-ray emissivities for both inner and
outer gas rings are consistent with the LIS value in both the spectra and
normalizations. For region II, the $\gamma$-ray emissivities’ spectral indices
in both rings are consistent with the LIS value, while the normalization in
the inner ring is higher. The $\gamma$-ray emissivities from H ii templates
are also significantly harder in both regions and are consistent with the
results derived from dust+H iimodels in both spectra and normalizations.

It should be noted that the normalization of the derived emissivities
associated with H ii component is uncertain. This is because the ionized gas
only contributes a small portion to the total gas column in our Galaxy.
Moreover, it is possible that the CRs associated with ionized gas also
illuminate other gas components. In Fig.6 and Fig.8, we assumed only H ii gas
are responsible for the $\gamma$ rays associated with the H ii templates.
Thus, the emissivities derived here can be regarded as the upper limits.

##  4 Discussion

In this work, we found a hard $\gamma$-ray component associated H ii in two
regions in the Galactic plane. The derived intensity and spectrum of such a
$\gamma$-ray component are stable when adopting different gas tracers and IC
models. Including such a component will also alter the $\gamma$-ray
emissivities associated with the gas in the Galactic plane. This can be seen
in both the dust and gas models in our analysis. In the dust model, the
inclusion of the H ii component will deduce a similar $\gamma$-ray
emissivities in the whole line of sight as the local value, which implies that
the average CR density and spectra in the line of sight are similar to the
local CR density and spectrum. Such a conclusion is different from those in
Acero et al. (2016) and Yang et al. (2016). Moreover, in the gas model, the
gases are split into two rings with galactocentric radii below and above 4.5
kpc. Without the H ii component, we reproduced the results in Acero et al.
(2016) and Yang et al. (2016); thus, the $\gamma$-ray emissivities per H atom
are higher in the inner ring and reveal a significantly harder spectrum.
Including the H ii template, we found the $\gamma$-ray emissivities in both
rings become consistent with the local value for region I. For region II, the
spectra are similar to the local value, while for the inner rings the
normalizations are still significantly higher. However, we note that the
normalization of $\gamma$-ray emissivities per H atom is inversely
proportional to the $X_{co}$ conversion factor, which is fixed in our
analysis. Indeed, $X_{co}$ depends on the metallicity and can be different in
different parts of the Galaxy (Acero et al. 2016).

Figure 9:  $\gamma$-ray flux associated with H ii gas obtained for dust+H ii
models (blue) and gas+H ii models (red) for region I (left panel) and region
II (right panel), respectively. Also plotted are the possible contributions
from IC and pion-decay processes. The details can be found in Sect.4 .

To conclude, we found a robust hard $\gamma$-ray emission component associated
with the H ii gas, and the $\gamma$-ray emissions associated with other gas
reveal a significantly softer spectrum with a similar spectral shape. In the
following discussion, we label the $\gamma$-ray emission and corresponding CR
population related to the dust and gas templates and H ii templates as ”soft”
and ”hard” components, respectively. As for the hard component, the absolute
normalization can not be determined decisively, as we discussed above. By
using the upper limit of the emissivities in the most extreme case, the hard
component is only produced in the H ii region, the derived CR density can be
100 times higher than the LIS value (see Fig.6 and Fig.8). Such enhancement of
CR densities is expected near CR sources such as SNRs and YMCs. The dimension
of the regions with enhanced CR density can be of more than 100 pc (Aharonian
et al. 2019). Thus in this case, the hard $\gamma$-ray component associated
with H ii gas can be interpreted as the enhanced $\gamma$-ray emission in the
vicinity of sources. To illustrate such a scenario, we plotted in Figs.9 the
$\gamma$-ray spectrum assuming a power-law parent proton spectrum with an
index of $-2.2$, which is similar to the CR spectrum derived in the vicinity
of YMCs (Aharonian et al. 2019). On the other hand, the hard $\gamma$-ray
emission can also have a leptonic origin, especially considering the enhanced
radiation fields in these H ii regions. The interstellar radiation fields
(ISRFs), which provide the low-energy photon targets for the IC process, can
be divided into three components in our Galaxy (Popescu et al. 2017),
including the optical and UV fields from the starlight, the infrared fields
from the dust emission, and the cosmic microwave backgrounds. In the H ii
regions, the ionizing massive stars will increase the optical and UV fields
significantly and thus produce additional IC emissions. We calculated such a
component assuming the electron spectrum is the same as the electron spectrum
measured in the solar neighborhood and from radio observations (Strong et al.
2011), and the optical and UV components are described by a grey-body spectrum
with a temperature of $5000\leavevmode\nobreak\ \rm K$. The derived IC
$\gamma$-ray fluxes are also shown in Fig.9. In the calculated $\gamma$-ray
flux from both pion decay and the IC component, the total normalization is
rather arbitrary, since both the average gas density and the energy densities
of the enhanced optical and UV radiation fields are unknown. But the derived
spectral shapes already reveal significant differences. The IC components, due
to the Klein-Ninisha effects, are softened above dozens of GeVs. The current
data cannot rule out the leptonic origin, but further observations above 100
GeV may distinguish two such possibilities. If there is another hard electron
spectrum in these H ii regions, the IC component can be harder than the curve
calculated here, but the Klein-Nishina effects for the optical and UV
components will inevitably play a role above several hundred GeV in the
$\gamma$-ray spectrum (Popescu et al. 2017). Thus, the observations above
several hundred GeV may distinguish such two scenarios.

Interestingly, a two-component CR model was used in Yang & Aharonian (2019) to
explain the hardening of CR spectrum above 200 GV, as well as the increasing
positron/electron and antiproton/proton ratios. In such a scenario, CRs are
injected from two populations of sources with different spectral properties
(more specifically, a 1st source population with soft spectra and a 2nd source
population with hard spectra). For the second source population, CRs
accumulate the grammage of about $1\leavevmode\nobreak\ \rm g/cm^{2}$ in the
vicinity of the sources. Accordingly, they should be effective $\gamma$-ray
emitters. Yang & Aharonian (2019) calculated diffuse $\gamma$-ray emissions
associated with the second (hard) source population and found it can dominate
over the first (soft) source population above about 100 GeV, which is in
agreement with our measurements shown in Fig.5. As mentioned in Yang &
Aharonian (2019), one possible candidate for the second source population can
be YMC. Indeed, the OB star distributions do have a strong peak at about 4 kpc
away from the GC (Bronfman et al. 2000). Therefore, we expect a strong
enhancement of the hard $\gamma$-ray emissions associated with these sources
in this region, which can naturally explain the high energy density and hard
spectrum of CRs derived from the 4-6 kpc rings in our Galaxy (Acero et al.
2016; Yang et al. 2016).

On the other hand, the similar density and spectral shape of the soft CR
component with the LIS reveal the possible CR sea all over the Galaxy. Such a
possibility is also unveiled in Aharonian et al. (2020), in which the
$\gamma$-ray observations on individual molecular clouds are used to measure
the CR density in different positions of our Galaxy. Such a uniform CR
distribution is not trivial, since if the CR propagation is dominated by
diffusion, the spatial distribution of the sources should always leave a trace
in the spatial distribution of CRs (Strong & Moskalenko 1998). However, all
the promising candidates of CR accelerators, such as SNRs (Green 2015),
pulsars (Yusifov & Küçük 2004), and OB stars (Bronfman et al. 2000), reveal a
significant inhomogeneous distribution in our Galaxy. Thus, a larger CR halo
or strong re-acceleration may be needed to explain the uniform distribution of
the CRs, which will lead to a significant modification to the current paradigm
of CR propagation in our Galaxy.

Finally, we note that although the fitting of $\gamma$-ray data is
significantly improved by introducing the H ii gas component, there are still
strong residuals (see Fig.10 and Fig.11). This is not unexpected if we assume
the hard $\gamma$-ray component is from the vicinity of CR sources. Since the
CR distribution should be significantly inhomogeneous near the source, we
expect no perfect spatial correlation between the gas and the $\gamma$-ray
emissions. In this regard, a detailed analysis in each compact H ii region may
be helpful. These regions are often extremely crowded and different sources
can be confused and overlap in the line of sight (Liu et al. 2019). However,
as shown in Aharonian et al. (2019) these sources can also be PeVatrons. If
this is the case, the corresponding $\gamma$-ray emission should be higher
than 10 TeV. In this energy band, the background due to the diffuse
$\gamma$-ray emission associated with soft CR components and individual
$\gamma$-ray sources such as pulsars, SNRs, and pulsar wind nebulae can be
strongly suppressed. The extensive shower arrays such as LHAASO (Sciascio
2016) and HAWC (Abeysekara et al. 2013), as well as the Cerenkov telescope
arrays such as H.E.S.S., MAGIC, VERITAS, and CTA would be ideal instruments to
explore these regions.

## Acknowledgements

Bing Liu is supported by the NSFC under grant 12103049. Ruizhi Yang is
supported by the NSFC under grants 11421303, 12041305 and the national youth
thousand talents program in China.

## References

  * Abeysekara et al. (2013) Abeysekara, A. U., Alfaro, R., Alvarez, C., et al. 2013, Astroparticle Physics, 50, 26
  * Abramowski et al. (2012) Abramowski, A., Acero, F., Aharonian, F., et al. 2012, A&A, 537, A114
  * Acero et al. (2016) Acero, F., Ackermann, M., Ajello, M., et al. 2016, ApJS, 223, 26
  * Ackermann et al. (2011) Ackermann, M., Ajello, M., Allafort, A., et al. 2011, Science, 334, 1103
  * Ackermann et al. (2012) Ackermann, M., Ajello, M., Atwood, W. B., et al. 2012, ApJ, 750, 3
  * Aharonian et al. (2020) Aharonian, F., Peron, G., Yang, R., Casanova, S., & Zanin, R. 2020, Phys. Rev. D, 101, 083018
  * Aharonian et al. (2019) Aharonian, F., Yang, R., & de Oña Wilhelmi, E. 2019, Nature Astronomy, 3, 561
  * Ballet et al. (2020) Ballet, J., Burnett, T. H., Digel, S. W., & Lott, B. 2020, arXiv e-prints, arXiv:2005.11208
  * Bronfman et al. (2000) Bronfman, L., Casassus, S., May, J., & Nyman, L. Å. 2000, A&A, 358, 521
  * Casandjian (2015) Casandjian, J.-M. 2015, arXiv e-prints, arXiv:1502.07210
  * Clemens (1985) Clemens, D. P. 1985, ApJ, 295, 422
  * Dame et al. (2001) Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792
  * Finkbeiner (2003) Finkbeiner, D. P. 2003, ApJS, 146, 407
  * Green (2015) Green, D. A. 2015, MNRAS, 454, 1517
  * Grenier et al. (2005) Grenier, I. A., Casandjian, J.-M., & Terrier, R. 2005, Science, 307, 1292
  * H.E.S.S. Collaboration et al. (2015) H.E.S.S. Collaboration, Abramowski, A., Aharonian, F., et al. 2015, Science, 347, 406
  * HI4PI Collaboration et al. (2016) HI4PI Collaboration, Ben Bekhti, N., Flöer, L., et al. 2016, A&A, 594, A116
  * Lande et al. (2012) Lande, J., Ackermann, M., Allafort, A., et al. 2012, ApJ, 756, 5
  * Liu et al. (2019) Liu, B., Yang, R.-z., Sun, X.-n., Aharonian, F., & Chen, Y. 2019, ApJ, 881, 94
  * Planck Collaboration et al. (2014) Planck Collaboration, Abergel, A., Ade, P. A. R., et al. 2014, A&A, 571, A11
  * Planck Collaboration et al. (2016) Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016, A&A, 594, A10
  * Planck Collaboration et al. (2011) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2011, A&A, 536, A19
  * Popescu et al. (2017) Popescu, C. C., Yang, R., Tuffs, R. J., et al. 2017, MNRAS, 470, 2539
  * Recchia et al. (2016) Recchia, S., Blasi, P., & Morlino, G. 2016, MNRAS, 462, L88
  * Sciascio (2016) Sciascio, G. 2016, Nuclear and Particle Physics Proceedings, 279-281, 166 , proceedings of the 9th Cosmic Ray International Seminar
  * Sodroski et al. (1997) Sodroski, T. J., Odegard, N., Arendt, R. G., et al. 1997, ApJ, 480, 173
  * Strong & Moskalenko (1998) Strong, A. W. & Moskalenko, I. V. 1998, ApJ, 509, 212
  * Strong et al. (2011) Strong, A. W., Orlando, E., & Jaffe, T. R. 2011, A&A, 534, A54
  * Sun et al. (2020a) Sun, X.-N., Yang, R.-Z., Liang, Y.-F., et al. 2020a, A&A, 639, A80
  * Sun et al. (2020b) Sun, X.-N., Yang, R.-Z., & Wang, X.-Y. 2020b, MNRAS, 494, 3405
  * Vladimirov et al. (2011) Vladimirov, A. E., Digel, S. W., Jóhannesson, G., et al. 2011, Computer Physics Communications, 182, 1156
  * Yang & Aharonian (2019) Yang, R. & Aharonian, F. 2019, Phys. Rev. D, 100, 063020
  * Yang et al. (2016) Yang, R., Aharonian, F., & Evoli, C. 2016, Phys. Rev. D, 93, 123007
  * Yang & Aharonian (2017) Yang, R.-z. & Aharonian, F. 2017, A&A, 600, A107
  * Yang et al. (2018) Yang, R.-z., de Oña Wilhelmi, E., & Aharonian, F. 2018, A&A, 611, A77
  * Yusifov & Küçük (2004) Yusifov, I. & Küçük, I. 2004, A&A, 422, 545

##  Appendix A Residual significance (signal-to-noise, S/N) maps for regions I
and II

Figure 10:  Comparison of the residual significance (signal-to-noise, S/N)
maps in the energy range of 0.2–200 GeV for models without an H ii component
(top) and with an H ii component (middle) of region I. The The left panels
show the results of Fermi and Fermi$+$H ii models, the middle panels show the
results of dust and dust$+$H ii models, and the right panels illustrate the
results of gas and gas$+$H ii models. The color scales are from $-5\sigma$ to
$5\sigma$. The bottom panel show the distribution of significance (bin size is
0.05 $\sigma$ of each model without or without the HII component, in which
each solid line illustrates the Gaussian function fitted to the corresponding
distribution and $\mu$ and $\sigma$ represent the mean and the standard
deviation of each fit.

Figure 11:  Same as Fig. 10, but for region II.


# Multi-Linear Kernel Regression
and Imputation in Data Manifolds

Duc Thien Nguyen and Konstantinos Slavakis1 1D. T. Nguyen and K. Slavakis are
with the Department of Information and Communications Engineering, Tokyo
Institute of Technology, Tokyo, Japan (e-mails: {nguyen.t.au,
slavakis.k.aa}@m.titech.ac.jp).

###### Abstract

This paper introduces an efficient multi-linear nonparametric (kernel-based)
approximation framework for data regression and imputation, and its
application to dynamic magnetic-resonance imaging (dMRI). Data features are
assumed to reside in or close to a smooth manifold embedded in a reproducing
kernel Hilbert space. Landmark points are identified to describe concisely the
point cloud of features by linear approximating patches which mimic the
concept of tangent spaces to smooth manifolds. The multi-linear model effects
dimensionality reduction, enables efficient computations, and extracts data
patterns and their geometry without any training data or additional
information. Numerical tests on dMRI data under severe under-sampling
demonstrate remarkable improvements in efficiency and accuracy of the proposed
approach over its predecessors, popular data modeling methods, as well as
recent tensor-based and deep-image-prior schemes.

###### Index Terms:

Imputation, kernel, manifold, MRI, regression.

##  I Introduction

Dynamic magnetic resonance imaging (dMRI) is a popular non-invasive imaging
modality for observing body organ movement, with rich potential in cardiac and
neurological diagnosis [1]. DMRI stands out as an application domain where
regression grapples with all of the archetypal data-analytic bottlenecks:
large dimensionality due to the image data, dynamic data patterns due to
dMRI’s time component, missing data due to under-sampling, and strong but
unknown spatio-temporal correlations since, often, dMRI monitors structured
movement; e.g., a beating heart [2].

It comes, thus, as no surprise that numerous data-modeling approaches have
been proposed for regression and imputation on dMRI data: compressed sensing
[3, 4, 5, 6], low-rank models [7, 8, 9], and learning strategies based on
dictionaries [10, 11, 12, 13, 14], transforms [15], manifolds [16, 17, 18, 19,
20], kernels [21, 22, 23], and tensors [24, 25, 26, 27]. Recent efforts on
imputation-by-regression revolve also around deep-learning (DeepL) approaches
[28, 29, 30, 31], which rely on time-consuming processes to learn from
training data prior to reconstructing test data. Notwithstanding, concerns
were raised in [32] via numerical tests which highlighted potential
instabilities of DeepL approaches. Motivated by deep image priors (DIP), DeepL
networks have also been used as implicit structural priors in regression for
dMRI to avoid the use of training data and potential over-fitting issues [33,
34].

Departing from all of the previous schemes, this paper offers an extension of
the novel nonparametric data-modeling approach of [35, 36], coined hereafter
multi-linear kernel regression and imputation in data manifolds (MultiL-KRIM).
MultiL-KRIM introduces a multi-linear matrix decomposition in data modeling to
offer two-pronged innovation over its bi-linear predecessors KRIM [36] and
BiLMDM [35]:  (i) unlike KRIM and BiLMDM, where dimensionality-reduction pre-
steps are disassociated from the regression task, MultiL-KRIM connects
dimensionality reduction directly with the regression task by enabling its
inverse-problem solution to identify the “optimal” dimensionality-reduced
rendition of a kernel matrix which contributes in minimum data-recovery error;
and  (ii) it exploits its multiple matrix factors to promote efficient
computations in its inverse-problem algorithmic solution.

MultiL-KRIM retains also the attributes which differentiate KRIM and BiLMDM
from state-of-the-art modeling approaches: unlike low-rank [7, 8, 9, 15],
dictionary-learning [10, 11, 12, 13, 14] and tensor [24, 25, 26, 27] models,
which promote a “blind decomposition” of the data matrix/tensor, MultiL-KRIM
incorporates the underlying data-manifold geometry directly into data
representations, but not via graph-Laplacian-matrix regularizers which are
widely used in manifold-learning approaches [16, 17, 18]. MultiL-KRIM adopts
instead a “collaborative-filtering” modeling approach to identify “optimal”
and manifold-cognizant combinations of the observed data features for
regression and imputation. MultiL-KRIM needs no training data to operate,
builds a nonparametric regression estimate to reduce the dependence of its
modeling assumptions on the probability distribution of the data [37], and
offers an explainable learning paradigm via simple geometric arguments, unlike
DeepL schemes, which are based, in general, on perplexed and cascading non-
linear function layers.

Numerical tests on synthetic dMRI data show that MultiL-KRIM outperforms
several state-of-the-art methods, including the total-variation tensor-based
scheme [27] and the deep-image-prior (DIP) propelled [34], while matching at
the same time the recovery-error performance of its predecessor KRIM [35, 36],
but with computational times which are lower than and can reach down to one-
third of those of KRIM.

##  II Data Collection and Formation in dMRI

(a) The (k,t)-space

(b) Cartesian sampling

(c) Radial sampling

(d) The image domain

Figure 1: (a) Complex-valued $N_{\text{f}}\times N_{\text{p}}\times
N_{\text{fr}}$ (k,t)-space dMRI data. The marked
$N_{\text{f}}\times\upsilon\times N_{\text{fr}}$ box corresponds to the
location of the faithful data, also known as “navigator/pilot” data in dMRI.
(b) 1D Cartesian and (c) radial sampling trajectories in k-space. (d) Complex-
valued $N_{\text{f}}\times N_{\text{p}}\times N_{\text{fr}}$ image-domain
data.

DMRI data take the form of a complex-valued ($\mathbb{C}$ is the set of all
complex numbers) three-way tensor $\mathbfscr{Y}$ defined on the
$(N_{\text{f}}\times N_{\text{p}}\times N_{\text{fr}})$-sized “(k,t)-space”
(Fig. 1a), with $N_{\text{f}},N_{\text{p}},N_{\text{fr}}\in\mathbb{N}_{*}$
($\mathbb{N}_{*}$ is the set of all positive integers). In the seldom case
where the k-space is “fully sampled,” the “slice/frame” $\mathbfcal{Y}_{t}$ of
$\mathbfscr{Y}$ ($t$ denotes discrete time with
$t\in\\{1,\ldots,N_{\text{fr}}\\}$) collects the $(N_{\text{f}}\times
N_{\text{p}})$-sized “k-space” measurements at $t$. In practice, it is often
the case that the k-space data is heavily under-sampled due to physical
limitations [2]. Popular sampling strategies are 1-D Cartesian (Fig. 1b) and
radial (Fig. 1c) sampling. Sampling is denoted by the entry-wise sampling
mapping
$\mathscr{S}(\cdot):\mathbfcal{Y}_{t}\mapsto\mathscr{S}(\mathbfcal{Y}_{t})$,
which nullifies the entry of $\mathbfcal{Y}_{t}$ when that entry is missing,
while retains the entry when that entry is successfully collected. Integers
$N_{\text{f}}$ and $N_{\text{p}}$ denote the numbers of frequency- and phase-
encoding lines, respectively [1], while $N_{\text{k}}\coloneqq
N_{\text{f}}N_{\text{p}}$ represents the number of entries of each k-space
frame. Typically, k-space is considered as the “frequency domain” of the
“image-data” domain (Fig. 1d), so that
$\mathbfcal{X}_{t}\coloneqq\mathscr{F}^{-1}(\mathbfcal{Y}_{t})$, where
$\mathscr{F}^{-1}(\cdot)$ is the 2D inverse DFT [1]. For convenience, the
columns of $\mathbfcal{Y}_{t}$ are stacked one below the other to create a
single $N_{\text{k}}\times 1$ vector
$\mathbf{y}_{t}\coloneqq\text{vec}(\mathbfcal{Y}_{t})$ and
$\mathbf{Y}\coloneqq[\mathbf{y}_{1},\ldots,\mathbf{y}_{N_{\text{fr}}}]\in\mathbb{C}^{N_{\text{k}}\times
N_{\text{fr}}}$.

It is often the case in an imputation framework for a subset of the sampled
data to be considered faithful. In the present context, the “low-frequency”
region, i.e., the central region of k-space, constitutes the faithful data
from which geometric information will be extracted [35, 36]. These data will
be called “navigator/pilot data.” The navigator data of the frame
$\mathbfcal{Y}_{t}$ (Fig. 1a) are gathered into the $\nu\times 1$ vector
$\check{\mathbf{y}}^{\text{f}}_{t}$, where $\nu\coloneqq\upsilon
N_{\text{f}}$. All of these vectors are finally stacked into the columns of
$\check{\mathbf{Y}}_{\text{f}}\coloneqq[\check{\mathbf{y}}_{1}^{\text{f}},\ldots,\check{\mathbf{y}}_{N_{\text{fr}}}^{\text{f}}]\in\mathbb{C}^{\nu\times
N_{\text{fr}}}$.

To offer an algorithmic scheme with manageable computational complexity as the
cardinality of the point-cloud
$\\{\check{\mathbf{y}}^{\text{f}}_{t}\\}_{t=1}^{N_{\text{fr}}}$ grows for
datasets with large $N_{\text{fr}}$, a subset
$\\{\mathbfit{l}_{k}\\}_{k=1}^{N_{\mathit{l}}}$, coined
landmark/representative points, with $N_{\mathit{l}}\leq N_{\text{fr}}$, is
selected from $\\{\check{\mathbf{y}}^{\text{f}}_{t}\\}_{t=1}^{N_{\text{fr}}}$.
Any selection strategy can be used to identify
$\\{\mathbfit{l}_{k}\\}_{k=1}^{N_{\mathit{l}}}$. Here, the min-max-distance
strategy of [38] is adopted along the lines of [35, 36]. For convenience, let
the $\nu\times N_{\mathit{l}}$ matrix
$\mathbf{L}\coloneqq[\mathbfit{l}_{1},\mathbfit{l}_{2},\ldots,\mathbfit{l}_{N_{\mathit{l}}}]$.

##  III Data Modeling

With $\varphi(\cdot)$ denoting the feature mapping which maps vector
$\mathbfit{l}_{k}$ to vector $\varphi(\mathbfit{l}_{k})$ in the feature space
$\mathscr{H}$, the crux of the modeling approach is that
$\\{\varphi(\mathbfit{l}_{k})\\}_{k=1}^{N_{\mathit{l}}}$ lie into or close to
an unknown-to-the-user smooth manifold $\mathscr{M}$ [39] embedded in
$\mathscr{H}$; see Fig. 2. To provide structured solutions, it is assumed that
$\mathscr{H}$ is a reproducing kernel Hilbert space (RKHS), equipped with a
reproducing kernel
$\kappa(\cdot,\cdot):\mathbb{C}^{\nu}\times\mathbb{C}^{\nu}\to\mathbb{C}$,
with well-documented merits in approximation theory [40]. To this end, the
feature mapping $\varphi(\cdot)$ is defined as the mapping induced by the
kernel $\kappa$ of $\mathscr{H}$; cf., [36]. A complex-valued $\mathscr{H}$ is
considered here [41, 36].

Figure 2: A “collaborative-filtering” approach: Points
$\\{\varphi(\mathbfit{l}_{k_{j}})\\}_{j=1}^{3}$, which lie into or close to
the unknown-to-the-user manifold $\mathscr{M}\subset\mathscr{H}$, collaborate
affinely to approximate $\varphi(\check{\bm{\mu}}_{t})$. All affine
combinations of $\\{\varphi(\mathbfit{l}_{k_{j}})\\}_{j=1}^{3}$ define the
approximating “linear patch” (gray-colored plane), which mimics the concept of
a tangent space to $\mathscr{M}$.

The $(i,t)$th entry $y_{it}$ of data $\mathbf{Y}$ is approximated as
$y_{it}\approx f_{i}(\check{\bm{\mu}}_{t})$, where
$f_{i}(\cdot):\mathbb{C}^{\nu}\to\mathbb{C}$ is an unknown non-linear function
that belongs to the functional space $\mathscr{H}$, and $\check{\bm{\mu}}_{t}$
is an unknown complex-valued $\nu\times 1$ vector. The well-known reproducing
property [40] of the RKHS $\mathscr{H}$ suggests that $y_{it}\approx
f_{i}(\check{\bm{\mu}}_{t})=\langle
f_{i}\mathop{}|\mathop{}\varphi(\check{\bm{\mu}}_{t})\rangle_{\mathscr{H}}$,
where $\langle\cdot\mathop{}|\mathop{}\cdot\rangle_{\mathscr{H}}$ stands for
the inner product of $\mathscr{H}$.

The concept of tangent spaces [39] (“linear patches” in simple words) to
manifolds will be used to model both $f_{i}$ and
$\varphi(\check{\bm{\mu}}_{t})$. Specifically, it is assumed that  (i)
$f_{i}$belongs to the linear span of
$\\{\varphi(\mathbfit{l}_{k})\\}_{k=1}^{N_{\mathit{l}}}$, i.e., there exists
$\mathbf{d}_{i}\coloneqq[d_{i1},\ldots,d_{iN_{\mathit{l}}}]^{\intercal}\in\mathbb{C}^{N_{\mathit{l}}}$
($\intercal$ denotes vector/matrix transposition) s.t.
$f_{i}=\sum_{k=1}^{N_{\mathit{l}}}d_{ik}\varphi(\mathbfit{l}_{k})=\bm{\Phi}(\mathbf{L})\mathbf{d}_{i}$,
where
$\bm{\Phi}(\mathbf{L})\coloneqq[\varphi(\mathbfit{l}_{1}),\ldots,\varphi(\mathbfit{l}_{N_{\mathit{l}}})]$;
and  (ii) $\varphi(\check{\bm{\mu}}_{t})$lies into or close to $\mathscr{M}$
and is approximated affinely by only a few members of
$\\{\varphi(\mathbfit{l}_{k})\\}_{k=1}^{N_{\mathit{l}}}$, i.e., there exists a
sparse $N_{\mathit{l}}\times 1$ vector $\mathbf{b}_{t}$ s.t.
$\varphi(\check{\bm{\mu}}_{t})=\bm{\Phi}(\mathbf{L})\mathbf{b}_{t}$, under the
affine constraint
$\mathbf{1}_{N_{\mathit{l}}}^{{\mathsf{H}}}\mathbf{b}_{t}=1$, where
$\mathbf{1}_{N_{\mathit{l}}}$ is the $N_{\mathit{l}}\times 1$ all-one vector
and $\mathsf{H}$ denotes complex conjugate vector/matrix transposition; cf.,
Fig. 2.

In other words, $y_{it}\approx f_{i}(\check{\bm{\mu}}_{t})=\langle
f_{i}\mathop{}|\mathop{}\varphi(\check{\bm{\mu}}_{t})\rangle_{\mathscr{H}}=\langle\bm{\Phi}(\mathbf{L})\mathbf{d}_{i}\mathop{}|\mathop{}\bm{\Phi}(\mathbf{L})\mathbf{b}_{t}\rangle_{\mathscr{H}}=\mathbf{d}_{i}^{\mathsf{H}}\mathbf{K}\mathbf{b}_{t}$,
where $\mathbf{K}$ is the complex-valued $N_{\mathit{l}}\times N_{\mathit{l}}$
matrix whose $(k,k^{\prime})$th entry is equal to
$\langle\varphi(\mathbfit{l}_{k})\mathop{}|\mathop{}\varphi(\mathbfit{l}_{k^{\prime}})\rangle_{\mathscr{H}}=\kappa(\mathbfit{l}_{k},\mathbfit{l}_{k^{\prime}})$.
To offer compact notations, if
$\mathbf{D}\coloneqq[\mathbf{d}_{1},\ldots,\mathbf{d}_{N_{\text{k}}}]^{{\mathsf{H}}}\in\mathbb{C}^{N_{\text{k}}\times
N_{\mathit{l}}}$ and
$\mathbf{B}\coloneqq[\mathbf{b}_{1},\ldots,\mathbf{b}_{N_{\text{fr}}}]\in\mathbb{C}^{N_{\mathit{l}}\times
N_{\text{fr}}}$, then data are modeled as
$\mathbf{Y}\approx\mathbf{D}\mathbf{K}\mathbf{B}$.

Choosing $\kappa(\cdot,\cdot)$ to define $\mathbf{K}$ entails the cumbersome
tasks of cross validation and fine tuning via extensive experimentation on
data sets. A popular way to surmount such cumbersome tasks, followed also in
[36], is via multiple kernels:
$\mathbf{Y}\approx\sum_{m=1}^{M}\mathbf{D}_{m}\mathbf{K}_{m}\mathbf{B}_{m}$,
with a dictionary of user-defined reproducing kernels
$\\{\kappa_{m}(\cdot,\cdot)\\}_{m=1}^{M}$, and thus kernel matrices
$\\{\mathbf{K}_{m}\\}_{m=1}^{M}$, complex-valued $N_{\text{k}}\times
N_{\mathit{l}}$ matrices $\\{\mathbf{D}_{m}\\}_{m=1}^{M}$, and
$N_{\mathit{l}}\times N_{\text{fr}}$ sparse matrices
$\\{\mathbf{B}_{m}\\}_{m=1}^{M}$ satisfying
$\mathbf{1}^{{\mathsf{H}}}_{N_{\mathit{l}}}\mathbf{B}_{m}=\mathbf{1}_{N_{\text{fr}}}^{{\mathsf{H}}}$.

To reduce the computational burden and effect low-rank constraints in the
resultant inverse problem, each $N_{\mathit{l}}\times N_{\mathit{l}}$ matrix
$\mathbf{K}_{m}$ is substituted in [36] by its low-dimensional $d\times
N_{\mathit{l}}$ rendition $\check{\mathbf{K}}_{m}$ in
$\mathbf{Y}\approx\sum_{m=1}^{M}\check{\mathbf{D}}_{m}\check{\mathbf{K}}_{m}\mathbf{B}_{m}$,
with $d\ll N_{\mathit{l}}$ and where $\\{\check{\mathbf{D}}_{m}\\}_{m=1}^{M}$
are low-rank $N_{\text{k}}\times d$ matrices. To this end,
$\check{\mathbf{K}}_{m}$ was computed from $\mathbf{K}_{m}$ via a
dimensionality-reduction module (pre-step) in [36]. However, such a
dimensionality-reduction pre-step introduces the following drawbacks:  (i)
when the numbers $M$ of kernels and $N_{\mathit{l}}$ of landmark points are
large, the dimensionality-reduction module inflicts heavy computations, while
fine-tuning its hyperparameters becomes a labor-intensive task; and  (ii) the
error from compressing $\mathbf{K}_{m}$ into $\check{\mathbf{K}}_{m}$ may
propagate to the next phase in the KRIM framework.

MultiL-KRIM avoids the previous drawbacks as follows:

$\displaystyle\mathbf{Y}\approx\sum\nolimits_{m=1}^{M}\mathbf{D}_{m}^{(1)}\mathbf{D}_{m}^{(2)}\cdots\mathbf{D}_{m}^{(Q)}\mathbf{K}_{m}\mathbf{B}_{m}\,,$
(1)

where $\mathbf{D}_{m}^{(q)}\in\mathbb{C}^{d_{q-1}\times d_{q}}$, the inner
matrix dimensions $\\{d_{q}\\}_{q=1}^{Q-1}$ are user-defined, with
$d_{0}\coloneqq N_{\text{k}}$ and $d_{Q}\coloneqq N_{\mathit{l}}$. Notice that
for $Q=2$, the term $\mathbf{D}_{m}^{(2)}\mathbf{K}_{m}$ may be considered as
the dimensionality-reduced $\check{\mathbf{K}}_{m}$ in [36]. Nonetheless,
$\\{\mathbf{D}_{m}^{(q)}\\}_{(q,m)}$ are identified during a single-stage
learning task, avoiding any pre-steps with their hyperparameter tuning and
errors. Note that $Q>1$ in (1) may offer considerable savings in computations
with respect to the $Q=1$ case. Indeed, the number of unknowns that need to be
identified in (1) for $Q>1$ is
$N_{Q>1}=M(\sum_{q=1}^{Q}d_{q-1}d_{q}+N_{\text{fr}}N_{\mathit{l}})$, as
opposed to $N_{Q=1}=M(N_{\text{k}}N_{\mathit{l}}+N_{\text{fr}}N_{\mathit{l}})$
in the $Q=1$ case. If $N_{\mathit{l}}$ is large, then
$\\{d_{q}\\}_{q=1}^{Q-1}$ can be chosen so that $N_{Q>1}\ll N_{Q=1}$.

##  IV Inverse Problem And its Iterative Solution

Letting
$\mathbfcal{D}_{1}\coloneqq[\mathbf{D}_{1}^{(1)},\mathbf{D}_{2}^{(1)},\ldots,\mathbf{D}_{M}^{(1)}]\in\mathbb{C}^{N_{\text{k}}\times
d_{1}M}$,
$\mathbfcal{B}\coloneqq[\mathbf{B}_{1}^{\mathsf{H}},\mathbf{B}_{2}^{\mathsf{H}},\ldots,\mathbf{B}_{M}^{\mathsf{H}}]^{\mathsf{H}}\in\mathbb{C}^{MN_{\mathit{l}}\times
N_{\text{fr}}}$, and the block diagonal matrices
$\mathbfcal{D}_{q}\coloneqq\operatorname{bdiag}{(\mathbf{D}_{1}^{(q)},\mathbf{D}_{2}^{(q)},\ldots,\mathbf{D}_{M}^{(q)})}\in\mathbb{C}^{d_{q-1}M\times
d_{q}M}$, $q\in\\{2,\ldots,Q\\}$,
$\mathbfcal{K}\coloneqq\operatorname{bdiag}{(\mathbf{K}_{1},\mathbf{K}_{2},\ldots,\mathbf{K}_{M})}\in\mathbb{C}^{MN_{\mathit{l}}\times
MN_{\mathit{l}}}$, (1) takes the form
$\mathbf{Y}\approx\mathbfcal{D}_{1}\mathbfcal{D}_{2}\cdots\mathbfcal{D}_{Q}\mathbfcal{K}\mathbfcal{B}$.

In dMRI, it is often the case that inverse problems are formulated in the
image domain, since, after all, this is the domain where the end-product lies
in. To this end, by using $\mathbf{X}$ to denote the image-domain data and by
defining $\mathbfcal{A}_{1}\coloneqq\mathscr{F}^{-1}(\mathbfcal{D}_{1})$ (let
also $\mathbfcal{A}_{q}\coloneqq\mathbfcal{D}_{q},\forall
q\in\\{2,\ldots,Q\\}$, for uniform notations), the following inverse problem
is postulated:

$\displaystyle\min_{(\mathbf{X},\mathbf{Z},\\{\mathbfcal{A}_{q}\\}_{q=1}^{Q},\mathbfcal{B})}{}\
$
$\displaystyle{}\tfrac{1}{2}\lVert\mathbf{X}-\mathbfcal{A}_{1}\mathbfcal{A}_{2}\cdots\mathbfcal{A}_{Q}\mathbfcal{K}\mathbfcal{B}\rVert^{2}_{\text{F}}$
$\displaystyle+\lambda_{1}\lVert\mathbfcal{B}\rVert_{1}+g(\mathbf{X},\mathbf{Z})+h(\\{\mathbfcal{A}_{q}\\}_{q=1}^{Q})$
(2a) s.to
$\displaystyle\mathscr{S}(\mathbf{Y})=\mathscr{S}\mathscr{F}(\mathbf{X})\,,$
(2b)
$\displaystyle\mathbf{1}_{N_{\mathit{l}}}^{{\mathsf{H}}}\mathbf{B}_{m}=\mathbf{1}_{N_{\text{fr}}}^{{\mathsf{H}}}\,,\forall
m\in\\{1,\ldots,M\\}\,,$ (2c) $\displaystyle\mathbfcal{A}_{q}\text{is block
diagonal},\forall q\in\\{2,\ldots,Q\\}\,,$ (2d)

where (2b) enforces consistency of the desired $\mathbf{X}$ with the data
collected in the (k,t) domain, $\lVert\cdot\rVert_{1}$ in (2a) is used to
impose sparsity on $\mathbfcal{B}$, and the convex regularizing functions
$g(\cdot)$ and $h(\cdot)$ are used to incorporate prior knowledge. More
specifically, $g(\cdot)$ employs not only variable $\mathbf{X}$ but also the
auxiliary $\mathbf{Z}$ to facilitate computations. For example, in the case
where the dMRI data capture a periodic organ movement over a static
background, then
$g(\mathbf{X},\mathbf{Z})\coloneqq(\lambda_{2}/2)\lVert\mathbf{Z}-\mathscr{F}_{\text{t}}(\mathbf{X})\rVert_{\text{F}}^{2}+\lambda_{3}\lVert\mathbf{Z}\rVert_{1}$,
where $\mathscr{F}_{\text{t}}(\cdot)$ stands for the temporal 1D DFT operator
which acts on rows of the matrix $\mathbf{X}$. Note that a row of $\mathbf{X}$
corresponds to the time series, of length $N_{\text{fr}}$, of a single pixel
in the image domain. Additionally, the designer can choose
$h(\\{\mathbfcal{A}_{q}\\}_{q=1}^{Q})\coloneqq(\lambda_{4}/2)\sum_{q=1}^{Q}\lVert\mathbfcal{A}_{q}\rVert_{\text{F}}^{2}$
to avoid unbounded solutions which may appear due to
$\mathbfcal{A}_{1}\mathbfcal{A}_{2}\cdots\mathbfcal{A}_{Q}$ in (2).

Algorithm 1 Solving MultiL-KRIM’s inverse problem

0: Limit point $\hat{\mathbf{X}}^{(*)}$ of sequence
$(\hat{\mathbf{X}}^{(n)})_{n\in\mathbb{N}}$.

1: Fix $\hat{\mathbfcal{S}}^{(0)}$, $\gamma_{0}\in(0,1]$, and $\zeta\in(0,1)$.

2: while $n\geq 0$ do

3: Available are $\hat{\mathbfcal{S}}^{(n)}$ and $\gamma_{n}$.

4: $\gamma_{n+1}\coloneqq\gamma_{n}(1-\zeta\gamma_{n})$.

5: Solve sub-tasks (3).

6:
$\hat{\mathbfcal{S}}^{(n+1)}\coloneqq\gamma_{n+1}\hat{\mathbfcal{S}}^{(n+1/2)}+(1-\gamma_{n+1})\hat{\mathbfcal{S}}^{(n)}$.

7: Set $n\leftarrow n+1$ and go to step 2.

8: end while

Alg. 1 sketches a solution to (2), where
$\hat{\mathbfcal{S}}^{(n+k/2)}\coloneqq(\hat{\mathbf{X}}^{(n+k/2)},\hat{\mathbf{Z}}^{(n+k/2)},\hat{\mathbfcal{A}}_{1}^{(n+k/2)},\ldots,\hat{\mathbfcal{A}}_{Q}^{(n+k/2)},\hat{\mathbfcal{B}}^{(n+k/2)})$,
$\forall n\in\mathbb{N}$, $\forall k\in\\{0,1\\}$. Alg. 1 is based on the
successive-convex-approximation framework of [42], which guarantees
convergence to a stationary point of the loss. The following convex sub-tasks
need to be solved per iteration: $\forall q\in\\{1,\ldots,Q\\}$,

$\displaystyle\hat{\mathbf{X}}^{(n+1/2)}$
$\displaystyle\in\arg\min_{{\mathbf{X}}}{}$ $\displaystyle{}\
\tfrac{1}{2}\lVert{\mathbf{X}}-\hat{\mathbfcal{A}}_{1}^{(n)}\cdots\hat{\mathbfcal{A}}_{Q}^{(n)}\mathbfcal{K}\hat{\mathbfcal{B}}^{(n)}\rVert_{\text{F}}^{2}$
$\displaystyle{}\
+\tfrac{\lambda_{2}}{2}\lVert\hat{\mathbf{Z}}^{(n)}-\mathscr{F}_{\text{t}}({\mathbf{X}})\rVert_{\text{F}}^{2}+\tfrac{\tau_{X}}{2}\lVert\mathbf{X}-\hat{\mathbf{X}}^{(n)}\rVert_{\text{F}}^{2}$
s.to $\displaystyle{}\
\mathscr{S}(\mathbf{Y})=\mathscr{S}\mathscr{F}({\mathbf{X}})\,,$ (3a)
$\displaystyle\hat{\mathbf{Z}}^{(n+1/2)}$
$\displaystyle\in\arg\min_{\mathbf{Z}}{}$ $\displaystyle{}\
\tfrac{\lambda_{2}}{2}\lVert\mathbf{Z}-\mathscr{F}_{\text{t}}(\hat{\mathbf{X}}^{(n)})\rVert^{2}_{\text{F}}+\lambda_{3}\lVert\mathbf{Z}\rVert_{1}$
$\displaystyle{}\
+\tfrac{\tau_{Z}}{2}\lVert\mathbf{Z}-\hat{\mathbf{Z}}^{(n)}\rVert_{\text{F}}^{2}\,,$
(3b) $\displaystyle\hat{\mathbfcal{A}}_{q}^{(n+1/2)}{}$
$\displaystyle{}\in\arg\min_{\mathbfcal{A}_{q}}{}$ $\displaystyle{}\
\tfrac{1}{2}\lVert\hat{\mathbf{X}}^{(n)}-\hat{\mathbfcal{A}}_{1}^{(n)}\cdots\mathbfcal{A}_{q}\cdots\hat{\mathbfcal{A}}_{Q}^{(n)}\mathbfcal{K}\hat{\mathbfcal{B}}^{(n)}\rVert_{\text{F}}^{2}$
$\displaystyle{}\
+\tfrac{\lambda_{4}}{2}\lVert\mathbfcal{A}_{q}\rVert_{\text{F}}^{2}+\tfrac{\tau_{A}}{2}\lVert\mathbfcal{A}_{q}-\hat{\mathbfcal{A}}_{q}^{(n)}\rVert_{\text{F}}^{2}$
s.to $\displaystyle{}\ \mathbfcal{A}_{q}\ \text{is block-diagonal}\,,\forall
q\in\\{2,\ldots,Q\\}\,,$ (3c) $\displaystyle\hat{\mathbfcal{B}}^{(n+1/2)}$
$\displaystyle\in\arg\min_{\mathbfcal{B}}{}$ $\displaystyle{}\
\tfrac{1}{2}\lVert\hat{\mathbf{X}}^{(n)}-\hat{\mathbfcal{A}}_{1}^{(n)}\cdots\hat{\mathbfcal{A}}_{Q}^{(n)}\mathbfcal{K}\mathbfcal{B}\rVert_{\text{F}}^{2}$
$\displaystyle{}\
+\lambda_{1}\lVert\mathbfcal{B}\rVert_{1}+\tfrac{\tau_{B}}{2}\lVert\mathbfcal{B}-\hat{\mathbfcal{B}}^{(n)}\rVert_{\text{F}}^{2}$
s.to $\displaystyle{}\
\mathbf{1}_{N_{\mathit{l}}}^{{\mathsf{H}}}\mathbf{B}_{m}=\mathbf{1}_{N_{\text{fr}}}^{{\mathsf{H}}}\,,\forall
m\in\\{1,\ldots,M\\}\,.$ (3d)

Sub-task (3d) is a composite convex minimization task under affine
constraints, hence it can be solved by [43], while (3a), (3b) and (3c) have
closed form solutions [36]. More specifically, the unique solution to (3b) is
provided by the well-known soft-thresholding operator [36].

##  V Numerical Tests

(a)

(b)

Figure 3: NRMSE values at different Cartesian LABEL:sub@fig:plot.cartesian and
radial LABEL:sub@fig:plot.radial acceleration rates. The NRMSE values are
plotted on a log-axis.

Following [9, 19, 20, 35], the proposed framework is validated on the magnetic
resonance extended cardiac-torso (MRXCAT) cine phantom dataset [44] under both
radial and Cartesian sampling. MultiL-KRIM[$M,Q$] is compared against its
predecessors KRIM [36] and BiLMDM [35], as well as against the low-rank tensor
factorization with total variation (LRTC-TV) [27] (designed also for dMRI
data), the DIP-based TDDIP [34] (designed especially for radial sampling), and
the popular PS-Sparse [9] and SToRM [18]. Tags KRIM[S] and KRIM[M] refer to
[36] for the case of a single ($M=1$) and multiple ($M>1$) kernels,
respectively. Comparisons of KRIM and BiLMDM against several other state-of-
the-art methods on the same data can be found in [36, 35]. All methods were
finely tuned to achieve best performance.

Parameter $Q\in\\{2,3,4\\}$, and for each $Q$ the inner matrix dimensions are
set as follows:  (i) if $Q=2$, then $d_{1}=6$;  (ii) if $Q=3$, then
$(d_{1},d_{2})=(2,6)$; and  (iii) if $Q=4$, then
$(d_{1},d_{2},d_{3})=(2,4,6)$.  The inner dimension parameter of KRIM is
$d\coloneqq 6$. Parameter $M\in\\{1,7\\}$ for both MultiL-KRIM and KRIM, with
choices of kernels as in [36]. Number of landmark points is
$N_{\mathit{l}}\coloneqq 100$. Since (2) is a non-convex task, and due to the
well-known fact that the limit point of any iterative algorithm which seeks a
stationary point of (2) depends on the starting point
$\hat{\mathbfcal{S}}^{(0)}$, Alg. 1 was run multiple times for each scenario,
with different $\hat{\mathbfcal{S}}^{(0)}$ per run, and all reported metric
values are the mean values of all those multiple runs. The software code for
(MultiL-)KRIM and BiLMDM was written in Julia [45]. All tests were run on an
8-core Intel(R) i7-11700 2.50GHz CPU with 32GB RAM.

The main evaluation metric is the normalized root mean square error
$\text{NRMSE}\coloneqq\lVert\mathbf{X}-\hat{\mathbf{X}}\rVert_{\text{F}}/\lVert\mathbf{X}\rVert_{\text{F}}$,
where $\mathbf{X}$ denotes the image data obtained from fully sampled
(k,t)-space data, and $\hat{\mathbf{X}}$ represents the estimate of
$\mathbf{X}$. Additionally, reconstructions of high-frequency regions are
evaluated by the high-frequency error norm (HFEN) and two sharpness measures
M1 (intensity-variance based) and M2 (energy of the image gradient) [46, (43)
and (46)]. Lastly, the structural similarity measure (SSIM) [47] captures
local similarities in pixel intensities after normalizing for luminance and
contrast.

It can be seen in Fig. 3, and more precisely in Tables I and II, that MultiL-
KRIM matches the performance of KRIM. Nevertheless, MultiL-KRIM shows
remarkable reduction in computational time by up to 65% for multiple kernels
and 30% for a single kernel. To showcase that the reduction in computational
time is because of the novel factorization approach, the reported times of
KRIM[M] and KRIM[S] do not include the time of KRIM’s dimensionality-reduction
pre-step. The tests also demonstrate the better performance of MultiL-KRIM
over the state-of-the-art LRTC-TV and TDDIP. Notice that TDDIP takes much
longer time to run on the CPU due to its deep CNN architecture. On the other
hand, LRTC-TV deteriorates quickly with high acceleration rates under
Cartesian sampling, with a high computational footprint.

TABLE I: Performance on Cartesian sampling (acceleration rate: 20x)

Methods $\setminus$ Metrics | NRMSE | SSIM | HFEN | M1 | M2 | Time
---|---|---|---|---|---|---
MultiL-KRIM[$M=7,Q=2$] | ${\mathbf{0.0443}}$ | $\mathbf{0.8698}$ | ${0.1147}$ | $37.26$ | $\mathbf{1.4\times 10^{6}}$ | 1.2hrs
MultiL-KRIM[$M=1,Q=2$] | ${0.0453}$ | ${0.8686}$ | ${0.1269}$ | $37.24$ | ${1.3\times 10^{6}}$ | 46min
MultiL-KRIM[$M=7,Q=3$] | ${{0.0444}}$ | ${0.8695}$ | ${0.1149}$ | $37.26$ | $\mathbf{1.4\times 10^{6}}$ | 1.3hrs
MultiL-KRIM[$M=1,Q=3$] | ${0.0454}$ | ${0.8684}$ | ${0.1270}$ | $37.24$ | ${1.3\times 10^{6}}$ | 42min
MultiL-KRIM[$M=7,Q=4$] | ${\mathbf{0.0443}}$ | ${0.8696}$ | ${0.1147}$ | $37.26$ | $\mathbf{1.4\times 10^{6}}$ | 1.4hrs
MultiL-KRIM[$M=1,Q=4$] | ${0.0453}$ | ${0.8686}$ | ${0.1269}$ | $37.24$ | ${1.3\times 10^{6}}$ | 43min
KRIM[M] [36] | $\mathbf{0.0443}$ | ${0.8696}$ | $\mathbf{0.1136}$ | $\mathbf{37.34}$ | $\mathbf{1.4\times 10^{6}}$ | 3.4hrs
KRIM[S] [36] | $0.0450$ | $0.8670$ | $0.1149$ | $37.25$ | $1.3\times 10^{6}$ | 55min
BiLMDM [35] | $0.0488$ | $0.8589$ | $0.1423$ | $37.08$ | $1.3\times 10^{6}$ | 40min
LRTC-TV [27] | $0.1645$ | $0.5942$ | $0.3834$ | $31.52$ | ${8.9\times 10^{5}}$ | 8.3hrs
SToRM [18] | $0.0850$ | $0.8110$ | $0.2504$ | $37.02$ | $1.2\times 10^{6}$ | 58min
PS-Sparse [9] | $0.0550$ | $0.8198$ | $0.1483$ | $37.18$ | $1.2\times 10^{6}$ | 15min

TABLE II: Performance of radial sampling (acceleration rate: 16x)

Methods $\setminus$ Metrics | NRMSE | SSIM | HFEN | M1 | M2 | Time
---|---|---|---|---|---|---
MultiL-KRIM[$M=7,Q=2$] | ${\mathbf{0.0448}}$ | $\mathbf{0.8680}$ | ${\mathbf{0.1023}}$ | $37.35$ | ${1.4\times 10^{6}}$ | 1.1hrs
MultiL-KRIM[$M=1,Q=2$] | ${0.0465}$ | ${0.8618}$ | ${0.1305}$ | $37.26$ | ${1.4\times 10^{6}}$ | 44min
MultiL-KRIM[$M=7,Q=3$] | ${{0.0450}}$ | ${0.8670}$ | ${{0.1130}}$ | $37.35$ | ${1.4\times 10^{6}}$ | 1.2hrs
MultiL-KRIM[$M=1,Q=3$] | ${0.0465}$ | ${0.8618}$ | ${0.1305}$ | $37.26$ | ${1.4\times 10^{6}}$ | 41min
MultiL-KRIM[$M=7,Q=4$] | ${\mathbf{0.0448}}$ | $\mathbf{0.8680}$ | ${\mathbf{0.1023}}$ | $37.35$ | ${1.4\times 10^{6}}$ | 1.2hrs
MultiL-KRIM[$M=1,Q=4$] | ${0.0465}$ | ${0.8616}$ | ${0.1307}$ | $37.26$ | ${1.4\times 10^{6}}$ | 42min
KRIM[M] [36] | ${0.0450}$ | ${0.8670}$ | ${0.1136}$ | $37.34$ | ${1.4\times 10^{6}}$ | 3hrs
KRIM[S] [36] | $0.0465$ | $0.8618$ | $0.1301$ | $37.26$ | $1.4\times 10^{6}$ | 58min
BiLMDM [35] | $0.0475$ | $0.8560$ | $0.1491$ | $37.30$ | $1.4\times 10^{6}$ | 1.4hrs
TDDIP [34] | $0.0579$ | $0.8167$ | $0.2195$ | $\mathbf{39.81}$ | $1.2\times 10^{6}$ | 20hrs
LRTC-TV [27] | $0.0738$ | $0.8063$ | $0.3725$ | $37.06$ | $\mathbf{1.7\times 10^{6}}$ | 8.6hrs
SToRM [18] | $0.0753$ | $0.8319$ | $0.3694$ | $37.38$ | $1.6\times 10^{6}$ | 30min
PS-Sparse [9] | $0.0496$ | $0.7908$ | $0.1733$ | $37.31$ | $1.4\times 10^{6}$ | 15min

##  VI Conclusions

This paper extends the KRIM framework [36] into a faster kernelized matrix
factorization framework which avoids KRIM’s dimensionality reduction pre-step.
Numerical tests demonstrate that the proposed data-modeling approach matches
the reconstruction performance of KRIM under both Cartesian and radial data
sampling, but with significant reduction in computational time, and
outperforms at the same time popular methods as well as state-of-the-art
tensor-based and deep-image-prior schemes.

## References

  * [1] Z.-P. Liang and P. C. Lauterbur, Principles of Magnetic Resonance Imaging: A Signal Processing Perspective.  IEEE Press, 2000.
  * [2] Z.-P. Liang and P. C. Lauterbur, “An efficient method for dynamic magnetic resonance imaging,” IEEE Trans. Medical Imaging, vol. 13, no. 4, pp. 677–686, 1994.
  * [3] M. Lustig, D. Donoho, and J. M. Pauly, “Sparse MRI: The application of compressed sensing for rapid MR imaging,” Magnetic Resonance in Medicine, vol. 58, no. 6, pp. 1182–1195, 2007.
  * [4] U. Gamper, P. Boesiger, and S. Kozerke, “Compressed sensing in dynamic MRI,” Magnetic Resonance in Medicine, vol. 59, no. 2, pp. 365–373, 2008\.
  * [5] D. Liang, E. V. R. DiBella, R.-R. Chen, and L. Ying, “k-t ISD: Dynamic cardiac MR imaging using compressed sensing with iterative support detection,” Magnetic Resonance in Medicine, vol. 68, no. 1, pp. 41–53, 2012.
  * [6] L. Feng, L. Axel, H. Chandarana, K. T. Block, D. K. Sodickson, and R. Otazo, “XD-GRASP: Golden-angle radial MRI with reconstruction of extra motion-state dimensions using compressed sensing,” Magnetic Resonance in Medicine, vol. 75, no. 2, pp. 775–788, 2016.
  * [7] R. Otazo, D. Kim, L. Axel, and D. K. Sodickson, “Combination of compressed sensing and parallel imaging for highly accelerated first-pass cardiac perfusion MRI,” Magnetic Resonance in Medicine, vol. 64, no. 3, pp. 767–776, 2010.
  * [8] S. G. Lingala, Y. Hu, E. V. R. DiBella, and M. Jacob, “Accelerated dynamic MRI exploiting sparsity and low-rank structure: k-t SLR,” IEEE Trans. Medical Imaging, vol. 30, no. 5, pp. 1042–1054, 2011.
  * [9] B. Zhao, J. P. Haldar, A. G. Christodoulou, and Z.-P. Liang, “Image reconstruction from highly undersampled (k,t)-space data with joint partial separability and sparsity constraints,” IEEE Trans. Medical Imaging, vol. 31, no. 9, pp. 1809–1820, 2012.
  * [10] S. P. Awate and E. V. R. DiBella, “Spatiotemporal dictionary learning for undersampled dynamic MRI reconstruction via joint frame-based and dictionary-based sparsity,” in Proceedings of ISBI, pp. 318–321, 2012\.
  * [11] Y. Wang and L. Ying, “Compressed sensing dynamic cardiac cine MRI using learned spatiotemporal dictionary,” IEEE Trans. Biomedical Eng., vol. 61, no. 4, pp. 1109–1120, 2014.
  * [12] J. Caballero, A. N. Price, D. Rueckert, and J. V. Hajnal, “Dictionary learning and time sparsity for dynamic MR data reconstruction,” IEEE Trans. Medical Imaging, vol. 33, no. 4, pp. 979–994, 2014.
  * [13] Y. Wang, N. Cao, Z. Liu, and Y. Zhang, “Real-time dynamic MRI using parallel dictionary learning and dynamic total variation,” Neurocomputing, vol. 238, pp. 410–419, 2017.
  * [14] S. Ravishankar, B. E. Moore, R. R. Nadakuditi, and J. A. Fessler, “Low-rank and adaptive sparse signal (LASSI) models for highly accelerated dynamic imaging,” IEEE Trans. Medical Imaging, vol. 36, no. 5, pp. 1116–1128, 2017.
  * [15] B. Wen, S. Ravishankar, and Y. Bresler, “FRIST—Flipping and rotation invariant sparsifying transform learning and applications,” Inverse Problems, vol. 33, no. 7, p. 074007, 2017.
  * [16] M. Usman, D. Atkinson, C. Kolbitsch, T. Schaeffter, and C. Prieto, “Manifold learning based ECG-free free-breathing cardiac CINE MRI,” J. Magnetic Resonance Imag., vol. 41, no. 6, pp. 1521–1527, 2015.
  * [17] X. Chen, M. Usman, C. F. Baumgartner, D. R. Balfour, P. K. Marsden, A. J. Reader, C. Prieto, and A. P. King, “High-resolution self-gated dynamic abdominal MRI using manifold alignment,” IEEE Trans. Medical Imaging, vol. 36, no. 4, pp. 960–971, 2017.
  * [18] S. Poddar and M. Jacob, “Dynamic MRI using smoothness regularization on manifolds (SToRM),” IEEE Trans. Medical Imaging, vol. 35, no. 4, pp. 1106–1115, 2016.
  * [19] U. Nakarmi, K. Slavakis, J. Lyu, and L. Ying, “M-MRI: A manifold-based framework to highly accelerated dynamic magnetic resonance imaging,” in Proceedings of ISBI, pp. 19–22, 2017.
  * [20] U. Nakarmi, K. Slavakis, and L. Ying, “MLS: Joint manifold-learning and sparsity-aware framework for highly accelerated dynamic magnetic resonance imaging,” in Proceedings of ISBI, pp. 1213–1216, 2018.
  * [21] U. Nakarmi, Y. Wang, J. Lyu, D. Liang, and L. Ying, “A kernel-based low-rank (KLR) model for low-dimensional manifold recovery in highly accelerated dynamic MRI,” IEEE Trans. Medical Imaging, vol. 36, no. 11, pp. 2297–2307, 2017.
  * [22] S. Poddar, Y. Q. Mohsin, D. Ansah, B. Thattaliyath, R. Ashwath, and M. Jacob, “Manifold recovery using kernel low-rank regularization: application to dynamic imaging,” IEEE Trans. Computational Imaging, vol. 5, no. 3, pp. 478–491, 2019.
  * [23] O. Arif, H. Afzal, H. Abbas, M. F. Amjad, J. Wan, and R. Nawaz, “Accelerated dynamic MRI using kernel-based low rank constraint,” J. Medical Systems, vol. 43, no. 8, p. 271, 2019.
  * [24] J. Liu, P. Musialski, P. Wonka, and J. Ye, “Tensor completion for estimating missing values in visual data,” IEEE Trans. Pattern Analysis and Machine Intelligence, vol. 35, no. 1, pp. 208–220, 2012.
  * [25] M. Signoretto, L. D. Lathauwer, and J. A. K. Suykens, “Learning tensors in reproducing kernel Hilbert spaces with multilinear spectral penalties,” arXiv:1310.4977, 2013.
  * [26] H. Kanagawa, T. Suzuki, H. Kobayashi, N. Shimizu, and Y. Tagami, “Gaussian process nonparametric tensor estimator and its minimax optimality,” in Proceedings of ICML, vol. 48, pp. 1632–1641, 2016.
  * [27] X. Li, Y. Ye, and X. Xu, “Low-rank tensor completion with total variation for visual data inpainting,” Proceedings of AAAI Conference on Artificial Intelligence, vol. 31, Feb. 2017.
  * [28] J. Schlemper, J. Caballero, J. V. Hajnal, A. N. Price, and D. Rueckert, “A deep cascade of convolutional neural networks for dynamic MR image reconstruction,” IEEE Trans. Medical Imaging, vol. 37, no. 2, pp. 491–503, 2018.
  * [29] S. Biswas, H. K. Aggarwal, and M. Jacob, “Dynamic MRI using model-based deep learning and SToRM priors: MoDL-SToRM,” Magnetic Resonance in Medicine, vol. 82, pp. 485–494, July 2019.
  * [30] D. Liang, J. Cheng, Z. Ke, and L. Ying, “Deep magnetic resonance image reconstruction: Inverse problems meet neural networks,” IEEE Signal Processing Magaz., vol. 37, no. 1, pp. 141–151, 2020.
  * [31] C. M. Sandino, J. Y. Cheng, F. Chen, M. Mardani, J. M. Pauly, and S. S. Vasanawala, “Compressed sensing: From research to clinical practice with deep neural networks: Shortening scan times for magnetic resonance imaging,” IEEE Signal Processing Magaz., vol. 37, no. 1, pp. 117–127, 2020.
  * [32] V. Antun, F. Renna, C. Poon, B. Adcock, and A. C. Hansen, “On instabilities of deep learning in image reconstruction and the potential costs of AI,” Proceedings of the National Academy of Sciences, vol. 117, no. 48, pp. 30088–30095, 2020.
  * [33] Q. Zou, A. H. Ahmed, P. Nagpal, S. Kruger, and M. Jacob, “Dynamic imaging using a deep generative SToRM (Gen-SToRM) model,” IEEE Trans. Medical Imaging, vol. 40, no. 11, pp. 3102–3112, 2021.
  * [34] J. Yoo, K. H. Jin, H. Gupta, J. Yerly, M. Stuber, and M. Unser, “Time-dependent deep image prior for dynamic MRI,” IEEE Trans. Medical Imaging, vol. 40, no. 12, pp. 3337–3348, 2021.
  * [35] G. N. Shetty, K. Slavakis, A. Bose, U. Nakarmi, G. Scutari, and L. Ying, “Bi-linear modeling of data manifolds for dynamic-MRI recovery,” IEEE Trans. Medical Imaging, vol. 39, no. 3, pp. 688–702, 2020.
  * [36] K. Slavakis, G. N. Shetty, L. Cannelli, G. Scutari, U. Nakarmi, and L. Ying, “Kernel regression imputation in manifolds via bi-linear modeling: The dynamic-MRI case,” IEEE Trans. Computational Imaging, vol. 8, pp. 133–147, 2022.
  * [37] L. Györfi, M. Kohler, A. Krzyżak, and H. Walk, A Distribution-Free Theory of Nonparametric Regression.  New York: Springer, 2010.
  * [38] V. De Silva and J. B. Tenenbaum, “Sparse multidimensional scaling using landmark points,” tech. rep., Stanford University, 2004.
  * [39] J. W. Robbin and D. A. Salamon, Introduction to Differential Geometry.  Berlin: Springer, 2022.
  * [40] N. Aronszajn, “Theory of reproducing kernels,” Trans. American Mathematical Society, vol. 68, no. 3, pp. 337–404, 1950.
  * [41] P. Bouboulis and S. Theodoridis, “Extension of Wirtinger’s calculus to reproducing kernel Hilbert spaces and the complex kernel LMS,” IEEE Trans. Signal Processing, vol. 59, no. 3, pp. 964–978, 2010.
  * [42] F. Facchinei, G. Scutari, and S. Sagratella, “Parallel selective algorithms for nonconvex big data optimization,” IEEE Trans. Signal Processing, vol. 63, no. 7, pp. 1874–1889, 2015.
  * [43] K. Slavakis and I. Yamada, “Fejér-monotone hybrid steepest descent method for affinely constrained and composite convex minimization tasks,” Optimization, vol. 67, no. 11, pp. 1963–2001, 2018.
  * [44] L. Wissmann, C. Santelli, W. P. Segars, and S. Kozerke, “MRXCAT: Realistic numerical phantoms for cardiovascular magnetic resonance,” J. Cardiovascular Magnetic Resonance, vol. 16, no. 1, pp. 1–11, 2014.
  * [45] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, “Julia: A fresh approach to numerical computing,” SIAM Review, vol. 59, no. 1, pp. 65–98, 2017.
  * [46] M. Subbarao, T.-S. Choi, and A. Nikzad, “Focusing techniques,” Optical Engineering, vol. 32, no. 11, pp. 2824–2836, 1993.
  * [47] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: From error visibility to structural similarity,” IEEE Trans. Medical Imaging, vol. 13, no. 4, pp. 600–612, 2004.