Daan Huybrechs

Haiyong Wang, Daan Huybrechs, Stefan Vandewalle

Barycentric interpolation is arguably the method of choice for numerical polynomial interpolation. The polynomial interpolant is expressed in terms of function values using the so-called barycentric weights, which depend on the interpolation points. Few explicit formulae for these barycentric weights are known. In [H. Wang and S. Xiang, Math. Comp., 81 (2012), 861--877], the authors have shown that the barycentric weights of the roots of Legendre polynomials can be expressed explicitly in terms of the weights of the corresponding Gaussian quadrature rule. This idea was subsequently implemented in the Chebfun package [L. N. Trefethen and others, The Chebfun Development Team, 2011] and in the process generalized by the Chebfun authors to the roots of Jacobi, Laguerre and Hermite polynomials. In this paper, we explore the generality of the link between barycentric weights and Gaussian quadrature and show that such relationships are related to the existence of lowering operators for orthogonal polynomials. We supply an exhaustive list of cases, in which all known formulae are recovered and also some new formulae are derived, including the barycentric weights for Gauss-Radau and Gauss-Lobatto points. Based on a fast ${\mathcal O}(n)$ algorithm for the computation of Gaussian quadrature, due to Hale and Townsend, this leads to an ${\mathcal O}(n)$ computational scheme for barycentric weights.

Ben Adcock, Daan Huybrechs, Jesus Martin-Vaquero

An effective means to approximate an analytic, nonperiodic function on a bounded interval is by using a Fourier series on a larger domain. When constructed appropriately, this so-called Fourier extension is known to converge geometrically fast in the truncation parameter. Unfortunately, computing a Fourier extension requires solving an ill-conditioned linear system, and hence one might expect such rapid convergence to be destroyed when carrying out computations in finite precision. The purpose of this paper is to show that this is not the case. Specifically, we show that Fourier extensions are actually numerically stable when implemented in finite arithmetic, and achieve a convergence rate that is at least superalgebraic. Thus, in this instance, ill-conditioning of the linear system does not prohibit a good approximation. In the second part of this paper we consider the issue of computing Fourier extensions from equispaced data. A result of Platte, Trefethen & Kuijlaars states that no method for this problem can be both numerically stable and exponentially convergent. We explain how Fourier extensions relate to this theoretical barrier, and demonstrate that they are particularly well suited for this problem: namely, they obtain at least superalgebraic convergence in a numerically stable manner.

Ben Adcock, Daan Huybrechs

Functions that are smooth but non-periodic on a certain interval possess Fourier series that lack uniform convergence and suffer from the Gibbs phenomenon. However, they can be represented accurately by a Fourier series that is periodic on a larger interval. This is commonly called a Fourier extension. When constructed in a particular manner, Fourier extensions share many of the same features of a standard Fourier series. In particular, one can compute Fourier extensions which converge spectrally fast whenever the function is smooth, and exponentially fast if the function is analytic, much the same as the Fourier series of a smooth/analytic and periodic function. With this in mind, the purpose of this paper is to describe, analyze and explain the observation that Fourier extensions, much like classical Fourier series, also have excellent resolution properties for representing oscillatory functions. The resolution power, or required number of degrees of freedom per wavelength, depends on a user-controlled parameter and, as we show, it varies between 2 and π. The former value is optimal and is achieved by classical Fourier series for periodic functions, for example. The latter value is the resolution power of algebraic polynomial approximations. Thus, Fourier extensions with an appropriate choice of parameter are eminently suitable for problems with moderate to high degrees of oscillation.

Haiyong Wang, Lun Zhang, Daan Huybrechs

In this paper, we study the asymptotics and fast computation of the one-sided oscillatory Hilbert transforms of the form $$H^{+}(f(t)e^{iωt})(x)=-int_{0}^{\infty}e^{iωt}\frac{f(t)}{t-x}dt,\qquad ω>0,\qquad x\geq 0,$$ where the bar indicates the Cauchy principal value and $f$ is a real-valued function with analytic continuation in the first quadrant, except possibly a branch point of algebraic type at the origin. When $x=0$, the integral is interpreted as a Hadamard finite-part integral, provided it is divergent. Asymptotic expansions in inverse powers of $ω$ are derived for each fixed $x\geq 0$, which clarify the large $ω$ behavior of this transform. We then present efficient and affordable approaches for numerical evaluation of such oscillatory transforms. Depending on the position of $x$, we classify our discussion into three regimes, namely, $x=\mathcal{O}(1)$ or $x\gg1$, $0<x\ll 1$ and $x=0$. Numerical experiments show that the convergence of the proposed methods greatly improve when the frequency $ω$ increases. Some extensions to oscillatory Hilbert transforms with Bessel oscillators are briefly discussed as well.

Andreas Asheim, Alfredo Deaño, Daan Huybrechs et al.

We investigate a Gaussian quadrature rule and the corresponding orthogonal polynomials for the oscillatory weight function $e^{iωx}$ on the interval $[-1,1]$. We show that such a rule attains high asymptotic order, in the sense that the quadrature error quickly decreases as a function of the frequency $ω$. However, accuracy is maintained for all values of $ω$ and in particular the rule elegantly reduces to the classical Gauss-Legendre rule as $ω\to 0$. The construction of such rules is briefly discussed, and though not all orthogonal polynomials exist, it is demonstrated numerically that rules with an even number of points are always well defined. We show that these rules are optimal both in terms of asymptotic order as well as in terms of polynomial order.

Daan Huybrechs, Anda-Elena Olteanu

The Wave Based Method (WBM) is a Trefftz method for the simulation of wave problems in vibroacoustics. Like other Trefftz methods, it employs a non-standard discretisation basis consisting of solutions of the partial differential equation (PDE) at hand. We analyse the convergence and numerical stability of the Wave Based Method for Helmholtz problems using tools from approximation theory. We show that the set of discretisation functions more closely resembles a frame, a redundant set of functions, than a basis. The redundancy of a frame typically leads to ill-conditioning, which indeed is common in Trefftz methods. Recent theoretical results on frames for function approximation suggest that the associated ill-conditioned system matrix can be successfully regularised, with error bounds available, when using a discrete least squares approach. While the original Wave Based Method is based on a weighted residual formulation, in this paper we pursue an oversampled collocation approach instead. We show that, for smooth scattering obstacles in two dimensions, the results closely follow the theory of frames. We identify cases where the method achieves very high accuracy whilst providing a solution with small norm coefficients, in spite of ill-conditioning. Moreover, the accurate results are reliably maintained even in parameter regimes associated with extremely high ill-conditioning.

Daan Huybrechs

We consider the computation of quadrature rules that are exact for a Chebyshev set of linearly independent functions on an interval $[a,b]$. A general theory of Chebyshev sets guarantees the existence of rules with a Gaussian property, in the sense that $2l$ basis functions can be integrated exactly with just $l$ points and weights. Moreover, all weights are positive and the points lie inside the interval $[a,b]$. However, the points are not the roots of an orthogonal polynomial or any other known special function as in the case of regular Gaussian quadrature. The rules are characterized by a nonlinear system of equations, and earlier numerical methods have mostly focused on finding suitable starting values for a Newton iteration to solve this system. In this paper we describe an alternative scheme that is robust and generally applicable for so-called complete Chebyshev sets. These are ordered Chebyshev sets where the first $k$ elements also form a Chebyshev set for each $k$. The points of the quadrature rule are computed one by one, increasing exactness of the rule in each step. Each step reduces to finding the unique root of a univariate and monotonic function. As such, the scheme of this paper is guaranteed to succeed. The quadrature rules are of interest for integrals with non-smooth integrands that are not well approximated by polynomials.

Kobe Bruyninckx, Daan Huybrechs, Karl Meerbergen

In parametrized linear systems $\mathsf{P}(μ)\mathsf{x}=\mathsf{b}$ the system matrix $\mathsf{P}$ depends nonlinearly on a parameter $μ$ and solutions are sought for many values of this parameter. We show that the compact rational Krylov (CORK) framework, originally introduced to solve nonlinear eigenvalue problems, can be used to efficiently produce approximate solutions to such a system for many values of the parameter at once. In this approach the parametrized system is first linearized, resulting in a large shifted linear system $(\boldsymbol{\mathsf{A}}-μ\boldsymbol{\mathsf{B}})\boldsymbol{\mathsf{y}}=\boldsymbol{\mathsf{d}}$. We formulate a left- and right-preconditioned rational Krylov GMRES method for shifted linear systems and show how the CORK framework can be used to speed up these methods in the setting of parametrized linear systems. Additionally, we show how to incorporate a right-hand side $\mathsf{b}(μ)$ that also depends on the parameter, how to choose the shifts to steer convergence and how to allow for inexact solves at these shifts throughout the iterations. As an application we consider the 'frequency sweeping' of Helmholtz scattering problems through the Boundary Element Method (BEM), enabled via an efficient representation of the dense but data-sparse wavenumber-dependent system matrix.

Daan Huybrechs

It is difficult to describe in mathematical terms what a neural network trained on data represents. On the other hand, there is a growing mathematical understanding of what neural networks are in principle capable of representing. Feedforward neural networks using the ReLU activation function represent continuous and piecewise linear functions and can approximate many others. The study of their expressivity addresses the question: which ones? Contributing to the available answers, we take the perspective of a neural network as an algorithm. In this analogy, a neural network is programmed constructively, rather than trained from data. An interesting example is a sorting algorithm: we explicitly construct a neural network that sorts its inputs exactly, not approximately, and that, in a sense, has optimal computational complexity if the input dimension is large. Such constructed networks may have several billion parameters. We construct and analyze several other examples, both existing and new. We find that, in these examples, neural networks as algorithms are typically recursive and parallel. Compared to conventional algorithms, ReLU networks are restricted by having to be continuous. Moreover, the depth of recursion is limited by the depth of the network, with deep networks having superior properties over shallow ones.

Daan Huybrechs, Peter Opsomer

Wave propagation and acoustic scattering problems require vast computational resources to be solved accurately at high frequencies. Asymptotic methods can make this cost potentially frequency independent by explicitly extracting the oscillatory properties of the solution. However, the high-frequency wave pattern becomes very complicated in the presence of multiple scattering obstacles. We consider a boundary integral equation formulation of the Helmholtz equation in two dimensions involving several obstacles, for which ray tracing schemes have been previously proposed. The existing analysis of ray tracing schemes focuses on periodic orbits between a subset of the obstacles. One observes that the densities on each of the obstacles converge to an equilibrium after a few iterations. In this paper we present an asymptotic approximation of the phases of those densities in equilibrium, in the form of a Taylor series. The densities represent a full cycle of reflections in a periodic orbit. We initially exploit symmetry in the case of two circular scatterers, but also provide an explicit algorithm for an arbitrary number of general 2D obstacles. The coefficients, as well as the time to compute them, are independent of the wavenumber and of the incident wave. The results may be used to accelerate ray tracing schemes after a small number of initial iterations.

Roel Matthysen, Daan Huybrechs

Fourier extension is an approximation scheme in which a function on an arbitary bounded domain is approximated using a classical Fourier series on a bounding box. On the smaller domain the Fourier series exhibits redundancy, and it has the mathematical structure of a frame rather than a basis. It is not trivial to construct approximations in this frame using function evaluations in points that belong to the domain only, but one way to do so is through a discrete least squares approximation. The corresponding system is extremely ill-conditioned, due to the redundancy in the frame, yet its solution via a regularized SVD is known to be accurate to very high (and nearly spectral) precision. Still, this computation requires ${\mathcal O}(N^3)$ operations. In this paper we describe an algorithm to compute such Fourier extension frame approximations in only ${\mathcal O}(N^2 \log^2 N)$ operations for general 2D domains. The cost improves to ${\mathcal O}(N \log^2N)$ operations for simpler tensor-product domains. The algorithm exploits a phenomenon called the plunge region in the analysis of time-frequency localization operators, which manifests itself here as a sudden drop in the singular values of the least squares matrix. It is known that the size of the plunge region scales like ${\mathcal O}(\log N)$ in one dimensional problems. In this paper we show that for most 2D domains in the fully discrete case the plunge region scales like ${\mathcal O}(N \log N)$, proving a discrete equivalent of a result that was conjectured by Widom for a related continuous problem. The complexity estimate depends on the Minkowski or box-counting dimension of the domain boundary, and as such it is larger than ${\mathcal O}(N \log N)$ for domains with fractal shape.

Daan Huybrechs, Peter Opsomer

Wave propagation and scattering problems in acoustics are often solved with boundary element methods. They lead to a discretization matrix that is typically dense and large: its size and condition number grow with increasing frequency. Yet, high frequency scattering problems are intrinsically local in nature, which is well represented by highly localized rays bouncing around. Asymptotic methods can be used to reduce the size of the linear system, even making it frequency independent, by explicitly extracting the oscillatory properties from the solution using ray tracing or analogous techniques. However, ray tracing becomes expensive or even intractable in the presence of (multiple) scattering obstacles with complicated geometries. In this paper, we start from the same discretization that constructs the fully resolved large and dense matrix, and achieve asymptotic compression by explicitly localizing the Green's function instead. This results in a large but sparse matrix, with a faster associated matrix-vector product and, as numerical experiments indicate, a much improved condition number. Though an appropriate localisation of the Green's function also depends on asymptotic information unavailable for general geometries, we can construct it adaptively in a frequency sweep from small to large frequencies in a way which automatically takes into account a general incident wave. We show that the approach is robust with respect to non-convex, multiple and even near-trapping domains, though the compression rate is clearly lower in the latter case. Furthermore, in spite of its asymptotic nature, the method is robust with respect to low-order discretizations such as piecewise constants, linears or cubics, commonly used in applications. On the other hand, we do not decrease the total number of degrees of freedom compared to a conventional classical discretization. The combination of the ...

Roel Matthysen, Daan Huybrechs

Fourier series of smooth, non-periodic functions on $[-1,1]$ are known to exhibit the Gibbs phenomenon, and exhibit overall slow convergence. One way of overcoming these problems is by using a Fourier series on a larger domain, say $[-T,T]$ with $T>1$, a technique called Fourier extension or Fourier continuation. When constructed as the discrete least squares minimizer in equidistant points, the Fourier extension has been shown shown to converge geometrically in the truncation parameter $N$. A fast ${\mathcal O}(N \log^2 N)$ algorithm has been described to compute Fourier extensions for the case where $T=2$, compared to ${\mathcal O}(N^3)$ for solving the dense discrete least squares problem. We present two ${\mathcal O}(N\log^2 N )$ algorithms for the computation of these approximations for the case of general $T$, made possible by exploiting the connection between Fourier extensions and Prolate Spheroidal Wave theory. The first algorithm is based on the explicit computation of so-called periodic discrete prolate spheroidal sequences, while the second algorithm is purely algebraic and only implicitly based on the theory.

Haiyong Wang, Daan Huybrechs

Chebyshev expansion coefficients can be computed efficiently by using the FFT, and for smooth functions the resulting approximation is close to optimal, with computations that are numerically stable. Given sufficiently accurate function samples, the Chebyshev expansion coefficients can be computed to machine precision accuracy. However, the accuracy is only with respect to absolute error, and this implies that very small expansion coefficients typically have very large relative error. Upon differentiating a Chebyshev expansion, this relative error in the small coefficients is magnified and accuracy may be lost, especially after repeated differentiation. At first sight, this seems unavoidable. Yet, in this paper, we focus on an alternative computation of Chebyshev expansion coefficients using contour integrals in the complex plane. The main result is that the coefficients can be computed with machine precision relative error, rather than absolute error. This implies that even very small coefficients can be computed with full floating point accuracy, even when they are themselves much smaller than machine precision. As a result, no accuracy is lost after differentiating the expansion, and even the 100th derivative of an analytic function can be computed with near machine precision accuracy using standard floating point arithmetic. In some cases, the contour integrals can be evaluated using the FFT, making the approach both highly accurate and fast.

A. Deano, D. Huybrechs, A. B. J. Kuijlaars

In this paper we study the asymptotic behavior of a family of polynomials which are orthogonal with respect to an exponential weight on certain contours of the complex plane. The zeros of these polynomials are the nodes for complex Gaussian quadrature of an oscillatory integral on the real axis with a high order stationary point, and their limit distribution is also analyzed. We show that the zeros accumulate along a contour in the complex plane that has the S-property in an external field. In addition, the strong asymptotics of the orthogonal polynomials is obtained by applying the nonlinear Deift--Zhou steepest descent method to the corresponding Riemann--Hilbert problem.