Wiener–Hopf Difference Equations and Semi-Cardinal Interpolation with Integrable Convolution Kernels

Abstract

Let $H\subset {\mathbb {Z}}^d$ be a half-space lattice, defined either relative to a fixed coordinate (e.g. $H = {\mathbb {Z}}^{d-1}\!\times \!{\mathbb {Z}}_+$), or relative to a linear order $\preceq$ on ${\mathbb {Z}}^d$, i.e. $H = \{j\in {\mathbb {Z}}^d: 0\preceq j\}$. We consider the problem of interpolation at the points of H from the space of series expansions in terms of the H-shifts of a decaying kernel $\phi$. Using the Wiener–Hopf factorization of the symbol for cardinal interpolation with $\phi$ on ${\mathbb {Z}}^d$, we derive some essential properties of semi-cardinal interpolation on H, such as existence and uniqueness, Lagrange series representation, variational characterization, and convergence to cardinal interpolation. Our main results prove that specific algebraic or exponential decay of the kernel $\phi$ is transferred to the Lagrange functions for interpolation on H, as in the case of cardinal interpolation. These results are shown to apply to a variety of examples, including the Gaussian, Matérn, generalized inverse multiquadric, box-spline, and polyharmonic B-spline kernels.

1 Introduction

We start by reviewing some background results concerning cardinal interpolation on the multi-integer lattice ${\mathbb {Z}}^d$, which represents a classical model in approximation theory.

Let $\phi :{\mathbb {R}}^d\rightarrow {\mathbb {R}}$ be a continuous and symmetric kernel, suitably decaying for a large norm of its argument, and denote by ${\mathcal {S}}(\phi )$ the shift-invariant space of pointwise convergent infinite linear combinations of ${\mathbb {Z}}^d$-translates of $\phi$ with bounded coefficients. Cardinal interpolation with the kernel $\phi$ is the problem of constructing a function $s\in {\mathcal {S}}(\phi )$ with prescribed values at the points of ${\mathbb {Z}}^d$, i.e. $s(j)=y_j$, $j\in {\mathbb {Z}}^d$, for a given bounded sequence $y=\{y_{j}\}_{j\in {\mathbb {Z}}^d}$.

This leads to a bi-infinite system of discrete convolution equations with Laurent matrix $L_\phi :=[\phi (j-k)]_{j,k\in {\mathbb {Z}}^d}$. The necessary and sufficient condition for the invertibility of this system is that the symbol $\sigma _\phi$, defined by the absolutely convergent series

$$\begin{aligned} \sigma _\phi ( \zeta ):= \sum _{k\in {\mathbb {Z}}^d} \phi (k) \zeta ^{k}, \quad \zeta \in T^d, \end{aligned}$$

(1.1)

has no zero on the d-dimensional unit torus $T^d$ (e.g. Chui et al. [19, Lemma 1.1]). If this condition is satisfied, then, by Wiener’s lemma, $1 / \sigma _\phi$ also admits an absolutely convergent Fourier representation, with Fourier coefficients $\{a_{k}\}_{k\in {\mathbb {Z}}^d}$, say. In turn, these define the inverse Laurent matrix $L_\phi ^{-1}=[a_{j-k}]_{j,k\in {\mathbb {Z}}^d}$, as well as the Lagrange (or fundamental) function

$$\begin{aligned} \chi ( x ):= \sum _{k\in {\mathbb {Z}}^d} a_k \phi (x-k),\quad x\in {\mathbb {R}}^d, \end{aligned}$$

(1.2)

that interpolates the Kronecker delta sequence:

$$\begin{aligned} \chi (j) = \delta _{j,0}, \quad j\in {\mathbb {Z}}^d. \end{aligned}$$

(1.3)

Hence, provided that $\chi$ also decays suitably, the unique solution s to the cardinal interpolation problem admits the pointwise convergent Lagrange representation in terms of the multi-integer shifts of $\chi$:

$$\begin{aligned} s ( x ) =\sum _{j\in {\mathbb {Z}}^d} y_{j} \chi ( x-j ),\quad x\in {\mathbb {R}}^d. \end{aligned}$$

(1.4)

Note that the localization quality of this representation is determined by the rate of dampening of $|\chi (x)|$ for large $\Vert x\Vert$. For instance, if $\phi$ is a univariate B-spline (Schoenberg [54, 55]) or a multivariate box-spline with non-vanishing symbol $\sigma _\phi$ (de Boor et al. [21, 22]), both of whom are compactly supported kernels, or if $\phi$ is the Gaussian kernel (Sivakumar [58]), then the Lagrange function $\chi$ decays exponentially. More generally, it follows via analyticity arguments that, if $|\phi (x)| = O (e^{-\alpha \Vert x\Vert })$, then $|\chi (x)| = O(e^{-\beta \Vert x\Vert })$, for some $\beta \in (0,\alpha )$ (see Sect. 2). On the other hand, if $\phi$ has an algebraic decay with rate $\alpha >d$, i.e. $|\phi (x)| = O ((1+\Vert x\Vert )^{-\alpha })$, then $\chi$ decays with exactly the same rate $|\chi (x)| = O((1+\Vert x\Vert )^{-\alpha })$, as proved by Bacchelli et al. [1], using a well-known result of Jaffard [38]. Recently, Fageot et al. [27] considered a weighted $\ell _1$-decay condition, for a weight function w:

$$\begin{aligned} \sup _{x\in [0,1]^d} \sum _{k\in {\mathbb {Z}}^d} | \phi (x-k) | w(k) < \infty . \end{aligned}$$

(1.5)

Assuming that w fulfills certain admissibility conditions and the cardinal symbol $\sigma _\phi$ does not vanish on the unit torus, they proved [27, Proposition 27] that property (1.5) is inherited by the corresponding Lagrange function $\chi$.

The above results on transferring the decay of the kernel $\phi$ to its associated Lagrange function $\chi$ for cardinal interpolation rely crucially on the fact that Wiener’s lemma holds not only in the classical Wiener algebra of symbols with absolutely convergent Fourier expansion, but also in subalgebras of symbols defined by stronger localization properties of Fourier coefficients, such as exponential or algebraic decay, or finiteness of some weighted $\ell _1$-norm. This fundamental principle of localization transfer under inversion, or inverse-closedness, has been extended in recent decades far beyond commutative algebras of symbols (or, equivalently, of Laurent matrices), to large classes of non-commutative algebras of matrices, integral operators, and pseudo-differential operators with off-diagonal decay (e.g. the surveys [33, 42, 57] and references therein).

Our paper deals with a new application of this principle to approximation theory, in the non-commutative setting of semi-infinite Toeplitz matrices generated by the problem of semi-cardinal kernel interpolation on a half-space lattice $H\subset {\mathbb {Z}}^d$. We will consider two types of half-space lattices, determined either relative to a fixed coordinate (e.g. $H = {\mathbb {Z}}^{d-1}\!\times \!{\mathbb {Z}}_+$), or relative to a linear order $\preceq$ on ${\mathbb {Z}}^d$, such that ${\mathbb {Z}}^d$ becomes an (additive) ordered group and $H = \{j\in {\mathbb {Z}}^d: 0\preceq j\}=:{\mathbb {Z}}^d_{\preceq ,+}$. In both cases, we formulate the problem of semi-cardinal interpolation on H with the kernel $\phi$, as follows:

Given a bounded sequence of real values $y=\{y_j\}_{j\in H}$, determine a bounded sequence of coefficients $c=\{c_k\}_{k\in H}$, such that the function

$$\begin{aligned} s ( x ) = \sum _{k\in H} c_k \phi (x-k),\quad x\in {\mathbb {R}}^d, \end{aligned}$$

(1.6)

satisfies the interpolation conditions $s(j) = y_j$, for all $j\in H$.

This problem is equivalent to solving the multi-index Wiener–Hopf system of difference equations

$$\begin{aligned} \sum _{k\in H} c_k \phi (j-k) = y_j,\quad j\in H, \end{aligned}$$

(1.7)

with the semi-infinite Toeplitz matrix $T_\phi :=[\phi (j-k)]_{j,k\in H}$. Based on methods developed in the case $d=1$ by Krein [41] and by Calderón et al. [18], the necessary/sufficient conditions for the invertibility of such systems, as well as the explicit form of their solution in the multi-index case ($d>1$), were given, for $H = {\mathbb {Z}}^{d-1}\!\times \!{\mathbb {Z}}_+$, by Goldenstein and Gohberg [29, 30], and, for $H={\mathbb {Z}}^d_{\preceq ,+}$, by van der Mee et al. [44] (see also Bötcher and Silbermann [16], Gohberg and Feldman [28]). However, these studies did not consider whether any specific off-diagonal decay of the Toeplitz matrix $T_\phi$, expressed as decay of the sequence $\{\phi (k)\}_{k\in {\mathbb {Z}}^d}$, may carry over to the entries of the inverse matrix $T_\phi ^{-1}$.

Our main results, Theorems 3.6 and 3.8, prove that, in analogy to cardinal interpolation, the algebraic or exponential decay of the kernel $\phi$ is converted into corresponding off-diagonal decay of the inverse matrix and further transferred to the Lagrange functions for semi-cardinal interpolation. We also establish, under minimal assumptions on $\phi$, some fundamental properties of semi-cardinal kernel interpolation that follow from the fact that the matrix $T_\phi ^{-1}$ satisfies the Schur condition (Theorem 3.4).

These results are obtained in Sect. 3, by constructing the Lagrange scheme for semi-cardinal interpolation on H. Whereas the cardinal interpolation scheme (1.4) uses the shifts of a single function $\chi$, this is no longer valid for the semi-cardinal Lagrange functions. Thus, for each $j\in H$, we construct a separate Lagrange function $\chi _j$, whose coefficient sequence in representation (1.6) is defined by means of the factors of the Wiener–Hopf (or canonical) factorization of the cardinal inverse symbol $\sigma _\phi ^{-1}$. To transfer the polynomial/exponential decay properties from the kernel $\phi$ to the family of Lagrange functions $\{\chi _j: j\in H\}$, with constants that are independent of j, we apply Jaffard’s well-known inverse-closedness results [38] for infinite matrices with off-diagonal decay.

We remark that, for $d>1$, interpolation on $H={\mathbb {Z}}^{d-1}\!\times \!{\mathbb {Z}}_+$ from the shifts of a decaying kernel has previously been considered by Bejancu and Sabin [7, 13] only for $d=2$ and $\phi = M_{2,2,2}$, the compactly supported three-direction box-spline of multiplicity 2 in each direction. In these references, the representation (1.6) is augmented with an extra layer of shifts of the kernel, which allows the formulation of boundary conditions in terms of kernel coefficients. As detailed in Example 5.4, the approximation results obtained for the stationary schemes treated in [7, 13] demonstrate that semi-cardinal interpolation can be an effective model for the analysis of boundary effects of kernel interpolation schemes under scaling. On the other hand, kernel interpolation on a half-space lattice $H={\mathbb {Z}}^d_{\preceq ,+}$ induced by a linear order on ${\mathbb {Z}}^d$ has not been studied before. As shown in Theorem 3.4(iii), the order structure leads to a Cholesky factorization of the inverse matrix $T_\phi ^{-1}$, via the notion of a multi-index triangular matrix.

The paper is laid out as follows. In Sect. 2, we review material pertaining to cardinal interpolation, including basic assumptions on the kernel $\phi$, as well as some results and concepts that are used in the next sections. The main results on semi-cardinal interpolation appear in Sect. 3, where we treat both types of half-space lattices in parallel, specifying their differences en route. Section 4 obtains two types of results, the first one dealing with the approximation relationship between the cardinal and semi-cardinal interpolation schemes, while the second proving a variational characterization that extends a classical spline property. In the last section, our results are applied to five classes of examples: the Gaussian, Matérn, generalized inverse multiquadric, box-spline, and polyharmonic B-spline kernels, leading not only to new results, but also to some improvements and extensions of current literature.

Notation remarks. In the sequel, we use $\ell _p:= \ell _p ({\mathbb {Z}}^d)$. Also, $T^d$ denotes the unit d-dimensional torus $T^d:= \{ \zeta =(\zeta _1,\dots ,\zeta _d) \in {\mathbb {C}}^d: |\zeta _1|=\ldots =|\zeta _d|=1 \}$, and $W:=W(T^d)$ is the Wiener algebra of continuous functions u on $T^d$ with Fourier coefficients $\{{\widehat{u}}_k\}_{k\in {\mathbb {Z}}^d}\in \ell _1$. The norm of $u\in W$ is $\Vert u \Vert _W = \sum _{k\in {\mathbb {Z}}^d} | {\widehat{u}}_k |$. For $\zeta =(\zeta _1,\dots ,\zeta _d) \in {\mathbb {C}}^d$ and $k=(k_1,\dots ,k_d)\in {\mathbb {Z}}^d$, we use $\zeta ^k = \zeta _1^{k_1}\cdots \zeta _d^{k_d}$.

2 Cardinal Interpolation on the Lattice ${\mathbb {Z}}^d$

In the first part of this section, we review basic questions on existence, uniqueness, and Lagrange representation for solutions of the cardinal interpolation problem, under mild assumptions on the kernel $\phi$. The second part shows that specific algebraic/exponential decay conditions on $\phi$ are transferred to the Lagrange function $\chi$ for cardinal interpolation. The presentation is a blueprint for the new results on semi-cardinal interpolation provided in the next section.

2.1 The Cardinal Interpolation Scheme

Let $\phi :{\mathbb {R}}^d\rightarrow {\mathbb {R}}$ be a continuous and symmetric function, i.e. $\phi (-x) = \phi (x)$, $x\in {\mathbb {R}}^d$. The main assumption we make on $\phi$ in this subsection is the condition

$$\begin{aligned} |\phi |_\infty := \sup _{x\in [0,1]^d} \sum _{k\in {\mathbb {Z}}^d} | \phi (x-k) | < \infty . \end{aligned}$$

(2.1)

This ensures that the shift-invariant space ${\mathcal {S}}(\phi )$ described in Sect. 1 consists of bounded functions s defined by a pointwise absolutely convergent series

$$\begin{aligned} s ( x ) = \sum _{k\in {\mathbb {Z}}^d} c_k \phi (x-k),\quad x\in {\mathbb {R}}^d, \end{aligned}$$

(2.2)

for an arbitrary bounded sequence $c=\{c_k\}_{k\in {\mathbb {Z}}^d}$.

Assumption (2.1) and related $L_p$-norm conditions have been employed by Jia and Micchelli [39] in their work on multiresolution analysis, and are also known to provide a natural setting for principal shift-invariant theory, e.g. de Boor and Ron [23], Johnson [40]. In particular, with the notation of [39], we have

$$\begin{aligned} \Vert \phi \Vert _{L_p ( {\mathbb {R}}^d )} \le |\phi |_p:= \left\| \sum _{k\in {\mathbb {Z}}^d} | \phi (\cdot -k) | \right\| _{L_p ( [0,1]^d )} \le |\phi |_\infty ,\quad p\in [1,\infty ], \end{aligned}$$

(2.3)

hence (2.1) implies $\phi \in L_p ({\mathbb {R}}^d)$ for $1\le p\le \infty$. In the context of cardinal interpolation, Fageot et al. [27] have recently used the stronger condition (1.5), which they relate to the general concept of Wiener amalgam spaces.

As stated in the previous section, for a bounded sequence $y=\{y_{j}\}_{j\in {\mathbb {Z}}^d}$ of real values, the problem of cardinal interpolation with $\phi$ is to find a function of the form (2.2), such that

$$\begin{aligned} s ( j ) =y_{j},\qquad \forall \ j\in {\mathbb {Z}}^d. \end{aligned}$$

(2.4)

This amounts to the infinite system of discrete convolution equations for the coefficient sequence c:

$$\begin{aligned} \sum _{k\in {\mathbb {Z}}^d} c_k \phi (j-k) = y_j,\quad j\in {\mathbb {Z}}^d. \end{aligned}$$

(2.5)

Recall the notation $L_\phi$ for the multi-index Laurent matrix $[\phi (j-k)]_{j,k\in {\mathbb {Z}}^d}$, identified with the convolution operator that maps a bounded sequence c to the sequence $L_\phi \, c$ whose j-th component is given by the left-hand side of (2.5). Under assumption (2.1), the associated symbol $\sigma _\phi$ defined by (1.1) belongs to the Wiener algebra W. Since $\sigma _\phi$ is continuous on $T^d$ and it takes only real values (due to the symmetry of $\phi$), the non-vanishing of $\sigma _\phi$ on $T^d$ may be replaced, without loss of generality, by the positivity condition

$$\begin{aligned} \sigma _\phi ( \zeta ) > 0,\quad \zeta \in T^d. \end{aligned}$$

(2.6)

The following result is well-known (e.g. Chui et al. [19, Lemma 1.1], Jia and Micchelli [39, Theorem 2.1]); for the sake of presentation, we will outline the main steps of the proof.

Theorem 2.1

Suppose that the continuous and symmetric kernel $\phi$ satisfies (2.1) and (2.6).

• (i) For any $p\in [1,\infty ]$, the Laurent operator $L_\phi$ is invertible in ${\mathcal {B}} (\ell _p)$, the Banach algebra of all bounded linear operators on $\ell _p$. Specifically, for any $y\in \ell _{\infty }$, the unique solution $c\in \ell _{\infty }$ of (2.5) is given explicitly by the formula

  $$\begin{aligned} c_k = \sum _{j\in {\mathbb {Z}}^d} a_{k-j} y_j,\quad k\in {\mathbb {Z}}^d, \end{aligned}$$

  (2.7)

  where $a=\{a_k\}_{k\in {\mathbb {Z}}^d}\in \ell _1$ is the sequence of Fourier coefficients of the inverse symbol $\omega _\phi := \sigma _\phi ^{-1}$. The corresponding interpolant s defined by (2.2) admits the pointwise absolutely convergent Lagrange representation (1.4).

• (ii) If $I_{\phi }$ denotes the interpolation operator assigning to each bounded sequence y its unique cardinal interpolant $I_{\phi }\, y:=s$ of the form (2.2), then $I_{\phi }$ is bounded from $\ell _p$ to $L_p ({\mathbb {R}}^d)$, and its operator norm satisfies

  $$\begin{aligned} \Vert I_{\phi } \Vert _p \le | \chi |_p \le \Vert \omega _\phi \Vert _W | \phi |_p, \quad 1\le p\le \infty , \end{aligned}$$

  (2.8)

  where $|\cdot |_p$ is defined in (2.3). For $p=\infty$, the operator norm (or Lebesgue constant) of $I_{\phi }$ further satisfies

  $$\begin{aligned} \Vert I_{\phi } \Vert _{\infty } = | \chi |_\infty . \end{aligned}$$

  (2.9)

Proof

(i) Note that $L_{\phi } \in {\mathcal {B}} (\ell _p)$ via Schur’s test (Gröchenig [32, Lemma 6.2.1]):

$$\begin{aligned} \Vert L_{\phi }\, c \Vert _p \le \Vert \sigma _\phi \Vert _W \Vert c \Vert _p, \quad c\in \ell _p,\ 1\le p\le \infty . \end{aligned}$$

By Wiener’s lemma (Rudin [52, Lemma 11.6]), we have $\omega _\phi = \sigma _\phi ^{-1}\in W$, hence, if $L_\phi ^{-1}:=[a_{j-k}]_{j,k\in {\mathbb {Z}}^d}$, another application of Schur’s test gives $L_\phi ^{-1} \in {\mathcal {B}} (\ell _p)$. The straightforward identities $L_{\phi } L_\phi ^{-1} y = y$ and $L_\phi ^{-1} L_\phi \, c =c$, valid for all $y,c \in \ell _\infty$, now imply the invertibility of $L_\phi$ in ${\mathcal {B}} (\ell _p)$, for $p\in [1,\infty ]$. The Lagrange representation (1.4) of the interpolant (2.2) follows via (1.2) and Fubini, since

$$\begin{aligned} \sum _{j,k\in {\mathbb {Z}}^d} |y_j a_{k-j} \phi (x-k)|= & {} \sum _{k\in {\mathbb {Z}}^d} |\phi (x-k)| \sum _{j\in {\mathbb {Z}}^d} |y_j| |a_{k-j}| \\\le & {} \Vert y\Vert _{\infty } \Vert \omega _\phi \Vert _W |\phi |_\infty \ < \ \infty . \end{aligned}$$

(ii) The two inequalities of (2.8) follow as special cases of Jia and Micchelli’s result [39, Theorem 2.1]) on periodizations of semi-discrete convolutions. When applied to the Lagrange representation (1.4), this result reads as:

$$\begin{aligned} \Vert I_{\phi }\, y \Vert _{L_p ( {\mathbb {R}}^d )} \le | \chi |_p \Vert y \Vert _p, \quad 1\le p\le \infty , \end{aligned}$$

implying the first inequality in (2.8). The second inequality is obtained directly by applying the cited result to the expression (1.2) of $\chi$. The equality (2.9) follows via a standard argument (de Boor and Ron [23, Proposition 2.2]). $\square$

Corollary 2.2

Under the hypotheses of Theorem 2.1, for any function s of the form (2.2), where c is a bounded sequence, the following reproduction (or sampling) formula holds:

$$\begin{aligned} s ( x ) =\sum _{j\in {\mathbb {Z}}^d} s(j) \chi ( x-j ),\quad x\in {\mathbb {R}}^d. \end{aligned}$$

(2.10)

In particular, $\phi =\sum _{j\in {\mathbb {Z}}^d} \phi (j) \chi ( \cdot -j )$.

2.2 Transfer of Kernel Decay to Lagrange Function

Here, we consider kernels with algebraic/exponential decay and we demonstrate how this carries over to the Lagrange function for cardinal interpolation. We begin by collecting the results needed on symbols $u\in W$, whose Fourier coefficients $\{{\widehat{u}}_k\}_{k\in {\mathbb {Z}}^d}$ decay algebraically/exponentially. For $x=(x_1,\ldots ,x_d)\in {\mathbb {R}}^d$, $\Vert x\Vert = (x_1^2 + \cdots + x_d^2)^{1/2}$ denotes the Euclidean norm and $|x| = |x_1| + \cdots + |x_d|$, the 1-norm.

Two inverse-closed classes of symbols

Definition 2.3

• (i) Let $\alpha >d$. A continuous symbol u on $T^d$ belongs to the class $J_\alpha :=J_\alpha (T^d)$ if there exists $C(u)>0$, such that

  $$\begin{aligned} |{\widehat{u}}_k| \le C(u) (1+ \Vert k\Vert )^{-\alpha }, \quad k\in {\mathbb {Z}}^d. \end{aligned}$$

  (2.11)

• (ii) For $\alpha >0$, a continuous symbol u on $T^d$ belongs to the class ${\mathcal {E}}_\alpha :={\mathcal {E}}_\alpha (T^d)$ if there exists $C(u)>0$, such that

  $$\begin{aligned} |{\widehat{u}}_k| \le C(u) \, e^{-\alpha |k|}, \quad k\in {\mathbb {Z}}^d. \end{aligned}$$

  (2.12)

  Also, define ${\mathcal {E}}:= \bigcup _{\alpha >0} {\mathcal {E}}_\alpha$.

The next lemma shows that $J_\alpha$ and ${\mathcal {E}}$ are subalgebras of the Wiener algebra W.

Lemma 2.4

• (i) If $u,v\in J_\alpha$, then $uv\in J_\alpha$, where the corresponding constant in (2.11) can be chosen as $C(uv):= 2^\alpha (C(u) \Vert v\Vert _W + C(v) \Vert u\Vert _W)$.

• (ii) If $u,v\in {\mathcal {E}}_\beta$, for some $\beta >0$, then $uv\in {\mathcal {E}}_{\beta '}$ for every $\beta ' \in (0,\beta )$, with a corresponding constant in (2.12) given by $C(uv):=C(u) C(v) \coth ^d ( \frac{\beta -\beta '}{2} )$.

Proof

This follows from a more general result of Jaffard [38, Proposition 1] concerning infinite matrices with off-diagonal decay. For the explicit decay constants provided in the lemma, we refer the reader to the first version of our paper, available as arXiv:2006.05282. Note that (i) was previously proved in Bacchelli et al. [1, Lemma 1], where it is equivalently formulated for sequences, endowed with the discrete convolution product. $\square$

The following result expresses the fact that $J_\alpha$ and ${\mathcal {E}}$ are inverse-closed subalgebras, i.e. Wiener’s lemma holds for these two algebras.

Lemma 2.5

• (i) If $u\in J_\alpha$ and u has no zeros on $T^d$, then $1/u\in J_\alpha$.

• (ii) If $u\in {\mathcal {E}}$ and u has no zeros on $T^d$, then $1/u\in {\mathcal {E}}$. More specifically, if $u\in {\mathcal {E}}_\alpha$ and $u(\zeta )>0$ for all $\zeta \in T^d$, then $1/u\in {\mathcal {E}}_\beta$, for some $\beta \in (0,\alpha )$.

Proof

Both parts can be deduced from the more general results of Jaffard [38, Propositions 3 and 2], concerning the class of invertible infinite matrices with off-diagonal algebraic or exponential decay, using the correspondence between symbols $u\in W$ and Laurent matrices $[ {\widehat{u}}_{j-k}]_{j,k\in {\mathbb {Z}}^d}$ generated with their Fourier coefficients. The proof of (i), formulated for sequences and their discrete convolution product, appeared in Bacchelli et al. [1, Lemma 2]. Fageot et al. [27, Theorem 13] gave an alternative proof of (ii) based on the fact that $u\in {\mathcal {E}}$ if and only if u is real analytic on $T^d$. For the second statement of (ii), a different proof, based on the complex analyticity arguments of Goodman et al. [31, Theorem 2.4], was given in the first version of our paper (arXiv:2006.05282). $\square$

Algebraic decay transfer. For the next result, we assume that $\phi$ has algebraic (or polynomial) decay, i.e. there exist constants $\alpha >d$, $C_0>0$, such that

$$\begin{aligned} |\phi ( x )| \le C_0 (1+ \Vert x\Vert )^{-\alpha }, \quad x\in {\mathbb {R}}^d. \end{aligned}$$

(2.13)

Clearly, a kernel $\phi$ with this property also satisfies (2.1).

Theorem 2.6

If the continuous symmetric kernel $\phi$ satisfies (2.13), with $\alpha >d$, and the cardinal symbol $\sigma _\phi$ satisfies the positivity condition (2.6), then both the sequence $\{a_k\}_{k\in {\mathbb {Z}}^d}$ of Fourier coefficients of the reciprocal symbol $\omega _\phi := 1 / \sigma _\phi$ and the associated Lagrange function $\chi$ defined by (1.2) decay algebraically with at least the same power as $\phi$, i.e. there exist positive constants $C_{1,2}=C_{1,2} (\alpha ,d,\phi )$, such that:

$$\begin{aligned}{} & {} |a_k| \le C_1 (1+ \Vert k\Vert )^{-\alpha }, \quad k\in {\mathbb {Z}}^d, \end{aligned}$$

(2.14)

$$\begin{aligned}{} & {} |\chi ( x )| \le C_2 (1+ \Vert x\Vert )^{-\alpha }, \quad x\in {\mathbb {R}}^d. \end{aligned}$$

(2.15)

Proof

Condition (2.13) ensures $\sigma \in J_{\alpha }$, therefore (2.6) and Lemma 2.5(i) imply that the sequence $\{a_k\}_{k\in {\mathbb {Z}}^d}$ of Fourier coefficients of $1/\sigma _\phi$ satisfies (2.14). The decay (2.15) of the Lagrange function follows from the result of Bacchelli et al. [1, Lemma 3] on semi-discrete convolutions with algebraic decay (for details, see the first version of our paper, at arXiv:2006.05282). $\square$

Remark

Theorem 2.6 can be adapted to the case of anisotropic decay with separate algebraic rates in each coordinate, using the inverse-closedness result proved by Gröchenig and Klotz [34, Theorem 1.1] for the corresponding subalgebra of matrices with anisotropic off-diagonal decay. Similarly, the transfer of the weighted $\ell _1$-condition (1.5) to the associated Lagrange function $\chi$ was proved by Fageot et al. [27, §7.1] via a weighted version of Wiener’s lemma (see Gröchenig [33, Theorem 5.24]).

Exponential decay transfer. Next, we consider the case of an exponentially decaying $\phi$, i.e. there exists a constant $C_0>0$ and a rate $\alpha >0$, such that

$$\begin{aligned} |\phi ( x )| \le C_0\, e^{-\alpha |x|}, \quad x\in {\mathbb {R}}^d. \end{aligned}$$

(2.16)

The use of the 1-norm in this condition is simply a matter of convenience.

Theorem 2.7

If the continuous symmetric kernel $\phi$ satisfies (2.16) and the associated cardinal symbol $\sigma _\phi$ satisfies (2.6), then there exist $\beta \in (0,\alpha )$, $C_1=C_1(\phi ,\beta ,d)$, $C_2=C_2(\phi ,\beta ,\alpha ,d)$, such that:

$$\begin{aligned}{} & {} |a_k| \le C_1\, e^{-\beta |k|}, \quad k\in {\mathbb {Z}}^d, \end{aligned}$$

(2.17)

$$\begin{aligned}{} & {} \vert \chi ( x ) \vert \le C_2\, e^{-\beta |x|}, \quad x\in {\mathbb {R}}^d. \end{aligned}$$

(2.18)

Proof

Since $\{ \phi (k) \}_{k\in {\mathbb {Z}}^d}$ is the sequence of Fourier coefficients of $\sigma _\phi$, we have $\sigma _\phi \in {\mathcal {E}}_\alpha$. Hence, by Lemma 2.5(ii), the Fourier coefficients $\{ a_k \}_{k\in {\mathbb {Z}}^d}$ of $1/\sigma _\phi$ satisfy (2.17), for some $\beta \in (0,\alpha )$. To establish the exponential decay (2.18) of the Lagrange function, we use (1.2), (2.16), and (2.17), which imply:

$$\begin{aligned} |\chi (x)|\le & {} \sum _{k\in {\mathbb {Z}}^d} |a_k|\, |\phi (x-k)| \\\le & {} C_1 C_0 \sum _{k\in {\mathbb {Z}}^d} \prod _{p=1}^d e^{-\beta |k_p|} e^{-\alpha |x_p - k_p|} \\= & {} C_1 C_0 \prod _{p=1}^d \sum _{k_p\in {\mathbb {Z}}} e^{-\beta |k_p|} e^{-\alpha |x_p - k_p|}, \quad x\in {\mathbb {R}}^d. \end{aligned}$$

The sum over $k_p$ is estimated using $-|k_p|\le |x_p-k_p|-|x_p|$:

$$\begin{aligned} \sum _{k_p\in {\mathbb {Z}}} e^{-\beta |k_p|} e^{-\alpha |x_p - k_p|} \le e^{-\beta |x_p|} \sum _{k_p\in {\mathbb {Z}}} e^{-(\alpha -\beta ) |x_p - k_p|}, \end{aligned}$$

and noting that the last sum admits the upper bound

$$\begin{aligned} \max _{x_p\in [0,1]} \left( e^{(\beta -\alpha ) x_p} + e^{(\beta -\alpha ) (1-x_p)} \right) \sum _{\nu =0}^\infty e^{(\beta -\alpha ) \nu } = \frac{1+e^{\beta -\alpha }}{1-e^{\beta -\alpha }} =: K. \end{aligned}$$

The conclusion follows, with $C_2:= C_1 C_0 K^d$. $\square$

3 Interpolation on a Half-Space Lattice

Half-space lattices. In the sequel, the expression ‘half-space lattice’ will be used to denote a subset H of ${\mathbb {Z}}^d$ belonging to one of two classes.

The first class is made of the following d lattice sets:

$$\begin{aligned} \{ (j_1,\ldots ,j_d)\in {\mathbb {Z}}^d: j_p \ge 0 \},\quad p\in \{1,\ldots ,d\}. \end{aligned}$$

Since the results obtained are formulated similarly for any of these d lattices, we will just consider $H = {\mathbb {Z}}^{d-1}\!\!\times \!{\mathbb {Z}}_+$ (for $p=d$) as a generic representative of this class of half-space lattices.

To describe a member of the second class, we assume that ${\mathbb {Z}}^d$ is endowed with a linear (or total) order relation $\preceq$ compatible with addition, in the sense that $j+l\preceq k+l$ whenever $j,k,l \in {\mathbb {Z}}^d$ and $j\preceq k$. We then let

$$\begin{aligned} H = {\mathbb {Z}}^{d}_{\preceq ,+}:= \{ j\in {\mathbb {Z}}^d :0\preceq j \}, \end{aligned}$$

the set of $\preceq$-nonnegative d-dimensional multi-integers. In this case, H satisfies the following properties:

• (a) $H+H\subset H$ (i.e. H is an additive semigroup);

• (b) $H\cap (-H) = \{0\}$;

• (c) $H\cup (-H) = {\mathbb {Z}}^d$.

Conversely, any set $H\subset {\mathbb {Z}}^d$ with these three properties induces a linear order $\preceq$ on ${\mathbb {Z}}^d$ compatible with addition, by simply defining $j\preceq k$ to mean that $k-j\in H$ (see Rudin [51, Chapter 8]).

Note that, if $d=1$, both types of half-space lattices reduce to $H = {\mathbb {Z}}_+$, while, if $d>1$, any lattice from the first class satisfies (a) and (c), but not (b). When $d>1$, the second class of lattices is infinite. The following two examples of order relations on ${\mathbb {Z}}^2$ have been used by Goodman et al. [31] in the context of Gram-Schmidt orthonormalization of the shifts of an integrable kernel:

• (i) The usual lexicographical order on ${\mathbb {Z}}^2$: $j\preceq k$ if either $j_1<k_1$, or $j_1 = k_1$ and $j_2\le k_2$;

• (ii) A different lexicographical order on ${\mathbb {Z}}^2$: $j\preceq k$ if either $j_1+j_2<k_1+k_2$, or $j_1+j_2 = k_1+k_2$ and $j_2\le k_2$.

Semi-cardinal interpolation. Let $H\subset {\mathbb {Z}}^d$ be a half-space lattice as above and suppose that the kernel $\phi$ satisfies the hypotheses of Theorem 2.1. The present section studies semi-cardinal interpolation on H with the kernel $\phi$, which is the problem of finding a bounded sequence of coefficients $c=\{c_k\}_{k\in H}$, such that the function s of the form (1.6) satisfies the interpolation conditions

$$\begin{aligned} s ( j ) =y_{j},\quad j\in H, \end{aligned}$$

(3.1)

for a given bounded sequence of real values $y=\{y_j\}_{j\in H}$. The solution, constructed in Sect. 3.2, will be expressed by means of the kernel expansion coefficients of Lagrange functions $\chi _j$, $j\in H$, defined by

$$\begin{aligned} \chi _{j} ( x ) = \sum _{k\in H} a_{k,j}\, \phi (x-k), \quad x\in {\mathbb {R}}^d, \end{aligned}$$

(3.2)

which satisfy the interpolation conditions

$$\begin{aligned} \chi _{j} ( l ) =\delta _{j,l} = \left\{ \begin{array}{ll} 1, &{}\quad l=j, \\ 0, &{}\quad l\in H \setminus \{j\}. \end{array} \right. \end{aligned}$$

(3.3)

This is equivalent to finding, for each $j\in H$, the sequence $\{a_{k,j}\}_{k\in H}$ satisfying the Wiener–Hopf (or semi-infinite Toeplitz) system of difference equations

$$\begin{aligned} \sum _{k\in H} a_{k,j} \phi (l-k) = \delta _{j,l}, \quad l\in H. \end{aligned}$$

(3.4)

In Sect. 3.1, the explicit solution of this system is obtained using the Wiener–Hopf factorization of the cardinal symbol $\sigma _\phi$ (equivalently, of its reciprocal) generated by $\phi$. This method is due, for $d=1$, to Krein [41] and Calderón et al. [18]. For $d>1$, the case $H = {\mathbb {Z}}^{d-1}\!\!\times \!{\mathbb {Z}}_+$ was treated by Goldenstein and Gohberg [30] and Goldenstein [29]. In the case $H = {\mathbb {Z}}^{d}_{\preceq ,+}$, the Wiener–Hopf factorization technique was studied by Goodman et al. [31] for obtaining Cholesky factorizations of bi-infinite Gram-Laurent matrices generated by shifts of multivariate kernels, while the application of this technique to the solution of more general multi-index semi-infinite block Toeplitz systems was considered in van der Mee et al. [44, p. 467].

The transfer of specific algebraic or exponential decay of the kernel $\phi$ to the family of Lagrange functions $\{\chi _j: j\in H\}$ is obtained in Sect. 3.3.

Notation remarks. For simplicity, in this section we remove the dependence on $\phi$ in most notation, e.g. we will use $\sigma :=\sigma _\phi$; nevertheless, this dependence remains implicit throughout the section. Also, for a d-dimensional vector x, we employ the partition notation $x=(x',x_d)\in {\mathbb {R}}^{d-1}\!\!\times \!{\mathbb {R}}$.

3.1 The Wiener–Hopf Factorization

The construction of the Wiener–Hopf (also known as ‘canonical’) factorization of $\omega =1/\sigma \in W$ is based on the fact that, under the hypotheses of Theorem 2.1, $\log \omega$ also belongs to the Wiener algebra. Indeed, this is ensured by the following multivariable version of Lévy’s extension of Wiener’s lemma, which is a special case (for the trivial weight) of a result proved by Goodman et. al. [31, Theorem 2.3] for certain weighted Wiener algebras.

Lemma 3.1

(Wiener–Lévy) If $\psi \in W$ and F is analytic in a neighborhood of the range of $\psi$, then $F\circ \psi \in W$.

Since $\omega =1/\sigma$ is positive on $T^d$ due to (2.6), letting $\psi :=\omega$ and $F:=\log$, it follows that $\log \omega \in W$, with an absolutely convergent Fourier expansion:

$$\begin{aligned} \log \omega ( \zeta ) =\sum _{k\in {\mathbb {Z}}^d} \lambda _{k} \zeta ^{k},\quad \zeta \in T^d. \end{aligned}$$

(3.5)

Note that $\lambda _{-k}=\lambda _{k}\in {\mathbb {R}}$, $\forall k\in {\mathbb {Z}}^d$, since the Fourier coefficients of $\omega$ are also real-valued and symmetric. The required Wiener–Hopf factorization of $\omega$ is now a consequence of the decomposition

$$\begin{aligned} \log \omega ( \zeta ) = \Lambda _+(\zeta ) + \Lambda _+(\zeta ^{-1}), \end{aligned}$$

which is obtained by selecting the function $\Lambda _+$ in the following way:

If $H = {\mathbb {Z}}^{d-1}\!\!\times \!{\mathbb {Z}}_+$, we set, as in Goldenstein and Gohberg [30],

$$\begin{aligned} \Lambda _{+} (\zeta ):= \frac{1}{2} \sum _{k'\in {\mathbb {Z}}^{d-1}} \lambda _{(k',0)} (\zeta ')^{k'} + \sum _{k_d=1}^\infty \sum _{k'\in {\mathbb {Z}}^{d-1}} \lambda _{k} \zeta ^{k},\quad \zeta \in T^d. \end{aligned}$$

(3.6)

If $H = {\mathbb {Z}}^{d}_{\preceq ,+}$, we define, as in Goodman et al. [31],

$$\begin{aligned} \Lambda _{+} (\zeta ):= \frac{1}{2} \lambda _0 + \sum _{k\succ 0} \lambda _{k} \zeta ^{k},\quad \zeta \in T^d. \end{aligned}$$

(3.7)

In both cases, letting

$$\begin{aligned} \omega _{+} (\zeta ):= \exp [ \Lambda _{+} (\zeta ) ], \end{aligned}$$

(3.8)

leads to the following Wiener–Hopf (or canonical) factorization:

$$\begin{aligned} \omega (\zeta ) = e^{\log \omega ( \zeta )} = e^{\Lambda _+(\zeta ) + \Lambda _+(\zeta ^{-1})} = \omega _+ (\zeta ) \omega _+ (\zeta ^{-1}), \quad \zeta \in T^d. \end{aligned}$$

(3.9)

Invoking again the Wiener–Lévy lemma, this time for the composite function (3.8), it follows that $\omega _+\in W$. Moreover, using the power series expansion of the exponential and the fact that H is closed under addition, we deduce that the Fourier coefficients of $\omega _+$, as those of $\Lambda _+$, are real-valued and supported on H:

$$\begin{aligned} \omega _{+} (\zeta ) = \sum _{k\in H} \gamma _{k} \zeta ^{k}, \quad \zeta \in T^d. \end{aligned}$$

(3.10)

The same arguments show that $\omega _+^{-1}=\exp (-\Lambda _+)\in W$ and its Fourier coefficients are real-valued and supported on H, as well.

Remark

The existence of a multi-index Wiener–Hopf factorization of a non-vanishing symbol $\omega$ is obtained in the literature under conditions more general than positivity, expressed in terms of the winding number $wind _p (\omega )$ of $\omega$ with respect to $\zeta _p$ about the origin. Specifically, if $H = {\mathbb {Z}}^{d-1}\!\!\times \!{\mathbb {Z}}_+$, then $wind _d (\omega ) = 0$ is required (Goldenstein and Gohberg [30]), while, if $H = {\mathbb {Z}}^{d}_{\preceq ,+}$, the extra condition is $wind _p (\omega ) = 0$, for all $p = 1,\ldots ,d$, which is seen to be independent of the underlying linear order (Ehrhardt and van der Mee [26]).

Semi-cardinal symbols. Let $P_+$ denote the truncation operator defined by

$$\begin{aligned} P_+ \left\{ \sum _{k\in {\mathbb {Z}}^{d}} {\widehat{\psi }}_{k} \zeta ^{k} \right\} = \sum _{k\in H} {\widehat{\psi }}_{k} \zeta ^{k}, \quad \zeta \in T^d, \end{aligned}$$

for any $\psi \in W$ with Fourier coefficients $\{{\hat{\psi }}_k\}_{k\in {\mathbb {Z}}^d}$.

With the factorization (3.9) in hand, we now follow Krein’s method [41] (for $d=1$) to express, for each $j\in H$, the solution $\{a_{k,j}\}_{k\in H}$ of the Wiener–Hopf system (1.7) as the sequence of Fourier coefficients of the symbol:

$$\begin{aligned} \omega _{j} (\zeta ):= \omega _{+} (\zeta ) P_{+}\left\{ \zeta ^{j} \omega _{+} (\zeta ^{-1}) \right\} , \quad \zeta \in T^d. \end{aligned}$$

(3.11)

The second factor can be seen as a twisted truncation of the symbol $\omega _{+} (\zeta ^{-1})$, which appears in the Wiener–Hopf factorization (3.9). Indeed, if $H = {\mathbb {Z}}^{d}_{\preceq ,+}$, then for each $j\succeq 0$, we have

$$\begin{aligned} P_+ \{ \zeta ^j \omega _+ (\zeta ^{-1}) \} = \zeta ^j \omega _{-}^{(j)} (\zeta ), \end{aligned}$$

(3.12)

where $\omega _{-}^{(j)}$ is the truncation of the symbol $\omega _+ (\zeta ^{-1})$ defined by

$$\begin{aligned} \omega _{-}^{(j)} (\zeta ):= \sum _{0\preceq l\preceq j} \gamma _l \zeta ^{-l}. \end{aligned}$$

(3.13)

Note that the last sum can have an infinite number of terms when $d>1$. More importantly, the following uniform bound holds:

$$\begin{aligned} \Vert \omega _{-}^{(j)} \Vert _W \le \Vert \omega _{+} \Vert _W, \quad \forall j\in H. \end{aligned}$$

(3.14)

Also, if $H = {\mathbb {Z}}^{d-1}\!\!\times \!{\mathbb {Z}}_+$, replacing (3.13) by the truncation with respect to the last coordinate:

$$\begin{aligned} \omega _{-}^{(j_d)} (\zeta ):= \sum _{l_d = 0} ^{j_d} \sum _{l'\in {\mathbb {Z}}^{d-1}} \gamma _l \zeta ^{-l}, \end{aligned}$$

(3.15)

it follows that (3.12) and (3.14) hold with $\omega _{-}^{(j_d)}$ in place of $\omega _{-}^{(j)}$.

Clearly, for both choices of H, we have $\omega _j \in W$ and the Fourier coefficients of $\omega _j$ are supported on H, since both factors of (3.11) possess this property. Hence, $\omega _j$ admits an absolutely convergent Fourier representation:

$$\begin{aligned} \omega _{j} (\zeta ) = \sum _{k\in H} a_{k,j}\, \zeta ^{k}, \quad \zeta \in T^d. \end{aligned}$$

(3.16)

However, it may not be immediately obvious that $\{a_{k,j}\}_{k\in H}$ is a solution of (3.4). This claim is settled next.

Proposition 3.2

For each $j\in H$, the sequence $\{a_{k,j}\}_{k\in H}$ of Fourier coefficients of $\omega _j$ satisfies the semi-infinite system (3.4).

Proof

Let ${\overline{H}}:={\mathbb {Z}}^{d}\setminus H$. Using the absolutely convergent Fourier expansion (3.10) of $\omega _+$ and the Wiener–Hopf factorization (3.9), definition (3.11) implies:

$$\begin{aligned} \omega _j (\zeta )= & {} \omega _+ (\zeta ) \, P_{+}\left\{ \zeta ^{j} \sum _{k\in H} \gamma _{k} \zeta ^{-k} \right\} \\= & {} \omega _+ (\zeta ) \left\{ \zeta ^{j} \omega _{+} (\zeta ^{-1}) - \sum _{l\in {\overline{H}}} \gamma _{j-l} \zeta ^{l} \right\} \\= & {} \zeta ^{j} \omega (\zeta ) - \omega _+ (\zeta ) \sum _{l\in {\overline{H}}} \gamma _{j-l} \zeta ^{l}, \end{aligned}$$

or, equivalently,

$$\begin{aligned} \omega _j (\zeta ) \sigma (\zeta ) = \zeta ^{j} - \omega _+^{-1} (\zeta ^{-1}) \sum _{l\in {\overline{H}}} \gamma _{j-l} \zeta ^{l}. \end{aligned}$$

Since $j\in H$ and the second term of the right-hand side admits an absolutely convergent expansion indexed over ${\overline{H}}$, it follows that

$$\begin{aligned} P_+ \left\{ \omega _j (\zeta ) \sigma (\zeta ) \right\} = \zeta ^{j}. \end{aligned}$$

On the other hand, using the Fourier expansions of $\omega _j$ and $\sigma$, we have, via Fubini:

$$\begin{aligned} P_+ \left\{ \omega _j (\zeta ) \sigma (\zeta ) \right\}= & {} \sum _{k\in H} a_{k,j}\, \zeta ^{k} \sum _{l\in H} \phi (l-k)\, \zeta ^{l-k} \\= & {} \sum _{l\in H} \sum _{k\in H} a_{k,j}\, \phi (l-k)\, \zeta ^{l}. \end{aligned}$$

The system (3.4) now follows by comparing the coefficient of $\zeta ^l$ in the right-hand sides of the last two displays. $\square$

Next, we obtain an explicit formula for the Fourier coefficients of the semi-cardinal symbols $\omega _j$, $j\in H$, in terms of those of the Wiener–Hopf factor $\omega _+$, from which we deduce the symmetry, as well as a crucial Schur-type property of these coefficients. As remarked by the referee, an alternative proof of the latter property can also be given via Wiener’s lemma for the Baskakov-Gohberg-Sjöstrand matrix algebra (see [42, 57]).

Proposition 3.3

Assume that $\phi$ satisfies the hypotheses of Theorem 2.1.

• (i) If $H = {\mathbb {Z}}^{d}_{\preceq ,+}$, then the coefficients of the symbol $\omega _j$ are given by

  $$\begin{aligned} a_{k,j} = \sum _{ 0\preceq l\preceq \min \{ j,k \} } \gamma _{k-l} \gamma _{j-l}, \quad j,k\in {\mathbb {Z}}^{d}_{\preceq ,+}. \end{aligned}$$

  (3.17)

  If $H = {\mathbb {Z}}^{d-1}\!\!\times \!{\mathbb {Z}}_{+}$, then, for all $j=(j',j_d)$ and $k=(k',k_d)$ in H,

  $$\begin{aligned} \omega _{(j',j_d)} (\zeta )= & {} (\zeta ')^{j'} \omega _{(0,j_d)} (\zeta ), \quad \zeta \in T^d, \end{aligned}$$

  (3.18)
  $$\begin{aligned} a_{k,j}= & {} a_{(k'-j',k_d),(0,j_d)}, \end{aligned}$$

  (3.19)
  $$\begin{aligned} a_{k,j}= & {} \sum _{l_d=0}^{\min \{j_d,k_d\}} \sum _{l'\in {\mathbb {Z}}^{d-1}} \gamma _{(k'-l',k_d-l_d)} \gamma _{(j'-l',j_d-l_d)}. \end{aligned}$$

  (3.20)

• (ii) For both types of half-space lattices H, we have $a_{k,j} = a_{j,k}$, for all $j,k \in H$, and the following Schur property holds:

  $$\begin{aligned} \sup _{j\in H} \sum _{k\in H} |a_{k,j}| < \infty . \end{aligned}$$

  (3.21)

Proof

• (i) Let $H = {\mathbb {Z}}^{d}_{\preceq ,+}$. Using (3.10) and (3.11), we have, for each $j\in {\mathbb {Z}}^{d}_{\preceq ,+}$,

  $$\begin{aligned} \omega _j (\zeta )= & {} \sum _{k\succeq 0} \gamma _k \zeta ^k \sum _{0\preceq l\preceq j} \gamma _{j-l} \zeta ^l \\= & {} \sum _{k\succeq 0} \zeta ^k \sum _{ 0\preceq l\preceq \min \{ j,k \} } \gamma _{k-l} \gamma _{j-l}, \end{aligned}$$

  which implies (3.17). For $H = {\mathbb {Z}}^{d-1}\!\!\times \!{\mathbb {Z}}_{+}$, relation (3.18) follows from (3.11), since, for any $p\in \{1,\ldots ,d-1\}$, multiplication by $\zeta _p$ commutes with the truncation operator $P_+$. Also, (3.18) implies (3.19). To obtain (3.20), we use re-arrangements of absolutely convergent d-dimensional Fourier series with respect to the last component of the summation multi-index, as detailed in the first version of the paper (arXiv:2006.05282).

• (ii) The symmetry of the coefficients follows directly from the explicit formulae (3.17) and (3.20). Note that (3.21) is equivalent to $\sup _{j\in H} \Vert \omega _j \Vert _W < \infty$. To prove this, we use the expression (3.11) of $\omega _j$ as a product, together with the Wiener algebra property $\Vert uv\Vert _W \le \Vert u\Vert _W \Vert v\Vert _W$, valid for any $u,v \in W$. Specifically, if $H = {\mathbb {Z}}^{d}_{\preceq ,+}$, then, due to (3.12) and the uniform bound (3.14) on the Wiener norm of the twisted truncation factor of (3.11), we have

  $$\begin{aligned} \Vert \omega _j\Vert _W \le \Vert \omega _+\Vert _W \Vert \omega _{-}^{(j)} \Vert _W \le \Vert \omega _+\Vert _W^2 < \infty , \quad j\in H. \end{aligned}$$

  (3.22)

  A similar proof applies to the case $H = {\mathbb {Z}}^{d-1}\!\!\times \!{\mathbb {Z}}_{+}$, since (3.12) and (3.14) hold with $\omega _{-}^{(j_d)}$ in place of $\omega _{-}^{(j)}$. $\square$

3.2 The Semi-Cardinal Interpolation Scheme

Semi-cardinal Lagrange functions. For each $j\in H$, the Fourier coefficients $\{a_{k,j}\}_{k\in H}$ of the symbol $\omega _j$ are now used to define $\chi _j$ via (3.2). That series is absolutely and uniformly convergent on ${\mathbb {R}}^d$, due to the boundedness of $\phi$ and the absolute summability of $\{a_{k,j}\}_{k\in H}$. Further, Proposition 3.2 shows that $\chi _j$, $j\in H$, are the expected Lagrange functions for interpolation on H, since they satisfy the interpolation conditions (3.3).

Remark

When $H = {\mathbb {Z}}^{d-1}\!\!\times \!{\mathbb {Z}}_{+}$, property (3.19) implies, via re-indexing, the translation symmetry of $\chi _j$ for shifts parallel to the boundary of the half-space domain $\Omega = {\mathbb {R}}^{d-1} \times [0,\infty )$, as follows:

$$\begin{aligned} \chi _{(j',j_d)} (x',x_d) = \chi _{(0,j_d)} (x' - j', x_d),\quad \forall x\in {\mathbb {R}}^d,\ j\in {\mathbb {Z}}^{d-1}\times {\mathbb {Z}}_{+}. \end{aligned}$$

(3.23)

Next, we consider the convolution operator represented by the semi-infinite multi-index Toeplitz matrix $T_{\phi }:= [\phi (j-k)]_{j,k\in H}$, i.e. for each bounded sequence $c=\{c_k\}_{k\in H}$, we set $T_{\phi } c:= y=\{y_j\}_{j\in H}$, where $y_j = \sum _{k\in H} \phi (j-k) c_k$, $j\in H$. The following result proves the basic properties of $T_\phi$ and of the semi-cardinal interpolation scheme under minimal conditions on the kernel $\phi$. We denote by ${\mathcal {B}} (\ell _p(H))$ the Banach algebra of all bounded linear operators on $\ell _p(H)$.

Theorem 3.4

Assume that $\phi$ satisfies the hypotheses of Theorem 2.1, and let the half-space lattice H be either ${\mathbb {Z}}^{d-1}\!\!\times \!{\mathbb {Z}}_+$ or ${\mathbb {Z}}^{d}_{\preceq ,+}$.

• (i) For any $p\in [1,\infty ]$, the Toeplitz operator $T_\phi$ is invertible in ${\mathcal {B}} (\ell _p(H))$. Specifically, for any $\{y_j\}_{j\in H}\in \ell _\infty (H)$, the unique solution $\{c_j\}_{j\in H}\in \ell _\infty (H)$ of the multi-index Wiener–Hopf system (1.7) is given explicitly by:

  $$\begin{aligned} c_k = \sum _{j\in H} a_{k,j} y_j,\quad k\in H, \end{aligned}$$

  (3.24)

  where $\{a_{k,j}\}_{k\in H}\in \ell _1(H)$ is the coefficient sequence of the Fourier representation (3.16) of the symbol $\omega _j$, $j\in H$. The corresponding bounded semi-cardinal interpolant s defined by (1.6) admits the pointwise absolutely convergent Lagrange representation

  $$\begin{aligned} s ( x ) =\sum _{j\in H} y_{j} \chi _{j} ( x ), \quad x\in {\mathbb {R}}^d. \end{aligned}$$

  (3.25)

• (ii) If $I_{\phi }^H: \ell _{\infty } (H)\rightarrow L_{\infty } ({\mathbb {R}}^d)$ denotes the linear operator assigning to each $y\in \ell _{\infty } (H)$ its unique semi-cardinal interpolant $I_{\phi }^H y:=s$ of the form (1.6), then $I_{\phi }^H$ is bounded from $\ell _p(H)$ to $L_p({\mathbb {R}}^d)$ and its operator norm satisfies

  $$\begin{aligned} \Vert I_{\phi }^H \Vert _p \le |\phi |_p \sup _{j\in H} \Vert \omega _j\Vert _W, \quad 1\le p \le \infty , \end{aligned}$$

  (3.26)

  where $|\phi |_p$ is defined in (2.3). For $p=\infty$, the operator norm (or Lebesgue constant) of $I_{\phi }^H$ can be expressed as

  $$\begin{aligned} \Vert I_{\phi }^H \Vert _\infty = \sup _{x\in {\mathbb {R}}^d} \sum _{j\in H} |\chi _j ( x )|. \end{aligned}$$

  (3.27)

• (iii) Let $G = [G_{kj}]_{k,j\in H}$ be the matrix defined by $G_{kj} = \gamma _{k-j}$ for $k,j\in H$, where $\{\gamma _k\}_{k\in H}$ provides the Fourier representation (3.10) of $\omega _+$, and $\gamma _k:= 0$ for $k\in {\mathbb {Z}}^d\setminus H$. Then the inverse matrix $T_{\phi }^{-1} = [a_{k,j}]_{k,j\in H}$ admits the factorization:

  $$\begin{aligned} T_{\phi }^{-1} = G G^{T}, \end{aligned}$$

  (3.28)

  where G is also invertible in ${\mathcal {B}}(\ell _p(H))$ and $G^T$ is the transpose of G.

Proof

We argue along the lines of the proof of Theorem 2.1, the difference being that the finiteness of the Wiener norm $\Vert \omega \Vert _W$ is replaced here by the uniform bound (3.22) on $\Vert \omega _j\Vert _W$, for $j\in H$, i.e. the Schur property (3.21).

• (i) Schur’s test for infinite matrices (Gröchenig [32, Lemma 6.2.1]) provides the estimate

  $$\begin{aligned} \Vert T_{\phi } c \Vert _p \le \Vert c \Vert _p\, \sup _{j\in H} \sum _{k\in H} |\phi (j-k)| \le \Vert c \Vert _p \sum _{k\in {\mathbb {Z}}^d} |\phi (k)|, \quad c\in \ell _p (H), \end{aligned}$$

  hence $T_{\phi } \in {\mathcal {B}} (\ell _p (H))$. Letting $T_{\phi }^{-1}:= [a_{k,j}]_{k,j\in H}$, Proposition 3.3(ii) enables another application of Schur’s test:

  $$\begin{aligned} \Vert T_{\phi }^{-1} y \Vert _p \le \Vert y \Vert _p\, \sup _{j\in H} \sum _{k\in H} |a_{k,j}| = \Vert y \Vert _p\, \sup _{j\in H} \Vert \omega _j\Vert _W, \quad y\in \ell _p (H), \qquad \end{aligned}$$

  (3.29)

  showing that $T_{\phi }^{-1} \in {\mathcal {B}} (\ell _p (H))$. It is now straightforward to verify the identities $T_{\phi } T_\phi ^{-1} y = y$ and $T_\phi ^{-1} T_\phi c =c$, for all $y,c \in \ell _\infty (H)$, which imply the invertibility of $T_\phi$ in ${\mathcal {B}} (\ell _p(H))$, for all $p\in [1,\infty ]$. The Lagrange representation (3.25) of the semi-cardinal interpolant (1.6) follows via (3.2) and Fubini, since

  $$\begin{aligned} \sum _{j,k\in H} |y_j a_{k,j} \phi (x-k)|= & {} \sum _{k\in H} |\phi (x-k)| \sum _{j\in H} |y_j| |a_{k,j}| \\\le & {} \Vert y\Vert _{\infty } \left( \sup _{k\in H} \sum _{j\in H} |a_{k,j}| \right) \sup _{x\in {\mathbb {R}}^d} \sum _{k\in H} |\phi (x-k)| \\\le & {} \Vert y\Vert _{\infty } \left( \sup _{k\in H} \Vert \omega _k\Vert _W \right) |\phi |_\infty , \end{aligned}$$

  where we have used the symmetry property $a_{k,j} = a_{j,k}$ and (3.21).

• (ii) For $p = \infty$, (3.26) follows directly from the last estimate, which implies

  $$\begin{aligned} \Vert I_{\phi }^H \Vert _{\infty } \le \sup _{x\in {\mathbb {R}}^d} \Lambda _H (x) \le |\phi |_\infty \sup _{k\in H} \Vert \omega _k\Vert _W, \end{aligned}$$

  where $\Lambda _H (x):=\sum _{j\in H} | \chi _j ( x ) |$ is the bounded ‘Lebesgue function’ for semi-cardinal interpolation on H. To prove the equality (3.27), let $\{x_n\}_{n=1}^\infty$ be a sequence of points in ${\mathbb {R}}^d$, such that $\Lambda _H (x_n) \rightarrow \sup _{x\in {\mathbb {R}}^d} \Lambda _H (x)$, as $n\rightarrow \infty$. For each n, choose a sequence $y^{(n)}=\{y^{(n)}_j\}_{j\in H}$ with terms $y^{(n)}_j = \pm 1$, $\forall j\in H$, such that $y^{(n)}_j \chi _j (x_n) = | \chi _j (x_n) |$, $\forall j\in H$. Since $\Vert y^{(n)} \Vert _\infty = 1$, we have

  $$\begin{aligned} \Lambda _H (x_n) = \sum _{j\in H} y^{(n)}_j \chi _j (x_n) = I_{\phi }^H y^{(n)} (x_n) \le \Vert I_{\phi }^H y^{(n)} \Vert _{\infty } \le \Vert I_{\phi }^H \Vert _{\infty }. \end{aligned}$$

  The proof of (3.27) is completed by letting $n\rightarrow \infty$. For $1\le p < \infty$, the estimate (3.26) is obtained by applying Jia and Micchelli’s result [39, Theorem 2.1]) to the kernel representation (1.6):

  $$\begin{aligned} \Vert I_{\phi }^H y \Vert _{ L_p({\mathbb {R}}^d) } \le | \phi |_p \Vert c \Vert _p, \end{aligned}$$

  and using the bound (3.29) on the $\ell _p (H)$-norm of $c=T_\phi ^{-1} y$.

• (iii) Note that, setting $\gamma _k:= 0$ for $k\not \in H$, the explicit formulae (3.17) and (3.20) both take the same form:

  $$\begin{aligned} a_{k,j} = \sum _{l\in H} \gamma _{k-l} \gamma _{j-l}, \quad k,j\in H. \end{aligned}$$

  Hence, the definition of the matrix G implies $a_{k,j} = \sum _{l\in H} G_{kl} G_{jl} = (GG^{T})_{kj}$ for all $k,j\in H$, which proves (3.28). Another application of the Schur test shows that $G\in {\mathcal {B}}(\ell _p (H))$. Further, one can verify that the bounded inverse of G is given by $G^{-1} = [{\widetilde{\gamma }}_{k-j}]_{k,j\in H}$, where $\{{\widetilde{\gamma }}_k\}_{k\in H}$ is the sequence of Fourier coefficients of the symbol $1/ \omega _+$ and ${\widetilde{\gamma }}_k:= 0$ for $k\in {\mathbb {Z}}^d\setminus H$. $\square$

Remark

• (i) As remarked by the referee, for $p=2$, the invertibility of the Toeplitz operator $T_\phi$ in ${\mathcal {B}} (\ell _2(H))$ can also be deduced from its positive definiteness as an operator on the Hilbert space $\ell _2(H)$. The latter property follows from the strict positivity condition (2.6) of the cardinal symbol $\sigma _\phi$.

• (ii) For $d=1$, the factorization (3.28) was given by Krein [41, (13.27)]. For $d>1$, there is an important distinction between the two types of half-space lattices. Indeed, if $H = {\mathbb {Z}}^{d}_{\preceq ,+}$, for a linear order $\preceq$, then the matrix G and its inverse $G^{-1}$ are lower triangular, in the sense that $G_{kj} = 0 = (G^{-1})_{kj}$ for $k\prec j$. In this case, (3.28) provides the Cholesky factorization of the inverse matrix $T_{\phi }^{-1}$ relative to $\preceq$. However, if $H = {\mathbb {Z}}^{d-1}\!\!\times \!{\mathbb {Z}}_+$, it can be readily seen that the lower triangular structure of G is lost.

The next result shows that the semi-cardinal interpolation scheme can be used to recover the kernel $\phi$ from its values sampled at a half-space lattice H.

Corollary 3.5

For a function s of the form (1.6), where $c\in \ell _{\infty } ( H)$, the following reproduction/sampling formula holds:

$$\begin{aligned} s ( x ) =\sum _{j\in H} s(j) \chi _j ( x ),\quad x\in {\mathbb {R}}^d. \end{aligned}$$

In particular, this applies to the kernel shifts $s:=\phi (\cdot - k)$, with $k\in H$, hence, for $k=0$, we have $\phi ( x ) =\sum _{j\in H} \phi (j) \chi _j ( x )$, for all $x\in {\mathbb {R}}^d$.

3.3 Decay Transfer

We now impose specific algebraic or exponential decay conditions on $\phi$ of the form considered in Sect. 2.2 and we prove that these carry over to all semi-cardinal Lagrange functions $\chi _j$, $j\in H$, with constants independent of j. To this aim, we will first obtain the requisite decay transfer to the Fourier coefficients of the semi-cardinal symbols $\omega _j$, $j\in H$, i.e. to the inverse matrix $T_{\phi }^{-1} = [a_{k,j}]_{k,j\in H}$.

Algebraic decay

Theorem 3.6

Assume that the continuous and symmetric kernel $\phi$ has algebraic decay of power $\alpha >d$, as in (2.13), that the cardinal symbol $\sigma _\phi$ satisfies the positivity condition (2.6), and that H is a half-space lattice of one of the two types considered in Sect. 3.1. Then:

• (i) The Fourier coefficients of the Wiener–Hopf factor $\omega _+$ and of the semi-cardinal symbols $\{\omega _j: j\in H\}$ decay with the same power $\alpha$, i.e. there exists $C_{1,2}=C_{1,2} (\alpha ,d,\phi )>0$, such that:

  $$\begin{aligned}{} & {} | \gamma _{k} | \le C_1 (1+ \Vert k\Vert )^{-\alpha }, \quad k\in H, \end{aligned}$$

  (3.30)
  $$\begin{aligned}{} & {} | a_{k,j} | \le C_2 (1+\Vert k-j\Vert )^{-\alpha },\quad j,k\in H. \end{aligned}$$

  (3.31)

• (ii) The semi-cardinal Lagrange functions $\{\chi _j:j\in H\}$ also decay with the same power $\alpha$, i.e. there exists $C_{3}=C_{3}(\alpha ,d,\phi )$, such that:

  $$\begin{aligned} | \chi _j ( x ) | \le C_3 (1+\Vert x-j\Vert )^{-\alpha }, \quad x\in {\mathbb {R}}^d,\ j\in H. \end{aligned}$$

  (3.32)

The key point about the decay estimates (3.31) and (3.32) is that $C_2$ and $C_3$ are independent of j. To prove (3.30), we need the next Wiener–Lévy type lemma for the algebra $J_\alpha$.

Lemma 3.7

If $\psi \in J_\alpha$, for $\alpha >d$, and F is analytic in a neighborhood of the range of $\psi$, then $F\circ \psi \in J_\alpha$.

Proof

This follows from Sun’s result [60, Theorem 5.3], specialized to the Jaffard algebra $Q_{\alpha }$ of matrices with polynomial off-diagonal decay of power $\alpha$ considered in [38]. The details are provided in the first version of our paper (arXiv:2006.05282). $\square$

Proof of Theorem 3.6

(i) By Theorem 2.6, the reciprocal symbol $\omega = 1 / \sigma \in J_{\alpha }$, therefore by Lemma 3.7 we obtain that $\log \omega \in J_{\alpha }$, as well. In turn, this implies that the function $\Lambda _+$ defined by either (3.6) or (3.7) satisfies $\Lambda _+ \in J_{\alpha }$. A new application of Lemma 3.7, with $F:=\exp$, shows that $\omega _+:= \exp (\Lambda _+) \in J_{\alpha }$, hence its Fourier coefficients satisfy (3.30), as stated.

For (3.31), the first version of this paper (arXiv:2006.05282) contains two proofs. The first is constructive, using the explicit product form (3.11) of $\omega _j$ and Lemma 2.4(i), i.e. the fact that $J_\alpha$ is a subalgebra of the Wiener algebra. For brevity, we include here only the second proof, based on Wiener’s lemma for Jaffard’s algebra $Q_\alpha (H)$ of matrices indexed by $H\times H$ with polynomial off-diagonal decay of power $\alpha$. Specifically, since (2.13) implies $\phi (j-k) = O ((1+\Vert j-k\Vert )^{-\alpha })$, for all $j,k\in H$, it follows that $T_\phi \in Q_\alpha (H)$. By Theorem 3.4, $T_\phi$ is also invertible as a bounded linear operator on $\ell _2 (H)$. Then [38, Proposition 3] (i.e. the inverse-closedness of $Q_\alpha (H)$) implies $T_\phi ^{-1}=[a_{k,j}]_{k,j\in H} \in Q_\alpha (H)$, which is equivalent to (3.31).

(ii) For $j\in H$, we use the kernel representation (3.2) of $\chi _j$, along with (2.13) and (3.31), to estimate

$$\begin{aligned} |\chi _j (x)|\le & {} \sum _{k\in H} |a_{k,j}|\, |\phi (x-k)| \\\le & {} C_2 C_0 \sum _{k\in H} (1+\Vert k-j\Vert )^{-\alpha } (1+\Vert x-k\Vert )^{-\alpha }. \end{aligned}$$

The above series admits the upper bound:

$$\begin{aligned} \sum _{k\in {\mathbb {Z}}^{d}} (1+\Vert k-j\Vert )^{-\alpha } (1+\Vert x-k\Vert )^{-\alpha }= & {} \sum _{k\in {\mathbb {Z}}^{d}} (1+\Vert k\Vert )^{-\alpha } (1+\Vert x-j-k\Vert )^{-\alpha } \\\le & {} C(\alpha ,d) (1+\Vert x-j\Vert )^{-\alpha }, \end{aligned}$$

where the last inequality follows by replacing x with $x-j$ in the corresponding estimate that proves (2.15) of Theorem 2.6, for cardinal interpolation. This establishes (3.32). $\square$

Remark

Theorem 3.6 can be adapted, as in the case of cardinal interpolation, to kernels $\phi$ with anisotropic algebraic decay, based on the version of Wiener’s lemma given in Gröchenig and Klotz [34, Theorem 1.1] for matrices indexed on $H\times H$. Similarly, for w satisfying a weak growth condition, the weighted $\ell _1$-condition (1.5) is transferred to the inverse Toeplitz matrix $T_\phi ^{-1}$ via the inverse-closed weighted Schur algebra ${\mathcal {A}}^1_w$ studied by Gröchenig and Leinert [35] (Remark 2 before Sect. 4) and by Sun [60, Theorem 4.1]. In this case, the decay induced on the semi-cardinal Lagrange functions is expressed via the finiteness of a weighted Lebesgue constant, namely:

$$\begin{aligned} \sup _{x\in {\mathbb {R}}^d} \sum _{j\in H} | \chi _j (x) | w(j) < \infty . \end{aligned}$$

Exponential decay

Theorem 3.8

Assume that the continuous and symmetric kernel $\phi$ decays exponentially with rate $\alpha >0$ as in (2.16), that the cardinal symbol $\sigma _\phi$ satisfies the positivity condition (2.6), and that H is a half-space lattice of one of the two types considered in Sect. 3.1. Then:

• (i) There exist constants $\beta \in (0,\alpha )$, $C_1=C_1(\beta ,\phi ,H,d)$, such that the Fourier coefficients of the Wiener–Hopf factor $\omega _+$ satisfy:

  $$\begin{aligned} \vert \gamma _{k} \vert \le C_1 \, e^{-\beta |k|}, \quad k\in H. \end{aligned}$$

  (3.33)

  Also, for each $\beta '\in (0,\beta )$, there exists $C_{2}=C_{2}(\beta ,\beta ',\phi ,H,d)$, such that the Fourier coefficients of the semi-cardinal symbols $\{\omega _j\}$ satisfy:

  $$\begin{aligned} \vert a_{k,j} \vert \le C_2 \, e^{-\beta ' |k-j|}, \quad j,k\in H, \end{aligned}$$

  (3.34)

  where the exponential factor can be replaced by $e^{-\beta ' |k'-j'|-\beta |k_d-j_d|}$, if $H={\mathbb {Z}}^{d-1}\times {\mathbb {Z}}_{+}$.

• (ii) For each $\beta '\in (0,\beta )$, there exists a constant $C_{3}=C_{3}(\beta ,\beta ',\phi ,H,d)$, such that the semi-cardinal Lagrange functions $\{\chi _j: j\in H\}$ satisfy:

  $$\begin{aligned} \vert \chi _j ( x ) \vert \le C_3 \, e^{-\beta ' |x-j|}, \quad x\in {\mathbb {R}}^d,\ j\in H. \end{aligned}$$

  (3.35)

  Further, if $H={\mathbb {Z}}^{d-1}\times {\mathbb {Z}}_{+}$, the last exponential factor can be replaced by $e^{-\beta ' |x'-j'|-\beta |x_d-j_d|}$.

Note that, crucially, $C_2$ and $C_3$ of (3.34) and (3.35) are independent of j.

Proof

(i) By reducing $\alpha$ if necessary, both the cardinal symbol $\sigma (\zeta )$ and its reciprocal $\omega (\zeta )=1/\sigma (\zeta )$, can be extended as analytic functions of d complex variables, with a positive real part in the neighborhood ${\mathcal {A}}^d_\alpha$ of the unit torus $T^d$, where

$$\begin{aligned} {\mathcal {A}}_\alpha :=\{z\in {\mathbb {C}}:e^{-\alpha }<|z|<e^{\alpha }\}. \end{aligned}$$

Hence $\log \omega (\zeta )$ is also analytic and admits an absolutely convergent Laurent-Fourier expansion (3.5) for $\zeta \in {\mathcal {A}}^d_\alpha$, with exponentially decaying coefficients $\{\lambda _k\}_{k\in {\mathbb {Z}}^d}$. It follows that $\Lambda _+$, as defined by (3.6) or (3.7), is also analytic on ${\mathcal {A}}^d_\alpha$, and the same property applies to $\omega _+ = \exp (\Lambda _+)$. We deduce that $\omega _+$ has an absolutely convergent Laurent-Fourier representation on ${\mathcal {A}}^d_\alpha$ with coefficients indexed by H, hence, by Cauchy’s estimates, (3.33) holds for some $\beta \in (0,\alpha )$ and $C_1:=\max \{ |\omega _+(\zeta )|: |\zeta _\tau |=e^{\pm \beta }, \tau = 1,\ldots ,d\}$. (The above arguments have previously been used by Goodman et al. [31, Theorem 2.4] in the context of Wiener–Hopf factorizations of symbols of Gram-Laurent matrices).

As in Theorem 3.6, the decay (3.34) of the semi-cardinal symbols’ coefficients admits two proofs: one is constructive, based on the explicit form (3.11) of $\omega _j$ and on Lemma 2.4(ii) (see arXiv:2006.05282v1), the other follows directly from Wiener’s lemma for the algebra of matrices indexed by $H\times H$ with exponential off-diagonal decay (i.e. Jaffard’s result [38, Proposition 2]) and the fact that, by Theorem 3.4, $T_\phi =[\phi (j-k)]_{j,k\in H}$ is invertible in ${\mathcal {B}}(\ell _2 (H))$. If $H={\mathbb {Z}}^{d-1}\times {\mathbb {Z}}_{+}$, the constructive proof also shows that the rate $\beta '$ can be improved to $\beta$ in the last component, as stated.

(ii) The estimate (3.35) follows from the decay properties (2.16), (3.34), and the representation (3.2):

$$\begin{aligned} |\chi _j (x)|\le & {} \sum _{k\in H} |a_{k,j}|\, |\phi (x-k)| \ \le \ C_2 C_0 \sum _{k\in H} e^{-\beta ' |k-j|} e^{-\alpha |x-k|} \\\le & {} C_2 C_0 \sum _{k\in {\mathbb {Z}}^d} e^{-\beta ' |k-j|} e^{-\alpha |x-k|} \ =\ C_2 C_0 \sum _{k\in {\mathbb {Z}}^d} e^{-\beta ' |k|} e^{-\alpha |x-j-k|} \\\le & {} C_3\, e^{-\beta ' |x-j|}, \quad x\in {\mathbb {R}}^d,\ j\in H. \end{aligned}$$

The last inequality is a consequence of the estimate obtained for the cardinal case in the proof of Theorem 2.7, with $C_3:=C_2 C_0 \coth ^d ( \frac{\alpha -\beta '}{2} )$.

It is now straightforward to see how, in case $H={\mathbb {Z}}^{d-1}\times {\mathbb {Z}}_{+}$, the improved decay property for the coefficients $a_{k,j}$ carries over to the Lagrange functions $\chi _j$, after making some necessary changes in the above estimates. $\square$

4 Further Properties

This section deals with two topics: firstly, we show that the Lagrange function for cardinal interpolation can be approximated with shifts of semi-cardinal Lagrange functions; secondly, we obtain the variational characterization of the cardinal and semi-cardinal Lagrange functions and their associated schemes.

4.1 Convergence of Semi-Cardinal Interpolation to Cardinal Interpolation

Let H be a half-plane lattice as in Sect. 3, and assume that $\phi$ satisfies the hypotheses of Theorem 2.1. Since any function of the form (1.6) satisfies the recovery formula (2.10) on the cardinal grid ${\mathbb {Z}}^d$, it follows that, for each $j\in H$, the Lagrange function $\chi _j$ of the corresponding semi-cardinal interpolation scheme on H admits the representation:

$$\begin{aligned} \chi _j ( x ) = \chi (x-j) + \sum _{k\in {\mathbb {Z}}^{d}\setminus H} \chi _j (k) \chi ( x-k ), \quad x\in {\mathbb {R}}^d. \end{aligned}$$

(4.1)

This identity suggests already that the difference $\chi _j - \chi (\cdot - j)$ may be made small as j becomes suitably large. To establish such a property, we employ a different representation, which is a more natural tool for the task.

We start by defining, for both types of half-space lattices, the function

$$\begin{aligned} \eta ( x ):= \sum _{k\in H} \gamma _k\, \phi ( x-k ), \quad x\in {\mathbb {R}}^d, \end{aligned}$$

(4.2)

where $\{\gamma _k\}_{k\in H}$ is the sequence of Fourier coefficients of the Wiener–Hopf factor $\omega _+$ appearing in (3.8)–(3.10). Since $\phi$ is bounded and $\{\gamma _k\}\in \ell _1(H)$, we deduce that (4.2) is an absolutely and uniformly convergent series, hence $\eta$ is continuous and bounded on ${\mathbb {R}}^d$. Moreover, in analogy with the decay properties of $\chi$ established in Sect. 2, we see that not only (2.1) transfers from $\phi$ to $\eta$, but also any additional algebraic or exponential decay that may apply to $\phi$.

Remark

The above function $\eta$ bears a close relationship with the so-called ‘limiting profile’ studied by Goodman et al. [31] in the different context of infinite Gram matrices generated by shifts of multivariate kernels.

Proposition 4.1

Let $\phi$ satisfy the hypotheses of Theorem 2.1. For both types of half-space lattices H, the following absolutely and uniformly convergent representation holds:

$$\begin{aligned} \chi (x) = \sum _{l\in H} \gamma _l\, \eta ( x+l ), \quad x\in {\mathbb {R}}^d. \end{aligned}$$

(4.3)

Further, if $H={\mathbb {Z}}^d_{\preceq ,+}$ and $j\in H$,

$$\begin{aligned} \chi _j (x) = \sum _{0\preceq l\preceq j} \gamma _{j-l}\, \eta ( x-l ), \quad x\in {\mathbb {R}}^d, \end{aligned}$$

(4.4)

while, for $H={\mathbb {Z}}^{d-1}\times {\mathbb {Z}}_+$ and $j=(j',j_d)\in H$, we have

$$\begin{aligned} \chi _j (x) = \sum _{l_d = 0}^{j_d} \sum _{l'\in {\mathbb {Z}}^{d-1}} \gamma _{(j'-l',j_d-l_d)}\, \eta ( x-l ), \quad x\in {\mathbb {R}}^d. \end{aligned}$$

(4.5)

Proof

For any half-space lattice H, we extend $\{\gamma _k\}_{k\in H}$ to a sequence indexed on ${\mathbb {Z}}^{d}$ by setting $\gamma _k:= 0$ for $k\not \in H$. With this convention, (3.9) and (3.10) imply

$$\begin{aligned} a_k = \sum _{l\in H} \gamma _{k+l} \gamma _{l}, \quad k\in {\mathbb {Z}}^{d}. \end{aligned}$$

Substituting $a_k$ in the kernel representation (1.2) of $\chi$ and rearranging,

$$\begin{aligned} \chi (x)= & {} \sum _{k\in {\mathbb {Z}}^{d}} \left( \sum _{l\in H} \gamma _{k+l} \gamma _{l} \right) \phi (x-k) \\= & {} \sum _{l\in H} \gamma _l \sum _{k+l\in H} \gamma _{k+l}\, \phi (x-k) \\= & {} \sum _{l\in H} \gamma _l \sum _{k\in H} \gamma _{k}\, \phi (x+l-k), \end{aligned}$$

which provides (4.3), since the sum of the last series is $\eta ( x+l )$.

Next, for $H={\mathbb {Z}}^d_{\preceq ,+}$, identity (4.4) is obtained by substituting the explicit formula (3.17) for $a_{k,j}$ into the kernel representation (3.2) of $\chi _j$:

$$\begin{aligned} \chi _j (x)= & {} \sum _{k\in H} \left( \sum _{0\preceq l\preceq j} \gamma _{k-l} \gamma _{j-l} \right) \phi (x-k) \\= & {} \sum _{0\preceq l\preceq j} \gamma _{j-l} \sum _{k\in H} \gamma _{k-l}\, \phi (x-k) \\= & {} \sum _{0\preceq l\preceq j} \gamma _{j-l} \sum _{k\in H} \gamma _{k}\, \phi (x-l-k). \end{aligned}$$

Also, (4.5) follows in a similar way if $H={\mathbb {Z}}^{d-1}\times {\mathbb {Z}}_+$, via re-arrangements with respect to the last coordinate of the summation index. $\square$

The next result establishes the uniform convergence of semi-cardinal Lagrange functions $\chi _j(\cdot + j)$ to $\chi$. Note that, if $H={\mathbb {Z}}^{d-1}\times {\mathbb {Z}}_+$ and $j=(j',j_d)\in H$, then, by (3.23), $\chi _{j} (x+j) = \chi _{(0,j_d)} (x',x_d+j_d)$, which is independent of $j'$.

Theorem 4.2

Assume that $\phi$ satisfies the hypotheses of Theorem 2.1.

• (i) If $H={\mathbb {Z}}^d_{\preceq ,+}$, define the nested sequence $\{H_n\}_{n=0}^\infty$ of finite sets

  $$\begin{aligned} H_n:= \{ k\in H: \Vert k\Vert _2 \le n \}, \end{aligned}$$

  and, for each $n\ge 0$, let $j^{(n)}$ denote the maximum element of $H_n$ with respect to the order relation $\preceq$. Then,

  $$\begin{aligned} \lim _{n\rightarrow \infty } \sup _{x\in {\mathbb {R}}^d} \left| \chi ( x ) - \chi _{j^{(n)}} \left( x+j^{(n)} \right) \right| = 0. \end{aligned}$$

  (4.6)

• (ii) If $H={\mathbb {Z}}^{d-1}\times {\mathbb {Z}}_+$, then

  $$\begin{aligned} \lim _{j_d\rightarrow \infty } \sup _{x\in {\mathbb {R}}^d} | \chi ( x ) - \chi _{(0,j_d)} (x',x_d+j_d) | = 0. \end{aligned}$$

  (4.7)

Proof

(i) If $H={\mathbb {Z}}^d_{\preceq ,+}$ and $j\in H$, expression (4.4) implies, on re-indexing and comparing with (4.3):

$$\begin{aligned} \chi (x) - \chi _j (x+j)= & {} \chi (x) - \sum _{0\preceq l\preceq j} \gamma _{l}\, \eta ( x+l ) \nonumber \\= & {} \sum _{l\succ j} \gamma _{l}\, \eta ( x+l ), \quad x\in {\mathbb {R}}^d. \end{aligned}$$

(4.8)

Hence, since $\eta$ is bounded on ${\mathbb {R}}^d$, we obtain the estimate

$$\begin{aligned} | \chi (x) - \chi _j (x+j) | \le M \sum _{l\succ j} | \gamma _{l} |, \quad x\in {\mathbb {R}}^d, \end{aligned}$$

(4.9)

where $M:=\sup _{x\in {\mathbb {R}}^d} |\eta (x)|$. The result follows from this, noting that the definition of $j^{(n)}$ implies $\{ l\in {\mathbb {Z}}^d: l \succ j^{(n)} \} \subset H{\setminus } H_n$, and that

$$\begin{aligned} \lim _{n\rightarrow \infty } \sum _{l\in H\setminus H_n} | \gamma _{l} | = 0, \end{aligned}$$

since $\{\gamma _l\}_{l\in H} \in \ell _1 (H)$ and $\cup H_n = H$.

(ii) If $H={\mathbb {Z}}^{d-1}\times {\mathbb {Z}}_+$, then (4.5) leads to a representation of $\chi - \chi _j(\cdot + j)$ similar to (4.8), in which the summation is taken over all multi-indices $l=(l',l_d)$ with $l'\in {\mathbb {Z}}^{d-1}$ and $l_d>j_d$. Since the absolute summability of $\{\gamma _l\}_{l\in H}$ implies

$$\begin{aligned} \lim _{j_d\rightarrow \infty } \sum _{l_d>j_d} \sum _{l'\in {\mathbb {Z}}^{d-1}} | \gamma _{(l',l_d)} | = 0, \end{aligned}$$

the conclusion is obtained via the estimate analogous to (4.9). $\square$

Remark

• (i) The above sequence $\{H_n\}_{n=0}^\infty$ may be replaced by more general nested sequences of finite sets with union H, as considered in [44, Proposition 3.3], where such sequences are used in connection with the projection (or finite section) method for semi-infinite linear systems.

• (ii) If the kernel $\phi$ satisfies algebraic or exponential decay, upper bounds on the rate of convergence in (4.6) and (4.7) can be obtained using standard tail estimates of the absolutely convergent series $\sum _{k\in H} |\gamma _k|$, based on the corresponding decay of the coefficients $\gamma _k$. Specifically, the uniform error in (4.7) is of magnitude $O ((1+j_d)^{d-\alpha })$ if $\phi$ decays with algebraic power $\alpha >d$, and of magnitude $O(e^{-\alpha 'j_d})$ if $\phi$ decays with exponential rate $\alpha >0$, where $\alpha '\in (0,\alpha )$. Similar rates apply to (4.6), replacing $j_d$ by n.

• (iii) Using Fourier transforms, a result similar to Theorem 4.2 was proved in [6, Theorem 8], for biharmonic splines and $H={\mathbb {Z}}\!\times \!{\mathbb {Z}}_+$. For the integrable kernels $\phi$ treated here, one can also establish, as in [6], the approximation of the cardinal scheme (1.4) by a sequence of semi-cardinal schemes of the type (3.25), suitably defined on shifted lattices $-k+H$, for $k\in H$.

4.2 Native Space Variational Property

The final result of the paper assumes that, in addition to the hypotheses of Theorem 2.1, the kernel $\phi$ is a positive definite function. In this case, it is known that its nonnegative Fourier transform ${\widehat{\phi }}$ belongs to $L_1 ({\mathbb {R}}^d)$; see Wendland [63, Corollary 6.12]. Further, by [63, Theorem 10.12], the translation invariant kernel $\phi (\cdot - \cdot )$ is the reproducing kernel of the ‘native space’ ${\mathcal {N}}_{\phi }$ identified with the real Hilbert space of functions $f\in L_2 ({\mathbb {R}}^d) \cap C ({\mathbb {R}}^d)$ for which there exists, say, ${\mathcal {Q}} f \in L_2 ({\mathbb {R}}^d)$, such that

$$\begin{aligned} {\widehat{f}} = ({\widehat{\phi }})^{1/2} {\mathcal {Q}} f. \end{aligned}$$

(4.10)

The inner product on ${\mathcal {N}}_{\phi }$ is then defined by

$$\begin{aligned} (f,g)_{\phi } = \frac{1}{(2\pi )^d} \int _{\textrm{supp}\,{\widehat{\phi }}} {\mathcal {Q}} f (t)\, \overline{{\mathcal {Q}} g (t)} \, dt. \end{aligned}$$

Theorem 4.3

Let $\phi$ be a positive definite function satisfying the hypotheses of Theorem 2.1.

• (i) The corresponding Lagrange function $\chi$ for cardinal interpolation belongs to the native space ${\mathcal {N}}_{\phi }$ and the following ‘fundamental identity’ holds:

  $$\begin{aligned} (f,\chi )_{\phi } = \sum _{k\in {\mathbb {Z}}^d} a_k f(k), \quad f\in {\mathcal {N}}_{\phi }, \end{aligned}$$

  (4.11)

  where $\{a_k\}_{k\in {\mathbb {Z}}^d}$ are the coefficients of the kernel representation (1.2) of $\chi$. Consequently, $\chi$ is the unique function of minimum-norm in ${\mathcal {N}}_{\phi }$ that satisfies the interpolation conditions (1.3).

• (ii) If H is a half-space lattice in ${\mathbb {R}}^d$, then, for each $j\in H$, the corresponding Lagrange function $\chi _j$ for semi-cardinal interpolation belongs to the native space ${\mathcal {N}}_{\phi }$ and the following ‘fundamental identity’ holds:

  $$\begin{aligned} (f,\chi _j)_{\phi } = \sum _{k\in H} a_{k,j} f(k), \quad f\in {\mathcal {N}}_{\phi }, \end{aligned}$$

  where $\{a_{k,j}\}_{k\in H}$ are the coefficients of the kernel representation (3.2) of $\chi _j$. Hence, $\chi _j$ is the unique function of minimum-norm in ${\mathcal {N}}_{\phi }$ that satisfies the interpolation conditions (3.3).

Proof

(i) The continuity of $\chi$ follows from the uniform convergence on compact sets of its kernel representation (1.2). Further, as shown in [39, section 2], assumption (2.1) implies $\phi \in L_p ({\mathbb {R}}^d)$, for any $p\in [1,\infty ]$. In particular, $\phi \in L_2 ({\mathbb {R}}^d)$ implies that (1.2) is an absolutely convergent series in $L_2 ({\mathbb {R}}^d)$, since $a\in \ell _1$ and $\Vert \chi (\cdot - k) \Vert _{L_2} = \Vert \chi \Vert _{L_2}$, for all $k\in {\mathbb {Z}}^d$. Therefore, the pointwise sum $\chi$ of this series is equal almost everywhere to its $L_2$-sum, hence $\chi \in L_2 ({\mathbb {R}}^d)$.

Similarly, $\phi \in L_1 ({\mathbb {R}}^d)$ implies that (1.2) is an absolutely convergent series in $L_1 ({\mathbb {R}}^d)$, and we may apply the Fourier transform termwise, to obtain

$$\begin{aligned} {\widehat{\chi }} (t) = \sum _{k\in {\mathbb {Z}}^d} a_k [\phi (\cdot -k)]^{\wedge } (t) = {\widehat{\phi }} (t) \sum _{k\in {\mathbb {Z}}^d} a_k e^{-ikt} = {\widehat{\phi }} (t) \, \omega (t), \quad t\in {\mathbb {R}}^d, \end{aligned}$$

(4.12)

where $\omega =1/\sigma _\phi$. Since the symbol $\omega$ is bounded and ${\widehat{\phi }} \in L_1 ({\mathbb {R}}^d)$, we can choose ${\mathcal {Q}}\chi := ({\widehat{\phi }})^{1/2} \omega \in L_2 ({\mathbb {R}}^d)$, hence (4.12) implies $\chi \in {\mathcal {N}}_{\phi }$.

Next, for any $f\in {\mathcal {N}}_{\phi }$, the factorization (4.10) implies ${\widehat{f}} \in L_1 ({\mathbb {R}}^d)$. This fact, together with $f\in C ({\mathbb {R}}^d)$ and the Fourier inversion formula $f=({\widehat{f}}\,)^\vee$ valid for $f\in L_2 ({\mathbb {R}}^d)$ (see Stein and Weiss [59, Theorem 2.4]), imply that the inversion formula holds everywhere as an identity between continuous functions:

$$\begin{aligned} f(x) = \frac{1}{(2\pi )^d} \int _{{\mathbb {R}}^d} {\widehat{f}} (t) \, e^{ixt} \, dt,\quad x\in {\mathbb {R}}^d. \end{aligned}$$

In particular, f is bounded on ${\mathbb {R}}^d$. Since $\{a_k\}_{k\in {\mathbb {Z}}^d}\in \ell _1$, the series (4.11) is absolutely convergent and we have

$$\begin{aligned} (2\pi )^d \, (f,\chi )_{\phi }= & {} \int _{\textrm{supp}\,{\widehat{\phi }}} {\mathcal {Q}}f\! (t)\, \overline{{\mathcal {Q}} \chi (t)} \, dt \ = \ \int _{\textrm{supp}\,{\widehat{\phi }}} {\mathcal {Q}}f\! (t)\, ({\widehat{\phi }})^{1/2}\! (t)\, \overline{\omega (t)} \, dt \\= & {} \int _{{\mathbb {R}}^d} {\widehat{f}} (t) \, \overline{\omega (t)} \, dt \ = \ \sum _{k\in {\mathbb {Z}}^d} a_k \int _{{\mathbb {R}}^d} {\widehat{f}}(t)\, e^{ikt} \,dt \\= & {} (2\pi )^d \sum _{k\in {\mathbb {Z}}^d} a_k f(k), \end{aligned}$$

by permitted interchange of the integration and summation operators in the fourth equality. The identity (4.11) implies that $\chi$ is orthogonal on any $f\in {\mathcal {N}}_{\phi }$, such that $f(k) = 0$ for all $k\in {\mathbb {Z}}^d$. A standard Pythagorean-type argument then establishes the last assertion of part (i).

(ii) This follows as in the cardinal case (i), by replacing (4.12) with:

$$\begin{aligned} \widehat{\chi _j} (t) = {\widehat{\phi }} (t)\, \omega _j (-t),\quad t\in {\mathbb {R}}^d, \end{aligned}$$

which is the Fourier transform of $\chi _j$. $\square$

Remark

A similar proof can be used to show that, for a data sequence $y\in \ell _1$, the cardinal interpolant $s=I_\phi y$ given by (2.2) belongs to the native space ${\mathcal {N}}_{\phi }$, it satisfies a fundamental identity similar to (4.11) (with coefficients $\{c_k\}$ replacing $\{a_k\}$), and it is the unique norm-minimizer among all functions in ${\mathcal {N}}_{\phi }$ that satisfy the interpolation conditions (2.4). The result also holds for the semi-cardinal interpolant $s=I^H_\phi y$ of (1.6) in place of $\chi _j$, if the data sequence $y\in \ell _1 (H)$.

5 Examples

In this section, we specialize our main results to five classes of integrable kernels that have been studied in the cardinal interpolation literature.

Example 5.1

(Gaussian) For a parameter $c>0$, the Gaussian kernel

$$\begin{aligned} \phi (x) = \exp (-c\Vert x\Vert ^2), \quad x\in {\mathbb {R}}^d, \end{aligned}$$

obviously verifies the exponential decay condition (2.16). Since

$$\begin{aligned} {\widehat{\phi }} (t) = \frac{1}{(4\pi c)^{d/2}} \exp \left( -\frac{\Vert t\Vert ^2}{4c} \right) > 0, \quad t\in {\mathbb {R}}^d, \end{aligned}$$

$\phi$ is positive definite and an application of the Poisson Summation Formula shows that the corresponding cardinal symbol satisfies

$$\begin{aligned} \sigma _\phi (e^{it}) = \sum _{l\in {\mathbb {Z}}^d} {\widehat{\phi }} (t+2\pi l) > 0. \end{aligned}$$

Cardinal interpolation on ${\mathbb {Z}}^d$ with the Gaussian kernel has been thoroughly studied by Riemenschneider and Sivakumar [49, 50]. However, semi-cardinal interpolation has so far been considered only in the univariate case by Baxter and Sivakumar [2, Theorem 2.9], who proved that, for $\tau \in {\mathbb {R}}$, the shifted Gaussian kernel $\phi (\cdot +\tau )$ generates an invertible semi-infinite Toeplitz matrix on $\ell _2 ({\mathbb {Z}}_+)$ if and only if $|\tau |<1/2$.

Our new results in this paper establish the main properties of semi-cardinal interpolation with the Gaussian kernel on any half-space lattice $H\subset {\mathbb {Z}}^d$ of one of the two types considered in Sect. 3.1, including the exponential decay of the associated Lagrange functions and of their kernel expansion coefficients (Theorem 3.8).

Example 5.2

(Matérn) The Matérn kernel is given by

$$\begin{aligned} \phi (x) = \Vert x\Vert ^{m-\frac{d}{2}} K_{m-\frac{d}{2}} (\Vert x\Vert ), \quad x\in {\mathbb {R}}^d, \end{aligned}$$

where $m\in {\mathbb {R}}$, $m>\frac{d}{2}$, and $K_{m-\frac{d}{2}}$ is a modified Bessel function. Since, for $\nu >0$, $K_{\nu } (r)$ is continuous on $(0,\infty )$ and $r^{\nu } K_{\nu } (r)$ can be extended by continuity at $r=0$, it follows that $\phi$ is continuous on ${\mathbb {R}}^d$. Also, the asymptotic behavior of $K_{\nu }$ for large argument [63, Lemma 5.13] implies $\phi (x) = O ( \Vert x\Vert ^{m-\frac{d+1}{2}} e^{ -\Vert x\Vert } )$, as $\Vert x\Vert \rightarrow \infty$, so the exponential decay condition (2.16) holds, for some $\alpha >0$. By [63, Theorem  6.13], ${\widehat{\phi }} (t) = \rho _{m,d} (1+\Vert t\Vert ^2)^{-m}$, $t\in {\mathbb {R}}^d$, where $\rho _{m,d} > 0$, hence $\phi$ is positive definite and its associated native space is the Bessel potential space

$$\begin{aligned} H_2^m ({\mathbb {R}}^d) = \{ f\in L_2 ({\mathbb {R}}^d): {\widehat{f}} = (1+\Vert \!\cdot \!\Vert ^2)^{-m/2}\,{\widehat{g}},\ g\in L_2 ({\mathbb {R}}^d) \}. \end{aligned}$$

This coincides with the usual Sobolev space when m is a positive integer; in this case, Matérn kernels are also known under the name ‘Sobolev splines’. Also, $\phi$ is the fundamental solution of the (pseudo)differential operator $(1-\Delta )^m$, where $\Delta$ is the Laplace operator in ${\mathbb {R}}^d$.

Since the condition $2m>d$ permits the application of Poisson’s Summation Formula, as in Example 5.1, the associated cardinal symbol $\sigma _\phi$ is seen to satisfy the positivity condition (2.6). Therefore Theorems 2.7 and 3.8 provide the exponential decay of the Lagrange functions for cardinal and semi-cardinal interpolation with any Matérn kernel $\phi$ satisfying $2m>d$. These results also extend to the class of ‘generalized Whittle-Matérn kernels’ introduced and studied by Bozzini et al. [15].

For integer values of m, the follow-up paper [10] by the author considers the scaled Matérn kernel $\phi _h:= \phi (h\cdot )$, proving the stronger result that the exponential decay of the Lagrange function for cardinal interpolation with $\phi _h$ on ${\mathbb {Z}}^d$ holds with constants independent of h, as $h\rightarrow 0$. As an application, the optimal $L_\infty$-convergence rate $O(h^{2m})$ for nonstationary ‘flat ladder’ interpolation with $\phi$ on the scaled grid $h{\mathbb {Z}}^d$ is established.

Note that, if $d=1$ and m is a positive integer, the Matérn kernel reduces to a piecewise exponential of class $C^{2m-2}({\mathbb {R}})$:

$$\begin{aligned} \phi (x) = |x|^{m-\frac{1}{2}} K_{m-\frac{1}{2}} (|x|) = e^{ -|x| } \sum _{r=0}^{m-1} c_{m,r} |x|^r, \end{aligned}$$

(5.1)

since the modified Bessel function of half-integer order is expressible in finite terms [62, Eq. 3$\cdot$71(12)], for some positive constants $\{c_{m,r}\}$. The main properties of cardinal interpolation with this exponential spline kernel on the real line follow from Micchelli’s general ‘cardinal L-spline’ theory [45]. Micchelli’s comprehensive treatment includes, among other results, the exponential decay of the Lagrange function $\chi$ for cardinal interpolation on ${\mathbb {Z}}$, obtained by different methods. Also, in [12] and [8], cardinal and semi-cardinal interpolation schemes for a scaled version of the univariate kernel (5.1) have played a prominent role in the construction, via separation of variables, of multivariate ‘polysplines’ interpolating continuous functions prescribed on equi-spaced parallel hyperplanes. These two references, as well as [10], employ Fourier transforms, rather than kernel representations, in order to obtain the exponential decay of the corresponding Lagrange functions.

For any dimension d, the choice $m=\frac{d+1}{2}$ provides the particularly simple kernel $\phi (x) = \Vert x\Vert ^{\frac{1}{2}} K_{\frac{1}{2}} (\Vert x\Vert ) = \sqrt{\frac{\pi }{2}} \, e^{ -\Vert x\Vert }$.

Example 5.3

(Generalized inverse multiquadric) This kernel is given by

$$\begin{aligned} \phi (x) = (c^2+\Vert x\Vert ^2)^{-m}, \quad x\in {\mathbb {R}}^d, \end{aligned}$$

where $2\,m > d$, $c>0$, so the algebraic decay condition (2.13) holds with the rate 2m. This example essentially exchanges $\phi$ with ${\widehat{\phi }}$ from Example 5.2, and by usual transform laws we have:

$$\begin{aligned} {\widehat{\phi }} (t) = \frac{1}{\rho _{m,d}} \left( \frac{\Vert t\Vert }{c} \right) ^{m-\frac{d}{2}} K_{m-\frac{d}{2}} (\Vert c\,t\Vert ), \quad t\in {\mathbb {R}}^d. \end{aligned}$$

Since $r^{\nu } K_{\nu } (r)$ takes only positive values for $\nu >0$ and $r\in [0,\infty )$, it follows that $\phi$ is positive definite and we may again conclude, via Poisson’s Summation Formula, that the associated cardinal symbol satisfies the positivity condition (2.6). Consequently, for cardinal and semi-cardinal interpolation with the above kernel $\phi$ satisfying $2m>d$, Theorems 2.6 and 3.6 provide the transfer of the algebraic decay rate 2m to the corresponding Lagrange functions and their kernel expansion coefficients.

In the case of cardinal interpolation, this result represents a substantial improvement over decay rates recently obtained by Hamm and Ledford in [36, 37], who derived their results from estimates based on the smoothness of Fourier transforms. Specifically, assuming $d=1$ and $m>1$ (with our notation), [36, Corollary 2] states the algebraic rate $\lceil 2m-2 \rceil$ for the decay of the Lagrange function $\chi$ for cardinal interpolation. Also, under the assumption $2m>2d+1$ (with our notation), [37, Theorem 3.1] proves the algebraic rate $\lfloor 2m-d \rfloor$ if $m\not \in {\mathbb {Z}}$, and $2\,m-d-1$ if $m\in {\mathbb {Z}}$, for the kernel expansion coefficients of $\chi$, while [37, Corollary 4.7 & Eq. (10)] imply the algebraic decay rate $d+1$ for $\chi$.

Example 5.4

(B-splines and box-splines) The study of univariate cardinal interpolation with the central B-spline $\phi =M_n$ of polynomial degree $n-1$, compactly supported on $[-n/2,n/2]$, and of class $C^{n-2}({\mathbb {R}})$, was initiated in Schoenberg’s 1946 paper [54]. As pointed in [54], the explicit piecewise polynomial form of $M_n$, as well as its Fourier transform

$$\begin{aligned} \widehat{M_n} (t) = \left( \frac{\sin (t/2)}{t/2}\right) ^n,\quad t\in {\mathbb {R}}, \end{aligned}$$

can be traced back to Laplace’s work on probability. Clearly, $M_n$ is positive definite iff n is an even integer. Schoenberg proved that the cardinal symbol $\sigma _\phi$ (expressed in terms of the so-called Euler–Frobenius polynomial) associated to $\phi =M_n$ satisfies the positivity condition (2.6), irrespective of the parity of n. This enabled his construction (see [55, §4.5]) of cardinal interpolation for the B-spline kernel $M_n$, with an exponentially decaying Lagrange function $\chi$.

Subsequently, in [55, 56], Schoenberg employed two methods to construct semi-cardinal schemes for interpolation on ${\mathbb {Z}}_+$, with splines of odd degree only. The existence and uniqueness properties of these schemes required the introduction of boundary conditions imposed at the origin. However, Schoenberg obtained the exponential decay of the corresponding semi-cardinal Lagrange functions without relying on B-spline representations in terms of the kernel $\phi =M_n$. For cubic and quintic splines, the construction of semi-cardinal interpolation with boundary conditions expressed in terms of finite differences (FD) of B-spline coefficients was given by Bejancu [7, 9] and Bejancu et al. [11].

The latter approach is briefly illustrated here in the cubic spline case $n=4$, to motivate the bivariate extension described further below. Since the support of $M_4$ is $[-2,2]$, any cubic spline defined on $[0,\infty )$, with knots at the nonnegative integers, admits a B-spline representation

$$\begin{aligned} s(x) = \sum _{k=-1}^\infty c_k M_4 (x-k), \quad x\in [0,\infty ), \end{aligned}$$

(5.2)

for some coefficients $c_k$, $k\ge -1$. After imposing the interpolation conditions $s(j) = y_j$, $j\in {\mathbb {Z}}_+$, there is still one ‘degree of freedom’ to be utilized by a boundary condition for s. The ‘natural’ condition $s'' (0) = 0$ used in [7] is equivalent to the second order FD condition

$$\begin{aligned} c_{-1} - 2c_0 + c_1 = 0, \end{aligned}$$

(5.3)

while the ‘not-a-knot’ condition $s'''(1-) = s'''(1+)$, which eliminates the spline knot at 1, can be expressed equivalently (see [53]) as:

$$\begin{aligned} c_{-1} - 4c_0 + 6c_1 - 4c_2 + c_3 = 0. \end{aligned}$$

(5.4)

In addition, for $\phi = M_4$ and $H = {\mathbb {Z}}_+$, representation (1.6) of the present paper is obtained by imposing the boundary condition $c_{-1} = 0$ in (5.2).

The extension of cardinal spline interpolation to several variables, in a non-tensor product setting, was initiated by de Boor, Höllig, and Riemenschneider [21], for certain bivariate ‘box spline’ kernels, which are compactly supported and piecewise polynomial over a three-direction partition of the plane. The ensuing theory of cardinal interpolation with box splines in arbitrary dimension is much more intricate than its univariate version; for a thorough introduction to the main results, see the monograph [22, Chapter IV]. Since the exponential decay condition (2.16) holds trivially if $\phi$ is a box spline in any number of variables, our Theorem 3.8 establishes the exponential decay of Lagrange functions for semi-cardinal interpolation whenever the box spline kernel is symmetric and has a positive symbol (2.6).

In [7], the author presented a semi-cardinal scheme for interpolation at the points of the half-plane lattice $H={\mathbb {Z}}\!\times \!{\mathbb {Z}}_+$ with the three-direction box-spline $\phi = M_{2,2,2}\in C^2 ({\mathbb {R}}^d)$, whose direction matrix has every multiplicity 2. This box-spline is an example of a symmetric non-radial positive definite kernel which may be regarded as a bivariate analog of $M_4$, and its only non-zero values on ${\mathbb {Z}}^2$ are:

$$\begin{aligned} \phi (0,0) = \frac{1}{2}, \ \ \phi (\pm (1,1)) = \phi (0,\pm 1) = \phi (\pm 1,0) = \frac{1}{12}. \end{aligned}$$

It follows that the associated cardinal symbol satisfies (2.6), since

$$\begin{aligned} \sigma _\phi (e^{it_1},e^{it_2}) = \frac{1}{6} (3 + \cos t_1 + \cos t_2 + \cos (t_1+t_2) ) \ge \frac{1}{4}, \quad t\in {\mathbb {R}}^2. \end{aligned}$$

The approach of [7] is based on the fact that, due to the size of the support of $M_{2,2,2}$, the restriction of a bivariate cardinal series $s = \sum _{k\in {\mathbb {Z}}^2} c_k M_{2,2,2} (\cdot -k)$ to the upper half-plane ${\mathbb {R}}\!\times \![0,\infty )$ has a representation analogous to (5.2):

$$\begin{aligned} s(x) = \sum _{k_2=-1}^\infty \sum _{k_1=-\infty }^\infty c_{k_1,k_2} M_{2,2,2} (x_1-k_1,x_2-k_2), \quad x\in {\mathbb {R}}\!\times \![0,\infty ). \end{aligned}$$

(5.5)

This permitted the formulation of a system of FD edge conditions in terms of the box-spline coefficients $c_{k_1,k_2}$, with the weights stencil (A) shown below, as the analog of the ‘natural’ univariate stencil (5.3). Further, Sabin and Bejancu [13, 53] introduced the stencil (B), which extends the ‘not-a-knot’ univariate stencil (5.4).

$$\begin{aligned} (A) \quad \begin{array}{rr} &{} 1 \\ -1 &{} -1 \\ 1 &{} \end{array} \qquad \qquad \qquad (B) \quad \begin{array}{rrr} &{} &{} 1 \\ &{} -2 &{} -2 \\ 1 &{} 4 &{} 1 \\ -2 &{} -2 &{} \\ 1 &{} &{} \end{array} \end{aligned}$$

(Note that the bottom weights of both stencils are applied to the coefficients $c_{k_1,k_2}$ corresponding to $k_2 = -1$.) Moreover, for $\phi = M_{2,2,2}$ and $H={\mathbb {Z}}\!\times \!{\mathbb {Z}}_+$, representation (1.6) of the present paper corresponds to simply imposing zero-order FD conditions in (5.5):

$$\begin{aligned} c_{k_1,k_2}=0,\quad (k_1,k_2)\in {\mathbb {Z}}\!\times \!\{-1\}. \end{aligned}$$

(5.6)

Using explicit expressions for box spline coefficients, [7, 13] established not only the exponential decay of the bivariate Lagrange functions, but also the polynomial reproduction properties of the semi-cardinal box-spline schemes, implying that the ‘stationary’ scaled scheme with box-spline edge conditions based on (A) achieves exactly half of the maximal approximation order 4 known to hold for scaled cardinal interpolation with $M_{2,2,2}$, while the edge conditions of type (B) have the effect of restoring the full approximation order.

On the other hand, as in the univariate case [9], it can be proved that the bivariate stationary scaled scheme based on representation (1.6) of this paper, i.e. on the boundary conditions (5.6), possesses approximation order zero for $\phi = M_{2,2,2}$. However, the advantage of this scheme over those in [7, 13] rests on the fact that it extends to arbitrary box-spline kernels that generate ‘correct’ cardinal interpolation, for any dimension and any half-space lattice H (Theorem 3.8). A similar extension for semi-cardinal schemes with FD boundary conditions of positive order poses a challenging problem in arbitrary dimension.

Example 5.5

(Polyharmonic B-splines) This example is related to cardinal interpolation with kernels $\phi (x)$ that actually grow, as $\Vert x\Vert \rightarrow \infty$, in which case the symbol $\sigma$ cannot be defined classically and the associated Laurent operator $L_\phi$ is not bounded anymore on $\ell _2$. Nevertheless, a distributional Fourier transform approach still allows the construction of a suitably decaying Lagrange function $\chi$ with an absolutely convergent kernel representation (1.2), as demonstrated, for polyharmonic kernels, by Madych and Nelson [43], and, for other non-decaying kernels including the multiquadrics, by Buhmann [17].

In particular, for each integer $m>d/2$, [43] considered cardinal interpolation with the m-harmonic kernel

$$\begin{aligned} \phi (x) = c_{m,d} \left\{ \begin{array}{ll} \Vert x\Vert ^{2m-d} \ln \Vert x\Vert , &{} \textrm{if}\ d\ \mathrm {is\ even}, \\ \Vert x\Vert ^{2m-d}, &{} \textrm{if}\ d\ \mathrm {is\ odd}, \end{array} \right. \quad x\in {\mathbb {R}}^d. \end{aligned}$$

For a suitable constant $c_{m,d}$, this satisfies $\Delta ^m \phi = \delta$, where $\Delta$ is the Laplace operator in ${\mathbb {R}}^d$ and $\delta$, the Dirac distribution. Based on the analyticity properties of the Fourier transform ${\widehat{\chi }}$, Madych and Nelson proved that the corresponding Lagrange function $\chi$ and its coefficients of representation (1.2) decay exponentially.

For each $i=1,\ldots ,d$, let $\nabla _i$ denote the second-order central difference operator of unit step in the i-th variable. In analogy with odd-degree polynomial B-splines, Rabut [46, 47] considered the ‘elementary m-harmonic cardinal B-spline’ defined by

$$\begin{aligned} \psi := \left( \sum _{i=1}^d \nabla _i \right) ^m \phi . \end{aligned}$$

In the thin plate spline case $m=d=2$, this bell-shaped kernel had been used before by Dyn and Levin [24, 25], for surface fitting in a bounded domain. Although $\psi$ does not have a compact support, it decays algebraically, with $\psi (x) = O (\Vert x\Vert ^{-(d+2)})$, for large $\Vert x\Vert$. Rabut showed that the cardinal interpolation Lagrange function $\chi$ generated by $\phi$ also admits a series representation of the type (1.2) with $\phi$ replaced by $\psi$ and with coefficients decaying at the same algebraic rate as $\psi$.

In [46, 48], Rabut replaced $\sum _{i=1}^d \nabla _i$ in the above definition of $\psi$ with a higher order discretization of the Laplace operator. This led him to ‘high level’ m-harmonic cardinal B-splines with faster decay, reaching $O (\Vert x\Vert ^{-(d+2\,m)})$, for large $\Vert x\Vert$. Also, Van De Ville et al. [61] introduced a related class of polyharmonic B-splines in two and three variables, using a more isotropic discretization of the Laplacean. Compared to Rabut’s ‘elementary’ B-splines, these possess some improved features, including a $O (\Vert x\Vert ^{-(d+4)})$ decay. Both the ‘high level’ and the ‘isotropic’ polyharmonic B-splines provide corresponding kernel representations for Madych and Nelson’s cardinal Lagrange function $\chi$.

For each choice of a polyharmonic B-spline kernel $\psi$ described above, our Theorem 3.6 establishes the algebraic decay of the Lagrange functions for semi-cardinal interpolation on a half-space lattice, with the same rate as that of $\psi$. Note that, for $H={\mathbb {Z}}^{d-1}\!\times \!{\mathbb {Z}}_+$, the semi-cardinal scheme generated here by such a kernel $\psi$ does not coincide with that constructed by Bejancu [5] via the distributional Fourier transform of $\phi$, due to different Wiener–Hopf factorizations of the corresponding inverse cardinal symbols.

The idea of using a pseudo B-spline kernel $\psi$ generated from differences of polyharmonic or of other non-decaying kernels, such as the multiquadrics, is a recurrent theme in the radial basis function literature, e.g., [3, 4, 14, 20, 24, 25]. In this context, it may be of interest to ask whether, in analogy to Example 5.4, the formulation of FD boundary conditions in terms of kernel coefficients can lead to improved practical schemes for interpolation in bounded domains.