Fast fitting of radial basis functions: Methods based on preconditioned GMRES iteration

Solving large radial basis function (RBF) interpolation problems with non‐customised methods is computationally expensive and the matrices that occur are typically badly conditioned. For example, using the usual direct methods to fit an RBF with N centres requires O(N ²) storage and O(N ³) flops. Thus such an approach is not viable for large problems with N ≥10,000. In this paper we present preconditioning strategies which, in combination with fast matrix–vector multiplication and GMRES iteration, make the solution of large RBF interpolation problems orders of magnitude less expensive in storage and operations. In numerical experiments with thin‐plate spline and multiquadric RBFs the preconditioning typically results in dramatic clustering of eigenvalues and improves the condition numbers of the interpolation problem by several orders of magnitude. As a result of the eigenvalue clustering the number of GMRES iterations required to solve the preconditioned problem is of the order of 10-20. Taken together, the combination of a suitable approximate cardinal function preconditioner, the GMRES iterative method, and existing fast matrix–vector algorithms for RBFs [4,5] reduce the computational cost of solving an RBF interpolation problem to O(N) storage, and O(N \log N) operations. This revised version was published online in June 2006 with corrections to the Cover Date.     

Spline approximation of a non-linear Riemann–Hilbert problem†

In this article we use linear spline approximation of a non-linear Riemann–Hilbert problem on the unit disk. The boundary condition for the holomorphic function is reformulated as a non-linear singular integral equation A(u) = 0, where A : H ¹(Γ) → H ¹(Γ) is defined via a Nemytski operator. We approximate A by A ⁿ : H ¹(Γ) → H ¹(Γ) using spline collocation and show that this defines a Fredholm quasi-ruled mapping. Following the results of (A.I. ?nirel'man, The degree of quasi-ruled mapping and a nonlinear Hilbert problem, Math. USSR-Sbornik 18 (1972), pp. 373–396; M.A. Efendiev, On a property of the conjugate integral and a nonlinear Hilbert problem, Soviet Math. Dokl. 35 (1987), pp. 535–539; M.A. Efendiev, W.L. Wendland, Nonlinear Riemann–Hilbert problems for multiply connected domains, Nonlinear Anal. 27 (1996), pp. 37–58; Nonlinear Riemann–Hilbert problems without transversality. Math. Nachr. 183 (1997), pp. 73–89; Nonlinear Riemann–Hilbert problems for doubly connected domains and closed boundary data, Topol. Methods Nonlinear Anal. 17 (2001), pp. 111–124; Nonlinear Riemann–Hilbert problems with Lipschitz, continuous boundary data without transversality, Nonlinear Anal. 47 (2001), pp. 457–466; Nonlinear Riemann–Hilbert problems with Lipschitz-continuous boundary data: Doubly connected domains, Proc. Roy. Soc. London Ser. A 459 (2003), pp. 945–955.), we define a degree of mapping and show the existence of the spline solutions of the fully discrete equations A ⁿ (u) = 0, for n large enough. We conclude this article by discussing the solvability of the non-linear collocation method, where we shall need an additional uniform strong ellipticity condition for employing the spline approximation.     

Cubic spline functions and initial value problems

A numerical process is presented which provides a cubic spline function approximation for the solution of initial value problems in ordinary differential equations. With interpolate cubic spline functions we can achieveO(h ⁴) convergence.     

Interior Capacities of Condensers in Locally Compact Spaces

The study deals with the theory of interior capacities of condensers in a locally compact space, a condenser being treated here as a finite collection of arbitrary sets with sign + 1 or − 1 prescribed such that the closures of oppositely signed sets are mutually disjoint. We are motivated by the known fact that, in the noncompact case, the main minimum-problem of the theory is in general unsolvable, and this occurs even under very natural assumptions (e.g., for the Newtonian, Green, or Riesz kernels in \mathbb Rⁿ\mathbb R^n and closed condensers). Therefore it was particularly interesting to find statements of variational problems dual to the main minimum-problem (and hence providing new equivalent definitions to the capacity), but now always solvable (e.g., even for nonclosed condensers). For all positive definite kernels satisfying Fuglede’s condition of consistency between the strong and vague (= weak*) topologies, problems with the desired properties are posed and solved. Their solutions provide a natural generalization of the well-known notion of interior equilibrium measures associated with a set. We describe those solutions and the corresponding equilibrium constants, analyze their uniqueness and continuity, and point out their characteristic properties. Such results are new even for classical kernels in \mathbb Rⁿ\mathbb R^n, which is important in applications.     

Fast generalized cross validation using Krylov subspace methods

The task of fitting smoothing spline surfaces to meteorological data such as temperature or rainfall observations is computationally intensive. The generalized cross validation (GCV) smoothing algorithm, if implemented using direct matrix techniques, is O(n ³) computationally, and memory requirements are O(n ²). Thus, for data sets larger than a few hundred observations, the algorithm is prohibitively slow. The core of the algorithm consists of solving series of shifted linear systems, and iterative techniques have been used to lower the computational complexity and facilitate implementation on a variety of supercomputer architectures. For large data sets though, the execution time is still quite high. In this paper we describe a Lanczos based approach that avoids explicitly solving the linear systems and dramatically reduces the amount of time required to fit surfaces to sets of data.      

Fast numerical solution of fredholm integral equations with stationary kernels

A fast recursive matrix method for the numerical solution of Fredholm integral equations with stationary kernels is derived. IfN denotes the number of nodal points, the complexity of the algorithm isO(N ²), which should be compared toO(N ³) for conventional algorithms for solving such problems. The method is related to fast algorithms for inverting Toeplitz matrices.Applications to equations of the first and second kind as well as miscellaneous problems are discussed and illustrated with numerical examples. These show that the theoretical improvement in efficiency is indeed obtained, and that no problems with numerical stability or accuracy are encountered.     

Hierarchies of higher order kernels

Summary Recent literature on functional estimation has shown the importance of kernels with vanishing moments although no general framework was given to build kernels of increasing order apart from some specific methods based on moment relationships. The purpose of the present paper is to develop such a framework and to show how to build higher order kernels with nice properties and to solve optimization problems about kernels. The proofs given here, unlike standard variational arguments, explain why some hierarchies of kernels do have optimality properties. Applications are given to functional estimation in a general context. In the last section special attention is paid to density estimates based on kernels of order (m, r), i.e., kernels of orderr for estimation of derivatives of orderm. Convergence theorems are easily derived from interpretation by means of projections inL ² spaces.     

The application of optimal control methodology to nonlinear programming problems

Dynamic programming techniques have proven to be more successful than alternative nonlinear programming algorithms for solving many discrete-time optimal control problems. The reason for this is that, because of the stagewise decomposition which characterizes dynamic programming, the computational burden grows approximately linearly with the numbern of decision times, whereas the burden for other methods tends to grow faster (e.g.,n ³ for Newton's method). The idea motivating the present study is that the advantages of dynamic programming can be brought to bear on classical nonlinear programming problems if only they can somehow be rephrased as optimal control problems.As shown herein, it is indeed the case that many prominent problems in the nonlinear programming literature can be viewed as optimal control problems, and for these problems, modern dynamic programming methodology is competitive with respect to processing time. The mechanism behind this success is that such methodology achieves quadratic convergence without requiring solution of large systems of linear equations.     

The Gibbs Phenomenon for n-Dimensional Kernels

We give conditions which ensure the existence in Rⁿ of the Gibbs phenomenon for a large class of kernels non necessarily rotation or translation invariant.     

Time dependent Lyapunov functions for some Kolmogorov semigroups perturbed by unbounded potentials

We study global regularity properties of transition kernels associated to second order differential operators in \mathbb R^N{\mathbb {R}^N} with unbounded drift and potential terms. Under suitable conditions, we prove pointwise upper bounds. We use time dependent Lyapunov function techniques allowing us to gain a better time behaviour of such kernels.     

Embedding problems with cyclic kernel of order 4

We consider Galois theoretical embedding problems with kernelC ₄, and prove that such an embedding problem can be ‘constructively’ reduced to two embedding problems, where the kernels are groups of roots of unity.     

Translation-Invariant Bilinear Operators with Positive Kernels

We study L ^r (or L ^{r, ∞}) boundedness for bilinear translation-invariant operators with nonnegative kernels acting on functions on \mathbb Rⁿ{\mathbb {R}^n}. We prove that if such operators are bounded on some products of Lebesgue spaces, then their kernels must necessarily be integrable functions on \mathbb R²ⁿ{\mathbb R^{2n}}, while via a counterexample we show that the converse statement is not valid. We provide certain necessary and some sufficient conditions on nonnegative kernels yielding boundedness for the corresponding operators on products of Lebesgue spaces. We also prove that, unlike the linear case where boundedness from L ¹ to L ¹ and from L ¹ to L ^{1, ∞} are equivalent properties, boundedness from L ¹ × L ¹ to L ^1/2 and from L ¹ × L ¹ to L ^{1/2, ∞} may not be equivalent properties for bilinear translation-invariant operators with nonnegative kernels.     

A two-dimensional minimum-derivative spline

We consider the construction of a C ^(1,1) interpolation parabolic spline function of two variables on a uniform rectangular grid, i.e., a function continuous in a given region together with its first partial derivatives which on every partial grid rectangle is a polynomial of second degree in x and second degree in y. The spline function is constructed as a minimum-derivative one-dimensional quadratic spline in one of the variables, and the spline coefficients themselves are minimum-derivative quadratic spline functions of the other variable.     

Finding an interior point in the optimal face of linear programs

We study the problem of finding a point in the relative interior of the optimal face of a linear program. We prove that in the worst case such a point can be obtained in O(n ³ L) arithmetic operations. This complexity is the same as the complexity for solving a linear program. We also show how to find such a point in practice. We report and discuss computational results obtained for the linear programming problems in the NETLIB test set.Research supported in part by NSF Grant CCR-8810107, CCR-9019469 and a grant from GTE Laboratories.Research supported in part by NSF Grant DDM-8922636 and NSF Coop. Agr. No. CCR-8809615 through Rice University.     

Tournament matrices and their generalizations,I.

If M is any complex matrix with rank (M + M ^* + I) = 1, we show that any eigenvalue of M that is not geometrically simple has 1/2 for its real part. This generalizes a recent finding of de Caen and Hoffman: the rank of any n × n tournament matrix is at least n ? 1. We extend several spectral properties of tournament matrices to this and related types of matrices. For example, we characterize the singular real matrices M with 0 diagonal for which rank (M + M^T + I) = 1 and we characterize the vectors that can be in the kernels of such matrices. We show that singular, irreducible n × n tournament matrices exist if and only n? {2,3,4,5} and exhibit many infinite families of such matrices. Connections with signed digraphs are explored and several open problems are presented.     

A Crank‐Nicolson Petrov‐Galerkin method with quadrature for semi‐linear parabolic problems

We propose and analyze an application of a fully discrete C² spline quadrature Petrov‐Galerkin method for spatial discretization of semi‐linear parabolic initial‐boundary value problems on rectangular domains. We prove second order in time and optimal order H¹ norm convergence in space for the extrapolated Crank‐Nicolson quadrature Petrov‐Galerkin scheme. We demonstrate numerically both L² and H¹ norm optimal order convergence of the scheme even if the nonlinear source term is not smooth. © 2005 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2005.     

Reproducing kernels of generalized Sobolev spaces via a Green function approach with distributional operators

In this paper we introduce a generalized Sobolev space by defining a semi-inner product formulated in terms of a vector distributional operator P consisting of finitely or countably many distributional operators P _n , which are defined on the dual space of the Schwartz space. The types of operators we consider include not only differential operators, but also more general distributional operators such as pseudo-differential operators. We deduce that a certain appropriate full-space Green function G with respect to L := P ^*T P now becomes a conditionally positive function. In order to support this claim we ensure that the distributional adjoint operator P ^* of P is well-defined in the distributional sense. Under sufficient conditions, the native space (reproducing-kernel Hilbert space) associated with the Green function G can be embedded into or even be equivalent to a generalized Sobolev space. As an application, we take linear combinations of translates of the Green function with possibly added polynomial terms and construct a multivariate minimum-norm interpolant s _f,X to data values sampled from an unknown generalized Sobolev function f at data sites located in some set X ì \mathbbR^d{X \subset \mathbb{R}^d}. We provide several examples, such as Matérn kernels or Gaussian kernels, that illustrate how many reproducing-kernel Hilbert spaces of well-known reproducing kernels are equivalent to a generalized Sobolev space. These examples further illustrate how we can rescale the Sobolev spaces by the vector distributional operator P. Introducing the notion of scale as part of the definition of a generalized Sobolev space may help us to choose the “best” kernel function for kernel-based approximation methods.     

Matrix decomposition algorithms for elliptic boundary value problems: a survey

We provide an overview of matrix decomposition algorithms (MDAs) for the solution of systems of linear equations arising when various discretization techniques are applied in the numerical solution of certain separable elliptic boundary value problems in the unit square. An MDA is a direct method which reduces the algebraic problem to one of solving a set of independent one-dimensional problems which are generally banded, block tridiagonal, or almost block diagonal. Often, fast Fourier transforms (FFTs) can be employed in an MDA with a resulting computational cost of O(N ² logN) on an N × N uniform partition of the unit square. To formulate MDAs, we require knowledge of the eigenvalues and eigenvectors of matrices arising in corresponding two–point boundary value problems in one space dimension. In many important cases, these eigensystems are known explicitly, while in others, they must be computed. The first MDAs were formulated almost fifty years ago, for finite difference methods. Herein, we discuss more recent developments in the formulation and application of MDAs in spline collocation, finite element Galerkin and spectral methods, and the method of fundamental solutions. For ease of exposition, we focus primarily on the Dirichlet problem for Poisson’s equation in the unit square, sketch extensions to other boundary conditions and to more involved elliptic problems, including the biharmonic Dirichlet problem, and report extensions to three dimensional problems in a cube. MDAs have also been used extensively as preconditioners in iterative methods for solving linear systems arising from discretizations of non-separable boundary value problems.     

Integral representation of slowly growing equidistant splines

In this paper we consider polynomial splines S(x) with equidistant nodes which may grow as O (|x|^s). We present an integral representation of such splines with a distribution kernel. This representation is related to the Fourier integral of slowly growing functions. The part of the Fourier exponentials herewith play the so called exponential splines by Schoenberg. The integral representation provides a flexible tool for dealing with the growing equidistant splines. First, it allows us to construct a rich library of splines possessing the property that translations of any such spline form a basis of corresponding spline space. It is shown that any such spline is associated with a dual spline whose translations form a biorthogonal basis. As examples we present solutions of the problems of projection of a growing function onto spline spaces and of spline interpolation of a growing function. We derive formulas for approximate evaluation of splines projecting a function onto the spline space and establish therewith exact estimations of the approximation errors.     

A Petrov‐Galerkin method with quadrature for semi‐linear second‐order hyperbolic problems

We propose and analyze a Crank–Nicolson quadrature Petrov–Galerkin (CNQPG) ‐spline method for solving semi‐linear second‐order hyperbolic initial‐boundary value problems. We prove second‐order convergence in time and optimal order H² norm convergence in space for the CNQPG scheme that requires only linear algebraic solvers. We demonstrate numerically optimal order H^k, k = 0,1,2, norm convergence of the scheme for some test problems with smooth and nonsmooth nonlinearities. © 2006 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2006