Geometric Variational Inference

Efficiently accessing the information contained in non-linear and high dimensional probability distributions remains a core challenge in modern statistics. Traditionally, estimators that go beyond point estimates are either categorized as Variational Inference (VI) or Markov-Chain Monte-Carlo (MCMC) techniques. While MCMC methods that utilize the geometric properties of continuous probability distributions to increase their efficiency have been proposed, VI methods rarely use the geometry. This work aims to fill this gap and proposes geometric Variational Inference (geoVI), a method based on Riemannian geometry and the Fisher information metric. It is used to construct a coordinate transformation that relates the Riemannian manifold associated with the metric to Euclidean space. The distribution, expressed in the coordinate system induced by the transformation, takes a particularly simple form that allows for an accurate variational approximation by a normal distribution. Furthermore, the algorithmic structure allows for an efficient implementation of geoVI which is demonstrated on multiple examples, ranging from low-dimensional illustrative ones to non-linear, hierarchical Bayesian inverse problems in thousands of dimensions.     

High-order time-splitting Hermite and Fourier spectral methods

In this paper, we are concerned with the numerical solution of the time-dependent Gross–Pitaevskii Equation (GPE) involving a quasi-harmonic potential. Primarily, we consider discretisations that are based on spectral methods in space and higher-order exponential operator splitting methods in time. The resulting methods are favourable in view of accuracy and efficiency; moreover, geometric properties of the equation such as particle number and energy conservation are well captured.Regarding the spatial discretisation of the GPE, we consider two approaches. In the unbounded domain, we employ a spectral decomposition of the solution into Hermite basis functions; on the other hand, restricting the equation to a sufficiently large bounded domain, Fourier techniques are applicable. For the time integration of the GPE, we study various exponential operator splitting methods of convergence orders two, four, and six.Our main objective is to provide accuracy and efficiency comparisons of exponential operator splitting Fourier and Hermite pseudospectral methods for the time evolution of the GPE. Furthermore, we illustrate the effectiveness of higher-order time-splitting methods compared to standard integrators in a long-term integration.     

Compact integration factor methods for complex domains and adaptive mesh refinement

Implicit integration factor (IIF) method, a class of efficient semi-implicit temporal scheme, was introduced recently for stiff reaction–diffusion equations. To reduce cost of IIF, compact implicit integration factor (cIIF) method was later developed for efficient storage and calculation of exponential matrices associated with the diffusion operators in two and three spatial dimensions for Cartesian coordinates with regular meshes. Unlike IIF, cIIF cannot be directly extended to other curvilinear coordinates, such as polar and spherical coordinates, due to the compact representation for the diffusion terms in cIIF. In this paper, we present a method to generalize cIIF for other curvilinear coordinates through examples of polar and spherical coordinates. The new cIIF method in polar and spherical coordinates has similar computational efficiency and stability properties as the cIIF in Cartesian coordinate. In addition, we present a method for integrating cIIF with adaptive mesh refinement (AMR) to take advantage of the excellent stability condition for cIIF. Because the second order cIIF is unconditionally stable, it allows large time steps for AMR, unlike a typical explicit temporal scheme whose time step is severely restricted by the smallest mesh size in the entire spatial domain. Finally, we apply those methods to simulating a cell signaling system described by a system of stiff reaction–diffusion equations in both two and three spatial dimensions using AMR, curvilinear and Cartesian coordinates. Excellent performance of the new methods is observed.     

Time-dependent generalized polynomial chaos

Generalized polynomial chaos (gPC) has non-uniform convergence and tends to break down for long-time integration. The reason is that the probability density distribution (PDF) of the solution evolves as a function of time. The set of orthogonal polynomials associated with the initial distribution will therefore not be optimal at later times, thus causing the reduced efficiency of the method for long-time integration. Adaptation of the set of orthogonal polynomials with respect to the changing PDF removes the error with respect to long-time integration. In this method new stochastic variables and orthogonal polynomials are constructed as time progresses. In the new stochastic variable the solution can be represented exactly by linear functions. This allows the method to use only low order polynomial approximations with high accuracy. The method is illustrated with a simple decay model for which an analytic solution is available and subsequently applied to the three mode Kraichnan–Orszag problem with favorable results.     

On the symplectic integration of the discrete nonlinear Schrödinger equation with disorder

We present several methods, which utilize symplectic integration techniques based on two and three part operator splitting, for numerically solving the equations of motion of the disordered, discrete nonlinear Schrödinger (DDNLS) equation, and compare their efficiency. Our results suggest that the most suitable methods for the very long time integration of this one-dimensional Hamiltonian lattice model with many degrees of freedom (of the order of a few hundreds) are the ones based on three part splits of the system’s Hamiltonian. Two part split techniques can be preferred for relatively small lattices having up to N?≈?70 sites. An advantage of the latter methods is the better conservation of the system’s second integral, i.e. the wave packet’s norm.     

On the performance of some unwrapping algorithms

This work attempts an objective comparison in terms of performance and execution speed of some of the best known phase unwrapping algorithms. First the algorithms chosen, grouped into 4 classes (sequential methods, residues methods, global least square integration methods, others) are described. Then the influence of the weighting function on each of them is investigated. So as to obtain quantitative results it is necessary to use synthetic images, so that the exact solution is known, the technique for generating these images is also illustrated. Lastly, the algorithms performance in terms of influence of the weight functions, robustness and execution time are discussed.     

Krylov implicit integration factor methods for spatial discretization on high dimensional unstructured meshes: Application to discontinuous Galerkin methods

Integration factor methods are a class of “exactly linear part” time discretization methods. In [Q. Nie, Y.-T. Zhang, R. Zhao, Efficient semi-implicit schemes for stiff systems, Journal of Computational Physics, 214 (2006) 521–537], a class of efficient implicit integration factor (IIF) methods were developed for solving systems with both stiff linear and nonlinear terms, arising from spatial discretization of time-dependent partial differential equations (PDEs) with linear high order terms and stiff lower order nonlinear terms. The tremendous challenge in applying IIF temporal discretization for PDEs on high spatial dimensions is how to evaluate the matrix exponential operator efficiently. For spatial discretization on unstructured meshes to solve PDEs on complex geometrical domains, how to efficiently apply the IIF temporal discretization was open. In this paper, we solve this problem by applying the Krylov subspace approximations to the matrix exponential operator. Then we apply this novel time discretization technique to discontinuous Galerkin (DG) methods on unstructured meshes for solving reaction–diffusion equations. Numerical examples are shown to demonstrate the accuracy, efficiency and robustness of the method in resolving the stiffness of the DG spatial operator for reaction–diffusion PDEs. Application of the method to a mathematical model in pattern formation during zebrafish embryo development shall be shown.     

A new class of exponential propagation iterative methods of Runge–Kutta type (EPIRK)

We propose a new class of the exponential propagation iterative methods of Runge–Kutta-type (EPIRK). The EPIRK schemes are exponential integrators that can be competitive with explicit and implicit methods for integration of large stiff systems of ODEs. Introducing the new, more general than previously proposed, ansatz for EPIRK schemes allows for more flexibility in deriving computationally efficient high-order integrators. Recent extension of the theory of B-series to exponential integrators [1] is used to derive classical order conditions for schemes up to order five. An algorithm to systematically solve these conditions is presented and several new fifth order schemes are constructed. Several numerical examples are used to verify the order of the methods and to illustrate the performance of the new schemes.     

A multirate time integrator for regularized Stokeslets

The method of regularized Stokeslets is a numerical approach to approximating solutions of fluid–structure interaction problems in the Stokes regime. Regularized Stokeslets are fundamental solutions to the Stokes equations with a regularized point-force term that are used to represent forces generated by a rigid or elastic object interacting with the fluid. Due to the linearity of the Stokes equations, the velocity at any point in the fluid can be computed by summing the contributions of regularized Stokeslets, and the time evolution of positions can be computed using standard methods for ordinary differential equations. Rigid or elastic objects in the flow are usually treated as immersed boundaries represented by a collection of regularized Stokeslets coupled together by virtual springs which determine the forces exerted by the boundary in the fluid. For problems with boundaries modeled by springs with large spring constants, the resulting ordinary differential equations become stiff, and hence the time step for explicit time integration methods is severely constrained. Unfortunately, the use of standard implicit time integration methods for the method of regularized Stokeslets requires the solution of dense nonlinear systems of equations for many relevant problems. Here, an alternate strategy using an explicit multirate time integration scheme based on spectral deferred corrections is incorporated that in many cases can significantly decrease the computational cost of the method. The multirate methods are higher-order methods that treat different portions of the ODE explicitly with different time steps depending on the stiffness of each component. Numerical examples on two nontrivial three-dimensional problems demonstrate the increased efficiency of the multi-explicit approach with no significant increase in numerical error.     

Coupling discontinuous Galerkin methods and retarded potentials for transient wave propagation on unbounded domains

This work deals with the numerical simulation of wave propagation on unbounded domains with localized heterogeneities. To do so, we propose to combine a discretization based on a discontinuous Galerkin method in space and explicit finite differences in time on the regions containing heterogeneities with the retarded potential method to account the unbounded nature of the computational domain. The coupling formula enforces a discrete energy identity ensuring the stability under the usual CFL condition in the interior. Moreover, the scheme allows to use a smaller time step in the interior domain yielding to quasi-optimal discretization parameters for both methods. The aliasing phenomena introduced by the local time stepping are reduced by a post-processing by averaging in time obtaining a stable and second order consistent (in time) coupling algorithm. We compute the numerical rate of convergence of the method for an academic problem. The numerical results show the feasibility of the whole discretization procedure.     

Fourier spectral and wavelet solvers for the incompressible Navier–Stokes equations with volume-penalization: Convergence of a dipole–wall collision

In this study, we use volume-penalization to mimic the presence of obstacles in a flow or a domain with no-slip boundaries. This allows in principle the use of fast Fourier spectral methods and coherent vortex simulation techniques (based on wavelet decomposition of the flow variables) to compute turbulent wall-bounded flow or flows around solid obstacles by simply adding one term in the equation. Convergence checks are reported using a recently revived, and unexpectedly difficult dipole–wall collision as a benchmark computation. Several quantities, like the vorticity isolines, truncation error, kinetic energy and enstrophy are inspected for a collision of a dipole with a no-slip wall and compared with available benchmark data obtained with a standard Chebyshev pseudospectral method. We quantify the possible deteriorating effects of the Gibbs phenomenon present in the Fourier based schemes due to continuity restrictions of the penalized Navier–Stokes equations on the wall. It is found that Gibbs oscillations have a negligible effect on the flow evolution allowing higher-order recovery of the accuracy on a Fourier basis by means of postprocessing. An advantage of coherent vortex simulations, on the other hand, is that the degrees of freedom of the flow computation can strongly be reduced. In this study, we quantify the possible reduction of degrees of freedom while keeping the accuracy. For an optimal convergence scenario the penalization parameter has to scale with the number of Fourier and wavelet modes. In addition, an implicit treatment of the Darcy drag term in the penalized Navier–Stokes equations is beneficial since this allows one to set the time step independent from the penalization parameter without additional computational or memory requirements.     

Operator Splitting Implicit Integration Factor Methods for Stiff Reaction-Diffusion-Advection Systems

For reaction-diffusion-advection equations, the stiffness from the reaction and diffusion terms often requires very restricted time step size, while the nonlinear advection term may lead to a sharp gradient in localized spatial regions. It is challenging to design numerical methods that can efficiently handle both difficulties. For reaction-diffusion systems with both stiff reaction and diffusion terms, implicit integration factor (IIF) method and its higher dimensional analog compact IIF (cIIF) serve as an efficient class of time-stepping methods, and their second order version is linearly unconditionally stable. For nonlinear hyperbolic equations, weighted essentially non-oscillatory (WENO) methods are a class of schemes with a uniformly high-order of accuracy in smooth regions of the solution, which can also resolve the sharp gradient in an accurate and essentially non-oscillatory fashion. In this paper, we couple IIF/cIIF with WENO methods using the operator splitting approach to solve reaction-diffusion-advection equations. In particular, we apply the IIF/cIIF method to the stiff reaction and diffusion terms and the WENO method to the advection term in two different splitting sequences. Calculation of local truncation error and direct numerical simulations for both splitting approaches show the second order accuracy of the splitting method, and linear stability analysis and direct comparison with other approaches reveals excellent efficiency and stability properties. Applications of the splitting approach to two biological systems demonstrate that the overall method is accurate and efficient, and the splitting sequence consisting of two reaction-diffusion steps is more desirable than the one consisting of two advection steps, because CWC exhibits better accuracy and stability.     

Dual timestepping methods for detailed combustion chemistry

In this work, we develop and study several dual time integration methods for the solution of stiff, explosive differential equations governing combustion chemistry. Dual time integration is an implicit method wherein the sub-iteration process of each timestep is performed as a steady-state integration process, rather than the commonly used Newton–Raphson method. This allows stabilisation when nonlinear ignition events are contained within a timestep, providing considerable freedom in the choice of resolved phenomena. Timesteps may be chosen so as to resolve relatively long process timescales accurately rather than fast chemical timescales, something not possible with the common Newton's method. We illustrate this method using several backward difference formula methods and demonstrate the efficacy of our method in resolving low-frequency solutions of continuous flow stirred-tank reactors with periodic ignition–extinction events. We are able to step over ignition–extinction events with our stable, adaptive dual time method, and we study numerical convergence and error scaling on process timescales.     

Construction of fiber-optic bundle light-collection systems and calculations of collection efficiency

This paper describes the position-sensitive light-collection system that we use in our fast-beam laser experiments. The collection system consists of fiber-optic bundles whose facets are arranged to accept light emitted from a beam of fluorescent atoms. The flexibility of the fiber bundles allows their use in scanning collection systems with precise position sensitivity. We describe calculations of geometrical collection efficiency using a numerical integration scheme and compare the results with measurements. We also compare the collection efficiencies of the different fiber bundle arrangements that we used as our apparatus evolved with the implementation of various improvements.     

High-order algorithms for compressible reacting flow with complex chemistry

In this paper we describe a numerical algorithm for integrating the multicomponent, reacting, compressible Navier–Stokes equations, targeted for direct numerical simulation of combustion phenomena. The algorithm addresses two shortcomings of previous methods. First, it incorporates an eighth-order narrow stencil approximation of diffusive terms that reduces the communication compared to existing methods and removes the need to use a filtering algorithm to remove Nyquist frequency oscillations that are not damped with traditional approaches. The methodology also incorporates a multirate temporal integration strategy that provides an efficient mechanism for treating chemical mechanisms that are stiff relative to fluid dynamical time-scales. The overall methodology is eighth order in space with options for fourth order to eighth order in time. The implementation uses a hybrid programming model designed for effective utilisation of many-core architectures. We present numerical results demonstrating the convergence properties of the algorithm with realistic chemical kinetics and illustrating its performance characteristics. We also present a validation example showing that the algorithm matches detailed results obtained with an established low Mach number solver.     

Total-variation-diminishing implicit–explicit Runge–Kutta methods for the simulation of double-diffusive convection in astrophysics

We put forward the use of total-variation-diminishing (or more generally, strong stability preserving) implicit–explicit Runge–Kutta methods for the time integration of the equations of motion associated with the semiconvection problem in the simulation of stellar convection. The fully compressible Navier–Stokes equation, augmented by continuity and total energy equations, and an equation of state describing the relation between the thermodynamic quantities, is semi-discretized in space by essentially non-oscillatory schemes and dissipative finite difference methods. It is subsequently integrated in time by Runge–Kutta methods which are constructed such as to preserve the total variation diminishing (or strong stability) property satisfied by the spatial discretization coupled with the forward Euler method. We analyse the stability, accuracy and dissipativity of the time integrators and demonstrate that the most successful methods yield a substantial gain in computational efficiency as compared to classical explicit Runge–Kutta methods.     

TIME INTEGRATION OF NON-LINEAR DYNAMIC EQUATIONS BY MEANS OF A DIRECT VARIATIONAL METHOD

Non-linear dynamic problems governed by ordinary (ODE) or partial differential equations (PDE) are herein approached by means of an alternative methodology. A generalized solution named WEM by the authors and previously developed for boundary value problems, is applied to linear and non-linear equations. A simple transformation after selecting an arbitrary interval of interest T allows using WEM in initial conditions problems and others with both initial and boundary conditions. When dealing with the time variable, the methodology may be seen as a time integration scheme. The application of WEM leads to arbitrary precision results. It is shown that it lacks neither numerical damping nor chaos which were found to be present with the application of some of the time integration schemes most commonly used in standard finite element codes (e.g., methods of central difference, Newmark, Wilson-θ, and so on). Illustrations include the solution of two non-linear ODEs which govern the dynamics of a single-degree-of-freedom (s.d.o.f.) system of a mass and a spring with two different non-linear laws (cubic and hyperbolic tangent respectively). The third example is the application of WEM to the dynamic problem of a beam with non-linear springs at its ends and subjected to a dynamic load varying both in space and time, even with discontinuities, governed by a PDE. To handle systems of non-linear equations iterative algorithms are employed. The convergence of the iteration is achieved by takingn partitions of T. However, the values of T/n are, in general, several times larger than the usual Δt in other time integration techniques. The maximum error (measured as a percentage of the energy) is calculated for the first example and it is shown that WEM yields an acceptable level of errors even when rather large time steps are used.     

Multitaper Techniques and Filter Diagonalization Methods—A Comparison

In the present contribution we compare the new Multitaper Filtering technique with the very popular Filter Diagonalization Method. The substitution of a time-independent problem, like the standard Schrödinger equation, by a time-dependent one from the Filter Diagonalization Method allows the employment of and comparison with standard signal processing filtration machinery. The use of zero-order prolate spheroidal tapers as filtering functions is here extended and exactly formulated using techniques originating from general investigations of prolate spheroidal wave functions. We investigate the modifications presented with respect to accuracy and general effectiveness. The approach may be useful in various branches of physics and engineering sciences including signal processing applications as well as possibly also in general time-dependent processes.     

Optical time-domain analog pattern correlator for high-speed real-time image recognition

The speed of image processing is limited by image acquisition circuitry. While optical pattern recognition techniques can reduce the computational burden on digital image processing, their image correlation rates are typically low due to the use of spatial optical elements. Here we report a method that overcomes this limitation and enables fast real-time analog image recognition at a record correlation rate of 36.7 MHz--1000 times higher rates than conventional methods. This technique seamlessly performs image acquisition, correlation, and signal integration all optically in the time domain before analog-to-digital conversion by virtue of optical space-to-time mapping.     

On the Fully Implicit Solution of a Phase-Field Model for Binary Alloy Solidification in Three Dimensions

A fully implicit numerical method, based upon a combination of adaptively refined hierarchical meshes and geometric multigrid, is presented for the simulation of binary alloy solidification in three space dimensions. The computational techniques are presented for a particular mathematical model, based upon the phase-field approach, however, their applicability is of greater generality than for the specific phase-field model used here. In particular, an implicit second order time discretization is combined with the use of second order spatial differences to yield a large nonlinear system of algebraic equations as each time step. It is demonstrated that these equations may be solved reliably and efficiently through the use of a nonlinear multigrid scheme for locally refined grids. In effect, this paper presents an extension of earlier research in two space dimensions (J. Comput. Phys., 225 (2007), pp. 1271-1287) to fully three-dimensional problems. This extension is validated against earlier two-dimensional results and against some of the limited results available in three dimensions, obtained using an explicit scheme. The efficiency of the implicit approach and the multigrid solver are then demonstrated and some sample computational results for the simulation of the growth of dendrite structures are presented.