11institutetext: University of Padova, Padua, Italy 11email:
<EMAIL_ADDRESS>22institutetext: Polytechnic
of Torino, Turin, Italy 22email<EMAIL_ADDRESS>33institutetext:
University of Genova, Genoa, Italy 33email:
<EMAIL_ADDRESS>
# Mapped Variably Scaled Kernels: Applications to Solar Imaging

Francesco Marchetti 11 0000-0003-1087-7589 Emma Perracchione 22
0000-0003-2663-7803 Anna Volpara 33 0000-0003-1479-0827 Anna Maria Massone 33
0000-0003-4966-8864 Stefano De Marchi 11 0000-0002-2832-8476 Michele Piana 33
0000-0003-1700-991X

###### Abstract

Variably scaled kernels and mapped bases constructed via the so-called fake
nodes approach are two different strategies to provide adaptive bases for
function interpolation. In this paper, we focus on kernel-based interpolation
and we present what we call mapped variably scaled kernels, which take
advantage of both strategies. We present some theoretical analysis and then we
show their efficacy via numerical experiments. Moreover, we test such a new
basis for image reconstruction tasks in the framework of hard X-ray
astronomical imaging.

###### Keywords:

Variably scaled kernels Mapped bases interpolation Hard X-ray imaging.

##  1 Introduction

Kernel-based interpolation is an effective approach to deal with the scattered
data interpolation problem, where data sites do not necessarily belong to some
particular structure or grid [17, 30]. Therefore, because of its flexibility
and the achievable accuracy, it finds application in many different contexts
[20], including image reconstruction and the numerical solution of partial
differential equations. The effectiveness of radial kernel-based
interpolation, which is also known as Radial Basis Function (RBF)
interpolation, very often relies on a good choice of the so-called shape
parameter, which rules the shape of basis functions. However, in many
situations, a fine tuning of this shape parameter is not sufficient to
construct basis functions that are tailored with respect to both the data
distribution and the target function to be recovered.

Therefore, in order to gain more adaptivity in the interpolation process,
Variably Scaled Kernels (VSKs) have been introduced in [6], and further
analyzed in [8] in the more general framework of kernel-based regression
networks, and in [16] in the context of persistent homology. The VSK setting
has also been employed in the reconstruction of functions presenting jumps,
leading to the definition of Variably Scaled Discontinuous Kernels (VSDKs)
[13], which turned out to be effective in medical image reconstruction tasks
in the field of Magnetic Particle Imaging (MPI) [12].

An alternative approach for constructing data-dependent or target-dependent
basis functions, the so-called Fake Nodes Approach (FNA), has been introduced
in [14] in the framework of univariate polynomial interpolation, and then
extended to rational barycentric approximation [4] and general multivariate
interpolation, including the RBF framework [15]. In particular, in latter
paper the authors showed that the VSDK setting and the FNA can lead to very
similar results when facing the Gibbs phenomenon in the reconstruction of
discontinuous functions. Moreover, in [11] an effective scheme for dealing
with both Gibbs and Runge’s phenomena has been proposed.

In this paper, our purpose is to design a unified kernel-based approach for
function interpolation by taking advantage of both the VSK setting and the
FNA, which are recalled in Section 2 and 3, respectively. The resulting
kernels that we call Mapped VSKs (MVSKs) are defined in Section 4 and provided
with an original theoretical contribution, which includes a focus on the
discontinuous case tested in Section 5. In Section 6 we present the
application of the proposed method to the inverse problems of solar hard X-ray
imaging and, specifically, we consider data from the on the ESA
Spectrometer/Telescope for Imaging X-rays (STIX) telescope, on board of Solar
Orbiter mission. Finally, in Section 7 we draw some conclusions.

##  2 Variably Scaled Kernel-Based Approximation

We refer to [17, 30] for the following introduction.

Let $\Omega\subseteq\mathbb{R}^{d}$, $d\in\mathbb{N}$, and
$X=X_{N}=\\{\bm{x}_{i},\;i=1,\ldots,N\\}\subset\Omega$ be a set of possible
scattered distinct nodes, $N\in\mathbb{N}$. Suppose that we wish to
reconstruct an unknown function $f\mathrel{\mathop{\mathchar
58\relax}}\Omega\longrightarrow\mathbb{R}$ from its values at $X_{N}$, i.e.
from the vector
$\bm{f}=(f(\bm{x}_{1}),\dots,f(\bm{x}_{n}))^{\intercal}=(f_{1},\dots,f_{n})^{\intercal}$.
In kernel-based interpolation, this is obtained by considering an
approximating function of the form

$R_{f,X}(\bm{x})=\sum_{i=1}^{n}c_{i}\kappa_{\varepsilon}(\bm{x},\bm{x}_{i}),\quad\bm{x}\in\Omega,$

where $\bm{c}=(c_{1},\dots,c_{n})^{\intercal}\in\mathbb{R}^{N}$ and
$\kappa_{\varepsilon}\mathrel{\mathop{\mathchar
58\relax}}\Omega\times\Omega\longrightarrow\mathbb{R}$ is a strictly positive
definite kernel, which depends on a shape parameter $\varepsilon>0$. In the
following, we may use the shortened notation $\kappa=\kappa_{\varepsilon}$.
The interpolation conditions are imposed by employing $\bm{c}$ so that

$\mathsf{K}\bm{c}=\bm{f},$ (1)

where $\mathsf{K}=(\mathsf{K}_{i,j})=\kappa(\bm{x}_{i},\bm{x}_{j})$,
$i,j=1,\dots,N$, is the so-called interpolation (or collocation or simply
kernel) matrix. The vector $\bm{c}$ that satisfies (1) is unique as long as
$\kappa$ is strictly positive definite. Moreover, we assume
$\kappa_{\varepsilon}$ to be radial, i.e. there exists a univariate function
$\varphi_{\varepsilon}\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}_{\geq
0}\longrightarrow\mathbb{R}$ such that
$\kappa_{\varepsilon}(\bm{x},\bm{y})=\varphi_{\varepsilon}(r)$, with
$r\coloneqq\lVert\bm{x}-\bm{y}\lVert_{2}$.

The interpolant $R_{f,\mathcal{X}}$ belongs to the dot-product space

$\mathcal{H}_{\kappa}=\mathrm{span}\left\\{\kappa(\cdot,\bm{x}),\;\bm{x}\in\Omega\right\\},$

whose completion with respect to the norm
$\lVert\cdot\lVert_{\mathcal{H}_{\kappa}}=\sqrt{(\cdot,\cdot)_{\mathcal{H}_{\kappa}}}$
induced by the bilinear form $(\cdot,\cdot)_{\mathcal{H}_{\kappa}}$ is the
native space $\mathcal{N}_{\kappa}$ associated to $\kappa$. The well-known
pointwise error bound [17, Theorem 14.2, p. 117]

$|f(\bm{x})-R_{f,X}(\bm{x})|\leq P_{\kappa,X}(\bm{x})\lVert
f\lVert_{\mathcal{N}_{\kappa}},\quad
f\in\mathcal{N}_{\kappa},\quad\bm{x}\in\Omega,$ (2)

involves the power function
$P_{\kappa,X}=\lVert\kappa(\cdot,x)-\bm{\kappa}(\cdot)^{\intercal}\mathsf{K}^{-1}\bm{\kappa}(\bm{x})\lVert_{\mathcal{H}_{\kappa}}$,
where
$\bm{\kappa}(\bm{x})\coloneqq(\kappa(\bm{x},\bm{x}_{1}),\dots,\kappa(\bm{x},\bm{x}_{N}))^{\intercal}$.
In the estimate (2), a noteworthy property is that the error is split into a
first term that only depends on the nodes and on the kernel, and a second term
that relies on the underlying function $f$.

A different perspective is provided by the following error bound, which takes
into account the so-called fill distance

$h_{X,\Omega}\coloneqq\sup_{\bm{x}\in\Omega}{\min_{\bm{x}_{k}\in{X}}\lVert\bm{x}-\bm{x}_{k}\lVert_{2}}.$

Assuming $\Omega$ to be bounded and satisfying an interior cone condition, and
$\kappa\in C^{2k}\left(\Omega\times\Omega\right)$, there exists
$C_{\kappa}(\bm{x})$ such that

$|f\left(\bm{x}\right)-R_{f,X}\left(\bm{x}\right)|\leq
C_{\kappa}(\bm{x})h^{k}_{X,\Omega}\lVert f\lVert_{{\cal N}_{\kappa}(\Omega)}.$
(3)

The factor $C_{\kappa}(\bm{x})$ depends on the maximum of kernel derivatives
of degree $2k$ in a neighborhood of $\bm{x}\in\Omega$, and $h^{k}_{X,\Omega}$
needs to be small enough; we refer to [17, Section 14.5] for a detailed
presentation of this result.

While the theoretically achievable convergence rate is influenced by the fill
distance and the smoothness of the kernel, two terms play an important role in
affecting the conditioning of the interpolation process:

  * •

The separation distance $q_{X}\coloneqq\dfrac{1}{2}\min_{i\neq
j}\lVert\bm{x}_{i}-\bm{x}_{j}\lVert_{2}.$

  * •

The value of the shape parameter $\varepsilon$ of the kernel
$\kappa_{\varepsilon}$.

Precisely, the interpolation process gets more ill-conditioned as the
separation distance becomes smaller, which is what usually happens in practice
when increasing the number of interpolation nodes, thus reducing the value of
the fill-distance (this is often denoted as a trade-off principle in RBF
literature). As far as the shape parameter is concerned, a large value
produces very localized basis functions that lead to a well-conditioned but
likely inaccurate approximation scheme, while by lowering such value we may
obtain a more accurate reconstruction at the price of an ill-conditioned
setting. As a consequence, the tuning of the shape parameter is a non trivial
problem in kernel-based interpolation.

In order to partially overcome such instability issues, in [6] the authors
introduced Variably Scaled Kernels (VSKs), which are defined as follows.
Letting $\psi\mathrel{\mathop{\mathchar
58\relax}}\Omega\longrightarrow\mathbb{R}$ be a scaling or shape function and
$\kappa_{\varepsilon}\mathrel{\mathop{\mathchar
58\relax}}(\Omega\times\mathbb{R})\times(\Omega\times\mathbb{R})\longrightarrow\mathbb{R}$
a kernel, a VSK $\kappa_{\psi}\mathrel{\mathop{\mathchar
58\relax}}\Omega\times\Omega\longrightarrow\mathbb{R}$ is defined as

$\kappa_{\psi}(\bm{{x}},\bm{{y}})\coloneqq\kappa_{\varepsilon}((\bm{{x}},\psi(\bm{x})),(\bm{y},\psi(\bm{y})),$

for $\bm{x},\bm{y}\in\Omega$. Note that we can consider the related native
space $\mathcal{N}_{\kappa_{\psi}}$ spanned by the functions
$\kappa_{\psi}(\cdot,\bm{{x}})$, $\bm{x}\in\Omega$. Among the properties of
VSKs, we recall that:

  1. A1.

The theoretical analysis of the VSK setting reduces to the analysis of the
classical framework in the augmented space $\Omega\times\mathbb{R}$.

  2. A2.

The spaces $\mathcal{N}_{\kappa}$ and $\mathcal{N}_{\kappa_{\psi}}$ are
isometrically isomorphic.

The shape parameter is often set to $\varepsilon=1$ in the VSK framework.
Therefore, the role played by the shape function is duplex. On the one hand,
the tuning of the shape parameter is substituted with the choice of the shape
function. On the other hand, $\psi$ can lead to an improvement in the
conditioning of the approximation process by increasing the value of the
separation distance in the augmented domain $\Omega\times\mathbb{R}$. In the
following, we do not necessarily restrict to a fixed value of the shape
parameter, as also done in the context of moving least squares [21]. Note that
$\kappa_{\psi}$ is (strictly) positive definite if so is $\kappa$.

##  3 Interpolation via Mapped Bases

In this section, we present the so called Fake Nodes Approach (FNA) for the
kernel-based approximation framework; for a more general and comprehensive
treatment we refer to [15].

Let $S\mathrel{\mathop{\mathchar
58\relax}}\Omega\longrightarrow\mathbb{R}^{d}$ be an injective map. Our
purpose is to construct an interpolant $R^{S}_{f,X}$ of the function $f$ in
the space

$\mathcal{H}_{\kappa^{S}}=\mathrm{span}\left\\{\kappa^{S}(\cdot,\bm{x}),\;\bm{x}\in\Omega\right\\},$

where $\kappa^{S}(\bm{x},\bm{y})\coloneqq\kappa(S(\bm{x}),S(\bm{y}))$. We have

$R^{S}_{f,X}(\bm{x})=\sum_{i=1}^{N}c^{S}_{i}\kappa^{S}(\bm{x},\bm{x}_{i})=\sum_{i=1}^{N}c^{S}_{i}\kappa(S(\bm{x}),S(\bm{x}_{i}))=R_{g,X}(S(\bm{x})),\quad\bm{x}\in\Omega,$

where $g_{|S(X)}=f_{|X}$. In other words, the construction of the interpolant
$R^{S}_{f,X}\in\mathcal{H}_{\kappa^{S}}$ is equivalent to the construction of
a classical interpolant $R_{g,S(X)}\in\mathcal{H}_{\kappa}$ at the _fake_
nodes $S(X)$. Similarly to the VSK setting, the FNA is provided with the
following properties (cf. A1 and A2):

  1. B1.

The theoretical analysis of the FNA interpolant reduces to the analysis of the
interpolant in the classical framework. This was proved in [15, Proposition
3.4] and it is linked to the inheritance property of the Lebesgue constant
[7].

  2. B2.

The spaces $\mathcal{N}_{\kappa}$ and $\mathcal{N}_{\kappa^{S}}$ are
isometrically isomorphic (a direct consequence of Theorem 4.1 and Proposition
1 stated below).

In previous works, the FNA has been mainly employed for two main purposes. The
first is using $S$ to obtain an interpolation design $S(X)$ that leads to a
more stable interpolation process with respect to the original set of nodes.
For example, this can be achieved in a polynomial-based framework by mapping
onto the set of (tensor-product) Chebyshev-Lobatto nodes in the (multi) one-
dimensional case. In addition, as we will discuss in Subsection 4, $S$ can be
constructed in a target-dependent fashion to emulate the possible
discontinuities of the underlying function and thus recover accuracy near jump
points.

##  4 Mapped VSKs

###  4.1 The general framework

In Sections 2 and 3, we outlined two different approaches, which however
present some similarities and are employed for analogous purposes. In the
following, we discuss how the VSK setting and the FNA can be merged in a
unified framework. We start by giving the following definition.

###### Definition 1

Let $S\mathrel{\mathop{\mathchar
58\relax}}\Omega\longrightarrow\mathbb{R}^{d}$ be an injective map,
$\psi\mathrel{\mathop{\mathchar 58\relax}}\Omega\longrightarrow\mathbb{R}$ be
a shape function and let $\kappa_{\varepsilon}\mathrel{\mathop{\mathchar
58\relax}}(S(\Omega)\times\mathbb{R})\times(S(\Omega)\times\mathbb{R})\longrightarrow\mathbb{R}$
be a kernel. Then, a Mapped VSK (MVSK)
$\kappa^{S}_{\psi}\mathrel{\mathop{\mathchar
58\relax}}\Omega\times\Omega\longrightarrow\mathbb{R}$ is defined as

$\kappa^{S}_{\psi}(\bm{{x}},\bm{{y}})\coloneqq\kappa_{\varepsilon}((S(\bm{{x}}),\psi(\bm{x})),(S(\bm{y}),\psi(\bm{y})),$

for $\bm{x},\bm{y}\in\Omega$.

We remark that $\kappa^{S}_{\psi}$ might be defined in different possible
equivalent manners. However, the advantage of Definition 1 lies in the
separation between the actions of $S$ and $\psi$: The function $S$ works in
the original dimension $\mathbb{R}^{d}$, while $\psi$ rules the coordinate of
the input in the augmented dimension. Under certain assumptions, a MVSK
reduces to a mapped or VSK.

###### Proposition 1

Let $\kappa^{S}_{\psi}$ be a MVSK on $\Omega\times\Omega$ built upon a radial
kernel $\kappa_{\varepsilon}$. Then:

  1. 1.

If $\psi(\bm{x})\equiv\alpha\in\mathbb{R}$, then
$\kappa^{S}_{\psi}=\kappa^{S}$.

  2. 2.

If $S(\bm{x})-S(\bm{y})=\bm{x}-\bm{y}$, (e.g., $S$ is the identity map), then
$\kappa^{S}_{\psi}=\kappa_{\psi}$.

###### Proof

By hypothesis, there exists $\varphi\mathrel{\mathop{\mathchar
58\relax}}[0,+\infty)\longrightarrow\mathbb{R}$ such that
$\kappa_{\varepsilon}(\bm{x},\bm{y})=\varphi(\lVert\bm{x}-\bm{y}\lVert_{2})$.
Therefore, we can write

$\begin{split}\kappa^{S}_{\psi}(\bm{{x}},\bm{{y}})&=\varphi(\lVert(S(\bm{{x}}),\psi(\bm{x}))-(S(\bm{{y}}),\psi(\bm{y}))\lVert_{2})\\\
&=\varphi(\sqrt{(S(\bm{x})-S(\bm{y}))^{2}+(\psi(\bm{x})-\psi(\bm{y}))^{2}}),\end{split}$

from which the two theses follow.

We also prove the following.

###### Theorem 4.1

The spaces $\mathcal{H}_{\kappa}$ and $\mathcal{H}_{\kappa^{S}_{\psi}}$ are
isometric.

###### Proof

By defining the map
$\Lambda_{\psi}^{S}(\bm{x})\coloneqq(S(\bm{x}),\psi(\bm{x}))$ we can see
$\kappa^{S}_{\psi}$ as the push-forward of $\kappa$ in the sense provided in
[29, Equation 2.52]. Then, since $S$ is injective, the proof follows from [29,
Theorem 2.9].

As a consequence of Theorem 4.1, the native spaces $\mathcal{N}_{\kappa}$ and
$\mathcal{N}_{\kappa^{S}_{\psi}}$ are isometrically isomorphic (cf. [6,
Section 3]).

Independently of the target function to be recovered, we recall that the
conditioning of the interpolation problem is related to the
$\ell_{2}$-conditioning of the kernel matrix
$\mathrm{cond}(\mathsf{K})=\lambda_{\max}/\lambda_{\min}$, being
$\lambda_{\max},\;\lambda_{\min}$ the maximum and minimum eigenvalue,
respectively. It is known that $\lambda_{\min}$ decays according to the
separation distance $q_{X}$, in a way that is influenced by the regularity of
the kernel (see [17, Chapter 16]). In this direction, MVSKs can be employed in
order to increase the separation distance and thus improve the conditioning of
the interpolation scheme, e.g., by separating clustered nodes. On the other
hand, diminishing the fill distance may improve the accuracy of the method
(see (3)). We will experiment on this in Section 5.

###  4.2 Working with mapped discontinuous kernels

In the following, we focus on the case of discontinuous functions. First, we
briefly review in which manners VSKs and the FNA were used in the
discontinuous setting.

####  4.2.1 Variably Scaled Discontinuous Kernels (VSDKs).

The idea proposed in [13] and further investigated in [12] was to define the
scaling function $\psi$ to be discontinuous at the jumps of the target
function $f\mathrel{\mathop{\mathchar
58\relax}}\Omega\longrightarrow\mathbb{R}$. In order to do so, we assume
$\Omega\subset\mathbb{R}^{d}$ to be a bounded set such that:

  * •

$\Omega$ is the union of $m$ pairwise disjoint sets $\Omega_{k}$,
$k\in\\{1,\ldots,m\\}$.

  * •

Each subset $\Omega_{k}$, $k=1,\dots,m$, has a Lipschitz boundary.

  * •

The discontinuity points of $f$ are contained in the union of the boundaries
of the subsets $\Omega_{k}$, $k=1,\dots,m$.

Then, letting $\bm{\alpha}=(\alpha_{1},\ldots,\alpha_{m})$,
$\alpha_{i}\in\mathbb{R}$, the scaling function $\psi$ is such that:

  * •

$\psi$ is piecewise constant, such that
${\left.\kern-1.2pt\psi(\bm{x})\right|_{\Omega_{k}}}=\alpha_{k}$.

  * •

$\alpha_{i}\neq\alpha_{j}$ if $\Omega_{i}$ and $\Omega_{j}$ are neighboring
sets.

The theoretical analysis carried out in the referring papers then focused on
radial kernels whose related univariate function $\varphi$ has the following
Fourier decay

$(\mathrm{F}{\varphi})(\bm{\omega})\sim(1+\lVert\bm{\omega}\lVert^{2})^{-s-\frac{1}{2}},\quad
s>\frac{d-1}{2}.$ (4)

The native space of the kernels that satisfy (4), e.g. Matérn and Wendland
kernels, is a Sobolev space [30, Chapter 10]. In order to present an error
bound in terms of Sobolev spaces norm for the VSDK setting, we introduce two
necessary ingredients:

  1. 1.

The regional fill-distance

$h_{k}\coloneqq
h_{X,\Omega_{k}}=\sup_{\bm{x}\in\Omega_{k}}\min_{\bm{x}_{k}\in{X}\cap\Omega_{k}}\lVert\bm{x}-\bm{x}_{k}\lVert_{2}$

and the global fill distance

$h\coloneqq\max_{k\in\\{1,\ldots,m\\}}h_{k}.$

  2. 2.

Letting $s\geq 0$ and $1\leq p\leq\infty$, we define the space

$\mathcal{W}\mathcal{P}_{p}^{s}(\Omega)\coloneqq\left\\{f\mathrel{\mathop{\mathchar
58\relax}}\Omega\longrightarrow\mathbb{R}\>|\>{\left.\kern-1.2ptf\right|_{\Omega_{k}}}\in\mathcal{W}_{p}^{s}(\Omega_{k}),\;k\in\\{1,\ldots,m\\}\right\\},$

which contains the piecewise smooth functions $f$ on $\Omega$ whose
restriction to any subregion $\Omega_{k}$, $i=k,\dots,m$, is contained in the
standard Sobolev space $\mathcal{W}_{p}^{s}(\Omega_{k})$. The space
$\mathcal{W}\mathcal{P}_{p}^{s}(\Omega)$ is endowed with the norm

$\lVert
f\lVert_{\mathcal{W}\mathcal{P}_{p}^{s}(\Omega)}^{p}=\sum_{k=1}^{m}\lVert{\left.\kern-1.2ptf\right|_{\Omega_{k}}}\lVert_{\mathcal{W}_{p}^{s}(\Omega_{k})}^{p}.$

In the outlined assumptions and letting $R^{\psi}_{f,X}$ be the kernel-based
interpolant of $f$ at $X$ built upon the VSDK $\kappa_{\psi}$, in [12, Theorem
3.4] the authors proved the following. Let $s>0$, $1\leq q\leq\infty$ and
$t\in\mathbb{N}_{0}$ such that $\lfloor s\rfloor>t+\frac{d}{2}$. Then, for
$f\in\mathcal{W}\mathcal{P}_{2}^{s}(\Omega)$ and sufficiently small $h$ we
have that

$\lVert f-R^{\psi}_{f,X}\lVert_{\mathcal{W}\mathcal{P}_{q}^{t}(\Omega)}\leq
Ch^{s-t-d(1/2-1/q)_{+}}\lVert
f\lVert_{\mathcal{W}\mathcal{P}_{2}^{s}(\Omega)},$ (5)

where the constant $C>0$ is independent of $h$.

####  4.2.2 The S-Gibbs map in the FNA.

In the mapped bases approach, similarly to the VSDK framework, the intuition
is to map the nodes in order to create some gaps in presence of
discontinuities. To do so, considering the collection of subsets
$\Omega_{1},\dots,\Omega_{m}$ employed to construct VSDKs, the so-called
S-Gibbs map is designed as

$S(\bm{x})=\bm{x}+\sum_{k=1}^{m}\bm{\beta}_{k}\mbox{\Large$\chi$}_{\Omega_{k}}(\bm{x}),$

where $\bm{\beta}_{k}=(k\beta,\dots,k\beta)\in\mathbb{R}^{d}$,
$\beta\in\mathbb{R}$, and $\mbox{\Large$\chi$}_{\Omega_{k}}$ is the
characteristic function corresponding to $\Omega_{k}$. In [15, Section 4.2],
the authors discussed the analogies between a VSDK $\kappa_{\psi}$ and a
mapped kernels $\kappa^{S}$ constructed via the S-Gibbs map. It turned out
that these two approaches can lead to similar results for certain values of
the vectors of parameters $\bm{\alpha}$ and $\bm{\beta}_{k}$. This is due to
the fact that the interpolation matrices are close being the kernel radial.

####  4.2.3 Mapped VSDKs (MVSDKs).

In the mixed approach with the mapped VSK kernel $\kappa_{\psi}^{S}$, we deal
with the jumps of the underlying function as follows.

  * •

We define $\psi$ as in the VSDKs framework. Therefore, the role of the shape
function is to mimic the jumps of the target function and thus to prevent the
appearance of the Gibbs phenomenon.

  * •

Since the discontinuities are already addressed by $\psi$, we employ $S$ to
map the set of nodes $X$ to obtain improved distributions locally on each
$\Omega_{k}$, meaning that we aim at diminishing the global fill distance and
increasing the global separation distance defined as

$q\coloneqq\min_{k\in\\{1,\ldots,m\\}}q_{k},$

being $q_{k}\coloneqq q_{X_{k}}=\dfrac{1}{2}\min_{i\neq
j}\lVert\bm{x}_{i}-\bm{x}_{j}\lVert_{2},\;\bm{x}_{i},\bm{x}_{j}\in\Omega_{k}$.

Recalling the error estimate in (5), this proposed construction for MVSDKs may
lead to better results, being $h$ smaller. Moreover, by increasing the
separation distance, an improvement in the conditioning of the scheme with
respect to classical VSDKs is likely to be obtained. We test these aspects in
some numerical examples in the next section.

##  5 Numerical tests with MVSDKs

The purpose of this section is to provide a numerical example to show the
benefits of the MVSDK framework in comparison to VSDKs and classical RBF
interpolation. To do so, we set $\Omega=[-1,1]^{2}$ and letting
$\bm{x}=(x_{1},x_{2})$ we consider the target function

$f\mathrel{\mathop{\mathchar 58\relax}}\Omega\longrightarrow\mathbb{R},\quad
f(\bm{x})=\left\\{\begin{array}[]{ll}x_{1}+x_{2},&\quad x_{1}<-0.3,\\\
\sin(x_{1}-2x_{2}),&\quad 0\leq x_{1}<0.5,\\\ 0,&\quad\textrm{otherwise.}\\\
\end{array}\right.$

To deal with the jumps of $f$, we define the shape function

$\psi\mathrel{\mathop{\mathchar
58\relax}}\Omega\longrightarrow\mathbb{R},\quad\psi(\bm{x})=\left\\{\begin{array}[]{ll}0,&\quad
x_{1}<-0.3,\\\ 1,&\quad-0.3\leq x_{1}<0,\\\ 2,&\quad 0\leq x_{1}<0.5,\\\
3,&\quad x_{1}\geq 0.5.\\\ \end{array}\right.$

As far as the nodes are concerned, we let $G_{N}$ be a set of $N$ nodes that
are sampled from a bivariate normal distribution with mean $\bm{\mu}=(0,0)$
and covariance matrix $\mathsf{\Sigma}=0.1\cdot\mathsf{I}$, being $\mathsf{I}$
the $2\times 2$ identity matrix. We can then consider the following mapping
function $S\mathrel{\mathop{\mathchar 58\relax}}\Omega\longrightarrow\Omega$
defined as

$S(x_{1},x_{2})=\bigg{(}1+\mathrm{erf}\bigg{(}\frac{x_{1}}{\sqrt{0.2}}\bigg{)},1+\mathrm{erf}\bigg{(}\frac{x_{2}}{\sqrt{0.2}}\bigg{)}\bigg{)}-1,$

where erf is the well-known error function. To clarify the idea behind the
construction of $S$, it is known from classical probability theory that if
$z_{1},\dots,z_{n}$ are sampled according to a normal distribution of mean
$\mu$ and standard deviation $\sigma$, then
$0.5(1+\textrm{erf}((z_{i}-\mu)/(\sqrt{2}\sigma)))$, $i=1,\dots,n$, are
distributed uniformly in $[0,1]$. Therefore, our set of nodes $G_{N}$ is
mapped to a uniform distribution in the square $[-1,1]^{2}$, as displayed in
Figure 1. We point out that possible nodes in $G_{N}$ that are not in $\Omega$
are removed from the set.

Figure 1: $N=400$. Left: $G_{N}$. Right: $S(G_{N})$.

Consequently, the mapped set $S(G_{N})$ is very likely to present smaller fill
distance and larger separation distance than $G_{N}$. The interpolation
results are evaluated on a finer $M\times M$ equispaced grid $\Xi_{M^{2}}$ in
$\Omega$. Precisely, we compute the Root Mean Square Error (RMSE)

$\displaystyle{\rm
RMSE}=\sqrt{\frac{1}{M^{2}}\sum_{k=1}^{M^{2}}\left(f(\bm{\xi}_{k})-\iota(\bm{\xi}_{k})\right)^{2}},$
(6)

where $\bm{\xi}_{k}\in\Xi_{M^{2}}$ and
$\iota=R_{f,G_{N}},\;R^{\psi}_{f,G_{N}},\;R^{\psi,S}_{f,G_{N}}$ are the
interpolants constructed via the classical, the VSDK and the MVSDK. The shape
parameter $\varepsilon$ is chosen between $200$ equispaced values in the
interval $[0.01,50]$ via Leave-One-Out Cross Validation (LOOCV) [23]. Finally,
we let $N$ vary between $10$ and $500$ and we test two radial kernels:

$\begin{split}&\varphi_{W}(r)=(1-\varepsilon r)_{+}^{2}\quad\textrm{Wendland
$C^{0}$,}\\\ &\varphi_{M}(r)=e^{-\varepsilon r}(15+15\varepsilon
r+6(\varepsilon r)^{2}+(\varepsilon r)^{3})\quad\textrm{Matérn
$C^{6}$.}\end{split}$

In Figure 2, we show the behavior of the separation and fill distances, while
in Figure 3 we show the RMSEs achieved with both $\varphi_{W}$ and
$\varphi_{M}$.

Figure 2: The separation distance (left) and the fill distance (right) varying
$N$.

Figure 3: The RMSE obtained using $\varphi_{W}$ (left) and $\varphi_{M}$
(right) varying $N$.

The plots in Figure 2 display the benefits in employing $S$ in the MVSDK in
terms of diminished fill distance and increased separation distance. In Figure
3, we observe that MVSDKs are more effective in the case of the chosen Matérn
kernel. This is due to the fact that this kernel is more regular than the
chosen Wendland kernel, therefore it is more prone to provide an ill-
conditioned interpolation process. Furthermore, we remark that the sets of
nodes $G_{N}$, $N=10,\dots,500$ are not nested, and they are clustered around
the origin. Therefore, we can not expect an accurate recovering of the
function $f$, nor convergence increasing $N$.

##  6 Applications to the STIX Imaging Framework

In order to test the proposed MVSKs in real applied sciences, we focus on
solar hard X-ray imaging and, specifically, on the ESA STIX telescope [28], on
board of Solar Orbiter mission (see Figure 4). Hard X-ray telescopes provide
experimental measurements, named visibilities, of the Fourier transform of the
incoming photon flux at specific points of the spatial frequency plane. In the
case of STIX, $30$ subcollimators relying on the Moiré pattern technology
provide $N=60$ visibilities on 10 circles of the frequency plane with
increasing radii from about $2.79\times 10^{-3}$ arcsec-1 to $7.02\times
10^{-2}$ arcsec-1 (see Figure 4). We observe that the visibilities lying in
the lower half plane are obtained by reflecting the visibilities in the upper
half with respect to the origin.

Figure 4: The Spectrometer/Telescope for Imaging X-rays (STIX).

Figure 5: STIX visibilities: reflecting with respect to the origin leads to a
total number of $60$ visibilities.

We denote by ${\bf{f}}$ the vector whose components are the discretized values
of the incoming flux, by ${\bf{F}}$ the discretized Fourier transform sampled
at the set of points $\\{{\bf u}_{i}=(u_{i},v_{i})\\}_{i=1}^{N}$ in the
$({\mbox{u}},{\mbox{v}})$-plane and with ${\bf{V}}$ the vector whose $n$
components are the observed visibilities. Then, the image formation model in
this framework can be approximated by

${\bf{V}}={\bf{F}}{\bf{f}}~{}.$ (7)

###  6.1 The imaging process

Many inversion methods have been formulated to express the STIX observations
as images, see e.g. [2, 3, 5, 9, 24, 25, 18]. The approach that we propose
consists of two steps: interpolation of the visibilities so that we obtain the
visibility surfaces and the inversion of the so generated surfaces with rather
standard techniques. This idea was already used for the dismissed NASA
telescope Reuven Ramaty High Energy Solar Spectroscopic Imager (RHESSI) [22],
and its IDL implementation called uv$\\_$smooth [25], can be find in the NASA
Solar SoftWare (SSW) tree.

The uv$\\_$smooth code addresses equation (7) by means of an interpolation and
extrapolation procedure in which the interpolation step is carried out via an
algorithm based on spline functions and the extrapolation step is realized by
means of a soft-thresholding scheme [10, 1]. As far as the second step is
concerned, we will rely on the well-established projected Landweber iterative
method [27]. The first step instead will be carried out with our MVSKs.

###  6.2 MVSKs for STIX

To employ MVSKs, we need to define a scaling function $\psi$ and a mapping
$S$.

####  6.2.1 The choice of $\psi$.

In order to define the scaling function, we take advantage of a first
approximation of the inverse problem obtained via a standard back-projection
algorithm [26] that computes the discretized inverse Fourier transform of the
visibilities by means of the IDL source code vis_bpmap available in the NASA
SSW tree. The so-constructed image is then forward Fourier transformed to
obtain $\psi$. Once the interpolated visibility surface $\overline{\bf{V}}$
has been computed, the image reconstruction problem reads as follows

${\overline{\bf{V}}}={\overline{\bf{F}}}{\overline{\bf{f}}}~{},$ (8)

where ${\overline{\bf{F}}}$ is the $N^{2}\times N^{2}$ discretized Fourier
transform and ${\overline{\bf{f}}}$ is the $N^{2}\times 1$ vector to
reconstruct. In the following we will point out the advantages of
interpolating the visibilities with our technique.

####  6.2.2 The choice of $S$.

To present in details the mapping $S$ chosen for the STIX imaging framework,
let us first deepen the definition of the visibilities represented in Figure 5
(left). We have

$u_{i}=(L_{1}+L_{2})\frac{\cos(\alpha_{i}^{f})}{\rho_{i}^{f}}-L_{2}\frac{\cos(\alpha_{i}^{r})}{\rho_{i}^{r}},\;v_{i}=(L_{1}+L_{2})\frac{\sin(\alpha_{i}^{f})}{\rho_{i}^{f}}-L_{2}\frac{\sin(\alpha_{i}^{r})}{\rho_{i}^{r}},$

where $L_{1}=550$, $L_{2}=47$ and
$\alpha_{i}^{f},\rho_{i}^{f},\alpha_{i}^{r},\rho_{i}^{r}$ are discussed in
[19]. We define the map

$S(u,v)=C\log(\lVert(u,v)\lVert_{2})(\cos(\arctan(v/u)),\sin(\arctan(v/u))),$

where $C$ is a normalizing factor used to retain the order of magnitude of the
original visibilities after the mapping. In Figure 6, we can observe that the
mapped visibilities are distributed in a circular crown with no clustering
around the origin.

Figure 6: Mapped STIX visibilities.

###  6.3 Numerical results

On Jul 2022 STIX recorded a flare during the during the time interval
23:41:15–23:41:41 UT. The energy range of the event is 15–25 keV. In Figure 7
we reported the reconstruction carried out with the
interpolation/extrapolation algorithm where we respectively interpolate with:
classical radial kernels, VSKs and MVSKs. Moreover, as a further comparison we
also consider mem$\\_$ge [24] which is a well established implementation of
the maximum entropy approach and it is used by the solar physics community. We
note that all methods present artifacts but the shape of the source
reconstructed by interpolating with MVSKs is more similar to the one computed
with mem$\\_$ge. To have a quantitative feedback on the accuracy, we show in
Figure 8 the visibility fits obtained with the four different approaches. We
further observe that the chi squares values of mem$\\_$ge and uv$\\_$smooth +
MVSKs are similar and are the lowest.

Figure 7: Left to right, top to bottom: reconstruction of the flaring source
with: uv$\\_$smooth + classical radial kernels, uv$\\_$smooth + VSKs,
uv$\\_$smooth + MVSKs and mem$\\_$ge.

Figure 8: Top to bottom: visibility fits of the flaring source with:
uv$\\_$smooth + classical radial kernels, uv$\\_$smooth + VSKs, uv$\\_$smooth
+ MVSKs and mem$\\_$ge.

##  7 Conclusions

In this work, we define MVSKs by mixing VSKs with the FNA. After providing
some theoretical details, we performed some numerical simulations that showed
the advantages of using MVSDKs with respect to classical RBF kernels and
VSDKs. Then, we applied our method in the framework of hard X-ray astronomical
imaging. The obtained results encourage further investigations on this
research line.

## Acknowledgments

EP and AM are supported by the project: “Physics-based AI for predicting
extreme weather and spaceweather events (AIxtreme)”, Founded by Fondazione
Compagnia di San Paolo. FM acknowledges the financial support of the Programma
Operativo Nazionale (PON) "Ricerca e Innovazione" 2014 - 2020. The authors
kindly acknowledge the financial contribution from the AI-FLARES ASI-INAF
n.2018-16-HH.0 agreement. This research has been accomplished within GNCS-
IN$\delta$AM.

## References

  * [1] Allavena, S., Piana, M., Benvenuto, F., Massone, A.M.: An interpolation/extrapolation approach to X-ray imaging of solar flares. Inverse Probl. Imag. 6, 147 (2012)
  * [2] Aschwanden, M.J., Schmal, E., the RHESSI Team: Reconstruction of RHESSI solar flare images with a forward-fitting method. Sol. Phys. 210, 193 – 211 (2002)
  * [3] Benvenuto, F., Schwartz, R., Piana, M., Massone, A.M.: Expectation maximization for hard X-ray count modulation profiles. Astron. Astrophys. 555, A61 (2013)
  * [4] Berrut, J.P., De Marchi, S., Elefante, G., Marchetti, F.: Treating the Gibbs phenomenon in barycentric rational interpolation and approximation via the $S$-Gibbs algorithm. Appl. Math. Lett. 103, 106196, 7 (2020)
  * [5] Bonettini, S., Anastasia, C., Prato, M.: A new semiblind deconvolution approach for Fourier-based image restoration: An application in astronomy. SIAM J. Imaging Sci. 6, 1736–1757 (2013)
  * [6] Bozzini, M., Lenarduzzi, L., Rossini, M., Schaback, R.: Interpolation with variably scaled kernels. IMA J. Numer. Anal. 35(1), 199–219 (2015)
  * [7] Brutman, L.: Lebesgue functions for polynomial interpolation - a survey. Ann. Numer. Math. 4(1/4), 111–128 (1996)
  * [8] Campi, C., Marchetti, F., Perracchione, E.: Learning via variably scaled kernels. Adv. Comput. Math. 47(4), Paper No. 51, 23 (2021). https://doi.org/10.1007/s10444-021-09875-6, https://doi.org/10.1007/s10444-021-09875-6
  * [9] Cornwell, T., Evans, K.F.: A simple maximum entropy deconvolution algorithm. Astron. Astrophys. 143(1), 77–83 (1985)
  * [10] Daubechies, I., Defrise, M., De Mol, C.: An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Commun. Pure Appl. Math. 57(11), 1413–1457 (2004)
  * [11] De Marchi, S., Elefante, G., Marchetti, F.: Stable discontinuous mapped bases: the Gibbs-Runge-avoiding stable polynomial approximation (GRASPA) method. Comput. Appl. Math. 40(8), Paper No. 299, 17 (2021)
  * [12] De Marchi, S., Erb, W., Marchetti, F., Perracchione, E., Rossini, M.: Shape-driven interpolation with discontinuous kernels: Error analysis, edge extraction, and applications in magnetic particle imaging. SIAM J. Sci. Comput. 42(2), B472–B491 (2020)
  * [13] De Marchi, S., Marchetti, F., Perracchione, E.: Jumping with Variably Scaled Discontinuous Kernels (VSDKs). BIT Numer. Math. 60, 441–463 (2020)
  * [14] De Marchi, S., Marchetti, F., Perracchione, E., Poggiali, D.: Polynomial interpolation via mapped bases without resampling. J. Comput. Appl. Math. 364, 112347, 12 (2020)
  * [15] De Marchi, S., Marchetti, F., Perracchione, E., Poggiali, D.: Multivariate approximation at fake nodes. Appl. Math. Comput. 391, Paper No. 125628, 17 (2021)
  * [16] De Marchi, S., Lot, F., Marchetti, F., Poggiali, D.: Variably scaled persistence kernels (VSPKs) for persistent homology applications. Journal of Computational Mathematics and Data Science 4, 100050 (2022)
  * [17] Fasshauer, G.E.: Meshfree Approximation Methods with MATLAB. World Scientific (2007)
  * [18] Felix, S., Bolzern, R., Battaglia, M.: A compressed sensing-based image reconstruction algorithm for solar flare x-ray observations. Astrophys. J. 849(1), 10 (2017)
  * [19] Giordano, S., Pinamonti, N., Piana, M., Massone, A.M.: The process of data formation for the Spectrometer/Telescope for Imaging X-rays (STIX) in Solar Orbiter. SIAM J. Imaging Sci. 8(2), 1315–1331 (2015)
  * [20] Guastavino, S., Benvenuto, F.: Convergence rates of spectral regularization methods: a comparison between ill-posed inverse problems and statistical kernel learning. SIAM J. Numer. Anal. 58(6), 3504–3529 (2020)
  * [21] Karimnejad Esfahani, M., De Marchı, S., Marchetti, F.: Moving least squares approximation using variably scaled discontinuous weight function. Constructive Mathematical Analysis 6(1), 38 – 54 (2023)
  * [22] Lin, R.P., Dennis, B.R., Hurford, G.J., Smith, D.M., Zehnder, A., Harvey, P.R., Curtis, D.W., Pankow, D., Turin, P., Bester, M., Csillaghy, A., Lewis, M., Madden, N., Beek, H.F.V., Appleby, M., Raudorf, T., McTiernan, J., Ramaty, R., Schmahl, E., Schwartz, R., Krucker, S., Abiad, R., Quinn, T., Berg, P., Hashii, M., Sterling, R., Jackson, R., Pratt, R., Campbell, R.D., Malone, D., Landis, D., Barrington-Leigh, C.P., Slassi-Sennou, S., Cork, C., Clark, D., Amato, D., Orwig, L., Boyle, R., Banks, I.S., Shirey, K., Tolbert, A.K., Zarro, D., Snow, F., Thomsen, K., Henneck, R., Mchedlishvili, A., Ming, P., Fivian, M., Jordan, J., Wanner, R., Crubb, J., Preble, J., Matranga, M., Benz, A., Hudson, H., Canfield, R.C., Holman, G.D., Crannell, C., Kosugi, T., Emslie, A.G., Vilmer, N., Brown, J.C., Johns-Krull, C., Aschwanden, M., Metcalf, T., Conway, A.: The Reuven Ramaty High-Energy Solar Spectroscopic Imager (RHESSI). Sol. Phy. 210(1-2), 3 –32 (2002)
  * [23] Ling, L., Marchetti, F.: A stochastic extended Rippa’s algorithm for LpOCV. Applied Mathematics Letters 129, 107955 (2022)
  * [24] Massa, P., Schwartz, R., Tolbert, A.K., Massone, A.M., Dennis, B.R., Piana, M., Benvenuto, F.: MEM_GE: A new maximum entropy method for image reconstruction from solar X-ray visibilities. Astrophys. J. 894(1), 46 (2020)
  * [25] Massone, A.M., Emslie, A.G., Hurford, G.J., Prato, M., Kontar, E.P., Piana, M.: Hard X-ray imaging of solar flares using interpolated visibilities. Astrophys. J. 703, 2004–2016 (2009)
  * [26] Mersereau, R., Oppenheim, A.: Digital reconstruction of multidimensional signals from their projections. Proc. of IEEE 62(10) (1974)
  * [27] Piana, M., Bertero, M.: Projected Landweber method and preconditioning. Inverse Probl. 13(2), 441–463 (1997)
  * [28] S. Krucker, G. J. Hurford, O. Grimm, S. Kögl, H. P. Gröbelbauer, L. Etesi, D. Casadei, A. Csillaghy, A. O. Benz, N. G. Arnold, F. Molendini, P. Orleanski, D. Schori, H. Xiao, M. Kuhar, N. Hochmuth, S. Felix, F. Schramka, S. Marcin, S. Kobler, L. Iseli, M. Dreier, H. J. Wiehl, L. Kleint, M. Battaglia, E. Lastufka, H. Sathiapal, K. Lapadula, M. Bednarzik, G. Birrer, St. Stutz, Ch. Wild, F. Marone, K. R. Skup, A. Cichocki, K. Ber, K. Rutkowski, W. Bujwan, G. Juchnikowski, M. Winkler, Darmetko, M. Michalska, K. Seweryn, A. Bialek, P. Osica, J. Sylwester, M. Kowalinski, D. ´Scislowski, M. Siarkowski, M. Ste´slicki, T. Mrozek, P. Podgórski, A. Meuris, O. Limousin, O. Gevin, I. Le Mer, S. Brun, A. Strugarek, N. Vilmer, S. Musset, M. Maksimovi´c, F. Fárník, Z. Kozácek, J. Kasparová, G. Mann, H. Önel, A. Warmuth, J. Rendtel, J. Anderson, S. Bauer, F. Dionies, J. Paschke, D. Plüschke, M. Woche, F. Schuller, A. M Veronig, E. C. M. Dickson, P. T. Gallagher, S. A. Maloney, D. S. Bloomfield, M. Piana, A. M. Massone, F. Benvenuto, P. Massa, R. A. Schwartz, B. R. Dennis, H. F. van Beek, J. Rodríguez-Pacheco, R. P. Lin: The Spectrometer/Telescope for Imaging X-rays (STIX). Astron. Astrophys. 642, A15 (2020)
  * [29] Saitoh, S., Sawano, Y.: Theory of reproducing kernels and applications, Developments in Mathematics, vol. 44. Springer, Singapore (2016)
  * [30] Wendland, H.: Scattered Data Approximation, Cambridge Monographs on Applied and Computational Mathematics, vol. 17. Cambridge University Press, Cambridge (2005)