An inexact Newton-Krylov algorithm for constrained diffeomorphic image registration

We propose numerical algorithms for solving large deformation diffeomorphic image registration problems. We formulate the nonrigid image registration problem as a problem of optimal control. This leads to an infinite-dimensional partial differential equation (PDE) constrained optimization problem. The PDE constraint consists, in its simplest form, of a hyperbolic transport equation for the evolution of the image intensity. The control variable is the velocity field. Tikhonov regularization on the control ensures well-posedness. We consider standard smoothness regularization based on $H^1$- or $H^2$-seminorms. We augment this regularization scheme with a constraint on the divergence of the velocity field rendering the deformation incompressible and thus ensuring that the determinant of the deformation gradient is equal to one, up to the numerical error. We use a Fourier pseudospectral discretization in space and a Chebyshev pseudospectral discretization in time. We use a preconditioned, globalized, matrix-free, inexact Newton-Krylov method for numerical optimization. A parameter continuation is designed to estimate an optimal regularization parameter. Regularity is ensured by controlling the geometric properties of the deformation field. Overall, we arrive at a black-box solver. We study spectral properties of the Hessian, grid convergence, numerical accuracy, computational efficiency, and deformation regularity of our scheme. We compare the designed Newton-Krylov methods with a globalized preconditioned gradient descent. We study the influence of a varying number of unknowns in time. The reported results demonstrate excellent numerical accuracy, guaranteed local deformation regularity, and computational efficiency with an optional control on local mass conservation. The Newton-Krylov methods clearly outperform the Picard method if high accuracy of the inversion is required.Comment: 32 pages; 10 figures; 9 table

A Perfect Match Condition for Point-Set Matching Problems Using the Optimal Mass Transport Approach

We study the performance of optimal mass transport--based methods applied to point-set matching problems. The present study, which is based on the L2 mass transport cost, states that perfect matches always occur when the product of the point-set cardinality and the norm of the curl of the nonrigid deformation field does not exceed some constant. This analytic result is justified by a numerical study of matching two sets of pulmonary vascular tree branch points whose displacement is caused by the lung volume changes in the same human subject. The nearly perfect match performance verifies the effectiveness of this mass transport--based approach.Read More: http://epubs.siam.org/doi/abs/10.1137/12086443

Learning and inference with Wasserstein metrics

Thesis: Ph. D., Massachusetts Institute of Technology, Department of Brain and Cognitive Sciences, 2018.Cataloged from PDF version of thesis.Includes bibliographical references (pages 131-143).This thesis develops new approaches for three problems in machine learning, using tools from the study of optimal transport (or Wasserstein) distances between probability distributions. Optimal transport distances capture an intuitive notion of similarity between distributions, by incorporating the underlying geometry of the domain of the distributions. Despite their intuitive appeal, optimal transport distances are often difficult to apply in practice, as computing them requires solving a costly optimization problem. In each setting studied here, we describe a numerical method that overcomes this computational bottleneck and enables scaling to real data. In the first part, we consider the problem of multi-output learning in the presence of a metric on the output domain. We develop a loss function that measures the Wasserstein distance between the prediction and ground truth, and describe an efficient learning algorithm based on entropic regularization of the optimal transport problem. We additionally propose a novel extension of the Wasserstein distance from probability measures to unnormalized measures, which is applicable in settings where the ground truth is not naturally expressed as a probability distribution. We show statistical learning bounds for both the Wasserstein loss and its unnormalized counterpart. The Wasserstein loss can encourage smoothness of the predictions with respect to a chosen metric on the output space. We demonstrate this property on a real-data image tagging problem, outperforming a baseline that doesn't use the metric. In the second part, we consider the probabilistic inference problem for diffusion processes. Such processes model a variety of stochastic phenomena and appear often in continuous-time state space models. Exact inference for diffusion processes is generally intractable. In this work, we describe a novel approximate inference method, which is based on a characterization of the diffusion as following a gradient flow in a space of probability densities endowed with a Wasserstein metric. Existing methods for computing this Wasserstein gradient flow rely on discretizing the underlying domain of the diffusion, prohibiting their application to problems in more than several dimensions. In the current work, we propose a novel algorithm for computing a Wasserstein gradient flow that operates directly in a space of continuous functions, free of any underlying mesh. We apply our approximate gradient flow to the problem of filtering a diffusion, showing superior performance where standard filters struggle. Finally, we study the ecological inference problem, which is that of reasoning from aggregate measurements of a population to inferences about the individual behaviors of its members. This problem arises often when dealing with data from economics and political sciences, such as when attempting to infer the demographic breakdown of votes for each political party, given only the aggregate demographic and vote counts separately. Ecological inference is generally ill-posed, and requires prior information to distinguish a unique solution. We propose a novel, general framework for ecological inference that allows for a variety of priors and enables efficient computation of the most probable solution. Unlike previous methods, which rely on Monte Carlo estimates of the posterior, our inference procedure uses an efficient fixed point iteration that is linearly convergent. Given suitable prior information, our method can achieve more accurate inferences than existing methods. We additionally explore a sampling algorithm for estimating credible regions.by Charles Frogner.Ph. D

Numerical Methods for Hamilton-Jacobi-Bellman Equations with Applications

Hamilton-Jacobi-Bellman (HJB) equations are nonlinear controlled partial differential equations (PDEs). In this thesis, we propose various numerical methods for HJB equations arising from three specific applications. First, we study numerical methods for the HJB equation coupled with a Kolmogorov-Fokker-Planck (KFP) equation arising from mean field games. In order to solve the nonlinear discretized systems efficiently, we propose a multigrid method. The main novelty of our approach is that we subtract artificial viscosity from the direct discretization coarse grid operators, such that the coarse grid error estimations are more accurate. The convergence rate of the proposed multigrid method is mesh-independent and faster than the existing methods in the literature. Next, we investigate numerical methods for the HJB formulation that arises from the mass transport image registration model. We convert the PDE of the model (a Monge-Ampère equation) to an equivalent HJB equation, propose a monotone mixed discretization, and prove that it is guaranteed to converge to the viscosity solution. Then we propose multigrid methods for the mixed discretization, where we set wide stencil points as coarse grid points, use injection at wide stencil points as the restriction, and achieve a mesh-independent convergence rate. Moreover, we propose a novel periodic boundary condition for the image registration PDE, such that when two images are related by a combination of a translation and a non-rigid deformation, the numerical scheme recovers the underlying transformation correctly. Finally, we propose a deep neural network framework for the HJB equations emerging from the study of American options in high dimensions. We convert the HJB equation to an equivalent Backward Stochastic Differential Equation (BSDE), introduce the least squares residual of the BSDE as the loss function, and propose a new neural network architecture that utilizes the domain knowledge of American options. Our proposed framework yields American option prices and deltas on the entire spacetime, not only at a given point. The computational cost of the proposed approach is quadratic in dimension, which addresses the curse of dimensionality issue that state-of-the-art approaches suffer

Coupling schemes and inexact Newton for multi-physics and coupled optimization problems

This work targets mathematical solutions and software for complex numerical simulation and optimization problems. Characteristics are the combination of different models and software modules and the need for massively parallel execution on supercomputers. We consider two different types of multi-component problems in Part I and Part II of the thesis: (i) Surface coupled fluid- structure interactions and (ii) analysis of medical MR imaging data of brain tumor patients. In (i), we establish highly accurate simulations by combining different aspects such as fluid flow and arterial wall deformation in hemodynamics simulations or fluid flow, heat transfer and mechanical stresses in cooling systems. For (ii), we focus on (a) facilitating the transfer of information such as functional brain regions from a statistical healthy atlas brain to the individual patient brain (which is topologically different due to the tumor), and (b) to allow for patient specific tumor progression simulations based on the estimation of biophysical parameters via inverse tumor growth simulation (given a single snapshot in time, only). Applications and specific characteristics of both problems are very distinct, yet both are hallmarked by strong inter-component relations and result in formidable, very large, coupled systems of partial differential equations. Part I targets robust and efficient quasi-Newton methods for black-box surface-coupling of parti- tioned fluid-structure interaction simulations. The partitioned approach allows for great flexibility and exchangeable of sub-components. However, breaking up multi-physics into single components requires advanced coupling strategies to ensure correct inter-component relations and effectively tackle instabilities. Due to the black-box paradigm, solver internals are hidden and information exchange is reduced to input/output relations. We develop advanced quasi-Newton methods that effectively establish the equation coupling of two (or more) solvers based on solving a non-linear fixed-point equation at the interface. Established state of the art methods fall short by either requiring costly tuning of problem dependent parameters, or becoming infeasible for large scale problems. In developing parameter-free, linear-complexity alternatives, we lift the robustness and parallel scalability of quasi-Newton methods for partitioned surface-coupled multi-physics simulations to a new level. The developed methods are implemented in the parallel, general purpose coupling tool preCICE. Part II targets MR image analysis of glioblastoma multiforme pathologies and patient specific simulation of brain tumor progression. We apply a joint medical image registration and biophysical inversion strategy, targeting at facilitating diagnosis, aiding and supporting surgical planning, and improving the efficacy of brain tumor therapy. We propose two problem formulations and decompose the resulting large-scale, highly non-linear and non-convex PDE-constrained optimization problem into two tightly coupled problems: inverse tumor simulation and medical image registration. We deduce a novel, modular Picard iteration-type solution strategy. We are the first to successfully solve the inverse tumor-growth problem based on a single patient snapshot with a gradient-based approach. We present the joint inversion framework SIBIA, which scales to very high image resolutions and parallel execution on tens of thousands of cores. We apply our methodology to synthetic and actual clinical data sets and achieve excellent normal-to-abnormal registration quality and present a proof of concept for a very promising strategy to obtain clinically relevant biophysical information. Advanced inexact-Newton methods are an essential tool for both parts. We connect the two parts by pointing out commonalities and differences of variants used in the two communities in unified notation