Olivier Verdier

Robert I McLachlan, Klas Modin, Hans Munthe-Kaas et al.

Butcher series appear when Runge-Kutta methods for ordinary differential equations are expanded in power series of the step size parameter. Each term in a Butcher series consists of a weighted elementary differential, and the set of all such differentials is isomorphic to the set of rooted trees, as noted by Cayley in the mid 19th century. A century later Butcher discovered that rooted trees can also be used to obtain the order conditions of Runge-Kutta methods, and he found a natural group structure, today known as the Butcher group. It is now known that many numerical methods also can be expanded in Butcher series; these are called B-series methods. A long-standing problem has been to characterize, in terms of qualitative features, all B-series methods. Here we tell the story of Butcher series, stretching from the early work of Cayley, to modern developments and connections to abstract algebra, and finally to the resolution of the characterization problem. This resolution introduces geometric tools and perspectives to an area traditionally explored using analysis and combinatorics.

Hans Munthe-Kaas, Olivier Verdier

We show that without other further assumption than affine equivariance and locality, a numerical integrator has an expansion in a generalized form of Butcher series (B-series) which we call aromatic B-series. We obtain an explicit description of aromatic B-series in terms of elementary differentials associated to aromatic trees, which are directed graphs generalizing trees. We also define a new class of integrators, the class of aromatic Runge-Kutta methods, that extends the class of Runge-Kutta methods, and have aromatic B-series expansion but are not B-series methods. Finally, those results are partially extended to the case of more general affine group equivariance.

Elena Celledoni, Bawfeh Kingsley Kometa, Olivier Verdier

We propose a class of semi-Lagrangian methods of high approximation order in space and time, based on spectral element space discretizations and exponential integrators of Runge-Kutta type. We discuss the extension of these methods to the Navier-Stokes equations, and their implementation using projections. Semi-Lagrangian methods up to order three are implemented and tested on various examples. The good performance of the methods for convection-dominated problems is demonstrated with numerical experiments.

Henrik Kalisch, Daulet Moldabayev, Olivier Verdier

In nonlinear dispersive evolution equations, the competing effects of nonlinearity and dispersion make a number of interesting phenomena possible. In the current work, the focus is on the numerical approximation of traveling-wave solutions of such equations. We describe our efforts to write a dedicated Python code which is able to compute traveling-wave solutions of nonlinear dispersive equations of the general form \begin{equation*} u_t + [f(u)]_{x} + \mathcal{L} u_x = 0, \end{equation*} where $\mathcal{L}$ is a self-adjoint operator, and $f$ is a real-valued function with $f(0) = 0$. The SpectraVVave code uses a continuation method coupled with a spectral projection to compute approximations of steady symmetric solutions of this equation. The code is used in a number of situations to gain an understanding of traveling-wave solutions. The first case is the Whitham equation, where numerical evidence points to the conclusion that the main bifurcation branch features three distinct points of interest, namely a turning point, a point of stability inversion, and a terminal point which corresponds to a cusped wave. The second case is the so-called modified Benjamin-Ono equation where the interaction of two solitary waves is investigated. It is found that is possible for two solitary waves to interact in such a way that the smaller wave is annihilated. The third case concerns the Benjamin equation which features two competing dispersive operators. In this case, it is found that bifurcation curves of periodic traveling-wave solutions may cross and connect high up on the branch in the nonlinear regime.

Hans Z. Munthe-Kaas, Ari Stern, Olivier Verdier

Motivated by numerical integration on manifolds, we relate the algebraic properties of invariant connections to their geometric properties. Using this perspective, we generalize some classical results of Cartan and Nomizu to invariant connections on algebroids. This has fundamental consequences for the theory of numerical integrators, giving a characterization of the spaces on which Butcher and Lie-Butcher series methods, which generalize Runge-Kutta methods, may be applied.

Hans Munthe-Kaas, Olivier Verdier

We consider numerical integrators of ODEs on homogeneous spaces (spheres, affine spaces, hyperbolic spaces). Homogeneous spaces are equipped with a built-in symmetry. A numerical integrator respects this symmetry if it is equivariant. One obtains homogeneous space integrators by combining a Lie group integrator with an isotropy choice. We show that equivariant isotropy choices combined with equivariant Lie group integrators produce equivariant homogeneous space integrators. Moreover, we show that the RKMK, Crouch--Grossman or commutator-free methods are equivariant. To show this, we give a novel description of Lie group integrators in terms of stage trees and motion maps, which unifies the known Lie group integrators. We then proceed to study the equivariant isotropy maps of order zero, which we call connections, and show that they can be identified with reductive structures and invariant principal connections. We give concrete formulas for connections in standard homogeneous spaces of interest, such as Stiefel, Grassmannian, isospectral, and polar decomposition manifolds. Finally, we show that the space of matrices of fixed rank possesses no connection.

James Benn, Stephen Marsland, Robert I McLachlan et al.

The nonlinear spaces of shapes (unparameterized immersed curves or submanifolds) are of interest for many applications in image analysis, such as the identification of shapes that are similar modulo the action of some group. In this paper we study a general representation of shapes that is based on linear spaces and is suitable for numerical discretization, being robust to noise. We develop the theory of currents for shape spaces by considering both the analytic and numerical aspects of the problem. In particular, we study the analytical properties of the current map and the $H^{-s}$ norm that it induces on shapes. We determine the conditions under which the current determines the shape. We then provide a finite element discretization of the currents that is a practical computational tool for shapes. Finally, we demonstrate this approach on a variety of examples.

Klas Modin, Olivier Verdier

We study a family of nonholonomic mechanical systems. These systems consist of harmonic oscillators coupled through nonholonomic constraints. In particular, the family includes the so called contact oscillator, which has been used as a test problem for numerical methods for nonholonomic mechanics. Furthermore, the systems under study constitute simple models for continuously variable transmission gearboxes. The main result is that each system in the family is integrable reversible with respect to the canonical reversibility map on the cotangent bundle. This result may explain previous numerical observations, that some discretisations of the contact oscillator have favourable structure preserving properties.

Robert I McLachlan, Klas Modin, Olivier Verdier et al.

We show that symplectic Runge-Kutta methods provide effective symplectic integrators for Hamiltonian systems with index one constraints. These include the Hamiltonian description of variational problems subject to position and velocity constraints nondegenerate in the velocities, such as those arising in sub-Riemannian geometry and control theory.

Geir Bogfjellmo, Klas Modin, Olivier Verdier

Cubic spline interpolation on Euclidean space is a standard topic in numerical analysis, with countless applications in science and technology. In several emerging fields, for example computer vision and quantum control, there is a growing need for spline interpolation on curved, non-Euclidean space. The generalization of cubic splines to manifolds is not self-evident, with several distinct approaches. One possibility is to mimic the acceleration minimizing property, which leads to Riemannian cubics. This, however, require the solution of a coupled set of non-linear boundary value problems that cannot be integrated explicitly, even if formulae for geodesics are available. Another possibility is to mimic De~Casteljau's algorithm, which leads to generalized Bézier curves. To construct C2-splines from such curves is a complicated non-linear problem, until now lacking numerical methods. Here we provide an iterative algorithm for C2-splines on Riemannian symmetric spaces, and we prove convergence of linear order. In terms of numerical tractability and computational efficiency, the new method surpasses those based on Riemannian cubics. Each iteration is parallel, thus suitable for multi-core implementation. We demonstrate the algorithm for three geometries of interest: the $n$-sphere, complex projective space, and the real Grassmannian.

David Cohen, Olivier Verdier

We present a new multi-symplectic formulation of constrained Hamiltonian partial differential equations, and we study the associated local conservation laws. A multi-symplectic discretisation based on this new formulation is exemplified by means of the Euler box scheme. When applied to the wave map equation, this numerical scheme is explicit, preserves the constraint and can be seen as a generalisation of the Shake algorithm for constrained mechanical systems. Furthermore, numerical experiments show excellent conservation properties of the numerical solutions.

Olivier Verdier

Many algorithms in numerical analysis are affine equivariant: they are immune to changes of affine coordinates. This is because those algorithms are defined using affine invariant constructions. There is, however, a crucial ingredient missing: most algorithms are in fact defined regardless of the underlying dimension. As a result, they are also invariant with respect to non-invertible affine transformation from spaces of different dimensions. We formulate this property precisely: these algorithms fall short of being natural transformations between affine functors. We give a precise definition of what we call a weak natural transformation between functors, and illustrate the point using examples coming from numerical analysis, in particular B-Series.

Justus Sagemüller, Olivier Verdier

Various types of saliency methods have been proposed for explaining black-box classification. In image applications, this means highlighting the part of the image that is most relevant for the current decision. Unfortunately, the different methods may disagree and it can be hard to quantify how representative and faithful the explanation really is. We observe however that several of these methods can be seen as edge cases of a single, more general procedure based on finding a particular path through the classifier's domain. This offers additional geometric interpretation to the existing methods. We demonstrate furthermore that ablation paths can be directly used as a technique of its own right. This is able to compete with literature methods on existing benchmarks, while giving more fine-grained information and better opportunities for validation of the explanations' faithfulness.

Camille Pouchol, Olivier Verdier

Linear inverse problems $A μ= δ$ with Poisson noise and non-negative unknown $μ\geq 0$ are ubiquitous in applications, for instance in Positron Emission Tomography (PET) in medical imaging. The associated maximum likelihood problem is routinely solved using an expectation-maximisation algorithm (ML-EM). This typically results in images which look spiky, even with early stopping. We give an explanation for this phenomenon. We first regard the image $μ$ as a measure. We prove that if the measurements $δ$ are not in the cone $\{A μ, μ\geq 0\}$, which is typical of short exposure times, likelihood maximisers as well as ML-EM cluster points must be sparse, i.e., typically a sum of point masses. On the other hand, in the long exposure regime, we prove that cluster points of ML-EM will be measures without singular part. Finally, we provide concentration bounds for the probability to be in the sparse case.

Ozan Öktem, Camille Pouchol, Olivier Verdier

Patient movement in emission tomography deteriorates reconstruction quality because of motion blur. Gating the data improves the situation somewhat: each gate contains a movement phase which is approximately stationary. A standard method is to use only the data from a few gates, with little movement between them. However, the corresponding loss of data entails an increase of noise. Motion correction algorithms have been implemented to take into account all the gated data, but they do not scale well, especially not in 3D. We propose a novel motion correction algorithm which addresses the scalability issue. Our approach is to combine an enhanced ML-EM algorithm with deep learning based movement registration. The training is unsupervised, and with artificial data. We expect this approach to scale very well to higher resolutions and to 3D, as the overall cost of our algorithm is only marginally greater than that of a standard ML-EM algorithm. We show that we can significantly decrease the noise corresponding to a limited number of gates.

Jonas Adler, Sebastian Lunz, Olivier Verdier et al.

The paper considers the problem of performing a task defined on a model parameter that is only observed indirectly through noisy data in an ill-posed inverse problem. A key aspect is to formalize the steps of reconstruction and task as appropriate estimators (non-randomized decision rules) in statistical estimation problems. The implementation makes use of (deep) neural networks to provide a differentiable parametrization of the family of estimators for both steps. These networks are combined and jointly trained against suitable supervised training data in order to minimize a joint differentiable loss function, resulting in an end-to-end task adapted reconstruction method. The suggested framework is generic, yet adaptable, with a plug-and-play structure for adjusting both the inverse problem and the task at hand. More precisely, the data model (forward operator and statistical model of the noise) associated with the inverse problem is exchangeable, e.g., by using neural network architecture given by a learned iterative method. Furthermore, any task that is encodable as a trainable neural network can be used. The approach is demonstrated on joint tomographic image reconstruction, classification and joint tomographic image reconstruction segmentation.

Robert I. McLachlan, Klas Modin, Olivier Verdier

Classical Hamiltonian spin systems are continuous dynamical systems on the symplectic phase space $(S^2)^n$. In this paper we investigate the underlying geometry of a time discretization scheme for classical Hamiltonian spin systems called the spherical midpoint method. As it turns out, this method displays a range of interesting geometrical features, that yield insights and sets out general strategies for geometric time discretizations of Hamiltonian systems on non-canonical symplectic manifolds. In particular, our study provides two new, completely geometric proofs that the discrete-time spin systems obtained by the spherical midpoint method preserve symplecticity. The study follows two paths. First, we introduce an extended version of the Hopf fibration to show that the spherical midpoint method can be seen as originating from the classical midpoint method on $T^*\mathbf{R}^{2n}$ for a collective Hamiltonian. Symplecticity is then a direct, geometric consequence. Second, we propose a new discretization scheme on Riemannian manifolds called the Riemannian midpoint method. We determine its properties with respect to isometries and Riemannian submersions and, as a special case, we show that the spherical midpoint method is of this type for a non-Euclidean metric. In combination with Kähler geometry, this provides another geometric proof of symplecticity.

Robert I. McLachlan, Klas Modin, Hans Munthe-Kaas et al.

Butcher series, also called B-series, are a type of expansion, fundamental in the analysis of numerical integration. Numerical methods that can be expanded in B-series are defined in all dimensions, so they correspond to \emph{sequences of maps}---one map for each dimension. A long-standing problem has been to characterise those sequences of maps that arise from B-series. This problem is solved here: we prove that a sequence of smooth maps between vector fields on affine spaces has a B-series expansion if and only if it is \emph{affine equivariant}, meaning it respects all affine maps between affine spaces.

Olivier Verdier, Huiyan Xue, Antonella Zanna

In earlier work, Lomeli and Meiss used a generalization of the symplectic approach to study volume preserving generating differential forms. In particular, for the $\mathbb{R}^3$ case, the first to differ from the symplectic case, they derived thirty-six one-forms that generate exact volume preserving maps. Xue and Zanna had studied these differential forms in connection with the numerical solution of divergence-free differential equations: can such forms be used to devise new volume preserving integrators or to further understand existing ones? As a partial answer to this question, Xue and Zanna showed how six of the generating volume form were naturally associated to consistent, first order, volume preserving numerical integrators. In this paper, we investigate and classify the remaining cases. The main result is the reduction of the thirty-six cases to five essentially different cases, up to variable relabeling and adjunction. We classify these five cases, identifying two novel classes and associating the other three to volume preserving vector fields under a Hamiltonian or Lagrangian representation. We demonstrate how these generating form lead to consistent volume preserving schemes for volume preserving vector fields in $\mathbb{R}^3$.

Robert I. McLachlan, Klas Modin, Olivier Verdier

We present a symplectic integrator, based on the canonical midpoint rule, for classical spin systems in which each spin is a unit vector in $\mathbb{R}^3$. Unlike splitting methods, it is defined for all Hamiltonians, and is $O(3)$-equivariant. It is a rare example of a generating function for symplectic maps of a noncanonical phase space. It yields an integrable discretization of the reduced motion of a free rigid body.

Olivier Verdier

We study problems associated with an operator pencil, i.e., a pair of operators on Banach spaces. Two natural problems to consider are linear constrained differential equations and the description of the generalized spectrum. The main tool to tackle either of those problems is the reduction of the pencil. There are two kinds of natural reduction operations associated to a pencil, which are conjugate to each other. Our main result is that those two kinds of reductions commute, under some mild assumptions that we investigate thoroughly. Each reduction exhibits moreover a pivot operator. The invertibility of all the pivot operators of all possible successive reductions corresponds to the notion of regular pencil in the finite dimensional case, and to the inf-sup condition for saddle point problems on Hilbert spaces. Finally, we show how to use the reduction and the pivot operators to describe the generalized spectrum of the pencil.

Robert I. McLachlan, Klas Modin, Olivier Verdier et al.

A geometric analysis of the Shake and Rattle methods for constrained Hamiltonian problems is carried out. The study reveals the underlying differential geometric foundation of the two methods, and the exact relation between them. In addition, the geometric insight naturally generalises Shake and Rattle to allow for a strictly larger class of constrained Hamiltonian systems than in the classical setting. In order for Shake and Rattle to be well defined, two basic assumptions are needed. First, a nondegeneracy assumption, which is a condition on the Hamiltonian, i.e., on the dynamics of the system. Second, a coisotropy assumption, which is a condition on the geometry of the constrained phase space. Non-trivial examples of systems fulfilling, and failing to fulfill, these assumptions are given.

Olivier Verdier

Matrix pencils, or pairs of matrices, may be used in a variety of applications. In particular, a pair of matrices (E,A) may be interpreted as the differential equation E x' + A x = 0. Such an equation is invariant by changes of variables, or linear combination of the equations. This change of variables or equations is associated to a group action. The invariants corresponding to this group action are well known, namely the Kronecker indices and divisors. Similarly, for another group action corresponding to the weak equivalence, a complete set of invariants is also known, among others the strangeness. We show how to define those invariants in a directly invariant fashion, i.e. without using a basis or an extra Euclidean structure. To this end, we will define a reduction process which produces a new system out of the original one. The various invariants may then be defined from operators related to the repeated application of the reduction process. We then show the relation between the invariants and the reduced subspace dimensions, and the relation with the regular pencil condition. This is all done using invariant tools only. Making special choices of basis then allows to construct the Kronecker canonical form. In a related manner, we construct the strangeness canonical form associated to weak equivalence.