# Random sampling and polynomial-free interpolation by Generalized
MultiQuadrics

A. Sommariva  University of Padova, Italy  M. Vianello  University of Padova,
Italy

###### Abstract

We prove that interpolation matrices for Generalized MultiQuadrics (GMQ) of
order greater than one are almost surely nonsingular without polynomial
addition, in any dimension and with any continuous random distribution of
sampling points. We also include a new class of generalized MultiQuadrics
recently proposed by Buhmann and Ortmann.

2020 AMS subject classification: 65D05,65D12. Keywords: multivariate
interpolation, random sampling, Radial Basis Functions, Generalized
MultiQuadrics, polynomial-free unisolvence.

##  1 Introduction

Looking at the relevant literature on RBF interpolation, the only results
ensuring unisolvence of CPD (Conditionally Positive Definite) RBF without
polynomial addition essentially date back to almost 40 years ago. These
concern distance functions and CPD RBF of order 1 like MultiQuadrics [9], and
univariate cubic splines [1] which are of order 3. Such a theoretical problem
has been clearly stated in the popular monograph [5]: “There is no result that
states that interpolation with Thin-Plate Splines (or any other strictly
conditionally positive definite function of order $m\geq 2$) without the
addition of an appropriate degree $m-1$ polynomial is well-posed”. With the
noteworthy exception of [1] which has been in some sense overlooked (probably
because the result is purely univariate), the problem remains substantailly
open in the deterministic sampling framework. On the other hand, the fact that
for example interpolation by Thin-Plate Splines without polynomial addition is
practically always unisolvent, has been substantially recognized in the
application literature, cf. e.g. [12] with the references therein.

Recently indeed, some advancement has been made in the random sampling
framework. In particular, it has been proved that polynomial-free
interpolation by Polyharmonic Splines of any order, in any dimension and with
any continuous distribution of random sampling points is a.s. (almost surely)
unisolvent. The proof is essentially based on a careful analysis of the
determinant of the interpolation matrix, together with the fact that
Polyharmonic Splines are real analytic functions up to their centers, where
they present a singularity of some derivative (with additional difficulties in
the case of Radial Powers of odd integer exponent); cf. [2, 13]. The key
underlying property is a famous basic result of measure theory, namely that
the zero set of a not identically zero real analytic function on an open
connected domain has null Lebesgue measure; cf. e.g. [10] with the references
therein. It is worth quoting that a similar probabilistic approach was also
recently used in the framework of other (almost everywhere) analytic functions
spaces [4, 14], for example polynomial spaces or spaces of RBF with fixed
center, while the results quoted above concern the more difficult case of RBF
whose centers coincide with the sampling points.

In the present note we extend the polynomial-free unisolvence result by random
sampling to Generalized MultiQuadrics (GMQ), that correspond to the radial
functions

$\phi(r)=(1+(\varepsilon
r)^{2k})^{\beta}\;,\;\;1<\beta\not\in\mathbb{N}\;,\;k\in\mathbb{N}\;.$ (1)

The scale $\varepsilon>0$ represents the so-called _shape parameter_
associated with RBF and is important to control the interpolation problem
conditioning, cf. e.g. [5, 8]. For $k=1$, it is well-known that with $\beta<0$
(Generalized Inverse Multiquadrics) they are Positive Definite and with
$\beta>0$ they are CPD of order $\lceil\beta\rceil$ (for $k=1$ and $\beta=1/2$
they are the classical Hardy’s MultiQuadrics); cf. [3, 5, 6]. For $k>1$ these
radial functions are a generalization of a new family recently proposed by
Buhmann and Ortmann in the case $\beta=1/2$, which are shown to be not
positive definite and where the parameter $k$ allows enhanced polynomial
reproduction in the framework of quasi-interpolation; cf. [11]. To our
knowledge, such a family with its generalization appears here for the first
time in the framework of scattered interpolation.

The main difficulty to prove polynomial-free unisolvence with respect to the
case of TPS, is that $\phi$ is here an analytic function with no real
singularities, so in the proof we have to embed the problem in the complex
framework to exploit the presence of complex singularities (namely, branch
points).

##  2 Polynomial-free random GMQ interpolation

Our main result is the following:

###### Theorem 1

Let $\Omega$ be an open connected subset of $\mathbb{R}^{d}$, $d\geq 1$, and
$\\{\mathbf{x}_{i}\\}_{i\geq 1}$ be a randomly distributed sequence on
$\Omega$ with respect to any given probability density $\sigma(\mathbf{x})$,
i.e. a point sequence produced by sampling a sequence of absolutely continuous
random variables $\\{X_{i}\\}_{i\geq 1}$ which are i.i.d. (independent and
identically distributed) in $\Omega$ with density $\sigma\in
L^{1}_{+}(\Omega)$. Moreover, let

$V_{n}=[\phi_{j}(\mathbf{x}_{i})]\;,\;\;\phi_{j}(\mathbf{x})=\phi(\|\mathbf{x}-\mathbf{x}_{j}\|_{2})\;,\;1\leq
i,j\leq n\;,$ (2)

be the interpolation matrix with respect to Generalized MultiQuadrics (1) of
order $\lceil\beta\rceil>1$. Then, for every $n\geq 1$ the points
$\mathbf{x}_{1},\dots,\mathbf{x}_{n}$ are a.s.(almost surely) distinct and the
matrix $V_{n}$ is a.s. nonsingular.

Proof. The proof proceeds by complete induction on $n$. The induction base
with $n=1$ is trivial, since $det(V_{1})=\phi_{1}(\mathbf{x}_{1})=\phi(0)=1$.
Let us assume as inductive hypothesis that $det(V_{k})$ is a.s. nonzero for
$1\leq k\leq n$. First, $\mathbf{x}_{n+1}$ is a.s. distinct from
$\mathbf{x}_{1},\dots,\mathbf{x}_{n}$, since the probability that a random
point falls in a finite set is clearly zero, as finite set having zero measure
with respect to any distribution with density.

Consider now the $(n+1)\times(n+1)$ matrix

$A(\mathbf{x})=\left(\begin{array}[]{cccccc}1&\phi_{2}(\mathbf{x}_{1})&\cdots&\phi_{n-1}(\mathbf{x}_{1})&\phi_{n}(\mathbf{x}_{1})&\phi_{1}(\mathbf{x})\\\
\\\
\phi_{1}(\mathbf{x}_{2})&1&\cdots&\phi_{n-1}(\mathbf{x}_{2})&\phi_{n}(\mathbf{x}_{2})&\phi_{2}(\mathbf{x})\\\
\\\ \vdots&\vdots&\ddots&\vdots&\vdots&\vdots\\\ \\\
\phi_{1}(\mathbf{x}_{n-1})&\phi_{2}(\mathbf{x}_{n-1})&\cdots&1&\phi_{n}(\mathbf{x}_{n-1})&\phi_{n-1}(\mathbf{x})\\\
\\\
\phi_{1}(\mathbf{x}_{n})&\phi_{2}(\mathbf{x}_{n})&\cdots&\phi_{n-1}(\mathbf{x}_{n})&1&\phi_{n}(\mathbf{x})\\\
\\\
\phi_{1}(\mathbf{x})&\phi_{2}(\mathbf{x})&\cdots&\phi_{n-1}(\mathbf{x})&\phi_{n}(\mathbf{x})&1\\\
\\\ \end{array}\right)$

and observe that $A(\mathbf{x}_{n+1})=V_{n+1}$, since
$\phi_{j}(\mathbf{x}_{n+1})=\phi_{n+1}(\mathbf{x}_{j})$. Developing the
determinant of $A(\mathbf{x})$ by Laplace rule along the last row, we get

$f(\mathbf{x})=det(A(\mathbf{x}))=-det(V_{n-1})(\phi_{n}(\mathbf{x}))^{2}+a(\mathbf{x})\phi_{n}(\mathbf{x})+b(\mathbf{x})\;,$
(3)

where $a,b\in span\\{\phi_{1},\dots,\phi_{n-1}\\}$. Notice that $f$ is an
analytic function in $\mathbb{R}^{d}$, since it is defined by sum and products
of functions in $span\\{\phi_{1},\dots,\phi_{n}\\}$ and thus belongs to the
function algebra generated by $\phi_{1},\dots,\phi_{n}$; cf. [7]. We claim
that $f$ is a.s. not identically zero in $\Omega$.

In fact, if $f$ were identically zero in $\Omega$, it would be identically
zero also in $\mathbb{R}^{d}$, by the already quoted fundamental result that
the zero set of a not identically zero real analytic function has null
Lebesgue measure [10], whereas $meas(\Omega)>0$. Then taking the line
$\mathbf{x}(t)=\mathbf{x}_{n}+t\mathbf{u}$ where
$\mathbf{u}=(u_{1},\dots,u_{d})$ is a given unit vector, we obtain that the
real univariate function $f(\mathbf{x}(t))$ would be identically zero for
$t\in\mathbb{R}$. Consequently, its analytic extension to the complex plane,
say $f(\mathbf{x}(z))$, would also be identically zero for $z\in\mathbb{C}$.
Observe that
$\|\mathbf{x}(z)-\mathbf{x}_{j}\|_{2}^{2}=\|\mathbf{x}_{n}+z\mathbf{u}-\mathbf{x}_{j}\|_{2}^{2}$
appearing in
$\phi_{j}(\mathbf{x}(z))=\phi(\|\mathbf{x}(z)-\mathbf{x}_{j}\|_{2})$, $1\leq
j\leq n$, has to be seen not as the square of the complex 2-norm, but as the
complex extension of the corresponding real function

$\|\mathbf{x}_{n}+t\mathbf{u}-\mathbf{x}_{j}\|_{2}^{2}=\sum_{\ell=1}^{d}{((\mathbf{x}_{n}-\mathbf{x}_{j})_{\ell}+tu_{\ell})^{2}}$
$=\sum_{\ell=1}^{d}{[(\mathbf{x}_{n}-\mathbf{x}_{j})_{\ell}^{2}+2tu_{\ell}(\mathbf{x}_{n}-\mathbf{x}_{j})_{\ell}+t^{2}u_{\ell}^{2}]}$
$=t^{2}+2t\langle\mathbf{u},\mathbf{x}_{n}-\mathbf{x}_{j}\rangle+\|\mathbf{x}_{n}-\mathbf{x}_{j}\|_{2}^{2}\;,$

that is

$\|\mathbf{x}(z)-\mathbf{x}_{j}\|_{2}^{2}=z^{2}+2z\langle\mathbf{u},\mathbf{x}_{n}-\mathbf{x}_{j}\rangle+\|\mathbf{x}_{n}-\mathbf{x}_{j}\|_{2}^{2}\;,$

where $\langle\cdot,\cdot\rangle$ denotes the Euclidean scalar product in
$\mathbb{R}^{d}$. Consequently

$\phi_{j}(\mathbf{x}(z))=\left(1+\varepsilon^{2k}\left(z^{2}+2z\langle\mathbf{u},\mathbf{x}_{n}-\mathbf{x}_{j}\rangle+\|\mathbf{x}_{n}-\mathbf{x}_{j}\|_{2}^{2}\right)^{k}\right)^{\beta}$
$=\exp\left(\beta\log\left(1+\varepsilon^{2k}\left(z^{2}+2z\langle\mathbf{u},\mathbf{x}_{n}-\mathbf{x}_{j}\rangle+\|\mathbf{x}_{n}-\mathbf{x}_{j}\|_{2}^{2}\right)^{k}\right)\right)\;,$

where $\log$ is the principal value of the complex logarithm.

In particular

$\phi_{n}(\mathbf{x}(z))=\left(1+(\varepsilon z)^{2k}\right)^{\beta}$

has a branch point at
$z_{\ast}=z_{\ast}(\varepsilon)=\varepsilon^{-1}e^{i\pi/(2k)}$, which is a
root of the complex polynomial ($\varepsilon z)^{2k}+1$. On the other hand,
the functions $\phi_{j}(\mathbf{x}(z))$, $j<n$, are a.s. analytic at
$z_{\ast}$, at least for a suitable choice of the unit vector $\mathbf{u}$. To
prove this fact, we have to distinguish the cases $k=1$ and $k>1$, and in the
latter the cases $d\geq 2$ and $d=1$. It is sufficient to check that the
complex numbers

$1+\varepsilon^{2k}\left(z_{\ast}^{2}+2z_{\ast}\langle\mathbf{u},\mathbf{x}_{n}-\mathbf{x}_{j}\rangle+\|\mathbf{x}_{n}-\mathbf{x}_{j}\|_{2}^{2}\right)^{k}\;,\;\;1\leq
j<n\;,$

have positive real part, or even nonzero imaginary part, so that they cannot
fall on the branch point or on the branch cut of the noninteger power
$(\cdot)^{\beta}$.

For $k=1$, it is sufficient to observe that for any $\mathbf{u}$ the complex
numbers

$1+\varepsilon^{2}\left((i/\varepsilon)^{2}+2i/\varepsilon\langle\mathbf{u},\mathbf{x}_{n}-\mathbf{x}_{j}\rangle+\|\mathbf{x}_{n}-\mathbf{x}_{j}\|_{2}^{2}\right)=\varepsilon^{2}\|\mathbf{x}_{n}-\mathbf{x}_{j}\|_{2}^{2}+2i\varepsilon\langle\mathbf{u},\mathbf{x}_{n}-\mathbf{x}_{j}\rangle$

have a.s. positive real part since a.s. $\mathbf{x}_{j}\neq\mathbf{x}_{n}$ for
$j<n$.

For $k>1$, consider the complex-valued polynomial in $\mathbf{u}$ on the unit
sphere

$p_{j}(\mathbf{u})=\varepsilon^{2k}\left(z_{\ast}^{2}+2z_{\ast}\langle\mathbf{u},\mathbf{x}_{n}-\mathbf{x}_{j}\rangle+\|\mathbf{x}_{n}-\mathbf{x}_{j}\|_{2}^{2}\right)^{k}\;.$

For $d\geq 2$, taking $\mathbf{u}=\mathbf{u}_{j}$ orthogonal to
$\mathbf{x}_{n}-\mathbf{x}_{j}$, we get

$p_{j}(\mathbf{u}_{j})=\left((\varepsilon
z_{\ast})^{2}+\varepsilon^{2}\|\mathbf{x}_{n}-\mathbf{x}_{j}\|_{2}^{2}\right)^{k}=\left(e^{i\pi/k}+\varepsilon^{2}\|\mathbf{x}_{n}-\mathbf{x}_{j}\|_{2}^{2}\right)^{k}\;,$

which has nonzero imaginary part, since the (principal) argument of
$e^{i\pi/k}+\varepsilon^{2}\|\mathbf{x}_{n}-\mathbf{x}_{j}\|_{2}^{2}$ is
positive and less than $\pi/k$ (this is geometrically seen at a glance by
adding up the corresponding vectors in the complex plane). Then, also
$p_{j}(\mathbf{u}_{j})$ and thus $1+p_{j}(\mathbf{u}_{j})$ have nonzero
imaginary part, which means that the real polynomial
$\mbox{Im}(1+p_{j}(\mathbf{u}))$ is not identically zero on the unit sphere
and hence its zeros have null surface measure. The latter assertion can be
proved for example by writing the unit vector $\mathbf{u}$ in spherical
coordinates and observing that a nonzero real polynomial on the unit sphere
becomes a multivariate trigonometric polynomial in the coordinate box, to
which we can apply the fundamental result [10] on the zero set of real
analytic functions. Consequently, choosing a unit vector $\mathbf{u}$ not
belonging to the union in $j$ of the zero sets of $\mbox{Im}(1+p_{j}(u))$, the
polynomial $1+p_{j}(\mathbf{u})$ has nonzero imaginary part and thus the
functions $\phi_{j}(\mathbf{x}(z))$ are analytic at $z_{\ast}$ for every
$j<n$.

For $d=1$, $\Omega$ is an open bounded interval and with no loss of generality
we can interchange $x_{n}$ with $\max\\{x_{j}:\,1\leq j\leq n\\}$, which
simply corresponds to interchange two rows and columns in the matrix $A(x)$.
Choosing $u=1$ we get that the complex number

$p_{j}(1)=\left(e^{i\pi/k}+2\varepsilon(x_{n}-x_{j})e^{i\pi/(2k)}+\varepsilon^{2}(x_{n}-x_{j})^{2}\right)^{k}$

has positive real part, since the (principal) argument of
$e^{i\pi/k}+2\varepsilon(x_{n}-x_{j})e^{i\pi/(2k)}+\varepsilon^{2}(x_{n}-x_{j})^{2}$
is positive and less than $\pi/k$. Then also $1+p_{j}(1)$ has positive real
part and again the functions $\phi_{j}(x(z))$ are analytic at $z_{\ast}$ for
every $j<n$.

Since for $k\geq 1$ the functions $\phi_{j}(\mathbf{x}(z))$ are a.s. analytic
at $z_{\ast}$ for every $j<n$, also the functions $a(\mathbf{x}(z))$ and
$b(\mathbf{x}(z))$ are both a.s. analytic at $z_{\ast}$, recalling that
$a,b\in span\\{\phi_{1},\dots,\phi_{n-1}\\}$. Now, from
$f(\mathbf{x}(z))\equiv 0$ we would get from (3)

$det(V_{n-1})\left(1+(\varepsilon z)^{2k}\right)^{2\beta}\equiv
a(\mathbf{x}(z))\left(1+(\varepsilon
z)^{2k}\right)^{\beta}+b(\mathbf{x}(z))\;.$

We show now that this identity leads to a contradiction for all real values of
$\beta>1$, $\beta\not\in\mathbb{N}$. Let us begin with rational values of
$\beta$, namely $\beta=p/q$ with $p$, $q$ relatively prime. If $q=2$, we would
get

$det(V_{n-1})\left(1+(\varepsilon z)^{2k}\right)^{p}-b(\mathbf{x}(z))\equiv
a(\mathbf{x}(z))\left(1+(\varepsilon z)^{2k}\right)^{p/2}$

which cannot hold, because the function on the left-hand side is analytic at
$z_{\ast}$, whereas that on the right-hand side has a branch point of the
complex square root there. If $q\neq 2$, let us restrict the $z$-variable to
the segment $z=yz_{\ast}$ with $1-\delta<y<1$, where $\delta$ corresponds to a
neighborhood $|z-z_{\ast}|<\delta$ where $a(\mathbf{x}(z))$ and
$b(\mathbf{x}(z))$ are analytic. Then we would have

$det(V_{n-1})\left(1-y^{2k}\right)^{2p/q}\equiv
a(\mathbf{x}(yz_{\ast}))\left(1-y^{2k}\right)^{p/q}+b(\mathbf{x}(yz_{\ast}))\;,$

where we can write $1-y^{2k}=(1-y)(1+y+\dots+y^{2k-1})$. Moreover, by
analiticilty and using the Taylor expansions of $a(\mathbf{x}(z))$ and
$b(\mathbf{x}(z))$ at $z=z_{\ast}$ we would get $a(\mathbf{x}(yz_{\ast}))\sim
c_{s}(1-y)^{s}$ and $b(\mathbf{x}(yz_{\ast}))\sim d_{\nu}(1-y)^{\nu}$ as $y\to
1^{-}$, where $c_{s},d_{\nu}\neq 0$ and $s,\nu$ are the orders of the first
respective nonvanishing derivatives.

Again, this cannot hold because if $\nu=0$ we would have distinct limits on
the two sides as $y\to 1^{-}$. If $\nu>0$, looking at the infinitesimal orders
for $y\to 1^{-}$, it would require that $2p/q=\min\\{s+p/q,\nu\\}$ i.e. either
$p/q=s$ (that cannot hold since $\beta=p/q\not\in\mathbb{N}$) or $2p/q=\nu$
(that cannot hold since $p$ and $q$ are relatively prime and $q\neq 2$).
Finally, consider the case of $\beta$ irrational. We would have

$det(V_{n-1})\left(1-y^{2k}\right)^{2\beta}\equiv
a(\mathbf{x}(yz_{\ast}))\left(1-y^{2k}\right)^{\beta}+b(\mathbf{x}(yz_{\ast}))$

and taking again the limits as $y\to 1^{-}$, this leads to a contradiction for
$\nu=0$ (distinct limits on the two sides) and for $\nu>0$ because it would
require that either $2\beta=s+\beta$ i.e. $\beta=s\in\mathbb{N}$, or
$2\beta=\nu$ i.e. $\beta=\nu/2\in\mathbb{Q}$.

We can now conclude the proof. If $f$ in (3) is not identically zero, being
analytic in $\Omega$ by [10], its zero set, say $Z_{f}$, has null Lebesgue
measure, and hence null measure with respect to any continuous probability
distribution with density. On the other hand,
$f(\mathbf{x}_{n+1})=det(V_{n+1})$. Then, we can write

$\mbox{prob}\\{\mbox{det}(V_{n+1})=0\\}=\mbox{prob}\\{f(\mathbf{x}_{n+1})=0\\}$
$=\mbox{prob}\\{f\equiv 0\\}+\mbox{prob}\\{f\not\equiv
0\;\&\;\mathbf{x}_{n+1}\in Z_{f}\\}=0+0=0\;,$

and the inductive step is completed. $\square$

## Acknowledgements

Work partially supported by the DOR funds of the University of Padova, and by
the INdAM-GNCS 2024 Project “Kernel and polynomial methods for approximation
and integration: theory and application software”. This research has been
accomplished within the RITA “Research ITalian network on Approximation”, the
SIMAI Activity Group ANA&A and the UMI Group TAA “Approximation Theory and
Applications”.

## References

  * [1] L.P. Bos, K. Salkauskas, On the matrix $[|x_{i}-x_{j}|^{3}]$ and the cubic spline continuity equations, J. Approx. Theory 51 (1987), 81–88.
  * [2] L. Bos, A. Sommariva, M. Vianello, A note on polynomial-free unisolvence of polyharmonic splines at random points, arXiv:2312.13710.
  * [3] M.E. Chenoweth, S.A. Sarra, A Numerical Study of Generalized Multiquadric Radial Basis Function Interpolation, SIAM Undergraduate Research Online 2(2) (2009).
  * [4] F. Dell’Accio, A. Sommariva, M. Vianello, Random sampling and unisolvent interpolation by almost everywhere analytic functions, Appl. Math. Lett. 145 (2023), 108734.
  * [5] G.E. Fasshauer, Meshfree Approximation Methods with Matlab, Interdisciplinary Mathematical Sciences, Vol. 6, World Scientific, 2007.
  * [6] R.L. Hardy, Multiquadric equations of topography and other irregular surfaces, J. Geophys. Res. 76 (1971), 1905–1915.
  * [7] S.G. Krantz and H.R. Parks, A Primer of Real Analytic Functions, Second Edition, Birkhäuser, Boston, 2002.
  * [8] E. Larsson, R. Schaback, Scaling of radial basis functions, IMA Journal of Numerical Analysis (2023), drad035.
  * [9] C.A. Micchelli, Interpolation of scattered data: distance matrices and conditionally positive definite functions, Constr. Approx. 2 (1986), 11–22.
  * [10] B.S. Mityagin, The Zero Set of a Real Analytic Function, Math. Notes 107 (2020), 529–530.
  * [11] M. Ortmann, M. Buhmann, High Accuracy Quasi-Interpolation using a new class of generalized Multiquadrics, J. Math. Anal. Appl., available online 29 March 2024, 128359.
  * [12] A. Pasioti, On the Constrained Solution of RBF Surface Approximation, MDPI Mathematics 10 (2022), 2582.
  * [13] A. Sommariva, M. Vianello, Polynomial-free unisolvence of polyharmonic splines with odd exponent by random sampling, arXiv:2401.13322.
  * [14] Y. Xu, A. Narayan, Randomized weakly admissible meshes, J. Approx. Theory 285 (2023), 105835.
  * [15] H. Wendland, Scattered Data Approximation, Cambridge Monogr. Appl. Comput. Math., vol. 17, Cambridge Univ. Press, Cambridge, 2005.