BZPEERPapers, people, evidence and research paths
Menu

Discover

SearchFind papers, people, topics, institutions and opportunities from one place.GraphSee how papers, authors, topics and evidence connect.

Research tools

HomeContinue from the papers, people and topics that matter now.CollectionsTurn saved papers into reading lists, evidence maps and shareable context.ResearchReconstruct the path that led you from search to paper to question.CollaborationMove from a paper, topic or expert request into a focused thread.PublishCreate full-text research articles with metadata, formulas and versions.

Network

ExpertsFind specialists whose work and reviews explain why they can help.CommunitiesJoin research groups organized around papers, topics and shared evidence.

Opportunities

ListingsBrowse roles, calls and collaborations connected to real research context.

Account

Log inReturn to your saved papers, graph views and active threads.Create accountStart saving papers, building a research trail and joining teams.
Log inCreate account

Source author record

George Biros

George Biros appears in the imported research catalog. Authorship, coauthor and topic links are available while profile ownership is still unclaimed.

ResearcherUnclaimed source record
math.NANumerical AnalysisQuantitative Methodsphysics.med-phComputer Visioneess.IVmath.OCComputational Engineering, Finance, and Sciencecond-mat.softData Structures and AlgorithmsDistributed, Parallel, and Cluster ComputingMachine Learningphysics.flu-dynTissues and Organsmath.DSMathematical Software

Catalog footprint

What is connected

26works
16topics
4close collaborators

Actions

Connect this record

Log in to claim
Open graphBrowse works

Research graph

See the researcher in context

Open full explorer

Inspect adjacent papers, topics, institutions and collaborators without losing the researcher page.

Building this map preview

BZPEER is loading the nearby papers, people, topics and institutions for this page.

Published work

26 published item(s)

preprint2026arXiv

Inverse problems for history-enriched linear model reduction

Standard projection-based model reduction for dynamical systems incurs closure error because it only accounts for instantaneous dependence on the resolved state. From the Mori-Zwanzig (MZ) perspective, projecting the full dynamics onto a low-dimensional resolved subspace induces additional noise and memory terms arising from the dynamics of the unresolved component in the orthogonal complement. The memory term makes the resolved dynamics explicitly history dependent. In this work, based on the MZ identity, we derive exact, history-enriched models for the resolved dynamics of linear driven dynamical systems and formulate inverse problems to learn model operators from discrete snapshot data via least-squares regression. We propose a greedy time-marching scheme to solve the inverse problems efficiently and analyze operator identifiability under full and partial observation data availability. For full observation data, we show that, under mild assumptions, the operators are identifiable even when the full-state dynamics are governed by a general time-varying linear operator, whereas with partial observation data the inverse problem has a unique solution only when the full-state operator is time-invariant. To address the resulting non-uniqueness in the time-varying case, we introduce a time-smoothing Tikhonov regularization. Numerical results demonstrate that the operators can be faithfully reconstructed from both full and partial observation data and that the learned history-enriched MZ models yield accurate trajectories of the resolved state.

0 citations0 reviewsNo code linkedOpen access
preprint2026arXiv

LNODE: latent dynamics reveal the shared spatiotemporal structure of amyloid-$β$ progression

We introduce LNODE, a mechanism-based phenomenological model for amyloid beta (A$β$) dynamics, calibrated using positron emission tomography (PET) imaging. A$β$ is a key biomarker of Alzheimer's disease. LNODE is designed to support the fusion, harmonization, quantitative analysis, and interpretation of Abeta PET scans. We evaluate LNODE on 1461 subjects in the ADNI cohort and 1070 subjects in the A4 Study, using MUSE and DKT anatomical atlases. LNODE is formulated as a regional neural ordinary differential equation (ODE) model that is jointly calibrated on all available scans within a cohort. The model captures the spatial propagation, proliferation, and clearance of A$β$ and incorporates a latent-state representation that modulates A$β$ dynamics. The temporal evolution of these latent states is governed by cohort-shared parameters, enabling LNODE to represent both population-level trajectories and subject-specific deviations. The proposed model demonstrates strong parameter identifiability and stability properties, supported by synthetic experiments and analytical analysis of the Hessian condition number. To mitigate overfitting and reduce spurious correlations, LNODE is intentionally underparameterized, employing approximately five to ten parameters per subject. Despite this parsimonious parameterization, LNODE achieves $R^2 > 0.99$ in both the ADNI and A4 datasets. LNODE exhibits strong predictive performance: in the A4 cohort, it accurately forecasts the A$β$ PET signal in previously unseen follow-up scans, including cases with inter-scan intervals exceeding four years. Clustering in the learned latent-state space reveals distinct subgroups, consistent with the existence of different subtypes of Alzheimer's disease progression.

0 citations0 reviewsNo code linkedOpen access
preprint2022arXiv

Overlapping Domain Decomposition Preconditioner for Integral Equations

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.

0 citations0 reviewsNo code linkedOpen access
preprint2022arXiv

Shape dynamics of a red blood cell in Poiseuille flow

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.

0 citations0 reviewsNo code linkedOpen access
preprint2021arXiv

Quantitative in vivo imaging to enable tumor forecasting and treatment optimization

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.

0 citations0 reviewsNo code linkedOpen access
preprint2020arXiv

An efficient method for modeling flow in porous media with immersed faults

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.

0 citations0 reviewsNo code linkedOpen access
preprint2020arXiv

Automatic MRI-Driven Model Calibration for Advanced Brain Tumor Progression Analysis

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.

0 citations0 reviewsNo code linkedOpen access
preprint2020arXiv

Integrated Biophysical Modeling and Image Analysis: Application to Neuro-Oncology

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.

0 citations0 reviewsNo code linkedOpen access
preprint2020arXiv

Multiatlas Calibration of Biophysical Brain Tumor Growth Models with Mass Effect

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.

0 citations0 reviewsNo code linkedOpen access
preprint2020arXiv

Stable shapes of three-dimensional vesicles in unconfined and confined Poiseuille flow

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.

0 citations0 reviewsNo code linkedOpen access
preprint2019arXiv

CLAIRE: A distributed-memory solver for constrained large deformation diffeomorphic image registration

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.

0 citations0 reviewsNo code linkedOpen access
preprint2019arXiv

Image-Driven Biophysical Tumor Growth Model Calibration

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.

0 citations0 reviewsNo code linkedOpen access
preprint2019arXiv

Where did the tumor start? An inverse solver with sparse localization for tumor growth models

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.

0 citations0 reviewsNo code linkedOpen access
preprint2016arXiv

AccFFT: A library for distributed-memory FFT on CPU and GPU architectures

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.

0 citations0 reviewsNo code linkedOpen access
preprint2016arXiv

Constrained $H^1$-regularization schemes for diffeomorphic image registration

We propose regularization schemes for deformable registration and efficient algorithms for their numerical approximation. We treat image registration as a variational optimal control problem. The deformation map is parametrized by its velocity. Tikhonov regularization ensures well-posedness. Our scheme augments standard smoothness regularization operators based on $H^1$- and $H^2$-seminorms with a constraint on the divergence of the velocity field, which resembles variational formulations for Stokes incompressible flows. In our formulation, we invert for a stationary velocity field and a mass source map. This allows us to explicitly control the compressibility of the deformation map and by that the determinant of the deformation gradient. We also introduce a new regularization scheme that allows us to control shear. We use a globalized, preconditioned, matrix-free, reduced space (Gauss--)Newton--Krylov scheme for numerical optimization. We exploit variable elimination techniques to reduce the number of unknowns of our system; we only iterate on the reduced space of the velocity field. Our current implementation is limited to the two-dimensional case. The numerical experiments demonstrate that we can control the determinant of the deformation gradient without compromising registration quality. This additional control allows us to avoid oversmoothing of the deformation map. We also demonstrate that we can promote or penalize shear while controlling the determinant of the deformation gradient.

0 citations0 reviewsNo code linkedOpen access
preprint2016arXiv

FFT, FMM, or Multigrid? A comparative Study of State-Of-the-Art Poisson Solvers for Uniform and Nonuniform Grids in the Unit Cube

In this work, we benchmark and discuss the performance of the scalable methods for the Poisson problem which are used widely in practice: the fast Fourier transform (FFT), the fast multipole method (FMM), the geometric multigrid (GMG), and algebraic multigrid (AMG). In total we compare five different codes, three of which are developed in our group. Our FFT, GMG, and FMM are parallel solvers that use high-order approximation schemes for Poisson problems with continuous forcing functions (the source or right-hand side). We examine and report results for weak scaling, strong scaling, and time to solution for uniform and highly refined grids. We present results on the Stampede system at the Texas Advanced Computing Center and on the Titan system at the Oak Ridge National Laboratory. In our largest test case, we solved a problem with 600 billion unknowns on 229,379 cores of Titan. Overall, all methods scale quite well to these problem sizes. We have tested all of the methods with different source functions (the right-hand side in the Poisson problem). Our results indicate that FFT is the method of choice for smooth source functions that require uniform resolution. However, FFT loses its performance advantage when the source function has highly localized features like internal sharp layers. FMM and GMG considerably outperform FFT for those cases. The distinction between FMM and GMG is less pronounced and is sensitive to the quality (from a performance point of view) of the underlying implementations. The high-order accurate versions of GMG and FMM significantly outperform their low-order accurate counterparts.

0 citations0 reviewsNo code linkedOpen access
preprint2016arXiv

Inv-ASKIT: A Parallel Fast Diret Solver for Kernel Matrices

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.

0 citations0 reviewsNo code linkedOpen access
preprint2015arXiv

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.

0 citations0 reviewsNo code linkedOpen access
preprint2015arXiv

An inverse problem formulation for parameter estimation of a reaction diffusion model of low grade gliomas

We present a numerical scheme for solving a parameter estimation problem for a model of low-grade glioma growth. Our goal is to estimate the spatial distribution of tumor concentration, as well as the magnitude of anisotropic tumor diffusion. We use a constrained optimization formulation with a reaction-diffusion model that results in a system of nonlinear partial differential equations (PDEs). In our formulation, we estimate the parameters using partially observed, noisy tumor concentration data at two different time instances, along with white matter fiber directions derived from diffusion tensor imaging (DTI). The optimization problem is solved with a Gauss-Newton reduced space algorithm. We present the formulation and outline the numerical algorithms for solving the resulting equations. We test the method using a synthetic dataset and compute the reconstruction error for different noise levels and detection thresholds for monofocal and multifocal test cases.

0 citations0 reviewsNo code linkedOpen access
preprint2015arXiv

ASKIT: Approximate Skeletonization Kernel-Independent Treecode in High Dimensions

We present a fast algorithm for kernel summation problems in high-dimensions. These problems appear in computational physics, numerical approximation, non-parametric statistics, and machine learning. In our context, the sums depend on a kernel function that is a pair potential defined on a dataset of points in a high-dimensional Euclidean space. A direct evaluation of the sum scales quadratically with the number of points. Fast kernel summation methods can reduce this cost to linear complexity, but the constants involved do not scale well with the dimensionality of the dataset. The main algorithmic components of fast kernel summation algorithms are the separation of the kernel sum between near and far field (which is the basis for pruning) and the efficient and accurate approximation of the far field. We introduce novel methods for pruning and approximating the far field. Our far field approximation requires only kernel evaluations and does not use analytic expansions. Pruning is not done using bounding boxes but rather combinatorially using a sparsified nearest-neighbor graph of the input. The time complexity of our algorithm depends linearly on the ambient dimension. The error in the algorithm depends on the low-rank approximability of the far field, which in turn depends on the kernel function and on the intrinsic dimensionality of the distribution of the points. The error of the far field approximation does not depend on the ambient dimension. We present the new algorithm along with experimental results that demonstrate its performance. We report results for Gaussian kernel sums for 100 million points in 64 dimensions, for one million points in 1000 dimensions, and for problems in which the Gaussian kernel has a variable bandwidth. To the best of our knowledge, all of these experiments are impossible or prohibitively expensive with existing fast kernel summation methods.

0 citations0 reviewsNo code linkedOpen access
preprint2015arXiv

Comparison of Multigrid Algorithms for High-order Continuous Finite Element Discretizations

We present a comparison of different multigrid approaches for the solution of systems arising from high-order continuous finite element discretizations of elliptic partial differential equations on complex geometries. We consider the pointwise Jacobi, the Chebyshev-accelerated Jacobi and the symmetric successive over-relaxation (SSOR) smoothers, as well as elementwise block Jacobi smoothing. Three approaches for the multigrid hierarchy are compared: 1) high-order $h$-multigrid, which uses high-order interpolation and restriction between geometrically coarsened meshes; 2) $p$-multigrid, in which the polynomial order is reduced while the mesh remains unchanged, and the interpolation and restriction incorporate the different-order basis functions; and 3), a first-order approximation multigrid preconditioner constructed using the nodes of the high-order discretization. This latter approach is often combined with algebraic multigrid for the low-order operator and is attractive for high-order discretizations on unstructured meshes, where geometric coarsening is difficult. Based on a simple performance model, we compare the computational cost of the different approaches. Using scalar test problems in two and three dimensions with constant and varying coefficients, we compare the performance of the different multigrid approaches for polynomial orders up to 16. Overall, both $h$- and $p$-multigrid work well; the first-order approximation is less efficient. For constant coefficients, all smoothers work well. For variable coefficients, Chebyshev and SSOR smoothing outperforms Jacobi smoothing. While all of the tested methods converge in a mesh-independent number of iterations, none of them behaves completely independent of the polynomial order. When multigrid is used as a preconditioner in a Krylov method, the iteration number decreases significantly compared to using multigrid as a solver.

0 citations0 reviewsNo code linkedOpen access
preprint2015arXiv

Far-Field Compression for Fast Kernel Summation Methods in High Dimensions

We consider fast kernel summations in high dimensions: given a large set of points in $d$ dimensions (with $d \gg 3$) and a pair-potential function (the {\em kernel} function), we compute a weighted sum of all pairwise kernel interactions for each point in the set. Direct summation is equivalent to a (dense) matrix-vector multiplication and scales quadratically with the number of points. Fast kernel summation algorithms reduce this cost to log-linear or linear complexity. Treecodes and Fast Multipole Methods (FMMs) deliver tremendous speedups by constructing approximate representations of interactions of points that are far from each other. In algebraic terms, these representations correspond to low-rank approximations of blocks of the overall interaction matrix. Existing approaches require an excessive number of kernel evaluations with increasing $d$ and number of points in the dataset. To address this issue, we use a randomized algebraic approach in which we first sample the rows of a block and then construct its approximate, low-rank interpolative decomposition. We examine the feasibility of this approach theoretically and experimentally. We provide a new theoretical result showing a tighter bound on the reconstruction error from uniformly sampling rows than the existing state-of-the-art. We demonstrate that our sampling approach is competitive with existing (but prohibitively expensive) methods from the literature. We also construct kernel matrices for the Laplacian, Gaussian, and polynomial kernels -- all commonly used in physics and data analysis. We explore the numerical properties of blocks of these matrices, and show that they are amenable to our approach. Depending on the data set, our randomized algorithm can successfully compute low rank approximations in high dimensions. We report results for data sets with ambient dimensions from four to 1,000.

0 citations0 reviewsNo code linkedOpen access
preprint2014arXiv

Adaptive Time Stepping for Vesicle Suspensions

We present an adaptive arbitrary-order accurate time-stepping numerical scheme for the flow of vesicles suspended in Stokesian fluids. Our scheme can be summarized as an approximate implicit spectral deferred correction (SDC) method. Applying a textbook fully implicit SDC scheme to vesicle flows is prohibitively expensive. For this reason we introduce several approximations. Our scheme is based on a semi-implicit linearized low-order time stepping method. (Our discretization is spectrally accurate in space.) We also use invariant properties of vesicle flows, constant area and boundary length in two dimensions, to reduce the computational cost of error estimation for adaptive time stepping. We present results in two dimensions for single-vesicle flows, constricted geometry flows, converging flows, and flows in a Couette apparatus. We experimentally demonstrate that the proposed scheme enables automatic selection of the step size and high-order accuracy.

0 citations0 reviewsNo code linkedOpen access
preprint2014arXiv

High-order adaptive time stepping for vesicle suspensions with viscosity contrast

We construct a high-order adaptive time stepping scheme for vesicle suspensions with viscosity contrast. The high-order accuracy is achieved using a spectral deferred correction (SDC) method, and adaptivity is achieved by estimating the local truncation error with the numerical error of physically constant values. Numerical examples demonstrate that our method can handle suspensions with vesicles that are tumbling, tank-treading, or both. Moreover, we demonstrate that a user-prescribed tolerance can be automatically achieved for simulations with long time horizons.

0 citations0 reviewsNo code linkedOpen access
preprint2014arXiv

High-volume fraction simulations of two-dimensional vesicle suspensions

We consider numerical algorithms for the simulation of the rheology of two-dimensional vesicles suspended in a viscous Stokesian fluid. The vesicle evolution dynamics is governed by hydrodynamic and elastic forces. The elastic forces are due to local inextensibility of the vesicle membrane and resistance to bending. Numerically resolving vesicle flows poses several challenges. For example, we need to resolve moving interfaces, address stiffness due to bending, enforce the inextensibility constraint, and efficiently compute the (non-negligible) long-range hydrodynamic interactions. Our method is based on the work of {\em Rahimian, Veerapaneni, and Biros, "Dynamic simulation of locally inextensible vesicles suspended in an arbitrary two-dimensional domain, a boundary integral method", Journal of Computational Physics, 229 (18), 2010}. It is a boundary integral formulation of the Stokes equations coupled to the interface mass continuity and force balance. We extend the algorithms presented in that paper to increase the robustness of the method and enable simulations with concentrated suspensions. In particular, we propose a scheme in which both intra-vesicle and inter-vesicle interactions are treated semi-implicitly. In addition we use special integration for near-singular integrals and we introduce a spectrally accurate collision detection scheme. We test the proposed methodologies on both unconfined and confined flows for vesicles whose internal fluid may have a viscosity contrast with the bulk medium. Our experiments demonstrate the importance of treating both intra-vesicle and inter-vesicle interactions accurately.

0 citations0 reviewsNo code linkedOpen access
preprint2014arXiv

On preconditioners for the Laplace double-layer in 2D

The discretization of the double-layer potential integral equation for the interior Dirichlet Laplace problem in a domain with smooth boundary results in a linear system that has a bounded condition number. Thus, the number of iterations required for the convergence of a Krylov method is, asymptotically, independent of the discretization size $N$. Using the Fast Multipole Method (FMM) to accelerate the matrix-vector products, we obtain an optimal $\mathcal{O}(N)$ solver. In practice, however, when the geometry is complicated, the number of Krylov iterations can be quite large---to the extend that necessitates the use of preconditioning. We summarize the different methodologies that have appeared in the literature (single-grid, multigrid, approximate sparse inverses) and we propose a new class of preconditioners based on an FMM-based spatial decomposition of the double-layer operator. We present an experimental study in which we compare the different approaches and we discuss the merits and shortcomings of our approach. Our method can be easily extended to other second-kind integral equations with non-oscillatory kernels in two and three dimensions.

0 citations0 reviewsNo code linkedOpen access