Using FFT-based techniques in polynomial and matrix computations: recent advances and applicatons

Computations with univariate polynomials, like the evaluation of product, quotient, remainder, greatest common divisor, etc, are closely related to linear algebra computations performed with structured matrices having the Toeplitz-like or the Hankel-like structures.

The discrete Fourier transform, and the FFT algorithms for its computation, constitute a powerful tool for the design and analysis of fast algorithms for solving problems involving polynomials and structured matrices.

We recall the main correlations between polynomial and matrix computations and present two recent results in this field: in particular, we show how Fourier methods can speed up the solution of a wide class of problems arising in queueing theory where certain Markov chains, defined by infinite block Toeplitz matrices in generalized Hessenberg form, must be solved. Moreover, we present a new method for the numerical factorization of polynomials based on a matrix generalization of Koenig's theorem. This method, that relies on the evaluation/interpolation technique at the Fourier points, reduces the problem of polynomial factorization to the computation of the LU decomposition of a banded Toeplitz matrix with its rows and columns suitably permuted. Numerical experiments that show the effectiveness of our algorithms are presented     

Fast Fourier-Galerkin methods for solving singular boundary integral equations: Numerical integration and precondition

We develop a fast fully discrete Fourier-Galerkin method for solving a class of singular boundary integral equations. We prove that the number of multiplications used in generating the compressed matrix is O(nlog³n), and the solution of the proposed method preserves the optimal convergence order O(n^−t), where n is the order of the Fourier basis functions used in the method and t denotes the degree of regularity of the exact solution. Moreover, we propose a preconditioning which ensures the numerical stability when solving the preconditioned linear system. Numerical examples are presented to confirm the theoretical estimates and to demonstrate the approximation accuracy and computational efficiency of the proposed algorithm.     

A fast Fourier-collocation method for second boundary integral equations

In this paper we develop a fast collocation method for second boundary integral equations by the trigonometric polynomials. We propose a convenient way to compress the dense matrix representation of a compact integral operator with a smooth kernel under the Fourier basis and the corresponding collocation functionals. The compression leads to a sparse matrix with only O(nlog²n) number of nonzero entries, where 2n+1 denotes the order of the matrix. Thus we develop a fast Fourier-collocation method. We prove that the fast Fourier-collocation method gives the optimal convergence order up to a logarithmic factor. Moreover, we design a fast scheme for solving the corresponding truncated linear system. We establish that this algorithm preserves the quasi-optimal convergence of the approximate solution with requiring a number of O(nlog³n) multiplications.     

Algorithms for cardinal interpolation using box splines and radial basis functions

Summary We describe an algorithm for (bivariate) cardinal interpolation which can be applied to translates of basis functions which include box splines or radial basis functions. The algorithm is based on a representation of the Fourier transform of the fundamental interpolant, hence Fast Fourier Transform methods are available. In numerical tests the 4-directional box spline (transformed to the characteristical submodule of ²), the thin plate spline, and the multiquadric case give comparably equal and good results.     

Fast enclosure for solutions in underdetermined systems

Fast algorithms for enclosing minimal 2-norm solutions in underdetermined systems are proposed. For developing these algorithms, theory for computing error bounds for numerical solutions is established. Moreover techniques for accelerating the enclosure and obtaining smaller error bounds are introduced. Numerical results show the properties of the proposed algorithms.     

Fast transforms for tridiagonal linear equations

In this paper we study the use of the Fourier, Sine and Cosine Transform for solving or preconditioning linear systems, which arise from the discretization of elliptic problems. Recently, R. Chan and T. Chan considered circulant matrices for solving such systems. Instead of using circulant matrices, which are based on the Fourier Transform, we apply the Fourier and the Sine Transform directly. It is shown that tridiagonal matrices arising from the discretization of an onedimensional elliptic PDE are connected with circulant matrices by congruence transformations with the Fourier or the Sine matrix. Therefore, we can solve such linear systems directly, using only Fast Fourier Transforms and the Sherman-Morrison-Woodbury formula. The Fast Fourier Transform is highly parallelizable, and thus such an algorithm is interesting on a parallel computer. Moreover, similar relations hold between block tridiagonal matrices and Block Toeplitz-plus-Hankel matrices of ordern ²×n ² in the 2D case. This can be used to define in some sense natural approximations to the given matrix which lead to preconditioners for solving such linear systems.     

Solution of augmented linear systems using orthogonal factorizations

We consider solvingx+Ay=b andA ^T x=c for givenb, c andm ×n A of rankn. This is called the augmented system formulation (ASF) of two standard optimization problems, which include as special cases the minimum 2-norm of a linear underdetermined system (b=0) and the linear least squares problem (c=0), as well as more general problems. We examine the numerical stability of methods (for the ASF) based on the QR factorization ofA, whether by Householder transformations, Givens rotations, or the modified Gram-Schmidt (MGS) algorithm, and consider methods which useQ andR, or onlyR. We discuss the meaning of stability of algorithms for the ASF in terms of stability of algorithms for the underlying optimization problems.We prove the backward stability of several methods for the ASF which useQ andR, inclusing a new one based on MGS, and also show under what circumstances they may be regarded as strongly stable. We show why previous methods usingQ from MGS were not backward stable, but illustrate that some of these methods may be acceptable-error stable. We point out that the numerical accuracy of methods that do not useQ does not depend to any significant extent on which of of the above three QR factorizations is used. We then show that the standard methods which do not useQ are not backward stable or even acceptable-error stable for the general ASF problem, and discuss how iterative refinement can be used to counteract these deficiencies.Dedicated to Carl-Eric Fröberg on the occasion of his 75th birthdayThis research was partially supported by NSERC of Canada Grant No. A9236.     

A parallel method for fast and practical high-order newton interpolation

We present parallel algorithms for the computation and evaluation of interpolating polynomials. The algorithms use parallel prefix techniques for the calculation of divided differences in the Newton representation of the interpolating polynomial. Forn+1 given input pairs, the proposed interpolation algorithm requires only 2 [log(n+1)]+2 parallel arithmetic steps and circuit sizeO(n ²), reducing the best known circuit size for parallel interpolation by a factor of logn. The algorithm for the computation of the divided differences is shown to be numerically stable and does not require equidistant points, precomputation, or the fast Fourier transform. We report on numerical experiments comparing this with other serial and parallel algorithms. The experiments indicate that the method can be very useful for very high-order interpolation, which is made possible for special sets of interpolation nodes.Supported in part by the National Science Foundation under Grant No. NSF DCR-8603722.Supported by the National Science Foundation under Grants No. US NSF MIP-8410110, US NSF DCR85-09970, and US NSF CCR-8717942 and AT&T under Grant AT&T AFFL67Sameh.     

On the quadrature error in operational quadrature methods for convolutions

Summary We analyze the quadrature error associated with operational quadrature methods for convolution equations. The assumptions are that the convolution kernel is inL ¹ and that its Laplace transform is analytic and bounded in an obtuse sector of the complex plane. Under these circumstances the Laplace transform has a slow variation property which admits a Fourier analysis of the quadrature error. We provide generalL ^p error estimates assuming suitable smoothness conditions on the function under convolution.     

A fast transform for spherical harmonics

Spherical harmonics arise on the sphere S² in the same way that the (Fourier) exponential functions {e^ik}k arise on the circle. Spherical harmonic series have many of the same wonderful properties as Fourier series, but have lacked one important thing: a numerically stable fast transform analogous to the Fast Fourier Transform (FFT). Without a fast transform, evaluating (or expanding in) spherical harmonic series on the computer is slow—for large computations probibitively slow. This paper provides a fast transform.For a grid ofO(N²) points on the sphere, a direct calculation has computational complexityO(N⁴), but a simple separation of variables and FFT reduce it toO(N³) time. Here we present algorithms with timesO(N^5/2 log N) andO(N²(log N)²). The problem quickly reduces to the fast application of matrices of associated Legendre functions of certain orders. The essential insight is that although these matrices are dense and oscillatory, locally they can be represented efficiently in trigonometric series.     

Perturbation results for nearly uncoupled Markov chains with applications to iterative methods

Summary The standard perturbation theory for linear equations states that nearly uncoupled Markov chains (NUMCs) are very sensitive to small changes in the elements. Indeed, some algorithms, such as standard Gaussian elimination, will obtain poor results for such problems. A structured perturbation theory is given that shows that NUMCs usually lead to well conditioned problems. It is shown that with appropriate stopping, criteria, iterative aggregation/disaggregation algorithms will achieve these structured error bounds. A variant of Gaussian elimination due to Grassman, Taksar and Heyman was recently shown by O'Cinneide to achieve such bounds.Supported by the National Science Foundation under grant CCR-9000526 and its renewal, grant CCR-9201692. This research was done in part, during the author's visit to the Institute for Mathematics and its Applications, 514 Vincent Hall, 206 Church St. S.E., University of Minnesota, Minneapolis, MN 55455, USA     

Time discretization and quantization methods for optimal multiple switching problem

In this paper, we study probabilistic numerical methods based on optimal quantization algorithms for computing the solution to optimal multiple switching problems with regime-dependent state process. We first consider a discrete-time approximation of the optimal switching problem, and analyse its rate of convergence. Given a time step h$h$, the error is in general of order (hlog(1/h))^1/2${(hlog(1/h))}^{1/2}$, and of order h^1/2$h^{1/2}$ when the switching costs do not depend on the state process. We next propose quantization numerical schemes for the space discretization of the discrete-time Euler state process. A Markovian quantization approach relying on the optimal quantization of the normal distribution arising in the Euler scheme is analysed. In the particular case of uncontrolled state process, we describe an alternative marginal quantization method, which extends the recursive algorithm for optimal stopping problems as in Bally (2003) [1]. A priori L^p$L^p$-error estimates are stated in terms of quantization errors. Finally, some numerical tests are performed for an optimal switching problem with two regimes.     

Local error control inSDIRK-methods

This paper describes some problems that are encountered in the implementation of a class of Singly Diagonally Implicit Runge-Kutta (SDIRK) methods. The contribution to the local error from the local truncation error and the residual error from the algebraic systems involved are analysed. A section describes a special interpolation formula. This is used as a prediction stage in the iterative solution of the algebraic equations. A strategy for computing a starting stepsize is presented. The techniques are applied to numerical examples.     

On the uniform convergence of a collocation method for a class of singular integral equations

We establish the uniform convergence of a collocation method for solving a class of singular integral equations. This method uses the Jacobi polynomials {P _n ^{(, )} } as basis elements and the zeros of a Chebyshev polynomial of the first kind as collocation points. Uniform convergence is shown to hold under the weak assumption that the kernel and the right-hand side are Hölder-continous functions. Convergence rates are also given.     

Error growth analysis via stability regions for discretizations of initial value problems

This paper deals with numerical methods for the solution of linear initial value problems. Two main theorems are presented on the stability of these methods. Both theorems give conditions guaranteeing a mild error growth, for one-step methods characterized by a rational function ϕ(z). The conditions are related to the stability regionS={z:z∈ℂ with |ϕ(z)|≤1}, and can be viewed as variants to the resolvent condition occurring in the reputed Kreiss matrix theorem. Stability estimates are presented in terms of the number of time stepsn and the dimensions of the space. The first theorem gives a stability estimate which implies that errors in the numerical process cannot grow faster than linearly withs orn. It improves previous results in the literature where various restrictions were imposed onS and ϕ(z), including ϕ′(z)≠0 forz∈σS andS be bounded. The new theorem is not subject to any of these restrictions. The second theorem gives a sharper stability result under additional assumptions regarding the differential equation. This result implies that errors cannot grow faster thann ^β, with fixed β<1. The theory is illustrated in the numerical solution of an initial-boundary value problem for a partial differential equation, where the error growth is measured in the maximum norm.     

A Levin-type algorithm for accelerating the convergence of Fourier series

A nonlinear sequence transformation is presented which is able to accelerate the convergence of Fourier series. It is tailored to be exact for a certain model sequence. As in the case of the Levin transformation and other transformations of Levin-type, in this model sequence the partial sum of the series is written as the sum of the limit (or antilimit) and a certain remainder, i.e., it is of Levin-type. The remainder is assumed to be the product of a remainder estimate and the sum of the first terms oftwo Poincaré-type expansions which are premultiplied by two different phase factors. This occurrence of two phase factors is the essential difference to the Levin transformation. The model sequence for the new transformation may also be regarded as a special case of a model sequence based on several remainder estimates leading to the generalized Richardson extrapolation process introduced by Sidi. An algorithm for the recursive computation of the new transformation is presented. This algorithm can be implemented using only two one-dimensional arrays. It is proved that the sequence transformation is exact for Fourier series of geometric type which have coefficients proportional to the powers of a numberq, |q|<1. It is shown that under certain conditions the algorithm indeed accelerates convergence, and the order of the convergence is estimated. Finally, numerical test data are presented which show that in many cases the new sequence transformation is more powerful than Wynn's epsilon algorithm if the remainder estimates are properly chosen. However, it should be noted that in the vicinity of singularities of the Fourier series the new sequence transformation shows a larger tendency to numerical instability than the epsilon algorithm.     

Penultimate limiting forms in extreme value theory

Summary Let {X _n}_n≧1 be a sequence of independent, identically distributed random variables. If the distribution function (d.f.) ofM _n=max (X ₁,…,X _n), suitably normalized with attraction coefficients {α_n}_n≧1(α_n>0) and {b _n}_n≧1, converges to a non-degenerate d.f.G(x), asn→∞, it is of interest to study the rate of convergence to that limit law and if the convergence is slow, to find other d.f.'s which better approximate the d.f. of(M _n−b_n)/a_n thanG(x), for moderaten. We thus consider differences of the formF ⁿ(a_nx+b_n)−G(x), whereG(x) is a type I d.f. of largest values, i.e.,G(x)≡Λ(x)=exp (-exp(−x)), and show that for a broad class of d.f.'sF in the domain of attraction of Λ, there is a penultimate form of approximation which is a type II [Ф _α(x)=exp (−x^−α), x>0] or a type III [Ψ _α(x)= exp (−(−x)^α), x<0] d.f. of largest values, much closer toF ⁿ(a_nx+b_n) than the ultimate itself.     

A communication-less parallel algorithm for tridiagonal Toeplitz systems

Diagonally dominant tridiagonal Toeplitz systems of linear equations arise in many application areas and have been well studied in the past. Modern interest in numerical linear algebra is often focusing on solving classic problems in parallel. In McNally [Fast parallel algorithms for tri-diagonal symmetric Toeplitz systems, MCS Thesis, University of New Brunswick, Saint John, 1999], an m processor Split & Correct algorithm was presented for approximating the solution to a symmetric tridiagonal Toeplitz linear system of equations. Nemani [Perturbation methods for circulant-banded systems and their parallel implementation, Ph.D. Thesis, University of New Brunswick, Saint John, 2001] and McNally (2003) adapted the works of Rojo [A new method for solving symmetric circulant tri-diagonal system of linear equations, Comput. Math. Appl. 20 (1990) 61–67], Yan and Chung [A fast algorithm for solving special tri-diagonal systems, Computing 52 (1994) 203–211] and McNally et al. [A split-correct parallel algorithm for solving tri-diagonal symmetric Toeplitz systems, Internat. J. Comput. Math. 75 (2000) 303–313] to the non-symmetric case. In this paper we present relevant background from these methods and then introduce an m processor scalable communication-less approximation algorithm for solving a diagonally dominant tridiagonal Toeplitz system of linear equations.     

Application of fast spherical Fourier transform to density estimation

This paper is on density estimation on the 2-sphere, S², using the orthogonal series estimator corresponding to spherical harmonics. In the standard approach of truncating the Fourier series of the empirical density, the Fourier transform is replaced with a version of the discrete fast spherical Fourier transform, as developed by Driscoll and Healy. The fast transform only applies to quantitative data on a regular grid. We will apply a kernel operator to the empirical density, to produce a function whose values at the vertices of such a grid will be the basis for the density estimation. The proposed estimation procedure also contains a deconvolution step, in order to reduce the bias introduced by the initial kernel operator. The main issue is to find necessary conditions on the involved discretization and the bandwidth of the kernel operator, to preserve the rate of convergence that can be achieved by the usual computationally intensive Fourier transform. Density estimation is considered in L²(S²) and more generally in Sobolev spaces H_v(S²), any v?0, with the regularity assumption that the probability density to be estimated belongs to H_s(S²) for some s>v. The proposed technique to estimate the Fourier transform of an unknown density keeps computing cost down to order O(n), where n denotes the sample size.