Search

We investigate theory and algorithms for pool-based active learning for multiclass classification using multinomial logistic regression. Using finite sample analysis, we prove that the Fisher Information Ratio (FIR) lower and upper bounds the excess risk. Based on our theoretical analysis, we propose an active learning algorithm that employs regret minimization to minimize the FIR. To verify our derived excess risk bounds, we conduct experiments on synthetic datasets. Furthermore, we compare FIRAL with five other methods and found that our scheme outperforms them: it consistently produces the smallest classification error in the multiclass logistic regression setting, as demonstrated through experiments on MNIST, CIFAR-10, and 50-class ImageNet., Comment: Accepted at the 37th Conference on Neural Information Processing Systems (NeurIPS 2023)

FIRAL is a recently proposed deterministic active learning algorithm for multiclass classification using logistic regression. It was shown to outperform the state-of-the-art in terms of accuracy and robustness and comes with theoretical performance guarantees. However, its scalability suffers when dealing with datasets featuring a large number of points $n$, dimensions $d$, and classes $c$, due to its $\mathcal{O}(c^2d^2+nc^2d)$ storage and $\mathcal{O}(c^3(nd^2 + bd^3 + bn))$ computational complexity where $b$ is the number of points to select in active learning. To address these challenges, we propose an approximate algorithm with storage requirements reduced to $\mathcal{O}(n(d+c) + cd^2)$ and a computational complexity of $\mathcal{O}(bncd^2)$. Additionally, we present a parallel implementation on GPUs. We demonstrate the accuracy and scalability of our approach using MNIST, CIFAR-10, Caltech101, and ImageNet. The accuracy tests reveal no deterioration in accuracy compared to FIRAL. We report strong and weak scaling tests on up to 12 GPUs, for three million point synthetic dataset., Comment: To be appeared at SC'24. Link: https://sc24.conference-program.com/presentation/?id=pap624&sess=sess397

We present a numerical method for the velocity-space, spatially homogeneous, collisional Boltzmann equation for electron transport in low-temperature plasma (LTP) conditions. Modeling LTP plasmas is useful in many applications, including advanced manufacturing, material processing, semiconductor processing, and hypersonics, to name a few. Most state-of-the-art methods for electron kinetics are based on Monte-Carlo sampling for collisions combined with Lagrangian particle-in-cell methods. We discuss an Eulerian solver that approximates the electron velocity distribution function using spherical harmonics (angular components) and B-splines (energy component). Our solver supports electron-heavy elastic and inelastic binary collisions, electron-electron Coulomb interactions, steady-state and transient dynamics, and an arbitrary nmber of angular terms in the electron distribution function. We report convergence results and compare our solver to two other codes: an in-house particle Monte-Carlo ethod; and Bolsig+, a state-of-the-art Eulerian solver for electron transport in LTPs. Furthermore, we use our solver to study the relaxation time scales of the higher-order anisotropic correction terms. Our code is open-source and provides an interface that allows coupling to multiphysics simulations of low-temperature plasmas.

We present a multi-species partial differential equation (PDE) model for tumor growth and a an algorithm for calibrating the model from magnetic resonance imaging (MRI) scans. The model is designed for glioblastoma (GBM) brain tumors. The modeled species correspond to proliferative, infiltrative, and necrotic tumor cells. The model calibration is formulated as an inverse problem and solved a PDE-constrained optimization method. The data that drives the calibration is derived by a single multi-parametric MRI image. This a typical clinical scenario for GBMs. The unknown parameters that need to be calibrated from data include ten scalar parameters and the infinite dimensional initial condition (IC) for proliferative tumor cells. This inverse problem is highly ill-posed as we try to calibrate a nonlinear dynamical system from data taken at a single time. To address this ill-posedness, we split the inversion into two stages. First we regularize the IC reconstruction by solving a single-species compressed sensing problem. Then, using the IC reconstruction, we invert for model parameters using a weighted regularization term. We construct the regularization term by using auxiliary 1D inverse problems. We apply our proposed scheme to clinical data. We compare our algorithm with single-species reconstruction and unregularized reconstructions. Our scheme enables the stable estimation of non-observable species and quantification of infiltrative tumor cells.

We propose a method that uses a two-species ordinary differential equation (ODE) model to characterize misfolded tau (or simply tau) protein spreading in Alzheimer's disease (AD) and calibrates it from clinical data. The unknown model parameters are the initial condition (IC) for tau and three scalar parameters representing the migration, proliferation, and clearance of tau proteins. Driven by imaging data, these parameters are estimated by formulating a constrained optimization problem with a sparsity regularization for the IC. This optimization problem is solved with a projection-based quasi-Newton algorithm. We investigate the sensitivity of our method to different algorithm parameters. We evaluate the performance of our method on both synthetic and clinical data. The latter comprises cases from the AD Neuroimaging Initiative (ADNI) datasets: 455 cognitively normal (CN), 212 mild cognitive impairment (MCI), and 45 AD subjects. We compare the performance of our approach to the commonly used Fisher-Kolmogorov (FK) model with a fixed IC at the entorhinal cortex (EC). Our method demonstrates an average improvement of 25.7% relative error compared to the FK model on the AD dataset. HFK also achieves an R-squared score of 0.664 for fitting AD data compared with 0.55 from FK model results under the same optimization scheme. Furthermore, for cases that have longitudinal data, we estimate a subject-specific AD onset time., Comment: 10 pages, 6 figures

We propose GrainGNN, a surrogate model for the evolution of polycrystalline grain structure under rapid solidification conditions in metal additive manufacturing. High fidelity simulations of solidification microstructures are typically performed using multicomponent partial differential equations (PDEs) with moving interfaces. The inherent randomness of the PDE initial conditions (grain seeds) necessitates ensemble simulations to predict microstructure statistics, e.g., grain size, aspect ratio, and crystallographic orientation. Currently such ensemble simulations are prohibitively expensive and surrogates are necessary. In GrainGNN, we use a dynamic graph to represent interface motion and topological changes due to grain coarsening. We use a reduced representation of the microstructure using hand-crafted features; we combine pattern finding and altering graph algorithms with two neural networks, a classifier (for topological changes) and a regressor (for interface motion). Both networks have an encoder-decoder architecture; the encoder has a multi-layer transformer long-short-term-memory architecture; the decoder is a single layer perceptron. We evaluate GrainGNN by comparing it to high-fidelity phase field simulations for in-distribution and out-of-distribution grain configurations for solidification under laser power bed fusion conditions. GrainGNN results in 80\%--90\% pointwise accuracy; and nearly identical distributions of scalar quantities of interest (QoI) between phase field and GrainGNN simulations compared using Kolmogorov-Smirnov test. GrainGNN's inference speedup (PyTorch on single x86 CPU) over a high-fidelity phase field simulation (CUDA on a single NVIDIA A100 GPU) is 150$\times$--2000$\times$ for 100-initial grain problem. Further, using GrainGNN, we model the formation of 11,600 grains in 220 seconds on a single CPU core.

Boundary value problems involving elliptic PDEs such as the Laplace and the Helmholtz equations are ubiquitous in mathematical physics and engineering. Many such problems can be alternatively formulated as integral equations that are mathematically more tractable. However, an integral-equation formulation poses a significant computational challenge: solving large dense linear systems that arise upon discretization. In cases where iterative methods converge rapidly, existing methods that draw on fast summation schemes such as the Fast Multipole Method are highly efficient and well-established. More recently, linear complexity direct solvers that sidestep convergence issues by directly computing an invertible factorization have been developed. However, storage and computation costs are high, which limits their ability to solve large-scale problems in practice. In this work, we introduce a distributed-memory parallel algorithm based on an existing direct solver named ``strong recursive skeletonization factorization.'' Specifically, we apply low-rank compression to certain off-diagonal matrix blocks in a way that minimizes computation and data movement. Compared to iterative algorithms, our method is particularly suitable for problems involving ill-conditioned matrices or multiple right-hand sides. Large-scale numerical experiments are presented to show the performance of our Julia implementation.

In this paper, we present a novel numerical scheme for simulating deformable and extensible capsules suspended in a Stokesian fluid. The main feature of our scheme is a partition-of-unity (POU) based representation of the surface that enables asymptotically faster computations compared to spherical-harmonics based representations. We use a boundary integral equation formulation to represent and discretize hydrodynamic interactions. The boundary integrals are weakly singular. We use the quadrature scheme based on the regularized Stokes kernels. We also use partition-of unity based finite differences that are required for the computational of interfacial forces. Given an N-point surface discretization, our numerical scheme has fourth-order accuracy and O(N) asymptotic complexity, which is an improvement over the O(N^2 log(N)) complexity of a spherical harmonics based spectral scheme that uses product-rule quadratures. We use GPU acceleration and demonstrate the ability of our code to simulate the complex shapes with high resolution. We study capsules that resist shear and tension and their dynamics in shear and Poiseuille flows. We demonstrate the convergence of the scheme and compare with the state of the art.

We study the performance of CLAIRE -- a diffeomorphic multi-node, multi-GPU image-registration algorithm, and software -- in large-scale biomedical imaging applications with billions of voxels. At such resolutions, most existing software packages for diffeomorphic image registration are prohibitively expensive. As a result, practitioners first significantly downsample the original images and then register them using existing tools. Our main contribution is an extensive analysis of the impact of downsampling on registration performance. We study this impact by comparing full-resolution registrations obtained with CLAIRE to lower-resolution registrations for synthetic and real-world imaging datasets. Our results suggest that registration at full resolution can yield a superior registration quality -- but not always. For example, downsampling a synthetic image from $1024^3$ to $256^3$ decreases the Dice coefficient from 92% to 79%. However, the differences are less pronounced for noisy or low-contrast high-resolution images. CLAIRE allows us not only to register images of clinically relevant size in a few seconds but also to register images at unprecedented resolution in a reasonable time. The highest resolution considered is CLARITY images of size $2816\times3016\times1162$. To the best of our knowledge, this is the first study on image registration quality at such resolutions., Comment: 32 pages, 9 tables, 8 figures

We use numerical simulations to study the dynamics of red blood cells (RBCs) in unconfined and confined Poiseuille flow. Previous numerical studies with 3D vesicles have indicated that the slipper shape observed in experiments at high capillary number can be attributed to the bistability due to the interplay of wall push and outward migration tendency at higher viscosity contrasts. In this paper, we study this outward migration and bistability using numerical simulations for 3D capsules and provide phase diagrams of RBC dynamics with and without viscosity contrast. We observe a bistability of slipper and croissants in confined Poiseuille flow with viscosity contrast as observed in experiments., Comment: 11 pages, 14 figures

The discretization of certain integral equations, e.g., the first-kind Fredholm equation of Laplace's equation, leads to symmetric positive-definite linear systems, where the coefficient matrix is dense and often ill-conditioned. We introduce a new preconditioner based on a novel overlapping domain decomposition that can be combined efficiently with fast direct solvers. Empirically, we observe that the condition number of the preconditioned system is $O(1)$, independent of the problem size. Our domain decomposition is designed so that we can construct approximate factorizations of the subproblems efficiently. In particular, we apply the recursive skeletonization algorithm to subproblems associated with every subdomain. We present numerical results on problem sizes up to $16\,384^2$ in 2D and $256^3$ in 3D, which were solved in less than 16 hours and three hours, respectively, on an Intel Xeon Platinum 8280M.

We propose a method for extracting physics-based biomarkers from a single multiparametric Magnetic Resonance Imaging (mpMRI) scan bearing a glioma tumor. We account for mass effect, the deformation of brain parenchyma due to the growing tumor, which on its own is an important radiographic feature but its automatic quantification remains an open problem. In particular, we calibrate a partial differential equation (PDE) tumor growth model that captures mass effect, parameterized by a single scalar parameter, tumor proliferation, migration, while localizing the tumor initiation site. The single-scan calibration problem is severely ill-posed because the precancerous, healthy, brain anatomy is unknown. To address the ill-posedness, we introduce an ensemble inversion scheme that uses a number of normal subject brain templates as proxies for the healthy precancer subject anatomy. We verify our solver on a synthetic dataset and perform a retrospective analysis on a clinical dataset of 216 glioblastoma (GBM) patients. We analyze the reconstructions using our calibrated biophysical model and demonstrate that our solver provides both global and local quantitative measures of tumor biophysics and mass effect. We further highlight the improved performance in model calibration through the inclusion of mass effect in tumor growth models -- including mass effect in the model leads to 10% increase in average dice coefficients for patients with significant mass effect. We further evaluate our model by introducing novel biophysics-based features and using them for survival analysis. Our preliminary analysis suggests that including such features can improve patient stratification and survival prediction., Comment: 10 pages, 3 supplementary pages

Current clinical decision-making in oncology relies on averages of large patient populations to both assess tumor status and treatment outcomes. However, cancers exhibit an inherent evolving heterogeneity that requires an individual approach based on rigorous and precise predictions of cancer growth and treatment response. To this end, we advocate the use of quantitative in vivo imaging data to calibrate mathematical models for the personalized forecasting of tumor development. In this chapter, we summarize the main data types available from both common and emerging in vivo medical imaging technologies, and how these data can be used to obtain patient-specific parameters for common mathematical models of cancer. We then outline computational methods designed to solve these models, thereby enabling their use for producing personalized tumor forecasts in silico, which, ultimately, can be used to not only predict response, but also optimize treatment. Finally, we discuss the main barriers to making the above paradigm a clinical reality.

We introduce a randomized algorithm, namely RCHOL, to construct an approximate Cholesky factorization for a given Laplacian matrix (a.k.a., graph Laplacian). From a graph perspective, the exact Cholesky factorization introduces a clique in the underlying graph after eliminating a row/column. By randomization, RCHOL only retains a sparse subset of the edges in the clique using a random sampling developed by Spielman and Kyng. We prove RCHOL is breakdown-free and apply it to solving large sparse linear systems with symmetric diagonally dominant matrices. In addition, we parallelize RCHOL based on the nested dissection ordering for shared-memory machines. We report numerical experiments that demonstrate the robustness and the scalability of RCHOL. For example, our parallel code scaled up to 64 threads on a single node for solving the 3D Poisson equation, discretized with the 7-point stencil on a $1024\times 1024 \times 1024$ grid, a problem that has one billion unknowns.

Modeling flow in geosystems with natural fault is a challenging problem due to low permeability of fault compared to its surrounding porous media. One way to predict the behavior of the flow while taking the effects of fault into account is to use the mixed finite element method. However, the mixed method could be time consuming due to large number of degree of freedom since both pressure and velocity are considered in the system. A new modeling method is presented in this paper. First, we introduce approximations of pressure based on the relation of pressure and velocity. We furthure decouple the approximated pressure from velocity so that it can be solved independently by continuous Galerkin finite element method. The new problem involves less degree of freedom than the mixed method for a given mesh . Moreover, local problem associated with a small subdomain around the fault is additionally solved to increase the accuracy of approximations around fault. Numerical experiments are conducted to examine the accuracy and efficiency of the new method. Results of three-dimensional tests show that our new method is up to 30$\times$ faster than the the mixed method at given $L^2$ pressure error., Comment: 24 pages

Clustering is a fundamental task in machine learning. One of the most successful and broadly used algorithms is DBSCAN, a density-based clustering algorithm. DBSCAN requires $\epsilon$-nearest neighbor graphs of the input dataset, which are computed with range-search algorithms and spatial data structures like KD-trees. Despite many efforts to design scalable implementations for DBSCAN, existing work is limited to low-dimensional datasets, as constructing $\epsilon$-nearest neighbor graphs is expensive in high-dimensions. In this paper, we modify DBSCAN to enable use of $\kappa$-nearest neighbor graphs of the input dataset. The $\kappa$-nearest neighbor graphs are constructed using approximate algorithms based on randomized projections. Although these algorithms can become inaccurate or expensive in high-dimensions, they possess a much lower memory overhead than constructing $\epsilon$-nearest neighbor graphs. We delineate the conditions under which $k$NN-DBSCAN produces the same clustering as DBSCAN. We also present an efficient parallel implementation of the overall algorithm using OpenMP for shared memory and MPI for distributed memory parallelism. We present results on up to 16 billion points in 20 dimensions, and perform weak and strong scaling studies using synthetic data. Our code is efficient in both low and high dimensions. We can cluster one billion points in 3D in less than one second on 28K cores on the Frontera system at the Texas Advanced Computing Center (TACC). In our largest run, we cluster 65 billion points in 20 dimensions in less than 40 seconds using 114,688 x86 cores on TACC's Frontera system. Also, we compare with a state of the art parallel DBSCAN code; on 20d/4M point dataset, our code is up to 37$\times$ faster., Comment: 26 pages, 10 figures

We present a Gauss-Newton-Krylov solver for large deformation diffeomorphic image registration. We extend the publicly available CLAIRE library to multi-node multi-graphics processing unit (GPUs) systems and introduce novel algorithmic modifications that significantly improve performance. Our contributions comprise ($i$) a new preconditioner for the reduced-space Gauss-Newton Hessian system, ($ii$) a highly-optimized multi-node multi-GPU implementation exploiting device direct communication for the main computational kernels (interpolation, high-order finite difference operators and Fast-Fourier-Transform), and ($iii$) a comparison with state-of-the-art CPU and GPU implementations. We solve a $256^3$-resolution image registration problem in five seconds on a single NVIDIA Tesla V100, with a performance speedup of 70% compared to the state-of-the-art. In our largest run, we register $2048^3$ resolution images (25 B unknowns; approximately 152$\times$ larger than the largest problem solved in state-of-the-art GPU implementations) on 64 nodes with 256 GPUs on TACC's Longhorn system., Comment: Proc ACM/IEEE Conference on Supercomputing 2020 (accepted for publication)

Aggregates of misfolded tau proteins (or just 'tau' for brevity) play a crucial role in the progression of Alzheimer's disease (AD) as they correlate with cell death and accelerated tissue atrophy. Longitudinal positron emission tomography (PET) scans can be used quantify the extend of abnormal tau spread. Such PET-based image biomarkers are a promising technology for AD diagnosis and prognosis. Here, we propose to calibrate an organ-scale biophysical mathematical model using longitudinal PET scans to extract characteristic growth patterns and spreading of tau. The biophysical model is a reaction-advection-diffusion partial differential equation (PDE) with only two scalar unknown parameters, one representing the spreading (the diffusion part of the PDE) and the other one the growth of tau (the reaction part of the PDE). The advection term captures tissue atrophy and is obtained from diffeomorphic registration of longitudinal magnetic resonance imaging (MRI) scans. We describe the method, present a numerical scheme for the calibration of the growth and spreading parameters, perform a sensitivity study using synthetic data, and we perform a preliminary evaluation on clinical scans from the ADNI dataset. We study whether such model calibration is possible and investigate the sensitivity of such calibration to the time between consecutive scans and the presence of atrophy. Our findings show that despite using only two calibration parameters, the model can reconstruct clinical scans quite accurately. We discovered that small time intervals between scans and the presence of background noise create difficulties. Our reconstructed model fits the data well, yet the study on clinical data also reveals shortcomings of the simplistic model. Interestingly, the parameters show significant variability across patients, an indication that these parameters could be useful biomarkers., Comment: 22 pages, 11 figures

We present a 3D fully-automatic method for the calibration of partial differential equation (PDE) models of glioblastoma (GBM) growth with mass effect, the deformation of brain tissue due to the tumor. We quantify the mass effect, tumor proliferation, tumor migration, and the localized tumor initial condition from a single multiparameteric Magnetic Resonance Imaging (mpMRI) patient scan. The PDE is a reaction-advection-diffusion partial differential equation coupled with linear elasticity equations to capture mass effect. The single-scan calibration model is notoriously difficult because the precancerous (healthy) brain anatomy is unknown. To solve this inherently ill-posed and ill-conditioned optimization problem, we introduce a novel inversion scheme that uses multiple brain atlases as proxies for the healthy precancer patient brain resulting in robust and reliable parameter estimation. We apply our method on both synthetic and clinical datasets representative of the heterogeneous spatial landscape typically observed in glioblastomas to demonstrate the validity and performance of our methods. In the synthetic data, we report calibration errors (due to the ill-posedness and our solution scheme) in the 10\%-20\% range. In the clinical data, we report good quantitative agreement with the observed tumor and qualitative agreement with the mass effect (for which we do not have a ground truth). Our method uses a minimal set of parameters and provides both global and local quantitative measures of tumor infiltration and mass effect., Comment: Provisionally accepted to MICCAI 2020

3D image registration is one of the most fundamental and computationally expensive operations in medical image analysis. Here, we present a mixed-precision, Gauss--Newton--Krylov solver for diffeomorphic registration of two images. Our work extends the publicly available CLAIRE library to GPU architectures. Despite the importance of image registration, only a few implementations of large deformation diffeomorphic registration packages support GPUs. Our contributions are new algorithms to significantly reduce the run time of the two main computational kernels in CLAIRE: calculation of derivatives and scattered-data interpolation. We deploy (i) highly-optimized, mixed-precision GPU-kernels for the evaluation of scattered-data interpolation, (ii) replace Fast-Fourier-Transform (FFT)-based first-order derivatives with optimized 8th-order finite differences, and (iii) compare with state-of-the-art CPU and GPU implementations. As a highlight, we demonstrate that we can register $256^3$ clinical images in less than 6 seconds on a single NVIDIA Tesla V100. This amounts to over 20$\times$ speed-up over the current version of CLAIRE and over 30$\times$ speed-up over existing GPU implementations., Comment: 20 pages, 6 figures, 8 tables

Central nervous system (CNS) tumors come with the vastly heterogeneous histologic, molecular and radiographic landscape, rendering their precise characterization challenging. The rapidly growing fields of biophysical modeling and radiomics have shown promise in better characterizing the molecular, spatial, and temporal heterogeneity of tumors. Integrative analysis of CNS tumors, including clinically-acquired multi-parametric magnetic resonance imaging (mpMRI) and the inverse problem of calibrating biophysical models to mpMRI data, assists in identifying macroscopic quantifiable tumor patterns of invasion and proliferation, potentially leading to improved (i) detection/segmentation of tumor sub-regions, and (ii) computer-aided diagnostic/prognostic/predictive modeling. This paper presents a summary of (i) biophysical growth modeling and simulation, (ii) inverse problems for model calibration, (iii) their integration with imaging workflows, and (iv) their application on clinically-relevant studies. We anticipate that such quantitative integrative analysis may even be beneficial in a future revision of the World Health Organization (WHO) classification for CNS tumors, ultimately improving patient survival prospects., Comment: 24 pages, 5 figures

Our objective is the calibration of mathematical tumor growth models from a single multiparametric scan. The target problem is the analysis of preoperative Glioblastoma (GBM) scans. To this end, we present a fully automatic tumor-growth calibration methodology that integrates a single-species reaction-diffusion partial differential equation (PDE) model for tumor progression with multiparametric Magnetic Resonance Imaging (mpMRI) scans to robustly extract patient specific biomarkers i.e., estimates for (i) the tumor cell proliferation rate, (ii) the tumor cell migration rate, and (iii) the original, localized site(s) of tumor initiation. Our method is based on a sparse reconstruction algorithm for the tumor initial location (TIL). This problem is particularly challenging due to nonlinearity, ill-posedeness, and ill conditioning. We propose a coarse-to-fine multi-resolution continuation scheme with parameter decomposition to stabilize the inversion. We demonstrate robustness and practicality of our method by applying the proposed method to clinical data of 206 GBM patients. We analyze the extracted biomarkers and relate tumor origin with patient overall survival by mapping the former into a common atlas space. We present preliminary results that suggest improved accuracy for prediction of patient overall survival when a set of imaging features is augmented with estimated biophysical parameters. All extracted features, tumor initial positions, and biophysical growth parameters are made publicly available for further analysis. To our knowledge, this is the first fully automatic scheme that can handle multifocal tumors and can localize the TIL to a few millimeters., Comment: 10 pages, 5 figures

In Part I (arXiv:1911.00619) of this article, we proposed an importance sampling algorithm to compute rare-event probabilities in forward uncertainty quantification problems. The algorithm, which we termed the "Bayesian Inverse Monte Carlo (BIMC) method", was shown to be optimal for problems in which the input-output operator is nearly linear. But applying the original BIMC to highly nonlinear systems can lead to several different failure modes. In this paper, we modify the BIMC method to extend its applicability to a wider class of systems. The modified algorithm, which we call "Adaptive-BIMC (A-BIMC)", has two stages. In the first stage, we solve a sequence of optimization problems to roughly identify those regions of parameter space which trigger the rare-event. In the second stage, we use the stage one results to construct a mixture of Gaussians that can be then used in an importance sampling algorithm to estimate rare event probability. We propose using a local surrogate that minimizes costly forward solves. The effectiveness of A-BIMC is demonstrated via several synthetic examples. Yet again, the modified algorithm is prone to failure. We systematically identify conditions under which it fails to lead to an effective importance sampling distribution., Comment: Added link to Part I in abstract

We consider the problem of estimating rare event probabilities, focusing on systems whose evolution is governed by differential equations with uncertain input parameters. If the system dynamics is expensive to compute, standard sampling algorithms such as the Monte Carlo method may require infeasible running times to accurately evaluate these probabilities. We propose an importance sampling scheme (which we call BIMC) that relies on solving an auxiliary, fictitious Bayesian inverse problem. The solution of the inverse problem yields a posterior PDF, a local Gaussian approximation to which serves as the importance sampling density. We apply BIMC to several problems and demonstrate that it can lead to computational savings of several orders of magnitude over the Monte Carlo method. We delineate conditions under which BIMC is optimal, as well as conditions when it can fail to yield an effective IS density.

We introduce a fast algorithm for entry-wise evaluation of the Gauss-Newton Hessian (GNH) matrix for the fully-connected feed-forward neural network. The algorithm has a precomputation step and a sampling step. While it generally requires $O(Nn)$ work to compute an entry (and the entire column) in the GNH matrix for a neural network with $N$ parameters and $n$ data points, our fast sampling algorithm reduces the cost to $O(n+d/\epsilon^2)$ work, where $d$ is the output dimension of the network and $\epsilon$ is a prescribed accuracy (independent of $N$). One application of our algorithm is constructing the hierarchical-matrix (H-matrix) approximation of the GNH matrix for solving linear systems and eigenvalue problems. It generally requires $O(N^2)$ memory and $O(N^3)$ work to store and factorize the GNH matrix, respectively. The H-matrix approximation requires only $O(N r_o)$ memory footprint and $O(N r_o^2)$ work to be factorized, where $r_o \ll N$ is the maximum rank of off-diagonal blocks in the GNH matrix. We demonstrate the performance of our fast algorithm and the H-matrix approximation on classification and autoencoder neural networks.

We use numerical simulations to study the dynamics of three dimensional vesicles in unconfined and confined Poiseuille flow. Previous numerical studies have shown that when the fluid viscosity inside and outside the vesicle is same (no viscosity contrast), a transition from asymmetric slippers to symmetric parachutes takes place as viscous forcing or capillary number is increased. At higher viscosity contrast, an outward migration tendency has also been observed in unconfined flow simulations. In this paper, we study how the presence of viscosity contrast and confining walls affect the dynamics of vesicles and present phase diagrams for confined Poiseuille flow with and without viscosity contrast. To our knowledge, this is the first study that provides a phase diagram for 3D vesicles with viscosity contrast in confined Poiseuille flow. The confining walls push the vesicle towards the center while the viscosity contrast has the opposite effect. This interplay leads to important differences in the dynamics like bistability at high capillary numbers., Comment: 11 pages, 20 figures

We present a novel formulation for the calibration of a biophysical tumor growth model from a single-time snapshot, MRI scan of a glioblastoma patient. Tumor growth models are typically nonlinear parabolic partial differential equations (PDEs). Thus, we have to generate a second snapshot to be able to extract significant information from a single patient snapshot. We create this two-snapshot scenario as follows. We use an atlas (an average of several scans of healthy individuals) as a substitute for an earlier, pretumor, MRI scan of the patient. Then, using the patient scan and the atlas, we combine image-registration algorithms and parameter estimation algorithms to achieve a better estimate of the healthy patient scan and the tumor growth parameters that are consistent with the data. Our scheme is based on our recent work (Scheufele et al, "Biophysically constrained diffeomorphic image registration, Tumor growth, Atlas registration, Adjoint-based methods, Parallel algorithms", CMAME, 2018), but apply a different and novel scheme where the tumor growth simulation in contrast to the previous work is executed in the patient brain domain and not in the atlas domain yielding more meaningful patient-specific results. As a basis, we use a PDE-constrained optimization framework. We derive a modified Picard-iteration-type solution strategy in which we alternate between registration and tumor parameter estimation in a new way. In addition, we consider an $\ell_1$ sparsity constraint on the initial condition for the tumor and integrate it with the new joint inversion scheme. We solve the subproblems with a reduced-space, inexact Gauss-Newton-Krylov/quasi-Newton methods. We present results using real brain data with synthetic tumor data that show that the new scheme reconstructs the tumor parameters in a more accurate and reliable way compared to our earlier scheme., Comment: 24 pages, 8 figures

We present a numerical scheme for solving an inverse problem for parameter estimation in tumor growth models for glioblastomas, a form of aggressive primary brain tumor. The growth model is a reaction-diffusion partial differential equation (PDE) for the tumor concentration. We use a PDE-constrained optimization formulation for the inverse problem. The unknown parameters are the reaction coefficient (proliferation), the diffusion coefficient (infiltration), and the initial condition field for the tumor PDE. Segmentation of Magnetic Resonance Imaging (MRI) scans from a single time snapshot drive the inverse problem where segmented tumor regions serve as partial observations of the tumor concentration. The precise time relative to tumor initiation is unknown, which poses an additional difficulty for inversion. We perform a frozen-coefficient spectral analysis and show that the inverse problem is severely ill-posed. We introduce a biophysically motivated regularization on the tumor initial condition. In particular, we assume that the tumor starts at a few locations (enforced with a sparsity constraint) and that the initial condition magnitude in the maximum norm equals one. We solve the resulting optimization problem using an inexact quasi-Newton method combined with a compressive sampling algorithm for the sparsity constraint. Our implementation uses PETSc and AccFFT libraries. We conduct numerical experiments on synthetic and clinical images to highlight the improved performance of our solver over an existing solver that uses a two-norm regularization for the calibration parameters. The existing solver is unable to localize the initial condition. Our new solver can localize the initial condition and recover infiltration and proliferation. In clinical datasets (for which the ground truth is unknown), our solver results in qualitatively different solutions compared to the existing solver., Comment: 24 pages, 8 figures

It has been observed that residual networks can be viewed as the explicit Euler discretization of an Ordinary Differential Equation (ODE). This observation motivated the introduction of so-called Neural ODEs, which allow more general discretization schemes with adaptive time stepping. Here, we propose ANODEV2, which is an extension of this approach that also allows evolution of the neural network parameters, in a coupled ODE-based formulation. The Neural ODE method introduced earlier is in fact a special case of this new more general framework. We present the formulation of ANODEV2, derive optimality conditions, and implement a coupled reaction-diffusion-advection version of this framework in PyTorch. We present empirical results using several different configurations of ANODEV2, testing them on multiple models on CIFAR-10. We report results showing that this coupled ODE-based framework is indeed trainable, and that it achieves higher accuracy, as compared to the baseline models as well as the recently-proposed Neural ODE approach.

Particulate Stokesian flows describe the hydrodynamics of rigid or deformable particles in Stokes flows. Due to highly nonlinear fluid-structure interaction dynamics, moving interfaces, and multiple scales, numerical simulations of such flows are challenging and expensive. In this Letter, we propose a generic machine-learning-augmented reduced model for these flows. Our model replaces expensive parts of a numerical scheme with multilayer perceptrons. Given the physical parameters of the particle, our model generalizes to arbitrary geometries and boundary conditions without the need to retrain the regression function. It is 10 times faster than a state-of-the-art numerical scheme having the same number of degrees of freedom and can reproduce several features of the flow quite accurately. We illustrate the performance of our model on integral equation formulation of vesicle suspensions in two dimensions.

Residual neural networks can be viewed as the forward Euler discretization of an Ordinary Differential Equation (ODE) with a unit time step. This has recently motivated researchers to explore other discretization approaches and train ODE based networks. However, an important challenge of neural ODEs is their prohibitive memory cost during gradient backpropogation. Recently a method proposed in [8], claimed that this memory overhead can be reduced from O(LN_t), where N_t is the number of time steps, down to O(L) by solving forward ODE backwards in time, where L is the depth of the network. However, we will show that this approach may lead to several problems: (i) it may be numerically unstable for ReLU/non-ReLU activations and general convolution operators, and (ii) the proposed optimize-then-discretize approach may lead to divergent training due to inconsistent gradients for small time step sizes. We discuss the underlying problems, and to address them we propose ANODE, an Adjoint based Neural ODE framework which avoids the numerical instability related problems noted above, and provides unconditionally accurate gradients. ANODE has a memory footprint of O(L) + O(N_t), with the same computational cost as reversing ODE solve. We furthermore, discuss a memory efficient algorithm which can further reduce this footprint with a trade-off of additional computational cost. We show results on Cifar-10/100 datasets using ResNet and SqueezeNext neural networks.

We propose domain decomposition preconditioners for the solution of an integral equation formulation of forward and inverse acoustic scattering problems with point scatterers. We study both forward and inverse problems and propose preconditioning techniques to accelerate the iterative solvers. For the forward scattering problem, we extend the domain decomposition based preconditioning techniques presented for partial differential equations in {\em "A restricted additive Schwarz preconditioner for general sparse linear systems", SIAM Journal on Scientific Computing, 21 (1999), pp. 792--797}, to integral equations. We combine this domain decomposition preconditioner with a low-rank correction, which is easy to construct, forming a new preconditioner. For the inverse scattering problem, we use the forward problem preconditioner as a building block for constructing a preconditioner for the Gauss-Newton Hessian. We present numerical results that demonstrate the performance of both preconditioning strategies.

Gliomas are the most common primary brain malignancies, with different degrees of aggressiveness, variable prognosis and various heterogeneous histologic sub-regions, i.e., peritumoral edematous/invaded tissue, necrotic core, active and non-enhancing core. This intrinsic heterogeneity is also portrayed in their radio-phenotype, as their sub-regions are depicted by varying intensity profiles disseminated across multi-parametric magnetic resonance imaging (mpMRI) scans, reflecting varying biological properties. Their heterogeneous shape, extent, and location are some of the factors that make these tumors difficult to resect, and in some cases inoperable. The amount of resected tumor is a factor also considered in longitudinal scans, when evaluating the apparent tumor for potential diagnosis of progression. Furthermore, there is mounting evidence that accurate segmentation of the various tumor sub-regions can offer the basis for quantitative image analysis towards prediction of patient overall survival. This study assesses the state-of-the-art machine learning (ML) methods used for brain tumor image analysis in mpMRI scans, during the last seven instances of the International Brain Tumor Segmentation (BraTS) challenge, i.e., 2012-2018. Specifically, we focus on i) evaluating segmentations of the various glioma sub-regions in pre-operative mpMRI scans, ii) assessing potential tumor progression by virtue of longitudinal growth of tumor sub-regions, beyond use of the RECIST/RANO criteria, and iii) predicting the overall survival from pre-operative mpMRI scans of patients that underwent gross total resection. Finally, we investigate the challenge of identifying the best ML algorithms for each of these tasks, considering that apart from being diverse on each instance of the challenge, the multi-institutional mpMRI BraTS dataset has also been a continuously evolving/growing dataset., Comment: The International Multimodal Brain Tumor Segmentation (BraTS) Challenge

In this article, we present a multispecies reaction-advection-diffusion partial differential equation (PDE) coupled with linear elasticity for modeling tumor growth. The model aims to capture the phenomenological features of glioblastoma multiforme observed in magnetic resonance imaging (MRI) scans. These include enhancing and necrotic tumor structures, brain edema and the so called "mass effect", that is, the deformation of brain tissue due to the presence of the tumor. The multispecies model accounts for proliferating, invasive and necrotic tumor cells as well as a simple model for nutrition consumption and tumor-induced brain edema. The coupling of the model with linear elasticity equations with variable coefficients allows us to capture the mechanical deformations due to the tumor growth on surrounding tissues. We present the overall formulation along with a novel operator-splitting scheme with components that include linearly-implicit preconditioned elliptic solvers, and semi-Lagrangian method for advection. Also, we present results showing simulated MRI images which highlight the capability of our method to capture the overall structure of glioblastomas in MRIs., Comment: J. Math. Biol. (2019)

We propose a segmentation framework that uses deep neural networks and introduce two innovations. First, we describe a biophysics-based domain adaptation method. Second, we propose an automatic method to segment white and gray matter, and cerebrospinal fluid, in addition to tumorous tissue. Regarding our first innovation, we use a domain adaptation framework that combines a novel multispecies biophysical tumor growth model with a generative adversarial model to create realistic looking synthetic multimodal MR images with known segmentation. Regarding our second innovation, we propose an automatic approach to enrich available segmentation data by computing the segmentation for healthy tissues. This segmentation, which is done using diffeomorphic image registration between the BraTS training data and a set of prelabeled atlases, provides more information for training and reduces the class imbalance problem. Our overall approach is not specific to any particular neural network and can be used in conjunction with existing solutions. We demonstrate the performance improvement using a 2D U-Net for the BraTS'18 segmentation challenge. Our biophysics based domain adaptation achieves better results, as compared to the existing state-of-the-art GAN model used to create synthetic data for training.

With this work, we release CLAIRE, a distributed-memory implementation of an effective solver for constrained large deformation diffeomorphic image registration problems in three dimensions. We consider an optimal control formulation. We invert for a stationary velocity field that parameterizes the deformation map. Our solver is based on a globalized, preconditioned, inexact reduced space Gauss--Newton--Krylov scheme. We exploit state-of-the-art techniques in scientific computing to develop an effective solver that scales to thousands of distributed memory nodes on high-end clusters. We present the formulation, discuss algorithmic features, describe the software package, and introduce an improved preconditioner for the reduced space Hessian to speed up the convergence of our solver. We test registration performance on synthetic and real data. We demonstrate registration accuracy on several neuroimaging datasets. We compare the performance of our scheme against different flavors of the Demons algorithm for diffeomorphic image registration. We study convergence of our preconditioner and our overall algorithm. We report scalability results on state-of-the-art supercomputing platforms. We demonstrate that we can solve registration problems for clinically relevant data sizes in two to four minutes on a standard compute node with 20 cores, attaining excellent data fidelity. With the present work we achieve a speedup of (on average) 5$\times$ with a peak performance of up to 17$\times$ compared to our former work., Comment: 37 pages

We solve a design optimization problem for deterministic lateral displacement (DLD) device to sort same-size biological cells by their deformability, in particular to sort red blood cells (RBCs) by their viscosity contrast between the fluid in the interior and the exterior of the cells. A DLD device optimized for efficient cell sorting enables rapid medical diagnoses of several diseases such as malaria since infected cells are stiffer than their healthy counterparts. The device consists of pillar arrays in which pillar rows are tilted and hence are not orthogonal to the columns. This arrangement leads cells to have different final vertical displacements depending on their deformability, therefore, it vertically separates the cells. Pillar cross section, tilt angle of the pillar rows and center-to-center distances between pillars define a unique device. For a given pair of viscosity contrast values of the cells we seek optimal DLD designs by fixing the tilt angle and the center-to-center distances. So the only design parameter is the pillar cross section which we parameterize with uniform 5th order B-splines. We propose an objective function to try to capture efficient cell sorting. The objective function is evaluated by simulating the cell flows through a device using our 2D model based on a boundary integral method (Kabacaoglu et al. Journal of Computational Physics, 357:43-77, 2018). We solve the optimization problem using the covariance matrix adaptation evolution strategy (CMA-ES), which is a stochastic, derivative-free algorithm. We present several scenarios where solving the optimization problem finds designs that can separate cells with similar viscosity contrast values. To the best of our knowledge, this is the first study which poses designing a DLD device as a constrained optimization problem and shows that solving this problem systematically discovers optimal designs., Comment: 28 pages

In this paper, we present algorithms for reconstructing an unknown compact scatterer embedded in a random noisy background medium, given measurements of the scattered field and information about the background medium and the sound profile. We present six different methods for the solution of this inverse problem using different amounts of scattered data and prior information about the random background medium and the scatterer. The different inversion algorithms are defined by a combination of stochastic programming methods and Bayesian formulation. Our basic results show that if we have data for just one instance of the random background medium the best strategy is to invert for both random medium and unknown scatterer with appropriate regularization. However, if we have data for multiple instances of the medium it may be worth solving a coupled set of multiple inverse problems. We present several numerical results for inverting for various scatterer geometries under different inversion scenarios. The main take-away of our study is that one should invert for both unknown scatterer and random medium, with appropriate, prior-information based regularization. Furthermore, if data from multiple realizations of the background medium is available, then it may be beneficial to combine results from multiple inversions., Comment: 40 pages

PDE-constrained optimization problems find many applications in medical image analysis, for example, neuroimaging, cardiovascular imaging, and oncological imaging. We review related literature and give examples on the formulation, discretization, and numerical solution of PDE-constrained optimization problems for medical imaging. We discuss three examples. The first one is image registration. The second one is data assimilation for brain tumor patients, and the third one data assimilation in cardiovascular imaging. The image registration problem is a classical task in medical image analysis and seeks to find pointwise correspondences between two or more images. The data assimilation problems use a PDE-constrained formulation to link a biophysical model to patient-specific data obtained from medical images. The associated optimality systems turn out to be sets of nonlinear, multicomponent PDEs that are challenging to solve in an efficient way. The ultimate goal of our work is the design of inversion methods that integrate complementary data, and rigorously follow mathematical and physical principles, in an attempt to support clinical decision making. This requires reliable, high-fidelity algorithms with a short time-to-solution. This task is complicated by model and data uncertainties, and by the fact that PDE-constrained optimization problems are ill-posed in nature, and in general yield high-dimensional, severely ill-conditioned systems after discretization. These features make regularization, effective preconditioners, and iterative solvers that, in many cases, have to be implemented on distributed-memory architectures to be practical, a prerequisite. We showcase state-of-the-art techniques in scientific computing to tackle these challenges.

We study deformability-based sorting of same-size RBCs via DLD using an in-house integral equation solver. Our goal is to quantitatively characterize the physical mechanisms that enable the cell separation. To this end, we systematically investigate the effects of the interior fluid viscosity and membrane elasticity of a cell on its behavior. In particular, we consider deep devices in which a cell can show rich dynamics such as taking a particular angular orientation depending on its mechanical property. We have found out that cells moving with a sufficiently high positive inclination angle with respect to the flow direction displace laterally while those with smaller angles travel with the flow streamlines. Thereby, deformability-based cell sorting is possible. The underlying mechanism here is cell migration due to the cell's positive inclination and the shear gradient. The higher the inclination is, the farther the cell can travel laterally. We also assess the efficiency of the technique for dense suspensions. It turns out that most of the cells in dense suspensions does not displace in the lateral direction no matter what their deformability is. As a result, separating cells using a DLD device becomes harder.

We present the SIBIA (Scalable Integrated Biophysics-based Image Analysis) framework for joint image registration and biophysical inversion and we apply it to analyse MR images of glioblastomas (primary brain tumors). In particular, we consider the following problem. Given the segmentation of a normal brain MRI and the segmentation of a cancer patient MRI, we wish to determine tumor growth parameters and a registration map so that if we "grow a tumor" (using our tumor model) in the normal segmented image and then register it to the patient segmented image, then the registration mismatch is as small as possible. We call this "the coupled problem" because it two-way couples the biophysical inversion and registration problems. In the image registration step we solve a large-deformation diffeomorphic registration problem parameterized by an Eulerian velocity field. In the biophysical inversion step we estimate parameters in a reaction-diffusion tumor growth model that is formulated as a partial differential equation. In this paper, our contributions are the introduction of a PDE-constrained optimization formulation of the coupled problem, the derivation of the optimality conditions, and the derivation of a Picard iterative scheme for the solution of the coupled problem. In addition, we perform several tests to experimentally assess the performance of our method on synthetic and clinical datasets. We demonstrate the convergence of the SIBIA optimization solver in different usage scenarios. We demonstrate that we can accurately solve the coupled problem in three dimensions ($256^3$ resolution) in a few minutes using 11 dual-x86 nodes. Also, we demonstrate that, with our coupled approach, we can successfully register normal MRI to tumor-bearing MRI while obtaining Dice coefficients that match those achieved when registering of normal-to-normal MRI., Comment: 38 pages, 19 figures

We present GOFMM (geometry-oblivious FMM), a novel method that creates a hierarchical low-rank approximation, "compression," of an arbitrary dense symmetric positive definite (SPD) matrix. For many applications, GOFMM enables an approximate matrix-vector multiplication in $N \log N$ or even $N$ time, where $N$ is the matrix size. Compression requires $N \log N$ storage and work. In general, our scheme belongs to the family of hierarchical matrix approximation methods. In particular, it generalizes the fast multipole method (FMM) to a purely algebraic setting by only requiring the ability to sample matrix entries. Neither geometric information (i.e., point coordinates) nor knowledge of how the matrix entries have been generated is required, thus the term "geometry-oblivious." Also, we introduce a shared-memory parallel scheme for hierarchical matrix computations that reduces synchronization barriers. We present results on the Intel Knights Landing and Haswell architectures, and on the NVIDIA Pascal architecture for a variety of matrices., Comment: 13 pages, accepted by SC'17

Vesicle suspensions appear in many biological and industrial applications. These suspensions are characterized by rich and complex dynamics of vesicles due to their interaction with the bulk fluid, and their large deformations and nonlinear elastic properties. Many existing state-of-the-art numerical schemes can resolve such complex vesicle flows. However, even when using provably optimal algorithms, these simulations can be computationally expensive, especially for suspensions with a large number of vesicles. In this paper, we investigate the effect of a number of algorithmic empirical fixes in an attempt to make low-resolution simulations more stable and more predictive. Based on our empirical studies for a number of flow configurations, we propose a scheme that attempts to integrate these fixes in a systematic way. This low-resolution scheme is an extension of our previous work Quaife and Biros (2014) and Quaife and Biros (2016). Our low-resolution correction algorithms (LRCA) include anti-aliasing and membrane reparametrization for avoiding spurious oscillations in vesicles' membranes, adaptive time stepping and a repulsion force for handling vesicle collisions and, correction of vesicles' area and arc-length for maintaining physical vesicle shapes. We perform a systematic error analysis by comparing the low-resolution simulations of dilute and dense suspensions with their high-fidelity, fully resolved, counterparts. We observe that the LRCA enables both efficient and statistically accurate low-resolution simulations of vesicle suspensions, while it can be 10x to 100x faster.

Kernel matrices appear in machine learning and non-parametric statistics. Given $N$ points in $d$ dimensions and a kernel function that requires $\mathcal{O}(d)$ work to evaluate, we present an $\mathcal{O}(dN\log N)$-work algorithm for the approximate factorization of a regularized kernel matrix, a common computational bottleneck in the training phase of a learning task. With this factorization, solving a linear system with a kernel matrix can be done with $\mathcal{O}(N\log N)$ work. Our algorithm only requires kernel evaluations and does not require that the kernel matrix admits an efficient global low rank approximation. Instead our factorization only assumes low-rank properties for the off-diagonal blocks under an appropriate row and column ordering. We also present a hybrid method that, when the factorization is prohibitively expensive, combines a partial factorization with iterative methods. As a highlight, we are able to approximately factorize a dense $11M\times11M$ kernel matrix in 2 minutes on 3,072 x86 "Haswell" cores and a $4.5M\times4.5M$ matrix in 1 minute using 4,352 "Knights Landing" cores., Comment: proceeding 31st IEEE International Parallel & Distributed Processing Symposium

Over the past two decades the field of computational science and engineering (CSE) has penetrated both basic and applied research in academia, industry, and laboratories to advance discovery, optimize systems, support decision-makers, and educate the scientific and engineering workforce. Informed by centuries of theory and experiment, CSE performs computational experiments to answer questions that neither theory nor experiment alone is equipped to answer. CSE provides scientists and engineers of all persuasions with algorithmic inventions and software systems that transcend disciplines and scales. Carried on a wave of digital technology, CSE brings the power of parallelism to bear on troves of data. Mathematics-based advanced computing has become a prevalent means of discovery and innovation in essentially all areas of science, engineering, technology, and society; and the CSE community is at the core of this transformation. However, a combination of disruptive developments---including the architectural complexity of extreme-scale computing, the data revolution that engulfs the planet, and the specialization required to follow the applications to new frontiers---is redefining the scope and reach of the CSE endeavor. This report describes the rapid expansion of CSE and the challenges to sustaining its bold advances. The report also presents strategies and directions for CSE research and education for the next decade., Comment: Major revision, to appear in SIAM Review

We study mixing in Stokesian vesicle suspensions in two dimensions on a cylindrical Couette apparatus using numerical simulations. The vesicle flow simulation is done using a boundary integral method and the advection-diffusion equation for the mixing of the solute is solved using a pseudo-spectral scheme. We study the effect of the area fraction, the viscosity contrast between the inside (the vesicles) and the outside (the bulk) fluid, the initial condition of the solute, and the mixing metric. We compare mixing in the suspension with mixing in the Couette apparatus without vesicles. On the one hand, the presence of vesicles in most cases, slightly suppresses mixing. This is because the solute can be only diffused across the vesicle interface and not advected. On the other hand, there exist spatial distributions of the solute for which the unperturbed Couette flow completely fails to mix whereas the presence of vesicles enables mixing. We derive a simple condition that relates the velocity and solute and can be used to characterize the cases in which the presence of vesicles promotes mixing., Comment: 25 pages, 14 figures

We present a parallel distributed-memory algorithm for large deformation diffeomorphic registration of volumetric images that produces large isochoric deformations (locally volume preserving). Image registration is a key technology in medical image analysis. Our algorithm uses a partial differential equation constrained optimal control formulation. Finding the optimal deformation map requires the solution of a highly nonlinear problem that involves pseudo-differential operators, biharmonic operators, and pure advection operators both forward and back- ward in time. A key issue is the time to solution, which poses the demand for efficient optimization methods as well as an effective utilization of high performance computing resources. To address this problem we use a preconditioned, inexact, Gauss-Newton- Krylov solver. Our algorithm integrates several components: a spectral discretization in space, a semi-Lagrangian formulation in time, analytic adjoints, different regularization functionals (including volume-preserving ones), a spectral preconditioner, a highly optimized distributed Fast Fourier Transform, and a cubic interpolation scheme for the semi-Lagrangian time-stepping. We demonstrate the scalability of our algorithm on images with resolution of up to $1024^3$ on the "Maverick" and "Stampede" systems at the Texas Advanced Computing Center (TACC). The critical problem in the medical imaging application domain is strong scaling, that is, solving registration problems of a moderate size of $256^3$---a typical resolution for medical images. We are able to solve the registration problem for images of this size in less than five seconds on 64 x86 nodes of TACC's "Maverick" system., Comment: accepted for publication at SC16 in Salt Lake City, Utah, USA; November 2016

We propose an efficient numerical algorithm for the solution of diffeomorphic image registration problems. We use a variational formulation constrained by a partial differential equation (PDE), where the constraints are a scalar transport equation. We use a pseudospectral discretization in space and second-order accurate semi-Lagrangian time stepping scheme for the transport equations. We solve for a stationary velocity field using a preconditioned, globalized, matrix-free Newton-Krylov scheme. We propose and test a two-level Hessian preconditioner. We consider two strategies for inverting the preconditioner on the coarse grid: a nested preconditioned conjugate gradient method (exact solve) and a nested Chebyshev iterative method (inexact solve) with a fixed number of iterations. We test the performance of our solver in different synthetic and real-world two-dimensional application scenarios. We study grid convergence and computational efficiency of our new scheme. We compare the performance of our solver against our initial implementation that uses the same spatial discretization but a standard, explicit, second-order Runge-Kutta scheme for the numerical time integration of the transport equations and a single-level preconditioner. Our improved scheme delivers significant speedups over our original implementation. As a highlight, we observe a 20$\times$ speedup for a two dimensional, real world multi-subject medical image registration problem.

We present a parallel algorithm for computing the approximate factorization of an $N$-by-$N$ kernel matrix. Once this factorization has been constructed (with $N \log^2 N $ work), we can solve linear systems with this matrix with $N \log N $ work. Kernel matrices represent pairwise interactions of points in metric spaces. They appear in machine learning, approximation theory, and computational physics. Kernel matrices are typically dense (matrix multiplication scales quadratically with $N$) and ill-conditioned (solves can require 100s of Krylov iterations). Thus, fast algorithms for matrix multiplication and factorization are critical for scalability. Recently we introduced ASKIT, a new method for approximating a kernel matrix that resembles N-body methods. Here we introduce INV-ASKIT, a factorization scheme based on ASKIT. We describe the new method, derive complexity estimates, and conduct an empirical study of its accuracy and scalability. We report results on real-world datasets including "COVTYPE" ($0.5$M points in 54 dimensions), "SUSY" ($4.5$M points in 8 dimensions) and "MNIST" (2M points in 784 dimensions) using shared and distributed memory parallelism. In our largest run we approximately factorize a dense matrix of size 32M $\times$ 32M (generated from points in 64 dimensions) on 4,096 Sandy-Bridge cores. To our knowledge these results improve the state of the art by several orders of magnitude., Comment: 11 pages, 2 figures, to appear in IPDPS 2016

We present a new library for parallel distributed Fast Fourier Transforms (FFT). The importance of FFT in science and engineering and the advances in high performance computing necessitate further improvements. AccFFT extends existing FFT libraries for CUDA-enabled Graphics Processing Units (GPUs) to distributed memory clusters. We use overlapping communication method to reduce the overhead of PCIe transfers from/to GPU. We present numerical results on the Maverick platform at the Texas Advanced Computing Center (TACC) and on the Titan system at the Oak Ridge National Laboratory (ORNL). We present the scaling of the library up to 4,096 K20 GPUs of Titan., Comment: Parallel FFT Library