A Lagrangian Gauss-Newton-Krylov Solver for Mass- and Intensity-Preserving Diffeomorphic Image Registration

We present an efficient solver for diffeomorphic image registration problems in the framework of Large Deformations Diffeomorphic Metric Mappings (LDDMM). We use an optimal control formulation, in which the velocity field of a hyperbolic PDE needs to be found such that the distance between the final state of the system (the transformed/transported template image) and the observation (the reference image) is minimized. Our solver supports both stationary and non-stationary (i.e., transient or time-dependent) velocity fields. As transformation models, we consider both the transport equation (assuming intensities are preserved during the deformation) and the continuity equation (assuming mass-preservation). We consider the reduced form of the optimal control problem and solve the resulting unconstrained optimization problem using a discretize-then-optimize approach. A key contribution is the elimination of the PDE constraint using a Lagrangian hyperbolic PDE solver. Lagrangian methods rely on the concept of characteristic curves that we approximate here using a fourth-order Runge-Kutta method. We also present an efficient algorithm for computing the derivatives of final state of the system with respect to the velocity field. This allows us to use fast Gauss-Newton based methods. We present quickly converging iterative linear solvers using spectral preconditioners that render the overall optimization efficient and scalable. Our method is embedded into the image registration framework FAIR and, thus, supports the most commonly used similarity measures and regularization functionals. We demonstrate the potential of our new approach using several synthetic and real world test problems with up to 14.7 million degrees of freedom.Comment: code available at: https://github.com/C4IR/FAIR.m/tree/master/add-ons/LagLDDM

A Multiscale Method for Model Order Reduction in PDE Parameter Estimation

Estimating parameters of Partial Differential Equations (PDEs) is of interest in a number of applications such as geophysical and medical imaging. Parameter estimation is commonly phrased as a PDE-constrained optimization problem that can be solved iteratively using gradient-based optimization. A computational bottleneck in such approaches is that the underlying PDEs needs to be solved numerous times before the model is reconstructed with sufficient accuracy. One way to reduce this computational burden is by using Model Order Reduction (MOR) techniques such as the Multiscale Finite Volume Method (MSFV). In this paper, we apply MSFV for solving high-dimensional parameter estimation problems. Given a finite volume discretization of the PDE on a fine mesh, the MSFV method reduces the problem size by computing a parameter-dependent projection onto a nested coarse mesh. A novelty in our work is the integration of MSFV into a PDE-constrained optimization framework, which updates the reduced space in each iteration. We also present a computationally tractable way of differentiating the MOR solution that acknowledges the change of basis. As we demonstrate in our numerical experiments, our method leads to computational savings particularly for large-scale parameter estimation problems and can benefit from parallelization.Comment: 22 pages, 4 figures, 3 table

Fireshape: a shape optimization toolbox for Firedrake

We introduce Fireshape, an open-source and automated shape optimization toolbox for the finite element software Firedrake. Fireshape is based on the moving mesh method and allows users with minimal shape optimization knowledge to tackle with ease challenging shape optimization problems constrained to partial differential equations (PDEs).Comment: 17 pages, 13 figures, 6 listing

Suitable Spaces for Shape Optimization

The differential-geometric structure of certain shape spaces is investigated and applied to the theory of shape optimization problems constrained by partial differential equations and variational inequalities. Furthermore, we define a diffeological structure on a new space of so-called $H^{1/2}$-shapes. This can be seen as a first step towards the formulation of optimization techniques on diffeological spaces. The $H^{1/2}$-shapes are a generalization of smooth shapes and arise naturally in shape optimization problems

Concave Quadratic Cuts for Mixed-Integer Quadratic Problems

The technique of semidefinite programming (SDP) relaxation can be used to obtain a nontrivial bound on the optimal value of a nonconvex quadratically constrained quadratic program (QCQP). We explore concave quadratic inequalities that hold for any vector in the integer lattice ${\bf Z}^n$, and show that adding these inequalities to a mixed-integer nonconvex QCQP can improve the SDP-based bound on the optimal value. This scheme is tested using several numerical problem instances of the max-cut problem and the integer least squares problem.Comment: 24 pages, 1 figur

Taylor approximation and variance reduction for PDE-constrained optimal control under uncertainty

In this work we develop a scalable computational framework for the solution of PDE-constrained optimal control under high-dimensional uncertainty. Specifically, we consider a mean-variance formulation of the control objective and employ a Taylor expansion with respect to the uncertain parameter either to directly approximate the control objective or as a control variate for variance reduction. The expressions for the mean and variance of the Taylor approximation are known analytically, although their evaluation requires efficient computation of the trace of the (preconditioned) Hessian of the control objective. We propose to estimate this trace by solving a generalized eigenvalue problem using a randomized algorithm that only requires the action of the Hessian on a small number of random directions. Then, the computational work does not depend on the nominal dimension of the uncertain parameter, but depends only on the effective dimension, thus ensuring scalability to high-dimensional problems. Moreover, to increase the estimation accuracy of the mean and variance of the control objective by the Taylor approximation, we use it as a control variate for variance reduction, which results in considerable computational savings (several orders of magnitude) compared to a plain Monte Carlo method. We demonstrate the accuracy, efficiency, and scalability of the proposed computational method for two examples with high-dimensional uncertain parameters: subsurface flow in a porous medium modeled as an elliptic PDE, and turbulent jet flow modeled by the Reynolds-averaged Navier--Stokes equations coupled with a nonlinear advection-diffusion equation characterizing model uncertainty. In particular, for the latter more challenging example we show scalability of our algorithm up to one million parameters resulting from discretization of the uncertain parameter field

Derivative-free optimization methods

In many optimization problems arising from scientific, engineering and artificial intelligence applications, objective and constraint functions are available only as the output of a black-box or simulation oracle that does not provide derivative information. Such settings necessitate the use of methods for derivative-free, or zeroth-order, optimization. We provide a review and perspectives on developments in these methods, with an emphasis on highlighting recent developments and on unifying treatment of such problems in the non-linear optimization and machine learning literature. We categorize methods based on assumed properties of the black-box functions, as well as features of the methods. We first overview the primary setting of deterministic methods applied to unconstrained, non-convex optimization problems where the objective function is defined by a deterministic black-box oracle. We then discuss developments in randomized methods, methods that assume some additional structure about the objective (including convexity, separability and general non-smooth compositions), methods for problems where the output of the black-box oracle is stochastic, and methods for handling different types of constraints

A shape optimization algorithm for cellular composites

We propose and investigate a mesh deformation technique for PDE constrained shape optimization. Introducing a gradient penalization to the inner product for linearized shape spaces, mesh degeneration can be prevented within the optimization iteration allowing for the scalability of employed solvers. We illustrate the approach by a shape optimization for cellular composites with respect to linear elastic energy under tension. The influence of the gradient penalization is evaluated and the parallel scalability of the approach demonstrated employing a geometric multigrid solver on hierarchically distributed meshes

On the optimal control of some nonsmooth distributed parameter systems arising in mechanics

Variational inequalities are an important mathematical tool for modelling free boundary problems that arise in different application areas. Due to the intricate nonsmooth structure of the resulting models, their analysis and optimization is a difficult task that has drawn the attention of researchers for several decades. In this paper we focus on a class of variational inequalities, called of the second kind, with a twofold purpose. First, we aim at giving a glance at some of the most prominent applications of these types of variational inequalities in mechanics, and the related analytical and numerical difficulties. Second, we consider optimal control problems constrained by these variational inequalities and provide a thorough discussion on the existence of Lagrange multipliers and the different types of optimality systems that can be derived for the characterization of local minima. The article ends with a discussion of the main challenges and future perspectives of this important problem class

Preconditioning Parametrized Linear Systems

Preconditioners are generally essential for fast convergence in the iterative solution of linear systems of equations. However, the computation of a good preconditioner can be expensive. So, while solving a sequence of many linear systems, it is advantageous to recycle preconditioners, that is, update a previous preconditioner and reuse the updated version. In this paper, we introduce a simple and effective method for doing this. Although our approach can be used for matrices changing slowly in any way, we focus on the important case of sequences of the type $(s_k\textbf{E}(\textbf{p}) + \textbf{A}(\textbf{p}))\textbf{x}_k = \textbf{b}_k$, where the right hand side may or may not change. More general changes in matrices will be discussed in a future paper. We update preconditioners by defining a map from a new matrix to a previous matrix, for example the first matrix in the sequence, and combine the preconditioner for this previous matrix with the map to define the new preconditioner. This approach has several advantages. The update is entirely independent from the original preconditioner, so it can be applied to any preconditioner. The possibly high cost of an initial preconditioner can be amortized over many linear solves. The cost of updating the preconditioner is more or less constant and independent of the original preconditioner. There is flexibility in balancing the quality of the map with the computational cost. In the numerical experiments section we demonstrate good results for several applications.Comment: V2 Model Reduction discussion. V3 Theoretical section replaced with experimental analysis of sparsity patterns. Top opt application added. Rail replaced with Flow. Only early shifts for THT. ILUTP, SAM implementations more efficient; results updated. ILUTP m-file included. New citations added. V4 New pattern added to top opt sparsity pattern analysis/results. ILUTP m-file minor update