Jumping with variably scaled discontinuous kernels (VSDKs)

Abstract

In this paper we address the problem of approximating functions with discontinuities via kernel-based methods. The main result is the construction of discontinuous kernel-based basis functions. The linear spaces spanned by these discontinuous kernels lead to a very flexible tool which sensibly or completely reduces the well-known Gibbs phenomenon in reconstructing functions with jumps. For the new basis we provide error bounds and numerical results that support our claims. The method is also effectively tested for approximating satellite images.

1 Introduction

Radial Basis Function (RBF) methods (refer e.g. to [13, 14, 29, 32]) have become one of the most popular tools for solving multidimensional scattered data problems. Thanks to their independence from the mesh and to their easy implementation, they apply in a variety of fields, such as population dynamics, machine (deep) learning, solution of PDEs and image registration.

Even if such meshfree approaches have been extensively studied in the recent years, especially focusing on the efficiency and stability of the interpolant (cf. e.g. [3, 4, 10, 16]) not much effort has been addressed to construct robust approximants for functions with jumps. Indeed, infinitely smooth RBFs, such as Gaussians and Multiquadrics, theoretically show spectral accuracy, which is no longer preserved when interpolating functions with discontinuities. This fact, first observed in the context of truncated Fourier expansions and later used to characterize non-physical oscillations in approximating discontinuous functions, is known as Gibbs phenomenon.

To mitigate this effect for kernel-based approximation, one can use linear RBFs (see e.g. [15] for a general overview). This has been done in [20], where the Multiquadric has been replaced by the linear spline in regions around discontinuities. Alternatively, post-processing techniques, such as Gegenbauer reconstruction procedure [18] or digital total variation [28], are well-established tools. Finally, we point out that also the so-called Variably Scaled Kernels (VSKs) [1] are truly performing when reconstructing functions with gradient discontinuities, as shown also in [27].

Based on the last mentioned paper, and on the considerations about the usage of discontinuous bases for approximation purposes discussed in [30], here we propose a novel method which uses discontinuous kernels. The associated basis, constructed by means of what we call Variably Scaled Discontinuous Kernels (VSDKs), enables us to naturally reconstruct jump discontinuities (even with the family of Gaussians). The only drawback of the procedure lies in the fact that the algorithm needs to know where the discontinuities occur. To this aim, we recall that in [21] a one-dimensional edge detection method based on RBFs has been constructed and subsequently extended in [22] to a multidimensional framework by considering one-dimensional slices. Moreover, very recently in [26] an effective edge detector that analyzes the behaviour of the coefficients of the RBF approximant has been proposed (a similar idea was previously investigated in [20]). In that paper, a kernel-based discontinuos interpolant is empirically applied in the univariate setting. In this work, for edge detection we consider widely used schemes, such as Canny or Sobel edge detection (cf. [2, 31]).

After providing a theoretical analysis of the scheme in the one dimensional case, we extend the idea to higher dimensions and provide very general error bounds in terms of the well-known power function. Extensive numerical experiments are devoted to show the effectiveness of the method. We also provide the Matlab software, freely available for the scientific community at

http://www.math.unipd.it/~demarchi/RBF/CAARBF.html

that can be used to reproduce the tests presented in the paper. To conclude, we also investigate an application to the reconstruction of satellite images, where we deal with edges of irregular shapes; refer e.g. to [24].

The paper is organized as follows. In Sect. 2, we briefly review the main theoretical aspects of kernel-based approximation methods and introduce the VSKs. Section 3 presents our method for constructing VSDKs and in Sect. 4 we provide extensive numerical experiments. Applications to real world data are reported in Sect. 5. Conclusions and future works are discussed in Sect. 6.

2 Kernel-based approximation methods

Let ${{\mathscr {X}}} = \{ {\varvec{x}}_i, \; i = 1, \ldots , N\} \subset \varOmega$ be a set of distinct data points (data sites or nodes) arbitrarily distributed on a domain $\varOmega \subseteq {\mathbb {R}}^{d}$ and let ${{\mathscr {F}}} = \{ f_i = f({\varvec{x}}_i) , \; i=1, \ldots , N \}$ be an associated set of data values (measurements or function values) obtained by sampling some (unknown) function $f: \varOmega \longrightarrow {\mathbb {R}}$ at the nodes ${\varvec{x}}_i$. We can reconstruct f by interpolation, that is by finding a function $\mathscr {R}: \varOmega \longrightarrow {\mathbb {R}}$ satisfying the conditions

$$\begin{aligned} \mathscr {R}({\varvec{x}}_i)= f_i, \quad i=1,\ldots ,N. \end{aligned}$$

(2.1)

The interpolation problem (2.1) has a unique solution if $\mathscr {R} \in \text {span} \{ \varPhi _{\varepsilon }(\cdot ,{\varvec{x}}_i), {\varvec{x}}_i \in {{\mathscr {X}}}\}$, where $\varPhi _{\varepsilon } : \varOmega \times \varOmega \longrightarrow {\mathbb {R}}$ is a strictly positive definite and symmetric kernel and $\varepsilon >0$ is the so-called shape parameter. The resulting kernel-based interpolant, denoted by $\mathscr {R}_{\varepsilon ,{\mathscr {X}}},$ assumes the form

$$\begin{aligned} \mathscr {R}_{\varepsilon ,{\mathscr {X}}}({\varvec{x}}) = \sum _{k = 1}^N c_k \varPhi _{\varepsilon }({\varvec{x}}, {\varvec{x}}_k), \quad {\varvec{x}}\in \varOmega . \end{aligned}$$

(2.2)

We can write (2.1) by using the matrix $A_{\varepsilon } \in {\mathbb {R}}^{N \times N}$ which has entries $(A_{\varepsilon })_{ik}= \varPhi _{\varepsilon } ({\varvec{x}}_i , {\varvec{x}}_k),$$i,k=1, \ldots , N$. Then, letting ${\varvec{f}}= (f_1, \ldots , f_N)^{\intercal }$ the vector of data values, we can find the coefficients ${\varvec{c}}= (c_1, \ldots , c_N)^{\intercal }$ by solving the linear system $A_{\varepsilon } {\varvec{c}}= {\varvec{f}}$. Since we consider strictly positive definite and symmetric kernels, the existence and uniqueness of the solution of the linear system is ensured.

More precisely, we are interested in the class of strictly positive definite and symmetric radial kernels$\varPhi _{\varepsilon }$ defined as follows.

Definition 2.1

$\varPhi _{\varepsilon }$ is called radial kernel if there exists a continuous function $\varphi _{\varepsilon }: [0,+\infty )\longrightarrow {\mathbb {R}}$, depending on the shape parameter $\varepsilon >0$, such that

$$\begin{aligned} \varPhi _{\varepsilon }({\varvec{x}},{\varvec{y}})=\varphi _{\varepsilon }(\Vert {\varvec{x}}-{\varvec{y}}\Vert _2), \end{aligned}$$

(2.3)

for all ${\varvec{x}}, {\varvec{y}} \in \varOmega$.

From (2.3) it follows that if $\varPhi _{\varepsilon }$ is radial, then it is completely identified by the function $\varphi _{\varepsilon }$ and we can indifferently use $\varPhi _{\varepsilon }$ or $\varphi _{\varepsilon }$ for denoting the interpolant in (2.2).

To $\varPhi _{\varepsilon }$ we associate a real pre-Hilbert space $H_{\varPhi _{\varepsilon }}(\varOmega )$ with reproducing kernel $\varPhi _{\varepsilon }$

$$\begin{aligned} H_{\varPhi _{\varepsilon }}(\varOmega )=\text {span} \{ \varPhi _{\varepsilon }(\cdot ,{\varvec{x}}),\;{\varvec{x}} \in \varOmega \}, \end{aligned}$$

equipped with the bilinear form $\left( \cdot ,\cdot \right) _{H_{\varPhi _{\varepsilon }}(\varOmega )}$. We then define the native space${{\mathscr {N}}}_{\varPhi _{\varepsilon }} (\varOmega )$ of $\varPhi _{\varepsilon }$ as the completion of $H_{\varPhi _{\varepsilon }}(\varOmega )$ with respect to the norm $||\cdot ||_{H_{\varPhi _{\varepsilon }}(\varOmega )}$, that is $||f||_{H_{\varPhi _{\varepsilon }}(\varOmega )} = ||f||_{{{\mathscr {N}}}_{\varPhi _{\varepsilon }}(\varOmega )}$ for all $f \in H_{\varPhi _{\varepsilon }}(\varOmega )$ (for details see the monographs [14, 32]).

The accuracy of the interpolation process is usually expressed in terms of the power function. Let $A_{\varepsilon }(\mathscr {X})$ be the interpolation matrix related to the set of nodes $\mathscr {X}$ and to the kernel $\varPhi _\varepsilon$. Also let $A_{\varepsilon }(\mathscr {Y})$ be the matrix associated to the augmented set $\mathscr {Y}:= \{\varvec{x}\}\cup \mathscr {X},\; \varvec{x} \in \varOmega$. The power function is a positive function obtained by the ratio of two determinants (cf. [5, 11])

$$\begin{aligned} P_{\varPhi _\varepsilon ,{{\mathscr {X}}}} (\varvec{x}):=\sqrt{\text {det}(A_{\varepsilon }(\mathscr {Y})) \over \text {det}(A_{\varepsilon }(\mathscr {X}))}\,. \end{aligned}$$

(2.4)

The following pointwise error bound, that uses the power function and the norm of f in the native space, holds (see e.g. [14, Th. 14.2, p.117]).

Theorem 2.1

Let $\varPhi _{\varepsilon }\in C(\varOmega \times \varOmega )$ be a strictly positive definite kernel and ${{\mathscr {X}}} = \{{\varvec{x}}_i, i=1, \ldots , N \} \subseteq \varOmega$ a set of distinct points. For all $f \in {{\mathscr {N}}}_{\varPhi _{\varepsilon }}(\varOmega )$

$$\begin{aligned} |f\left( {\varvec{x}}\right) -\mathscr {R}_{\varepsilon ,{\mathscr {X}}}\left( {\varvec{x}}\right) | \le P_{\varPhi _{\varepsilon },{{\mathscr {X}}}}({\varvec{x}}) ||f||_{{{\mathscr {N}}}_{\varPhi _{\varepsilon }}(\varOmega )}, \quad {\varvec{x}} \in \varOmega . \end{aligned}$$

Note that this theorem bounds the pointwise error in terms of the power function which depends on the kernel and on the data points but is independent of the function values.

2.1 From RBF to VSK interpolation

As well-known, the choice of the shape parameter $\varepsilon$ is a crucial computational issue for RBF interpolation. If it is not properly chosen we might see instability effects. To overcome such problems for the Gaussian kernel, a solution is provided by the so-called RBF-QR method which is truly effective (see e.g. [17, 23]). An alternative, which can be applied to any kernel, is the use of Variably Scaled Kernels (or VSKs), introduced in [1]. We also notice, that VSK have been successfully used to reconstruct functions with gradient discontinuities in [27].

Starting from this idea, we analyse the case of jump discontinuities and we introduce a new family of kernels that we call Variably Scaled Discontinuous Kernels or simply VSDKs.

Considering VSKs, the classical tuning strategy of finding the optimal shape parameter might be substituted by the choice of a scale function which plays the role of a density function. More precisely (cf. [1, Def. 2.1]):

Definition 2.2

Letting $\mathscr {I}\subseteq (0,+\infty )$ and $\varPhi _{\varepsilon }$ a positive definite radial kernel on $\varOmega \times \mathscr {I}$ depending on the shape parameter $\varepsilon >0$. Given a scale function $\psi :\varOmega \longrightarrow \mathscr {I}$, we define a VSK $\varPhi _{\psi }$ on $\varOmega$ as

$$\begin{aligned} \varPhi _{\psi }(\varvec{x},\varvec{y}):= \varPhi _1((\varvec{x},\psi (\varvec{x})),(\varvec{y},\psi (\varvec{y}))), \end{aligned}$$

(2.5)

for $\varvec{x},\varvec{y}\in \varOmega$.

Defining then the map $\varPsi (\varvec{x})=(\varvec{x},\psi (\varvec{x}))$ on $\varOmega$, the interpolant on the set of nodes $\varPsi (\mathscr {X}):=\{(\varvec{x}_k,\psi (\varvec{x}_k)),\;\varvec{x}_k\in \mathscr {X}\}$ with fixed shape parameter $\varepsilon =1$ takes the form

$$\begin{aligned} \mathscr {R}_{1,\varPsi (\mathscr {X})}(\varPsi ({\varvec{x}})) = \sum _{k = 1}^N {c_k \varPhi _{1}(\varPsi ({\varvec{x}}),\varPsi ( {\varvec{x}}_k))}, \end{aligned}$$

(2.6)

with ${\varvec{x}}\in \varOmega ,\;{\varvec{x}}_k\in \mathscr {X}$.

From now on we take $\varepsilon =1$, as in [1, Def. 2.1]. Hence, if otherwise state, we omit to indicate the shape parameter.

By analogy with the interpolant in (2.2), the vector of coefficients ${\varvec{c}}= (c_1, \ldots , c_N)^{\intercal }$ in (2.6) is determined by solving the linear system $A_{\psi }{\varvec{c}}= {\varvec{f}}$, where $(A_{\psi })_{ik}= \varPhi _{\psi } ({\varvec{x}}_i , {\varvec{x}}_k)$ and $\varvec{f}$ is the vector of data values.

Once we have the interpolant $\mathscr {R}_{\varPsi (\mathscr {X})}$ on $\varOmega \times \mathscr {I}$, we can project back on $\varOmega$ the points $(\varvec{x},\psi (\varvec{x}))\in \varOmega \times \mathscr {I}$. In this way, we obtain a VSK interpolant $\mathscr {V}_{\psi }$ on $\varOmega$ that is, using (2.5),

$$\begin{aligned} \mathscr {V}_{\psi }(\varvec{x}):=\sum _{k = 1}^N {c_k \varPhi _{\psi }({\varvec{x}}, {\varvec{x}}_k)}= \sum _{k = 1}^N {c_k \varPhi (\varPsi ({\varvec{x}}),\varPsi ( {\varvec{x}}_k))}=\mathscr {R}_{\varPsi (\mathscr {X})}(\varPsi ({\varvec{x}})).\qquad \end{aligned}$$

(2.7)

The error and stability analysis of this varying scale process on $\varOmega$ coincides with the analysis of a fixed scale kernel on $\varOmega \times \mathscr {I}$. For details and analysis of these continuous scale functions, we refer the reader to [1].

3 Variably scaled discontinuous kernels

We introduce the VSDKs by considering the one dimensional case and observing that the extension to the multidimensional case is almost straightforward, as we will show later in Sect. 3.2.

Let $\varOmega =(a,b)\subset {{\mathbb {R}}}$ be an open interval and let $\xi \in \varOmega$. We consider the discontinuous function $f:\varOmega \longrightarrow {{\mathbb {R}}}$

$$\begin{aligned} f(x):= \left\{ \begin{array}{ll} f_1(x), &{} \quad a< x<\xi ,\\ f_2(x), &{} \quad \xi \le x<b,\\ \end{array}\right. \end{aligned}$$

where $f_1, f_2$ are real valued smooth functions such that $\displaystyle {\lim _{x\rightarrow a^+}{f_1(x)}}$ and $\displaystyle {\lim _{x\rightarrow b^-}{f_2(x)}}$ exist finite and

$$\begin{aligned} f_2(\xi )\ne \lim _{x\rightarrow \xi }{f_1(x)}\,. \end{aligned}$$

Our aim consists in approximating the function f on the set of nodes ${{\mathscr {X}}} \subset \varOmega$. Unfortunately the presence of jumps is the cause of oscillations in the reconstructing process.

To approximate f on ${\mathscr {X}}$ we take interpolants of the form (2.7) and we consider discontinuous scale functions in the interpolation process.

Let $\alpha ,\beta \in {{\mathbb {R}}},\;\alpha \ne \beta$ and $\mathscr {S}= \left\{ \alpha ,\beta \right\}$. We propose the following scale function $\psi :\varOmega \longrightarrow \mathscr {S}$ defined as:

$$\begin{aligned} \psi (x):= \left\{ \begin{array}{ll} \alpha ,\quad &{}\quad x<\xi ,\\ \beta ,\quad &{}\quad x\ge \xi .\\ \end{array} \right. \end{aligned}$$

The function $\psi$ is piecewise constant, having a jump discontinuity at $\xi$ as the function f. Let $\varPhi _{\varepsilon }$ be a positive definite radial kernel on $\varOmega \times \mathscr {S}$, possibly depending on a shape parameter $\varepsilon >0$ or alternatively a variably scaled kernel $\varPhi _{\psi }$ on $\varOmega$ as in (2.5). We now analyze the function $\varphi$ related to the kernel $\varPhi$ with the fixed shape parameter $\varepsilon =1$, that is

$$\begin{aligned} \begin{array}{ll} \varphi (\Vert \varPsi (x)-\varPsi (y)\Vert _2) &{} = \varphi (\Vert (x,\psi (x))-(y,\psi (y))\Vert _2)\\ &{} = \varphi \left( \sqrt{(x-y)^2+(\psi (x)-\psi (y))^2}\right) \,. \end{array} \end{aligned}$$

This implies

$$\begin{aligned} \varphi (\Vert \varPsi (x)-\varPsi (y)\Vert _2)= \left\{ \begin{array}{llll} \varphi (|x-y|),\quad &{}\quad x,y<\xi &{} \text { or } &{} x,y\ge \xi ,\\ \varphi (\Vert (x,\alpha )-(y,\beta )\Vert _2),\quad &{}\quad x<\xi \le y &{} \text { or } &{} y<\xi \le x,\\ \end{array} \right. \end{aligned}$$

since $\varphi (\Vert (x,\alpha )-(y,\beta )\Vert _2)=\varphi (\Vert (x,\beta )-(y,\alpha )\Vert _2).$

The so-constructed interpolant $\mathscr {V}_{\psi }:\varOmega \longrightarrow {{\mathbb {R}}}$ on the set $\mathscr {X}=\{ x_k,\;k=1,\dots ,N\}$ is a discontinuous linear combination of functions $\varPhi _{\psi }(\cdot , x_k)$ so defined:

• if $a<x_k<\xi$

  $$\begin{aligned} \varPhi _{\psi }(x, x_k)= \left\{ \begin{array}{ll} \varphi (|x-x_k|), &{}\quad \quad x<\xi ,\\ \varphi (\Vert (x,\alpha )-(x_k,\beta )\Vert _2), &{}\quad \quad x\ge \xi ,\\ \end{array} \right. \end{aligned}$$

• if $\xi \le x_k<b$

  $$\begin{aligned} \varPhi _{\psi }(x, x_k)= \left\{ \begin{array}{ll} \varphi (|x-x_k|), &{}\quad \quad x\ge \xi ,\\ \varphi (\Vert (x,\alpha )-(x_k,\beta )\Vert _2), &{}\quad \quad x<\xi .\\ \end{array} \right. \end{aligned}$$

Therefore, the interpolant $\mathscr {V}_{\psi }$ is a linear combination of functions having a discontinuity at $\xi$. We can easily generalize this procedure for a set of distinct discontinuity points on $\varOmega$; see the next section.

3.1 VSDKs: one dimensional case

To generalize the discussion carried out above, we need the following definition.

Definition 3.1

Let $\varOmega =(a,b)\subset {{\mathbb {R}}}$ be an open interval, $\mathscr {S}=\{\alpha ,\beta \}$ with $\alpha ,\beta \in {{\mathbb {R}}}_{>0},\;\alpha \ne \beta$ and let $\mathscr {D}=\{\xi _j,\;j=1,\dots ,\ell \}\subset \varOmega$ be a set of distinct points such that $\xi _j<\xi _{j+1}$ for every j. Let $\psi :\varOmega \longrightarrow \mathscr {S}$ be defined as

$$\begin{aligned} \psi (x):= \left\{ \begin{array}{lll} \alpha , &{} \quad x\in (a,\xi _1)\text { or }x\in [\xi _j,\xi _{j+1}), &{} \text{ where } j \text { is even},\\ \beta , &{} \quad x\in [\xi _j,\xi _{j+1}), &{} \text{ where } j \text { is odd},\\ \end{array} \right. \end{aligned}$$

and

$$\begin{aligned} {\psi (x)}|_{[\xi _{\ell },b)}:= \left\{ \begin{array}{ll} \alpha , &{} \quad \ell \text{ is } \text{ even }, \\ \beta , &{} \quad \ell \text{ is } \text{ odd }. \end{array} \right. \end{aligned}$$

With this choice of the scale function $\psi$ and similarly to (2.5), we call the kernel $\varPhi _{\psi }$ a VSDK on $\varOmega$.

For the analysis of the VSDKs introduced in Definition 3.1 we cannot rely on some important and well-known results of RBF interpolation. Therefore, before stating upper bounds for the VSDK interpolants in terms of the power function, we give a preliminary analysis.

Let $\varOmega$ and $\mathscr {D}$ be as in Definition 3.1 and $n\in {{\mathbb {N}}}$. We define $\psi _n:\varOmega \longrightarrow \mathscr {I} \subseteq (0,+ \infty )$ as

$$\begin{aligned}&\psi _n(x):= \left\{ \begin{array}{lll} \alpha ,\quad &{} x\in (a,\xi _1-1/n), \text { or }x\in [\xi _j+1/n,\xi _{j+1}-1/n), &{} j \text{ is } \text{ even },\\ \beta ,\quad &{} x\in [\xi _j+1/n,\xi _{j+1}-1/n) &{} j \text{ is } \text{ odd },\\ \gamma _1(x),\quad &{} x\in [\xi _j-1/n,\xi _{j}+1/n), &{} j \text{ is } \text{ odd },\\ \gamma _2(x),\quad &{} x\in [\xi _j-1/n,\xi _{j}+1/n), &{} j \text{ is } \text{ even }, \end{array} \right. \nonumber \\&{\psi _n(x)}|_{[\xi _{\ell }+1/n,b)}:= \left\{ \begin{array}{ll} \alpha ,\quad &{} \ell \text{ is } \text{ even },\\ \beta ,\quad &{} \ell \text{ is } \text{ odd }, \end{array} \right. \end{aligned}$$

(3.1)

where $\gamma _1,\gamma _2$ are continuous, strictly monotone functions so that

$$\begin{aligned} \lim _{x\rightarrow \xi _{j+1}+1/n}{\gamma _1(x)}=\gamma _2(\xi _j-1/n)=\beta ,\quad \lim _{x\rightarrow \xi _{j+1}+1/n}{\gamma _2(x)}=\gamma _1(\xi _j-1/n)=\alpha . \end{aligned}$$

From Definition 3.1, it is straightforward to verify that $\forall x\in \varOmega$ the following pointwise convergence result holds

$$\begin{aligned} \lim _{n\rightarrow \infty }{\psi _n(x)}=\psi (x). \end{aligned}$$

We point out that for every fixed $n\in {{\mathbb {N}}}$ the kernel $\varPhi _{\psi _n}$ is a continuous VSK, hence it satisfies the error bound of Theorem 2.1. For VSDKs instead we have the following result.

Theorem 3.1

For every $x,y\in \varOmega$, we have

$$\begin{aligned} \lim _{n\rightarrow \infty }{\varPhi _{\psi _n}(x,y)}=\varPhi _{\psi }(x,y), \end{aligned}$$

where $\varPhi _{\psi }$ is the kernel considered in Definition 3.1.

Proof

Let us consider the map $\varPsi _n(x)=(x,\psi _n(x))$ on $\varOmega$. We can write

$$\begin{aligned} \lim _{n\rightarrow \infty }{\varPhi _{\psi _n}(x,y)}=\lim _{n\rightarrow \infty }{\varPhi (\varPsi _n(x),\varPsi _n(y))}=\lim _{n\rightarrow \infty }{\varphi (\Vert \varPsi _n(x)-\varPsi _n(y)\Vert _2)}. \end{aligned}$$

Recalling (3.1), we get

$$\begin{aligned} \begin{aligned} \lim _{n\rightarrow \infty }{\varphi (\Vert \varPsi _n(x)-\varPsi _n(y)\Vert _2)}&=\varphi \big (\lim _{n\rightarrow \infty }{\Vert \varPsi _n(x)-\varPsi _n(y)\Vert _2}\big )\\&=\varphi (\Vert \varPsi (x)-\varPsi (y)\Vert _2)\\&=\varPhi _{\psi }(x,y). \end{aligned} \end{aligned}$$

This concludes the proof. $\square$

Corollary 3.1

Let $H_{\varPhi _{\psi _n}}(\varOmega )=\text {span} \{ \varPhi _{\psi _n}(\cdot ,x),\;x\in \varOmega \}$ be equipped with the bilinear form $\left( \cdot ,\cdot \right) _{H_{\varPhi _{\psi _n}}(\varOmega )}$ and let ${{\mathscr {N}}}_{\varPhi _{\psi _n}} (\varOmega )$ be the related native space. Then, taking the limit of the basis functions, we obtain the space $H_{\varPhi _{\psi }}(\varOmega )=\text {span} \{ \varPhi _{\psi }(\cdot ,x),\;x\in \varOmega \}$ equipped with the bilinear form $\left( \cdot ,\cdot \right) _{H_{\varPhi _{\psi }}(\varOmega )}$ and the related native space ${{\mathscr {N}}}_{\varPhi _{\psi }}(\varOmega )$.

Proof

If $f\in H_{\varPhi _{\psi }}(\varOmega )$, then it can be expressed as a linear combination of basis functions $\varPhi _{\psi }(\cdot ,x),\;x\in \varOmega$. From Theorem 3.1, we get that for every $x\in \varOmega$

$$\begin{aligned} \lim _{n\rightarrow \infty }{\varPhi _{\psi _n}(\cdot ,x)}=\varPhi _{\psi }(\cdot ,x), \end{aligned}$$

and so f is also a linear combination of the functions $\lim _{n\rightarrow \infty }{\varPhi _{\psi _n}(\cdot ,x)},\;x\in \varOmega$, as required. $\square$

We get an immediate consequence for the interpolant $\mathscr {V}_{\psi }$ too.

Corollary 3.2

Let $\varOmega$, $\mathscr {S}$ and $\mathscr {D}$ be as in Definition 3.1. Let $f:\varOmega \longrightarrow {{\mathbb {R}}}$ be a discontinuous function whose step discontinuities are located at the points belonging to $\mathscr {D}$. Moreover, let $\psi _n$ and $\psi$ be as in Theorem 3.1. Then, considering the interpolation problem with nodes $\mathscr {X}=\{ x_k,\;k=1,\dots ,N\}$ on $\varOmega$, we have

$$\begin{aligned} \lim _{n\rightarrow \infty }{\mathscr {V}_{\psi _n}(x)}=\mathscr {V}_{\psi }(x), \end{aligned}$$

for every $x\in \varOmega$.

Proof

Since $\mathscr {V}_{\psi }$ is a linear combination of the basis functions, the thesis follows from Theorem 3.1 and Corollary 3.1. $\square$

To provide error bounds, we now only need to introduce the power function for a VSDK $\varPhi _{\psi }$ on the set of nodes $\mathscr {X}$. From (2.4), we know that it is defined as

$$\begin{aligned} P_{\varPhi _{\psi },{{\mathscr {X}}}}(x)=\sqrt{\frac{\text {det}(A_{\psi }(\mathscr {Y}))}{\text {det}(A_{\psi }(\mathscr {X}))}}. \end{aligned}$$

From Theorem 3.1 and Corollary 3.1, it easily follows that $\forall x\in \varOmega$

$$\begin{aligned} P_{\varPhi _{\psi },{{\mathscr {X}}}}(x)=\lim _{n\rightarrow \infty }{P_{\varPhi _{\psi _n},{{\mathscr {X}}}}(x)}. \end{aligned}$$

These results allow to state an error bound for interpolation via VSDKs.

Proposition 3.1

Let $\varPhi _{\psi }$ be a VSDK on $\varOmega =(a,b)\subset {{\mathbb {R}}}$. Suppose that ${{\mathscr {X}}} = \{x_i, i=1, \ldots , N \} \subseteq \varOmega$ have distinct points. For all $f \in {{\mathscr {N}}}_{\varPhi _{\psi }}(\varOmega )$ we have

$$\begin{aligned} |f(x)-\mathscr {V}_{\psi }(x)| \le P_{\varPhi _{\psi },{{\mathscr {X}}}}(x) \Vert f\Vert _{{{\mathscr {N}}}_{\varPhi _{\psi }}(\varOmega )}, \quad x \in \varOmega . \end{aligned}$$

Proof

For every $n\in {{\mathbb {N}}}$ and $x\in \varOmega$, since the VSK $\varPhi _{\psi _n}$ is continuous, we know that (see Theorem 2.1)

$$\begin{aligned} |f(x)-\mathscr {V}_{\psi _n}(x)| \le P_{\varPhi _{\psi _n},{{\mathscr {X}}}}(x) \Vert f\Vert _{{{\mathscr {N}}}_{\varPhi _{\psi _n}}(\varOmega )}. \end{aligned}$$

Then, taking the limit $n\rightarrow \infty$ and recalling the results of this subsection, the thesis follows. $\square$

Proposition 3.1, as the classical bound for the RBF interpolants, limits the error in terms of the power function and consequently takes into account both the kernel and data.

3.2 VSDKs: multidimensional case

The VSDKs rely upon the classical RBF bases and therefore in principle they are suitable to be implemented in any dimension. However, since the geometry is more complex than in 1D, we need to carefully define the scale function $\psi$.

Let $\varOmega \subset {{\mathbb {R}}}^d$ be an open subset with Lipschitz boundary. In our discussion, we consider step functions $f:\varOmega \longrightarrow {{\mathbb {R}}}$ such that there exists a disjoint partition $\mathscr {P}=\{R_1,\ldots ,R_m\}$ of regions having Lipschitz boundaries. That is all the jumps of f lie along $(d-1)$-dimensional manifolds $\gamma _1,\dots ,\gamma _p$ such that

$$\begin{aligned} \gamma _i\subseteq \bigcup _{i=1}^{m} {\partial R_{i}}{\setminus }\partial \varOmega , \;\; \forall i=1,\dots ,p. \end{aligned}$$

Then, a suitable scale function $\psi$ for interpolating f via VSDKs can be defined as follows.

Definition 3.2

Let $\varOmega \subset {{\mathbb {R}}}^d$ be an open subset with Lipschitz boundary, $\mathscr {S}=\{\alpha _1,\ldots ,\alpha _m\}$ real distinct values and $\mathscr {P}=\{R_1, \ldots , R_m\}$ a partition of $\varOmega$ whose elements are regions having Lipschitz boundaries. Define $\psi :\varOmega \longrightarrow \mathscr {S}$ as

$$\begin{aligned} {\psi (\varvec{x})}|_{R_i}:= \alpha _i. \end{aligned}$$

With this choice of the scale function $\psi$ and referring to (2.5), we call again the kernel $\varPhi _{\psi }$ a VSDK on $\varOmega$.

Remark 3.1

In Definition 3.2 we choose a scale function which emulates the properties of the one-dimensional function of Definition 3.1. The difference is that the multidimensional $\psi$ could be discontinuous not exclusively at the same points as f, but also at other nodes. Precisely, if we are able to choose $\mathscr {P}$ so that

$$\begin{aligned} \bigcup _{i=1}^{p}{\gamma _i}=\bigcup _{i=1}^{m} {\partial R_{i}}{{\setminus }}\partial \varOmega , \end{aligned}$$

then f and $\psi$ have the same discontinuities. Otherwise, if

$$\begin{aligned} \bigcup _{i=1}^{p}{\gamma _i}\subset \bigcup _{i=1}^{m} {\partial R_{i}}{\setminus }\partial \varOmega , \end{aligned}$$

then $\psi$ is discontinuous along $\bigcup _{i=1}^{m} {\partial R_{i}}{\setminus }\big (\partial \varOmega \cup \bigcup _{i=1}^{p}{\gamma _i}\big )$, while f is not.

The theoretical analysis in the multidimensional case is done along the same path of the one-dimensional setting. Indeed, we consider continuous scale functions $\psi _n:\varOmega \longrightarrow \mathscr {I} \subseteq (0,+ \infty )$ such that $\forall \varvec{x}\in \varOmega$,

$$\begin{aligned} \lim _{n\rightarrow \infty }{\psi _n(\varvec{x})}=\psi (\varvec{x}), \end{aligned}$$

and

$$\begin{aligned} \lim _{n\rightarrow \infty }{\mathscr {V}_{\psi _n}(\varvec{x})}=\mathscr {V}_{\psi }(\varvec{x}), \end{aligned}$$

for every $\varvec{x}\in \varOmega$.

We omit the easy extension of all results discussed in Sect. 3.1 and we state directly the error bound.

Proposition 3.2

Let $\varPhi _{\psi }$ be a VSDK as in Definition 3.2. Suppose that ${{\mathscr {X}}} = \{\varvec{x}_i, i=1, \ldots , N \} \subseteq \varOmega$ have distinct points. For all $f \in {{\mathscr {N}}}_{\varPhi _{\psi }}(\varOmega )$ we have

$$\begin{aligned} |f(\varvec{x})-\mathscr {V}_{\psi }(\varvec{x})| \le P_{\varPhi _{\psi },{{\mathscr {X}}}}(\varvec{x}) \Vert f\Vert _{{{\mathscr {N}}}_{\varPhi _{\psi }}(\varOmega )}, \quad \varvec{x} \in \varOmega . \end{aligned}$$

Proof

Refer to Proposition 3.1 and to the remarks made in this section. $\square$

4 Numerical experiments

We consider three strictly positive definite RBFs having different regularities

$$\begin{aligned} \begin{array}{ll} \varphi ^1_{\varepsilon }(r)=e^{-\varepsilon r},\quad &{}\textit{Matern } C^0,\\ \varphi ^2_{\varepsilon }(r)=e^{-\varepsilon r} \left( 15+15\varepsilon r+6\varepsilon ^2 r^2+\varepsilon ^2 r^3 \right) ,\quad &{} \textit{Matern } C^6,\\ \varphi ^3_{\varepsilon }(r)=e^{-\varepsilon ^2 r^2},\quad &{} \textit{Gaussian }C^{\infty }.\\ \end{array} \end{aligned}$$

(4.1)

To point out the accuracy of our tests, both in 1D and 2D cases, we refer to the Maximum Absolute Error (MAE) and to the Root Mean Square Error (RMSE). Once the interpolant is constructed, we evaluate it on a grid of S evaluation points $\{\varvec{z}_1, \ldots , \varvec{z}_S\}$ so that

$$\begin{aligned} \mathrm{MAE} = \max _{1\le i \le S} |f(\varvec{z}_i)-\mathscr {V}_{\psi }(\varvec{z}_i)|, \qquad \mathrm{RMSE} = \sqrt{\frac{1}{S}\sum _{i=1}^{S} \left( f(\varvec{z}_i)-\mathscr {V}_{\psi }(\varvec{z}_i)\right) ^2}. \end{aligned}$$

Further, we also estimate the upper bound for the error by computing the Maximum of the Power Function (MPF)

$$\begin{aligned} \mathrm{MPF} = \max _{1\le i \le S} P_{\varPhi _{\psi },{{\mathscr {X}}}}(\varvec{z}_i). \end{aligned}$$

In the 2D case we also deal with grey-scale images and thus we consider the Structural Similarity Index (SSIM), which is a well-known parameter that indicates the similarity between two images. The SSIM index lies in the interval [0, 1] (the value 1 corresponds to identical images).

The numerical experiments have been carried out with Matlab on an Intel(R) Core(TM) i5-4200U CPU @ 2.30 GHz processor.

4.1 A toy example

Let $\varOmega = (-1,1)$,

$$\begin{aligned} f_1(x)= \left\{ \begin{array}{ll} e^{-x}, &{} \quad -1< x<-0.5,\\ x^3, &{} \quad -0.5\le x<0.5,\\ 1, &{} \quad \; \; \; 0.5\le x<1,\\ \end{array} \right. \end{aligned}$$

and

$$\begin{aligned} \mathscr {X}=\left\{ x_j=-1+\dfrac{(j-1)}{39},\;j=1,\dots ,79 \right\} . \end{aligned}$$

We evaluate then the interpolant on a grid of equispaced points on $\varOmega$ with step size 5.00E−4.

First, as described in (2.2), we interpolate the function $f_1$ via classical RBF interpolation on $\mathscr {X}$, using the kernel functions reported in (4.1). For such RBFs we select the optimal shape parameter $\varepsilon ^*$ via Leave One Out Cross Validation (LOOCV) (refer to [13] for a general overview). The resulting interpolants are plotted in Fig. 1.

Fig. 1

The function $f_1$ and the classical interpolants on $\mathscr {X}$ obtained using different kernel functions

In Table 1 we report the values of MAE and RMSE with respect to $f_1$ and the maximum values of the power function evaluated at the same set of nodes. The errors are similar except for $\varphi ^1_{\varepsilon }$ that turns out to be more stable. Indeed, as expected, the corresponding reconstruction is less affected by the Gibbs phenomenon, due to the poor regularity of $\varphi ^1_{\varepsilon }$.

The function $f_1$ presents two discontinuities at $\xi _1 = -0.5$ and $\xi _2=0.5$. Taking into account (3.1), we define the scale function for VSKs as

$$\begin{aligned} \psi _n(x):= \left\{ \begin{array}{ll} 1, \quad &{} x\in (-1,\xi _1-1/n)\text { or }x\in [\xi _2+1/n,1),\\ 2, \quad &{} x\in [\xi _1+1/n,\xi _{2}-1/n),\\ (nx-\xi _1n+3)/2, \quad &{} x \in [\xi _1-1/n,\xi _1+1/n),\\ (-nx+\xi _2n+3)/2, \quad &{} x\in [\xi _2-1/n,\xi _2+1/n). \end{array} \right. \end{aligned}$$

(4.2)

Moreover, considering the pointwise limit as $n\rightarrow \infty$ of $\psi _n(x)$ we consider for VSDKs the discontinuous scale function:

$$\begin{aligned} \psi (x):= \left\{ \begin{array}{ll} 1, \quad &{} x\in (-1,\xi _1)\text { or }x\in [\xi _2,1),\\ 2, \quad &{} x\in [\xi _1,\xi _{2}). \end{array} \right. \end{aligned}$$

(4.3)

Table 1 Results of classical RBF reconstruction for $f_1$

In Fig. 2, we plot the scale functions $\psi _n$ for some values of n and $\psi$.

Fig. 2

The continuous scale functions $\psi _{10},\;\psi _{50},\;\psi _{500}$ and the discontinuous one $\psi$; see (4.2) and (4.3)

Finally, for each of the RBFs considered in (4.1) we first take the scale function $\psi _n$ and compute the VSK interpolant for $f_1$ by considering increasing values of n ($n=10,50,500$). Then we approximate $f_1$ using the VSDK determined by $\psi$ (see formulae (4.2) and (4.3)). The graphical results for the three variably scaled (discontinuous) kernels are reported in Figs. 3, 4 and 5, while the accuracy indicators are shown in Tables 2, 3 and 4.

From the figures, we can graphically note how the reconstruction via VSDKs is indeed the limit of the continuous case. As expected the $C^0$ RBF combined with the VSKs is not truly affected by the Gibbs phenomenon. Using the other kernels, we note that such oscillations are progressively reduced as n increases and graphically disappear when using VSDKs.

Fig. 3

VSK and VSDK reconstructions of $f_1$ on $\mathscr {X}$ using $\varphi ^1$

Fig. 4

VSK and VSDK reconstructions of $f_1$ on $\mathscr {X}$ using $\varphi ^2$

Fig. 5

VSK and VSDK reconstructions of $f_1$ on $\mathscr {X}$ using $\varphi ^3$

From the results reported in this subsection it is evident that the maximum value of the power function usually decreases as n increases. However, there are no theoretical results about the fact that the power function for VSDKs assumes smaller values than the one of VSKs and this is confirmed numerically. Indeed, in some cases the maximum value attained by the power function is sensibly higher for VSDKs, even compared to the classical RBF reconstruction.

4.2 Tests with artificial data

Let $\varOmega = (-1,1)^2$. We consider two test functions,

$$\begin{aligned} f_2(x,y)= & {} \left\{ \begin{array}{ll} e^{-(x^2+y^2)}, &{} \quad x^2+y^2\le 0.6,\\ x+y, &{} \quad x^2+y^2> 0.6,\\ \end{array} \right. \\ f_3(x,y)= & {} \left\{ \begin{array}{ll} 2(1-e^{-(y+0.5)^2}), &{} \quad \;\;\; |x|\le 0.5,\;|y|\le 0.5,\\ 4(x+0.8), &{} \quad \quad -0.8\le x\le -0.65,\;|y|\le 0.8,\\ 0.5, &{} \quad \quad 0.65\le x\le 0.8,\;|y|\le 0.2,\\ 0, &{} \quad \quad \textit{otherwise.}\\ \end{array} \right. \end{aligned}$$

We take 1089 Halton points on $\varOmega$ as interpolation nodes and we evaluate the approximant on equispaced points with mesh size 1.00E−2.

Similarly to the one-dimensional case in Sect. 4.1, we interpolate the functions $f_2$ and $f_3$ via classical RBF interpolation on the set of nodes $\mathscr {X}$, using the kernel functions in (4.1) and selecting the optimal shape parameter $\varepsilon$ via LOOCV. Finally, we apply VSDKs and we evaluate the final results.

We start our discussion with the function $f_2$. The resulting standard RBF interpolants for $f_2$ are plotted in Fig. 6. As expected, the infinitely smooth Gaussian RBF introduces huge oscillations, while with functions with low regularity, the Gibbs phenomenon is less evident.

In Table 5 we report the accuracy indicators. We can observe that $\varphi ^1_{\varepsilon }$ outperforms the other two kernels in terms of SSIM index. Indeed, the related reconstruction is less affected by the Gibbs phenomenon, as graphically visible.

Considering the function $f_2$, for VSDKs we need a suitable scale function satisfying the Definition 3.2. For this purpose, we consider the function $\psi ^2$ defined as

$$\begin{aligned} \psi ^2(x,y)= \left\{ \begin{array}{ll} 1, &{} \quad x^2+y^2\le 0.6,\\ 2, &{} \quad x^2+y^2> 0.6.\\ \end{array} \right. \end{aligned}$$

(4.4)

We show the final reconstructions using VSDKs with different kernels in Fig. 7, while in Table 6 we report the values of the considered errors and parameters. We recover also for the 2D case the pattern already discovered about the fact that VSDKs reconstruct the jumps without graphically introducing oscillations, also with $C^\infty$ RBFs.

Table 2 Results of VSK and VSDK reconstruction for $f_1$ using $\varphi ^1$

Table 3 Results of VSK and VSDK reconstruction for $f_1$ using $\varphi ^2$

Table 4 Results of VSK and VSDK reconstruction for $f_1$ using $\varphi ^1_3$

Fig. 6

The function $f_2$ and the classical RBF interpolants on $\mathscr {X}$ obtained using different kernel functions

Table 5 Results of classical RBF reconstruction for $f_2$

Fig. 7

VSDK reconstructions of $f_2$ on $\mathscr {X}$ using the scale function $\psi ^2$

Table 6 Results of VSDK reconstructions for $f_2$

Considering now the function $f_3$, we show in Fig. 8 and in Table 7 the results obtained via classical RBF interpolation.

We can observe a behavior that is similar to what we obtained for $f_2$. Switching to VSDKs, we consider the scale function $\psi ^3$ defined as

$$\begin{aligned} \psi ^3(x,y)= \left\{ \begin{array}{ll} 1, &{} \quad \;\;\; |x|\le 0.5,\;|y|\le 0.5,\\ 2, &{} \quad \quad -0.8\le x\le -0.65,\;|y|\le 0.8,\\ 3, &{} \quad \quad 0.65\le x\le 0.8,\;|y|\le 0.2,\\ 0, &{} \quad \quad \textit{otherwise.}\\ \end{array} \right. \end{aligned}$$

(4.5)

We point out that the set of discontinuity points of $f_3$ is strictly contained in the set of discontinuity points of $\psi ^3$, which is the case considered in Remark 3.1. We present the final results in Fig. 9 and in Table 8.

Fig. 8

The function $f_3$ and the classical interpolants on $\mathscr {X}$ obtained using different kernel functions

We can observe that the VSDK reconstructions using $\varphi ^2_{1}$ and $\varphi ^3_{1}$ are not affected by the Gibbs phenomenon as in the classical RBF reconstructions and they outperform the reconstruction obtained using $\varphi ^1_{1}$.

Furthermore, for both the test functions considered in this section, the SSIM with VSDKs is about 1, which means that graphically there is a high similarity between the original and reconstructed image.

Table 7 Results of classical RBF reconstruction for $f_3$

Fig. 9

VSDK reconstructions of $f_3$ on $\mathscr {X}$ using the scale function $\psi ^3$

Concerning the maximum value of the power function for VSDKs, we can note that for both $f_2$ and $f_3$ it is comparable to the one obtained via standard RBFs. However, the reader should note that usually it reaches lower values than the ones achieved via the classical schemes. This reflects directly on the error (as theoretically proved in Proposition 3.2).

In general, for both 1D and 2D, the most promising results are the ones obtained via VSDKs and the Gaussian function. Indeed, it is well known that $C^{\infty }$ RBFs introduce the Gibbs phenomenon, which is sensibly reduced via VSDKs.

5 Application to satellite images

The modeling and analysis of data, for instance, coming from distributed measurements of physical quantities and satellite images is a challenging computational issue. Because of the huge size that some of these data sets achieve, reduced models such as the one presented in [24] are strongly advised. Nevertheless, the Gibbs phenomenon might affect also in this case the accuracy of the approximation. Thus, in this example, we show how VSDKs can intervene in this direction, sensibly reducing the oscillations. We consider the satellite image reported in Fig. 10, consisting of soil moisture data taken by NASA Soil Moisture Active Passive (SMAP) mission in April 2015 [12]. It is composed by a grid of $3856 \times 1624$ pixels. For dealing with the whole image, one needs to use reduced models, such as the one investigated in [25]. Moreover, if one only concentrates on a small portion of the image, e.g. a portion of Europe, the high resolution is lost (trivially due to zooming). In this case, a reconstruction scheme is necessary.

Table 8 Results of VSDK reconstructions for $f_3$

Focusing on a portion of of Europe, we obtain an image composed by $N = 144 \times 191$ pixels. After using these data to reconstruct the image, we evaluate it on a finer grid of evaluation points, composed by $216 \times 286$ pixels. Such a computation can be seen as a standard zoom, which might introduce Gibbs oscillations. They are indeed visible if, for instance, we reconstruct the image with the Wendland’s $C^2$ RBF defined by:

$$\begin{aligned} \varphi ^4_{\varepsilon }(r)= (1-\varepsilon r)_+^4(4 \varepsilon r +1) ,\quad \mathrm{Wendland }\,C^2, \end{aligned}$$

where $(\cdot )_+$ denotes the truncated power function. We consider the Wendland’s compactly supported $C^2$ RBF because it is well-known that by properly scaling the support of the basis function, it might lead to sparse interpolation systems and thus it gains in terms of stability, reducing the usual high condition number of the interpolation matrix. Despite this ability, the reconstruction via the classical method suffers from the Gibbs phenomenon, see Fig. 11 (left). Such oscillations are removed by VSDKs; refer to Fig. 11 (right). This example with real data is also devoted to show that VSDKs are performing also when the discontinuity is analytically unknown. In this case the curve defining the discontinuity has been approximated by means of Sobel detection scheme; see e.g. [31].

Fig. 10

Example of satellite image measuring the soil moisture

Fig. 11

The approximation of soil moisture over a portion of Europe via classical RBF reconstruction (left) and VSDKs (right)

6 Conclusions

In this paper we have presented a robust method for sensibly reducing non-physical oscillations due to the Gibbs phenomenon. The accuracy of the proposed method has been studied theoretically and many numerical experiments confirm its effectiveness, showing that the reconstruction via VSDKs outperforms the standard techniques when jumps occur.

Future work consists in investigating on the detection of discontinuities when dealing with real data and to extend the current work to other bases [8]. We stress that the current study might be useful for image reconstruction in the context of Magnetic Particle Imaging [6, 7]. Finally, for smooth RBFs, we should study the behaviour of VSDKs when rational RBF interpolants are used [9, 19].