NANov 22, 2017
Fast and Stable Pascal Matrix AlgorithmsSamuel F. Potter, Ramani Duraiswami
In this paper, we derive a family of fast and stable algorithms for multiplying and inverting $n \times n$ Pascal matrices that run in $O(n log^2 n)$ time and are closely related to De Casteljau's algorithm for Bézier curve evaluation. These algorithms use a recursive factorization of the triangular Pascal matrices and improve upon the cripplingly unstable $O(n log n)$ fast Fourier transform-based algorithms which involve a Toeplitz matrix factorization. We conduct numerical experiments which establish the speed and stability of our algorithm, as well as the poor performance of the Toeplitz factorization algorithm. As an example, we show how our formulation relates to Bézier curve evaluation.
NAAug 29, 2018
Computing the quasipotential for nongradient SDEs in 3DShuo Yang, Samuel F. Potter, Maria K. Cameron
Nongradient SDEs with small white noise often arise when modeling biological and ecological time-irreversible processes. If the governing SDE were gradient, the maximum likelihood transition paths, transition rates, expected exit times, and the invariant probability distribution would be given in terms of its potential function. The quasipotential plays a similar role for nongradient SDEs. Unfortunately, the quasipotential is the solution of a functional minimization problem that can be obtained analytically only in some special cases. We propose a Dijkstra-like solver for computing the quasipotential on regular rectangular meshes in 3D. This solver results from a promotion and an upgrade of the previously introduced ordered line integral method with the midpoint quadrature rule for 2D SDEs. The key innovations that have allowed us to keep the CPU times reasonable while maintaining good accuracy are $(i)$ a new hierarchical update strategy, $(ii)$ the use of Karush-Kuhn-Tucker theory for rejecting unnecessary simplex updates, and $(iii)$ pruning the number of admissible simplexes and a fast search for them. An extensive numerical study is conducted on a series of linear and nonlinear examples where the quasipotential is analytically available or can be found at transition states by other methods. In particular, the proposed solver is applied to Tao's examples where the transition states are hyperbolic periodic orbits, and to a genetic switch model by Lv et al. (2014). The C source code implementing the proposed algorithm is available at M. Cameron's web page.
8.3NAMay 20
A Butterfly-Accelerated Manifold Harmonic TransformPaul G. Beckman, Samuel F. Potter, Michael O'Neil
The eigenfunctions of the Laplacian are a natural basis of functions for many tasks in computational mathematics. On the circle and sphere, the eigenfunctions are given by complex periodic exponentials and spherical harmonics, respectively, and much work has been done to develop fast algorithms for analyzing and synthesizing data in these bases. In this work, we generalize these special-case transforms to Laplace-Beltrami eigenfunctions of arbitrary surfaces, referred to as manifold harmonics. The resulting fast algorithm for computing linear combinations of the manifold harmonics is based on a butterfly factorization, which hierarchically compresses the transform matrix by constructing nested low-rank approximations of carefully selected submatrices. Several numerical examples are provided which demonstrate the speedups and reduction in memory requirements achieved by our algorithm for a variety of geometries, discretizations, and applications. In addition, a detailed analysis of the algorithm is given in the case that the underlying manifold is the flat periodic square.