Numerical computation of periodic orbits and isochrons for state-dependent delay perturbation of an ODE in the plane

We present algorithms and their implementation to compute limit cycles and their isochrons for state-dependent delay equations (SDDE's) which are perturbed from a planar differential equation with a limit cycle. Note that the space of solutions of an SDDE is infinite dimensional. We compute a two parameter family of solutions of the SDDE which converge to the solutions of the ODE as the perturbation goes to zero in a neighborhood of the limit cycle. The method we use formulates functional equations among periodic functions (or functions converging exponentially to periodic). The functional equations express that the functions solve the SDDE. Therefore, rather than evolving initial data and finding solutions of a certain shape, we consider spaces of functions with the desired shape and require that they are solutions. The mathematical theory of these invariance equations is developed in a companion paper, which develops "a posteriori" theorems. They show that, if there is a sufficiently approximate solution (with respect to some explicit condition numbers), then there is a true solution close to the approximate one. Since the numerical methods produce an approximate solution, and provide estimates of the condition numbers, we can make sure that the numerical solutions we consider approximate true solutions. In this paper, we choose a systematic way to approximate functions by a finite set of numbers (Taylor-Fourier series) and develop a toolkit of algorithms that implement the operators -- notably composition -- that enter into the theory. We also present several implementation results and present the results of running the algorithms and their implementation in some representative cases.Comment: 32 page

Metric Gaussian Variational Inference

Solving Bayesian inference problems approximately with variational approaches can provide fast and accurate results. Capturing correlation within the approximation requires an explicit parametrization. This intrinsically limits this approach to either moderately dimensional problems, or requiring the strongly simplifying mean-field approach. We propose Metric Gaussian Variational Inference (MGVI) as a method that goes beyond mean-field. Here correlations between all model parameters are taken into account, while still scaling linearly in computational time and memory. With this method we achieve higher accuracy and in many cases a significant speedup compared to traditional methods. MGVI is an iterative method that performs a series of Gaussian approximations to the posterior. We alternate between approximating the covariance with the inverse Fisher information metric evaluated at an intermediate mean estimate and optimizing the KL-divergence for the given covariance with respect to the mean. This procedure is iterated until the uncertainty estimate is self-consistent with the mean parameter. We achieve linear scaling by avoiding to store the covariance explicitly at any time. Instead we draw samples from the approximating distribution relying on an implicit representation and numerical schemes to approximately solve linear equations. Those samples are used to approximate the KL-divergence and its gradient. The usage of natural gradient descent allows for rapid convergence. Formulating the Bayesian model in standardized coordinates makes MGVI applicable to any inference problem with continuous parameters. We demonstrate the high accuracy of MGVI by comparing it to HMC and its fast convergence relative to other established methods in several examples. We investigate real-data applications, as well as synthetic examples of varying size and complexity and up to a million model parameters.Comment: Code is part of NIFTy5 release at https://gitlab.mpcdf.mpg.de/ift/NIFT

Semantics, Representations and Grammars for Deep Learning

Deep learning is currently the subject of intensive study. However, fundamental concepts such as representations are not formally defined -- researchers "know them when they see them" -- and there is no common language for describing and analyzing algorithms. This essay proposes an abstract framework that identifies the essential features of current practice and may provide a foundation for future developments. The backbone of almost all deep learning algorithms is backpropagation, which is simply a gradient computation distributed over a neural network. The main ingredients of the framework are thus, unsurprisingly: (i) game theory, to formalize distributed optimization; and (ii) communication protocols, to track the flow of zeroth and first-order information. The framework allows natural definitions of semantics (as the meaning encoded in functions), representations (as functions whose semantics is chosen to optimized a criterion) and grammars (as communication protocols equipped with first-order convergence guarantees). Much of the essay is spent discussing examples taken from the literature. The ultimate aim is to develop a graphical language for describing the structure of deep learning algorithms that backgrounds the details of the optimization procedure and foregrounds how the components interact. Inspiration is taken from probabilistic graphical models and factor graphs, which capture the essential structural features of multivariate distributions.Comment: 20 pages, many diagram

Operator learning approach for the limited view problem in photoacoustic tomography

In photoacoustic tomography, one is interested to recover the initial pressure distribution inside a tissue from the corresponding measurements of the induced acoustic wave on the boundary of a region enclosing the tissue. In the limited view problem, the wave boundary measurements are given on the part of the boundary, whereas in the full view problem, the measurements are known on the whole boundary. For the full view problem, there exist various fast and robust reconstruction methods. These methods give severe reconstruction artifacts when they are applied directly to the limited view data. One approach for reducing such artefacts is trying to extend the limited view data to the whole region boundary, and then use existing reconstruction methods for the full view data. In this paper, we propose an operator learning approach for constructing an operator that gives an approximate extension of the limited view data. We consider the behavior of a reconstruction formula on the extended limited view data that is given by our proposed approach. Approximation errors of our approach are analyzed. We also present numerical results with the proposed extension approach supporting our theoretical analysis

Non-Archimedean Ergodic Theory and Pseudorandom Generators

The paper develops techniques in order to construct computer programs, pseudorandom number generators (PRNG), that produce uniformly distributed sequences. The paper exploits an approach that treats standard processor instructions (arithmetic and bitwise logical ones) as continuous functions on the space of 2-adic integers. Within this approach, a PRNG is considered as a dynamical system and is studied by means of the non-Archimedean ergodic theory.Comment: Submitted to The Computer Journa

Projecting to a Slow Manifold: Singularly Perturbed Systems and Legacy Codes

The long-term dynamics of many dynamical systems evolve on an attracting, invariant "slow manifold" that can be parameterized by a few observable variables. Yet a simulation using the full model of the problem requires initial values for all variables. Given a set of values for the observables parameterizing the slow manifold, one needs a procedure for finding the additional values such that the state is close to the slow manifold to some desired accuracy. We consider problems whose solution has a singular perturbation expansion, although we do not know what it is nor have any way to compute it. We show in this paper that, under some conditions, computing the values of the remaining variables so that their (m+1)-st time derivatives are zero provides an estimate of the unknown variables that is an mth-order approximation to a point on the slow manifold in sense to be defined. We then show how this criterion can be applied approximately when the system is defined by a legacy code rather than directly through closed form equations.Comment: 21 pages, 1 figur

PBBFMM3D: a parallel black-box algorithm for kernel matrix-vector multiplication

We introduce \texttt{PBBFMM3D}, a parallel black-box method for computing kernel matrix-vector multiplication, where the underlying kernel is a non-oscillatory function in three dimensions. While a naive method requires \O(N^2) computation, \texttt{PBBFMM3D} reduces the cost to \O(N) work. In particular, our algorithm requires only the ability to evaluate the kernel function, and is thus a black-box method. To further accelerate the computation on shared-memory machines, a parallel algorithm is presented and implemented using \verb|OpenMP|, which achieved at most $19\times$ speedup on 32 cores in our numerical experiments. A real-world application in geostatistics is also presented, where \texttt{PBBFMM3D} is used in computing the truncated eigen-decomposition (a.k.a., principle component analysis) of a covariance matrix (a.k.a., graph Laplacian)

Discrete Symbol Calculus

This paper deals with efficient numerical representation and manipulation of differential and integral operators as symbols in phase-space, i.e., functions of space $x$ and frequency $\xi$. The symbol smoothness conditions obeyed by many operators in connection to smooth linear partial differential equations allow to write fast-converging, non-asymptotic expansions in adequate systems of rational Chebyshev functions or hierarchical splines. The classical results of closedness of such symbol classes under multiplication, inversion and taking the square root translate into practical iterative algorithms for realizing these operations directly in the proposed expansions. Because symbol-based numerical methods handle operators and not functions, their complexity depends on the desired resolution $N$ very weakly, typically only through $\log N$ factors. We present three applications to computational problems related to wave propagation: 1) preconditioning the Helmholtz equation, 2) decomposing wavefields into one-way components and 3) depth-stepping in reflection seismology.Comment: 32 page

Symbolic-numeric interface: A review

A survey of the use of a combination of symbolic and numerical calculations is presented. Symbolic calculations primarily refer to the computer processing of procedures from classical algebra, analysis, and calculus. Numerical calculations refer to both numerical mathematics research and scientific computation. This survey is intended to point out a large number of problem areas where a cooperation of symbolic and numerical methods is likely to bear many fruits. These areas include such classical operations as differentiation and integration, such diverse activities as function approximations and qualitative analysis, and such contemporary topics as finite element calculations and computation complexity. It is contended that other less obvious topics such as the fast Fourier transform, linear algebra, nonlinear analysis and error analysis would also benefit from a synergistic approach

Time-Limited Toeplitz Operators on Abelian Groups: Applications in Information Theory and Subspace Approximation

Toeplitz operators are fundamental and ubiquitous in signal processing and information theory as models for linear, time-invariant (LTI) systems. Due to the fact that any practical system can access only signals of finite duration, time-limited restrictions of Toeplitz operators are naturally of interest. To provide a unifying treatment of such systems working on different signal domains, we consider time-limited Toeplitz operators on locally compact abelian groups with the aid of the Fourier transform on these groups. In particular, we survey existing results concerning the relationship between the spectrum of a time-limited Toeplitz operator and the spectrum of the corresponding non-time-limited Toeplitz operator. We also develop new results specifically concerning the eigenvalues of time-frequency limiting operators on locally compact abelian groups. Applications of our unifying treatment are discussed in relation to channel capacity and in relation to representation and approximation of signals