# Cardinal Interpolation With General Multiquadrics

Keaton Hamm Department of Mathematics, Vanderbilt University,Nashville, TN, 37240, USA  and  Jeff Ledford Department of Mathematics, Virginia Commonwealth University, Richmond, VA, 23284, USA
###### Abstract.

This paper studies the cardinal interpolation operators associated with the general multiquadrics, , . These operators take the form

 Iα,cy(x)=∑j∈ZdyjLα,c(x−j),y=(yj)j∈Zd,x∈Rd,

where is a fundamental function formed by integer translates of which satisfies the interpolatory condition .

We consider recovery results for interpolation of bandlimited functions in higher dimensions by limiting the parameter . In the univariate case, we consider the norm of the operator acting on spaces as well as prove decay rates for using a detailed analysis of the derivatives of its Fourier transform, .

###### Key words and phrases:
Cardinal Interpolation and General Multiquadric and Cardinal Functions and Paley-Wiener Functions
###### 2010 Mathematics Subject Classification:
41A05 and 41A30 and 41A63
The first author was partially supported by National Science Foundation grant DMS 1160633. The second author was partially supported by the 2014 Workshop in Analysis and Probability at Texas A& M University.

## 1. Introduction

The term cardinal interpolation refers to interpolation of data given at the integer lattice (or multi-integer lattice in higher dimensions). It was I. J. Schoenberg’s work on cardinal spline interpolation that brought about an intense study of the subject. Many avenues of study have been explored, including forming interpolation operators from translates of certain radial basis functions (RBFs). Works by Buhmann, Baxter, Riemenschneider, and Sivakumar [3, 7, 29] (see also [28] and references therein) have explored many cardinal interpolation operators of this type. Some of the radial basis functions that have been considered are the Gaussian kernel, the thin plate spline, the Hardy multiquadric, and the inverse multiquadric.

Radial basis cardinal interpolation also enjoys a strong connection with classical sampling theory, as evidenced by much of the aforementioned literature. This connection continues to be explored in recent developments by the second author [23, 24], and by parts of this article. As this is one of the most interesting aspects of cardinal interpolation, we provide some of the motivation. Recall the one-dimensional Whittaker–Kotelnikov–Shannon Sampling Theorem, which states that a bandlimited function, , (say with band size ) can be recovered uniformly by the series , where the sinc function is suitably defined so that it takes value 1 at the origin, and 0 at all the other integers. One observation about this series is that the sinc function decays slowly (as ), and so to approximate the series by truncation for example, one would have to use quite a lot of data points of to get a reasonable degree of accuracy.

However, there is a way of approximating the sinc series above: we seek to replace the sinc function with a so-called fundamental function (or cardinal function), , that preserves the property that takes value 1 at the origin and 0 at all other integers. We then form a function

 If(x)=∑j∈Zf(j)L(x−j).

The trade-off here is that while is not pointwise equal to the function , it does interpolate at the integer lattice, and moreover, the fundamental function may be constructed so that it decays much more rapidly than the sinc function. Precisely, one may construct a fundamental function from a given radial basis function, , which has the form

 Lϕ(x)=∑j∈Zajϕ(x−j).

In the case of the Gaussian kernel, , the fundamental function decays exponentially, whereas the fundamental function for the Hardy multiquadric, , decays as . So we give up the pointwise equality of the WKS Sampling Theorem in exchange for a series that converges more rapidly, while also ensuring that is close to in the norm.

This paper primarily considers the fundamental functions and cardinal interpolation operators associated with general multiquadrics, , which have thus far only been considered for certain exponents . Interpolation with fundamental functions has too long a history to recount here; however, [9] offers a good introduction using radial basis functions. Using (6) below as a starting point is especially popular since it allows one to solve problems in the Fourier transform domain. Many authors have used similar techniques for various radial basis functions. In [28, 29, 30], Riemenschneider and Sivakumar proved several results pertaining to the Gaussian. Multiquadric cardinal interpolation has been studied in a similar way by Buhmann and Micchelli [10], Baxter [3], Baxter and Sivakumar [4], Riemenschneider and Sivakumar [29], among others. Compactly supported radial basis functions have been studied by Buhmann [8] and Wendland [33].

The rest of the paper is laid out as follows. Section 2 provides the necessary preliminaries and a discussion of applications and calculations of the fundamental functions; Section 3 shows recovery results for cardinal interpolation of bandlimited functions in any dimension via interpolants of the form discussed above. Section 4 contains decay rates and information about the univariate fundamental functions associated with the general multiquadrics for a broad range of exponents. In Section 5, we consider the cardinal interpolation operators acting on data in traditional sequence spaces and calculate decay rates, bounds on the operator norms, and also explore some convergence properties in terms of the parameter . Section 6 provides some interesting concrete examples based on the theoretical results from the previous section. Finally, Section 7 provides the technical proofs of the statements in Section 4.

## 2. Basic Notions

If is an interval, then let , , be the usual Lebesgue space over with its usual norm. If no set is specified, we mean . Similarly, denote by the usual sequence spaces indexed by the set ; if no index set is given, we refer to . We will use to denote the natural numbers including 0. Let be the space of Schwartz functions on , that is the collection of infinitely differentiable functions such that for all multi-indices and ,

 supx∈Rd∣∣xαDβϕ(x)∣∣<∞.

The Fourier transform of a Schwartz function is given by

 (1)

Thus the inversion formula is

 (2) ϕ(x)=1(2π)d∫Rdˆϕ(ξ)ei⟨x,ξ⟩dξ,x∈Rd.

In the event that these formulas do not hold, then the Fourier transform should be interpreted in the sense of tempered distributions. Recall that if is a tempered distribution, then its Fourier transform is the tempered distribution defined by , .

Let and be fixed; then define the -dimensional general multiquadric by

 (3) ϕα,c(x):=(∥x∥2+c2)α,x∈Rd,

where denotes the Euclidean distance on .

If , the generalized Fourier transform of is given by the following (see, for example, [34, Theorem 8.15]):

 (4) ˆϕα,c(ξ)=21+αΓ(−α)(c∥ξ∥)α+d2Kα+d2(c∥ξ∥),ξ∈Rd∖{0},

where [32, p.185]

 (5) Kν(r)=Γ(12)rν2νΓ(ν+12)∫∞1e−rx(x2−1)ν−12dx,ν≥0,r>0.

is called the modified Bessel function of the second kind. We note that both and its Fourier transform are radial. It is also clear from the definition that is symmetric in its order; that is, for any . If , then the generalized Fourier transform of involves a measure and so cannot be expressed as a function.

### 2.1. Fundamental Functions

Now suppose that is fixed. To define the fundamental function associated with the general multiquadric, we first define the following function

 (6) ˆLα,c(ξ):=ˆϕα,c(ξ)∑j∈Zdˆϕα,c(ξ+2πj),ξ∈Rd.

We will see that , and so the function

 (7) Lα,c(x):=1(2π)d∫RdˆLα,c(ξ)ei⟨x,ξ⟩dξ,x∈Rd,

is well-defined and continuous. Furthermore, we will show that is a fundamental function, also called a cardinal function, which means that

 (8) Lα,c(j)=δ0,j,j∈Zd,

where is the Kronecker delta.

 (9) Lα,c(x)=∑j∈Zdcjϕα,c(x−j),x∈Rd.

Throughout the paper, we will use to denote an absolute constant due to the use of as the shape parameter of the multiquadric. The value of the particular constant may change from line to line, and we use subscripts to denote dependence on certain parameters when needed.

### 2.2. Evaluation of Fundamental Functions and Applications

Interpolation schemes involving fundamental functions as in (9) have been studied for quite some time, and there are many aspects to the theory. For example, such methods enjoy applications to geoscience [13] and sampling theory [23]. Recently, investigations have considered interpolation via radial kernels on manifolds [20, 21]. For a Galerkin type method for solving PDEs using meshless interpolation on the sphere, see [27].

Given the widespread applications of radial basis function approximation, it is of import to the computational community to determine stable ways of evaluating the approximants. Consequently, there is a substantial literature dealing with accuracy and stability of different computational methodologies for radial basis function approximation. We do not claim to list all of these methods, but at least a sampling is in order. We note that approximating (9) is typically very difficult, especially if one dilates the lattice. One way around this is the use of indirect computational methods to approximate the RBF interpolant [14, 15] (for a discussion specifically related to multiquadrics, see [16]). Another technique involves a change of basis method [5], while work by Fasshauer and McCourt [12] uses an eigenfunction decomposition to provide stable reconstruction using Gaussians.

Another quite promising method has recently been considered in which so-called local Lagrange functions are used to approximate rather than the global ones [17, 19]. Many of these results revolve around the situation of interpolation at finitely many data sites, which is of a somewhat different nature than we are considering here. For the interested reader, we also mention that these methods of cardinal interpolation have, at their core, a deep connection to the classical spline theory instigated by Schoenberg (see [31] and references therein) and continued by many followers. The underlying principle is that many of the results in spline interpolation theory have natural analogues via RBFs, and the problem at hand may determine which method is more useful.

### 2.3. Examples

Here, we provide some brief illustrations of the fundamental functions we have mentioned above. Figure 1 shows the univariate fundamental function associated with the inverse multiquadric and the sinc function for comparison.

The fundamental function was calculated by truncating the series in (6) and using a fast Fourier transform (FFT) method to approximate . The same method may be applied for the multivariate version, as Figure 2 shows.

## 3. Recovery of Multivariate Bandlimited Functions

When an interpolation scheme depends on a parameter, questions of convergence naturally arise. This question has been addressed by several authors. In [3], Baxter examines the Hardy multiquadric, while the Gaussian is studied in [30] by Riemenschneider and Sivakumar and in [18] by Hangelbroek, Narcowich, Madych, and Ward. In a more general context, ‘increasingly flat’ radial basis functions are the focus of Driscoll and Fornberg in [11], while the second author worked with ‘regular families’ of cardinal interpolants in [23].

In this section, we show that the result obtained by Baxter [3] holds not only for the traditional Hardy multiquadric (corresponding to ) but rather for all . We consider interpolation of bandlimited (or Paley-Wiener) functions in any dimension, and show that the cardinal interpolant converges to the function as the shape parameter tends to infinity. General multiquadrics were not considered for quite some time in this setting, but in [23] convergence results for cardinal interpolation of bandlimited functions were shown for a restricted range of exponents. However, the analysis there was of a more general nature, so here we show that a more specific analysis yields convergence results for the full range

We note that the results of this section have very close ties to classical sampling theory, which studies the reconstruction of certain classes of signals from their discrete samples. As mentioned above, these considerations lead to alternative methods for approximating the sampling series given by the WKS Sampling Theorem for bandlimited signals.

Let be the dimension, and and be fixed. It is evident from (4) that does not change sign. Therefore, for all . From (6), it is also evident that . To show that we begin with the following lemma.

###### Lemma 1.

Let , , and . Then

 |ˆϕα,c(R)|≤e−c(R−r)|ˆϕα,c(r)|.
###### Proof.

Defining , equations (4) and (5) yield the following series of estimates:

 |ˆϕα,c(R)|=|λ|∫∞1e−cRt(t2−1)α+d−12dt=|λ|∫∞1e−c(R−r)te−crt(t2−1)α+d−12dt≤|λ|e−c(R−r)∫∞1e−crt(t2−1)α+d−12dt=e−c(R−r)|ˆϕα,c(r)|.

We note that if , then , and so we have a purely exponential upper bound.

Let and . Then .

###### Proof.

First, choose an large. Then since for all , we have that

 ∫[−M,M]d|ˆLα,c(ξ)|dξ≤(2M)d.

Now we need to estimate

 I:=∫Rd∖[−M,M]d|ˆLα,c(ξ)|dξ.

To do this, we establish a pointwise estimate for . Let be fixed. Since is large, there exists some such that . Additionally, there is some constant which depends on and such that if , (see, for example, [34, Corollary 5.12]). Therefore, choose large enough so that for , we have . Then if is the constant from Lemma 1,

 ∣∣ ∣∣∑k∈Zdˆϕα,c(ξ+2πk)∣∣ ∣∣ ≥|ˆϕα,c(ξ+2πkξ)| ≥γ|λ|∥ξ+2πkξ∥−α−d2e−c∥ξ+2πkξ∥(c∥ξ+2πkξ∥)−12.

Now depending on the sign of , the above expression is minimized by plugging in or for in the appropriate places. Consequently, there is a positive constant such that

 ∣∣ ∣∣∑k∈Zdˆϕα,c(ξ+2πk)∣∣ ∣∣≥De−4πc.

We also find from [34, Lemma 5.13] that for every , Consequently, by adjusting if need be so that for , we find that there is a positive constant such that . We conclude that

 I ≤D−1e4πc∫Rd∖[−M,M]d|ˆϕα,c(ξ)|dξ ≤βD−1|λ|e4πc∫Rd∖[−M,M]d∥ξ∥−α−d2e−c∥ξ∥dξ.

The integral on the right hand side above is convergent, and the constants outside are all finite, so . ∎

Now we turn our attention to the function .

###### Proposition 3.

Let and . Then the function

 Lα,c(x)=1(2π)d∫RdˆLα,c(ξ)ei⟨x,ξ⟩dξ

is continuous, square-integrable, and satisfies the interpolatory condition , for every .

###### Proof.

Proposition 2 implies that is continuous and square-integrable, and indeed that is its Fourier transform. To see the interpolatory condition, first define . Then we have via the substitution that

The interchange of sum and integral in the third line is justified by the Dominated Convergence Theorem, for example. ∎

It is an important observation that much of the cardinal interpolation theory for bandlimited functions revolves around the fact that the fundamental functions converge to the function , which is equivalent to the statement that the Fourier transform of the fundamental function converges almost everywhere to the indicator function of the cube . The story is no different here. Defining to be the function that takes value 1 whenever , and 0 elsewhere, the following holds.

###### Proposition 4.

Let . Then

 limc→∞ˆLα,c(ξ)=I(ξ)

for all such that .

###### Proof.

First suppose that . Then there exists some such that . Therefore by Lemma 1,

 |ˆϕα,c(ξ)|≤e−c(∥ξ∥−∥ξ+2πk0∥)|ˆϕα,c(ξ+2πk0)|≤e−c(∥ξ∥−∥ξ+2πk0∥)∑k∈Zd|ˆϕα,c(ξ+2πk)|.

Consequently, since is of one sign,

 0≤ˆLα,c(ξ)≤e−c(∥ξ|−∥ξ+2πk0∥).

Since the exponent is negative, the limit of the right hand side as is 0. Therefore, for , .

Now suppose that . Then for all , . Due to (6), we may write

 ˆLα,c(ξ)=⎛⎝1+∑k≠0ˆϕα,c(ξ+2πk)ˆϕα,c(ξ)⎞⎠−1,

and therefore it suffices to show that

 limc→∞∑k≠0ˆϕα,c(ξ+2πk)ˆϕα,c(ξ)=0.

By Lemma 1,

 0≤∑k≠0ˆϕα,c(ξ+2πk)ˆϕα,c(ξ)≤∑k≠0e−c(∥ξ+2πk∥−∥ξ∥).

The series on the right hand side is convergent and dominated by the convergent series where is replaced by 1, so

 limc→∞∑k≠0ˆϕα,c(ξ+2πk)ˆϕα,c(ξ)=0

as desired. Convergence of the series stems from the fact that the series is majorized by , which converges. ∎

We now consider interpolation of bandlimited functions at the lattice by translates of the function . Define the -dimensional Paley-Wiener space by

 PW(d)π:={f∈L2(Rd):ˆf=0 a.e. % outside [−π,π]d}.

We begin our analysis with an version of the Poisson Summation Formula:

###### Lemma 5 (cf. [3] Lemma 3.2).

If , then

 (10) ∑j∈Zdˆf(ξ+2πj)=∑j∈Zdf(j)e−i⟨j,ξ⟩,

where the second series is convergent in .

###### Lemma 6.

Let . For , define

 ˆImα,cf(ξ):=⎛⎝∑∥k∥1≤mf(k)e−i⟨k,ξ⟩⎞⎠ˆLα,c(ξ),ξ∈Rd,

where for . Then forms a Cauchy sequence in .

###### Proof.

Define via

 Qm(ξ)=∑∥k∥1≤mf(k)e−i⟨k,ξ⟩.

Thus, . From Lemma 5, it is clear that is a Cauchy sequence in . So

 ∥ˆImα,cf−ˆIℓα,cf∥2L2(Rd)≤∫Rd|Qm(ξ)−Qℓ(ξ)|2(ˆLα,c(ξ))2dξ=∑k∈Zd∫[−π,π]d|Qm(ξ+2πk)−Qℓ(ξ+2πk)|2×(ˆLα,c(ξ+2πk))2dξ=∫[−π,π]d|Qm(ξ)−Qℓ(ξ)|2∑k∈Zd(ˆLα,c(ξ+2πk))2dξ≤∫[−π,π]d|Qm(ξ)−Qℓ(ξ)|2dξ.

The interchange of sum and integral is valid by Tonelli’s Theorem, and the last inequality follows from the fact that

 (11) ∑k∈Zd(ˆLα,c(ξ+2πk))2=∑k∈Zdˆϕα,c2(ξ+2πk)(l∈Zd∑ˆϕα,c(ξ+2πl))2≤1.

We also used the fact that for , . We conclude that is a Cauchy sequence in because , and the latter is Cauchy. ∎

Lemmas 5 and 6 allow us to define

 (12) ˆIα,cf(ξ):=ˆLα,c(ξ)∑k∈Zdf(k)e−i⟨k,ξ⟩,

where the series is convergent in . By a periodization argument similar to that in the proof of Lemma 6, one can show that . Thus applying the Fourier inversion formula we see that

 Iα,cf(x)=∑k∈Zdf(k)Lα,c(x−k),x∈Rd.
###### Theorem 7.

Let . If , then

 limc→∞∥Iα,cf−f∥L2(Rd)=0,

and uniformly on .

###### Proof.

We will first prove uniform convergence. The proof is the same as in [3]. Again let be the characteristic function of the cube. Then we see by the inversion formula and the oft-exploited periodization argument, that

 Iα,cf(x)−f(x)=1(2π)d∫Rd∑k∈Zdˆf(ξ+2πk)(ˆLα,c(ξ)−I(ξ))e−i⟨x,ξ⟩dξ=1(2π)d∫[−π,π]dˆf(ξ)∑k∈Zd(ˆLα,c(ξ+2πk)−I(ξ+2πk))×e−i⟨x,ξ+2πk⟩dξ.

Therefore, we find that

 |Iα,cf(x)−f(x)|≤1(2π)d∫[−π,π]d|ˆf(ξ)|∑k∈Zd∣∣ˆLα,c(ξ+2πk)−I(ξ+2πk)∣∣dξ=1(2π)d∫[−π,π]d|ˆf(ξ)|⎛⎝1−ˆLα,c(ξ)+∑k≠0ˆLα,c(ξ+2πk)⎞⎠dξ.

But then by definition,

 ∑k≠0ˆLα,c(ξ+2πk)=∑k∈Zdˆϕα,c(ξ+2πk)−ˆϕα,c(ξ)∑l∈Zdˆϕα,c(ξ+2πl)=1−ˆLα,c(ξ).

Therefore,

 |Iα,cf(x)−f(x)|≤21(2π)d∫[−π,π]d|ˆf(ξ)|(1−ˆLα,c(ξ))dξ.

As the integrand is non-negative and bounded by , and , the Dominated Convergence Theorem implies that

 limc→∞|Iα,cf(x)−f(x)|=0,x∈Rd.

The upper bound is independent of , hence the convergence is uniform.

We now turn to the proof of convergence. By Parseval’s Identity, it suffices to show that . This breaks up into two estimates. We first show this for the cube . Recall that since is an orthonormal basis for , we may write . Moreover,

 ∥ˆf∥L2[−π,π]d=∥f(k)∥ℓ2(Zd).

Thus using (12),

 ∥ˆIα,cf−ˆf∥2L2[−π,π]d =∫[−π,π]d∣∣ ∣∣∑k∈Zdf(k)(ˆLα,c(ξ)−1)e−i⟨k,ξ⟩∣∣ ∣∣2dξ

The right hand side is bounded by , and so by the Dominated Convergence Theorem and Proposition 4, .

Now for the rest of the space, for define . Then we see that since is bandlimited,

 ∥ˆIα,cf−ˆf∥2L2(Rd∖[−π,π]d)=∥ˆIα,cf∥2L2(Rd∖[−π,π]d)=∑l≠0∥ˆIα,cf∥2L2(Ql).

Consequently,

 ∫Rd∖[−π,π]d|ˆIα,cf(ξ)|2dξ=∑l≠0∫Ql∣∣ ∣∣ˆLα,c(ξ)∑k∈Zdf(k)e−i⟨k,ξ⟩∣∣ ∣∣2dξ=∫[−π,π]d∑l≠0|ˆLα,c(ξ+2πl)|2∣∣ ∣∣∑k∈Zdf(k)e−i⟨k,ξ⟩∣∣ ∣∣2dξ,

by the Monotone Convergence Theorem.

Recall that , so , and as calculated above, . Consequently, the integrand is bounded by

 |1−ˆLα,c(ξ)|∣∣ ∣∣∑k∈Zdf(k)e−i⟨k,ξ⟩∣∣ ∣∣2≤2∥f(k)∥2ℓ2(Zd).

Therefore, the Dominated Convergence Theorem and Proposition 4 imply that , and the proof is complete. ∎

To illustrate the convergence given by Theorem 7 above, the following figure shows the inverse multiquadric interpolant of the function whose Fourier transform is in dimension 1.

## 4. Properties of the Fundamental Function

For the rest of the paper, we turn our attentions to the one-dimensional cardinal interpolation operator associated with the general multiquadric. This section is devoted to the one-dimensional fundamental function , whose Fourier transform can be rewritten as

 (13) ˆLα,c(ξ)=⎡⎣1+∑j≠0ˆϕα,c(ξ+2πj)ˆϕα,c(ξ)⎤⎦−1.

The proofs of the results in this section are quite technical, so we postpone them until Section 7 and simply state our conclusions here.

To determine decay rates for , we determine how many derivatives has in , which we accomplish by establishing pointwise estimates. We begin by fixing , so that our estimates fall into three ranges: , , and the -length blocks for . Due to the differing behavior of for positive and negative values of , we must make corresponding distinctions in our calculations.

Following the insightful techniques of Riemenschneider and Sivakumar found in [29], we begin by defining some auxiliary functions to aid in the analysis of the fundamental function. We abbreviate (13) as , where

 (14) sα,c(ξ):=∑j≠0ˆϕα,c(ξ+2πj)/ˆϕα,c(ξ)=:∑j≠0aj(ξ),

and study the properties of .

###### Proposition 8.

Suppose that , , , and . If and , then there exists a constant such that

 (15) |a(k)j(ξ)|≤Aα,k(ε)c(k+1)(2α−⌊α⌋)+ke−2πcεe−2πc(|j|−1),

where as .

This estimate leads to the following bounds on and its derivatives.

###### Proposition 9.

Suppose that , , , and . If , then there exist constants such that

1. whenever ,

2. whenever , and

3. whenever and ,

where