Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

Fast multi-dimensional scattered data approximation with Neumann boundary conditions

Denis Grishin, Thomas Strohmer

TL;DR

The paper addresses stable reconstruction of a function from nonuniform samples by imposing Neumann boundary conditions and using cosine polynomials, mitigating boundary artifacts that arise with exponential bases. It shows that the resulting least-squares system has a scaled Toeplitz+Hankel structure A = D(T+H)D with a_k = 0.5 * sum_j w_j cos(pi k x_j) and derives a condition-number bound kappa(A) ≤ ((1+delta M)^2)/((1-delta M)^2) under a max-gap delta < 1/M; A is invertible when M < r. A fast solution is developed using conjugate gradient with fast A-vec products via the 1D DCT-I and, for nonuniform sampling, nonuniform DCT (NDCT) to compute initial quantities, achieving O(M log M) per CG iteration; the method generalizes to higher dimensions through block Toeplitz+Hankel structures and 2D DCT-I. A multilevel extension automatically selects the degree M by enforcing a discrepancy principle. The approach is validated on a geophysical scattered-data problem, showing reduced boundary artifacts and improved accuracy over standard methods such as ACT and cubic splines.

Abstract

An important problem in applications is the approximation of a function $f$ from a finite set of randomly scattered data $f(x_j)$. A common and powerful approach is to construct a trigonometric least squares approximation based on the set of exponentials $\{e^{2πi kx}\}$. This leads to fast numerical algorithms, but suffers from disturbing boundary effects due to the underlying periodicity assumption on the data, an assumption that is rarely satisfied in practice. To overcome this drawback we impose Neumann boundary conditions on the data. This implies the use of cosine polynomials $\cos (πkx)$ as basis functions. We show that scattered data approximation using cosine polynomials leads to a least squares problem involving certain Toeplitz+Hankel matrices. We derive estimates on the condition number of these matrices. Unlike other Toeplitz+Hankel matrices, the Toeplitz+Hankel matrices arising in our context cannot be diagonalized by the discrete cosine transform, but they still allow a fast matrix-vector multiplication via DCT which gives rise to fast conjugate gradient type algorithms. We show how the results can be generalized to higher dimensions. Finally we demonstrate the performance of the proposed method by applying it to a two-dimensional geophysical scattered data problem.

Fast multi-dimensional scattered data approximation with Neumann boundary conditions

TL;DR

The paper addresses stable reconstruction of a function from nonuniform samples by imposing Neumann boundary conditions and using cosine polynomials, mitigating boundary artifacts that arise with exponential bases. It shows that the resulting least-squares system has a scaled Toeplitz+Hankel structure A = D(T+H)D with a_k = 0.5 * sum_j w_j cos(pi k x_j) and derives a condition-number bound kappa(A) ≤ ((1+delta M)^2)/((1-delta M)^2) under a max-gap delta < 1/M; A is invertible when M < r. A fast solution is developed using conjugate gradient with fast A-vec products via the 1D DCT-I and, for nonuniform sampling, nonuniform DCT (NDCT) to compute initial quantities, achieving O(M log M) per CG iteration; the method generalizes to higher dimensions through block Toeplitz+Hankel structures and 2D DCT-I. A multilevel extension automatically selects the degree M by enforcing a discrepancy principle. The approach is validated on a geophysical scattered-data problem, showing reduced boundary artifacts and improved accuracy over standard methods such as ACT and cubic splines.

Abstract

An important problem in applications is the approximation of a function from a finite set of randomly scattered data . A common and powerful approach is to construct a trigonometric least squares approximation based on the set of exponentials . This leads to fast numerical algorithms, but suffers from disturbing boundary effects due to the underlying periodicity assumption on the data, an assumption that is rarely satisfied in practice. To overcome this drawback we impose Neumann boundary conditions on the data. This implies the use of cosine polynomials as basis functions. We show that scattered data approximation using cosine polynomials leads to a least squares problem involving certain Toeplitz+Hankel matrices. We derive estimates on the condition number of these matrices. Unlike other Toeplitz+Hankel matrices, the Toeplitz+Hankel matrices arising in our context cannot be diagonalized by the discrete cosine transform, but they still allow a fast matrix-vector multiplication via DCT which gives rise to fast conjugate gradient type algorithms. We show how the results can be generalized to higher dimensions. Finally we demonstrate the performance of the proposed method by applying it to a two-dimensional geophysical scattered data problem.

Paper Structure

This paper contains 6 sections, 3 theorems, 65 equations, 1 figure, 1 algorithm.

Table of Contents

  1. Introduction
  2. Nonuniform sampling, cosine polynomials, and Toeplitz+Hankel matrices
  3. Fast computation of the least squares approximation
  4. Multilevel scattered data approximation
  5. Two-dimensional scattered data approximation
  6. Numerical experiments: An example from geophysics

Key Result

theorem 1

Assume we are given nonuniformly spaced sampling points $\{x_j\}_{j=1}^r \in [0,1]$, sampling values $s=\{s_j\}_{j=1}^r$ and positive weights $\{w_j\}_{j=1}^r$. Define $A:= V^{T} V$, where $V$ is as in vander, and set $b=V^{T}{s^{(w)}}$. There holds: (i) The matrix $A$ is a scaled Toeplitz+Hankel ma where with and $D=\operatorname{diag} (\frac{1}{\sqrt{2}},1,\dots,1)$. (ii) If $M < r$ then $A$ i

Figures (1)

  • Figure 1: Approximation of gravitational anomaly from noisy scattered data (5% noise) by proposed method and comparison to standard algorithms.

Theorems & Definitions (5)

  • theorem 1
  • definition 1
  • theorem 2
  • definition 2
  • theorem 3