George Biros

Numerical simulation of deformable particle suspensions in Stokes flow is computationally expensive due to nonlinear fluid-structure interactions, evolving interfaces, and multiscale hydrodynamics. We present VesNet, a hybrid framework that accelerates two-dimensional vesicle suspension simulations by approximating vesicle self interactions, including background flow coupling and short-range lubrication forces, while retaining conventional modules for boundary reparameterization and far-field hydrodynamics. A GPU-accelerated implementation achieves over 100x speedup compared to a multithreaded MATLAB CPU boundary integral solver and about 5x relative to its GPU counterpart. VesNet accurately captures key dynamics, including single-vesicle phase behavior, pair interactions, and large-scale suspensions in Taylor-Green and Poiseuille flows, enabling efficient simulations of thousands of vesicles on modest computational resources.

Numerical simulation of deformable particle suspensions in Stokes flow is computationally expensive due to nonlinear fluid-structure interactions, evolving interfaces, and multiscale hydrodynamics. We present VesNet, a hybrid framework that accelerates two-dimensional vesicle suspension simulations by approximating vesicle self interactions, including background flow coupling and short-range lubrication forces, while retaining conventional modules for boundary reparameterization and far-field hydrodynamics. A GPU-accelerated implementation achieves over 100x speedup compared to a multithreaded MATLAB CPU boundary integral solver and about 5x relative to its GPU counterpart. VesNet accurately captures key dynamics, including single-vesicle phase behavior, pair interactions, and large-scale suspensions in Taylor-Green and Poiseuille flows, enabling efficient simulations of thousands of vesicles on modest computational resources.

We present GPU algorithms for Ewald summation methods for accelerating N-body Stokes flow problems in periodic domains. Like most N-body codes, Ewald sums use a near-field/far-field decomposition. The near field involves particle-to-particle (P2P) interactions. The far field primarily involves particle-to-grid (P2G) and grid-to-particle (G2P) interactions, as well as Fast Fourier Transforms. For each interaction, we investigate several algorithmic variants. Our implementation uses PyKokkos, a Python interface for the Kokkos C++ parallel programming framework, which supports portability to AMD/NVIDIA GPU and ARM/x86 CPU architectures. Double and single-precision numerical results, alongside analytical performance models, confirm the efficiency of our algorithms on AMD and NVIDIA GPU and on ARM and AMD CPU architectures. The P2P interaction achieves around 73% compute efficiency on NVIDIA H200, 84% on NVIDIA A100, 60% on AMD MI300, 52% on Grace CPU, and 68% on AMD Epyc CPU. A straightforward implementation of the P2G kernel can become a computational bottleneck. We introduce a novel P2G algorithm that achieves up to 16$\times$ speedup compared to a baseline GPU implementation. The overall Ewald sum code processes approximately 8 million particles per second on a H200 GPU, and about a half-million particles per second on a Grace CPU, for nine digits of accuracy. We also perform a multi-GPU weak scaling test on up to 256 million particles (64 GPUs) that shows bounded communication cost for all stages except the all-to-all particle sorting, which can be reduced to neighbor communication in the relevant time-stepping regime.

We present GPU algorithms for Ewald summation methods for accelerating N-body Stokes flow problems in periodic domains. Like most N-body codes, Ewald sums use a near-field/far-field decomposition. The near field involves particle-to-particle (P2P) interactions. The far field primarily involves particle-to-grid (P2G) and grid-to-particle (G2P) interactions, as well as Fast Fourier Transforms. For each interaction, we investigate several algorithmic variants. Our implementation uses PyKokkos, a Python interface for the Kokkos C++ parallel programming framework, which supports portability to AMD/NVIDIA GPU and ARM/x86 CPU architectures. Double and single-precision numerical results, alongside analytical performance models, confirm the efficiency of our algorithms on AMD and NVIDIA GPU and on ARM and AMD CPU architectures. The P2P interaction achieves around 73% compute efficiency on NVIDIA H200, 84% on NVIDIA A100, 60% on AMD MI300, 52% on Grace CPU, and 68% on AMD Epyc CPU. A straightforward implementation of the P2G kernel can become a computational bottleneck. We introduce a novel P2G algorithm that achieves up to 16$\times$ speedup compared to a baseline GPU implementation. The overall Ewald sum code processes approximately 8 million particles per second on a H200 GPU, and about a half-million particles per second on a Grace CPU, for nine digits of accuracy. We also perform a multi-GPU weak scaling test on up to 256 million particles (64 GPUs) that shows bounded communication cost for all stages except the all-to-all particle sorting, which can be reduced to neighbor communication in the relevant time-stepping regime.

We introduce a novel neural operator architecture designed to approximate solutions of linear elliptic partial differential equations with high-contrast, spatially varying coefficients. The network, termed the Iterated V-shaped Net (IV-Net), realizes a mapping from the input coefficients and righthand side to the corresponding solution field. The architecture of IV-Net is informed by, and closely resembles, a V-cycle multigrid solver. The IV-Net model is parameterized via convolutional layers defined in the physical domain. For coercive problems with highly heterogeneous coefficients, the proposed network exhibits superior performance relative to a proper orthogonal decomposition (POD) approach and several existing neural operator architectures. For low-frequency oscillatory Helmholtz problems with smooth coefficients, its performance is similar to that of a Fourier neural operator. We analyze the approximation error and convergence behavior of IV-Net, its data efficiency, and its dependence on the underlying discretization mesh. Furthermore, we demonstrate the practical effectiveness of the architecture through a series of numerical experiments, including applications to uncertainty quantification, inverse problems, and prediction of quantities of interest.

We introduce a novel neural operator architecture designed to approximate solutions of linear elliptic partial differential equations with high-contrast, spatially varying coefficients. The network, termed the Iterated V-shaped Net (IV-Net), realizes a mapping from the input coefficients and righthand side to the corresponding solution field. The architecture of IV-Net is informed by, and closely resembles, a V-cycle multigrid solver. The IV-Net model is parameterized via convolutional layers defined in the physical domain. For coercive problems with highly heterogeneous coefficients, the proposed network exhibits superior performance relative to a proper orthogonal decomposition (POD) approach and several existing neural operator architectures. For low-frequency oscillatory Helmholtz problems with smooth coefficients, its performance is similar to that of a Fourier neural operator. We analyze the approximation error and convergence behavior of IV-Net, its data efficiency, and its dependence on the underlying discretization mesh. Furthermore, we demonstrate the practical effectiveness of the architecture through a series of numerical experiments, including applications to uncertainty quantification, inverse problems, and prediction of quantities of interest.

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.

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.

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.

We consider the problem of sampling from a posterior distribution arising in Bayesian inverse problems in science, engineering, and imaging. Our method belongs to the family of independence Metropolis-Hastings (IMH) sampling algorithms, which are common in Bayesian inference. Relying on the existence of an approximate posterior distribution that is cheaper to sample from but may have significant bias, we introduce Proximal-IMH, a scheme that removes this bias by correcting samples from the approximate posterior through an auxiliary optimization problem. This yields a local adjustment that trades off adherence to the exact model against stability around the approximate reference point. For idealized settings, we prove that the proximal correction tightens the match between approximate and exact posteriors, thereby improving acceptance rates and mixing. The method applies to both linear and nonlinear input-output operators and is particularly suitable for inverse problems where exact posterior sampling is too expensive. We present numerical experiments including multimodal and data-driven priors with nonlinear input-output operators. The results show that Proximal-IMH reliably outperforms existing IMH variants.

We consider the problem of sampling from a posterior distribution arising in Bayesian inverse problems in science, engineering, and imaging. Our method belongs to the family of independence Metropolis-Hastings (IMH) sampling algorithms, which are common in Bayesian inference. Relying on the existence of an approximate posterior distribution that is cheaper to sample from but may have significant bias, we introduce Proximal-IMH, a scheme that removes this bias by correcting samples from the approximate posterior through an auxiliary optimization problem. This yields a local adjustment that trades off adherence to the exact model against stability around the approximate reference point. For idealized settings, we prove that the proximal correction tightens the match between approximate and exact posteriors, thereby improving acceptance rates and mixing. The method applies to both linear and nonlinear input-output operators and is particularly suitable for inverse problems where exact posterior sampling is too expensive. We present numerical experiments including multimodal and data-driven priors with nonlinear input-output operators. The results show that Proximal-IMH reliably outperforms existing IMH variants.

We study sampling from posterior distributions in Bayesian linear inverse problems where $A$, the parameters to observables operator, is computationally expensive. In many applications, $A$ can be factored in a manner that facilitates the construction of a cost-effective approximation $\tilde{A}$. In this framework, we introduce Latent-IMH, a sampling method based on the Metropolis-Hastings independence (IMH) sampler. Latent-IMH first generates intermediate latent variables using the approximate $\tilde{A}$, and then refines them using the exact $A$. Its primary benefit is that it shifts the computational cost to an offline phase. We theoretically analyze the performance of Latent-IMH using KL divergence and mixing time bounds. Using numerical experiments on several model problems, we show that, under reasonable assumptions, it outperforms state-of-the-art methods such as the No-U-Turn sampler (NUTS) in computational efficiency. In some cases, Latent-IMH can be orders of magnitude faster than existing schemes.

We study sampling from posterior distributions in Bayesian linear inverse problems where $A$, the parameters to observables operator, is computationally expensive. In many applications, $A$ can be factored in a manner that facilitates the construction of a cost-effective approximation $\tilde{A}$. In this framework, we introduce Latent-IMH, a sampling method based on the Metropolis-Hastings independence (IMH) sampler. Latent-IMH first generates intermediate latent variables using the approximate $\tilde{A}$, and then refines them using the exact $A$. Its primary benefit is that it shifts the computational cost to an offline phase. We theoretically analyze the performance of Latent-IMH using KL divergence and mixing time bounds. Using numerical experiments on several model problems, we show that, under reasonable assumptions, it outperforms state-of-the-art methods such as the No-U-Turn sampler (NUTS) in computational efficiency. In some cases, Latent-IMH can be orders of magnitude faster than existing schemes.

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.

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.

As the complexity and scale of modern parallel machines continue to grow, programmers increasingly rely on composition of software libraries to encapsulate and exploit parallelism. However, many libraries are not designed with composition in mind and assume they have exclusive access to all resources. Using such libraries concurrently can result in contention and degraded performance. Prior solutions involve modifying the libraries or the OS, which is often infeasible.

As the complexity and scale of modern parallel machines continue to grow, programmers increasingly rely on composition of software libraries to encapsulate and exploit parallelism. However, many libraries are not designed with composition in mind and assume they have exclusive access to all resources. Using such libraries concurrently can result in contention and degraded performance. Prior solutions involve modifying the libraries or the OS, which is often infeasible.

Aβ Positron Emission Tomography (PET) is often used to manage Alzheimer’s disease (AD). To better understand Aβ progression, we introduce and evaluate a mathematical model that couples Aβ at parcellated gray matter regions. We term this model LNODE for “latent network ordinary differential equations”. At each region, we track normal Aβ, abnormal Aβ, and m latent states that intend to capture unobservable mechanisms coupled to Aβ progression. LNODE is parameterized by subject-specific parameters and cohort parameters. We jointly invert for these parameters by fitting the model to Aβ-PET data from 585 subjects from the ADNI dataset. Although underparameterized, our model achieves population R2≥98% compared to R2≤60% when fitting without latent states. Furthermore, these preliminary results suggest the existence of different subtypes of Aβ progression.

Abnormal tau and amyloid beta are two primary imaging biomarkers used to assist in the diagnosis of Alzheimer's disease (AD). Recent efforts have focused on developing mechanism-based biophysical models to explain the spatiotemporal dynamics of these biomarkers. In this study, we adopt a connectome-based ODE model to capture the dynamics of tau and amyloid beta (Aβ), aiming to predict personalized future values of these biomarkers. The ODE model includes diffusion, reaction, and clearance terms, and accounts for tau-Aβ interactions. Additionally, it assumes a sparse initial condition (IC) of abnormalities, based on the assumption of localized initiation. Besides tau and Aβ, brain atrophy is used as a marker of neurodegeneration. We discuss the mathematical model of atrophy integrated into the tau-Aβ model. A common limitation in popular models is the use of chronological age as the time axis, which prevents the unification of subject trajectories onto a common time scale and hinders comprehensive cohort analysis. To address this issue, we use a normalized disease age that relates chronological age to biomarker values. In the ODE model, we use the disease age to track time and the biomarker dynamics. Furthermore, our analysis of region-of-interest-wise tau-Aβ temporal correlation reveals that different regions of interest (ROIs) play distinct roles across various disease stages.

We study parallel particle-in-cell (PIC) methods for low-temperature plasmas (LTPs), which discretize kinetic formulations that capture the time evolution of the probability density function of particles as a function of position and velocity. We use a kinetic description for electrons and a fluid approximation for heavy species. In this paper, we focus on GPU acceleration of algorithms for velocity-space interactions and in particular, collisions of electrons with neutrals, ions, and electrons. Our work has two thrusts. The first is algorithmic exploration and analysis. The second is examining the viability of rapid-prototyping implementations using Python-based HPC tools, in particular PyKokkos. We discuss several common PIC kernels and present performance results on NVIDIA Volta V100 and AMD MI250X GPUs. Overall, the MI250X is slightly faster for most kernels but shows more sensitivity to register pressure. We also report scaling results for a distributed memory implementation on up to 16 MPI ranks.

We propose a method that uses a two-species ordinary differential equation (ODE) model for subject-specific misfolded tau protein spreading in Alzheimer's disease (AD) and calibrates it from magnetic resonance imaging (MRI) and positron emission tomography (PET) scans. The ODE model is a variant of the heterodimer Fisher-Kolmogorov (HFK) model. 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 evaluate the performance of our method on both synthetic and clinical data. Subjects are 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 19.6% relative error compared to the FK model on the AD dataset. HFK also achieves an R-squared score of 0.591 for fitting AD data compared with 0.256 from FK model results with IC fixing at EC. The inverted IC from our scheme indicates that the EC is the most likely initial seeding region if subcortical regions are excluded from the analysis. However, other regions also have probability to be the IC seeding regions. Furthermore, for cases that have longitudinal data, we estimate a subject-specific AD onset time.

We present a multi-species partial differential equation (PDE) model for tumor growth and an algorithm for calibrating the model from magnetic resonance imaging (MRI) scans. The model is designed for glioblastoma multiforme (GBM) a fast-growing type of brain cancer. The modeled species correspond to proliferative, infiltrative, and necrotic tumor cells. The model calibration is formulated as an inverse problem and solved by a PDE-constrained optimization method. The data that drives the calibration is derived by a single multi-parametric MRI image. This is a typical clinical scenario for GBMs. The unknown parameters that need to be calibrated from data include 10 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. Our regularization improves the tumor Dice score by 5%-10% compared to single-species model reconstruction. Also, our regularization reduces model parameter reconstruction errors by 4%-80% in cases with known initial condition and brain anatomy compared to cases without regularization. Importantly, our model can estimate infiltrative tumor cells using observable tumor species.

Recent trends in high-performance computing show an increase in the adoption of performance portable frameworks such as Kokkos and interpreted languages such as Python. PyKokkos follows these trends and enables programmers to write performance-portable kernels in Python which greatly increases productivity. One issue that programmers still face is how to organize parallel code, as splitting code into separate kernels simplifies testing and debugging but may result in suboptimal performance. To enable programmers to organize kernels in any way they prefer while ensuring good performance, we present PyFuser, a program analysis framework for automatic fusion of performance portable PyKokkos kernels. PyFuser dynamically traces kernel calls and lazily fuses them once the result is requested by the application. PyFuser generates fused kernels that execute faster due to better reuse of data, improved compiler optimizations, and reduced kernel launch overhead, while not requiring any changes to existing PyKokkos code. We also introduce automated code transformations that further optimize the fused kernels generated by PyFuser. Our experiments show that on average PyFuser achieves speedups compared to unfused kernels of 3.8x on NVIDIA and AMD GPUs, as well as Intel and AMD CPUs.

Constructing k-nearest neighbor (kNN) graphs is a fundamental component in many machine learning and scientific computing applications. Despite its prevalence, efficiently building all-nearest-neighbor graphs at scale on distributed heterogeneous HPC systems remains challenging, especially for large sparse non-integer datasets. We introduce optimizations for algorithms based on forests of random projection trees. Our novel GPU kernels for batched, within leaf, exact searches achieve 1.18× speedup over sparse reference kernels with less peak memory, and up to 19× speedup over CPU for memory-intensive problems. Our library, PyRKNN , implements distributed randomized projection forests for approximate kNN search. Optimizations to reduce and hide communication overhead allow us to achieve 5× speedup, in per iteration performance, relative to GOFMM (another projection tree, MPI-based kNN library), for a 64M 128d dataset on 1,024 processes. On a single-node we achieve speedup over FAISS-GPU for dense datasets and up to 10× speedup over CPU-only libraries. PyRKNN uniquely supports distributed memory kNN graph construction for both dense and sparse coordinates on CPU and GPU accelerators.

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.

We present an algorithm for solving the one-dimensional space collisional Boltzmann transport equation (BTE) for electrons in low-temperature plasmas (LTPs). Modeling LTPs is useful in many applications, including advanced manufacturing, material processing, and hypersonic flows, to name a few. The proposed BTE solver is based on an Eulerian formulation. It uses Chebyshev collocation method in physical space and a combination of Galerkin and discrete ordinates in velocity space. We present self-convergence results and cross-code verification studies compared to an in-house particle-in-cell (PIC) direct simulation Monte Carlo (DSMC) code. Boltzsim is our open source implementation of the solver. Furthermore, we use Boltzsim to simulate radio-frequency glow discharge plasmas (RF-GDPs) and compare with an existing methodology that approximates the electron BTE. We compare these two approaches and quantify their differences as a function of the discharge pressure. The two approaches show an 80x, 3x, 1.6x, and 0.98x difference between cycle-averaged time periodic electron number density profiles at 0.1 Torr, 0.5 Torr, 1 Torr, and 2 Torr discharge pressures, respectively. As expected, these differences are significant at low pressures, for example less than 1 Torr.

C++ remains one of the most widely used languages in various computing fields, from embedded programming to high-performance computing. While new features are constantly being added to C++, an important aspect of the language that is often overlooked is its compilation time. Merely including a few header files can cause compilation time to increase significantly. An alternative to including header files is using forward declarations; however, the rules for forward declaring classes and functions are non obvious and confusing to most developers. Additionally, forward declaring methods, as well as functions that accept lambdas as arguments, is not possible. In this paper, we present a novel technique, termed Header Substitution, to automatically detect opportunities for forward declarations with the goal of replacing includes of header files and improving compilation time. Header Substitution also introduces function wrappers as an alternative to forward declaring methods and functions with lambda arguments. We implemented Header Substitution in a tool, dubbed Yalla, and applied it to various C++ projects in order to speed up the development cycle, i.e., the debugging, editing, compiling, and rerunning loop, achieving up to a 24.5x speedup when compiling C++ files and a 4.68x speedup of the development cycle.

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}\left(c^{2} d^{2}+n c^{2} d\right)$ storage and $\mathcal{O}\left(c^{3}\left(n d^{2}+b d^{3}+b n\right)\right)$ 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}\left(n(d+c)+c d^{2}\right)$ and a computational complexity of $\mathcal{O}\left(b n c d^{2}\right)$. 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.

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 ϵ-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 ϵ-nearest neighbor graphs can be expensive in high-dimensions. This article introduces a modified DBSCAN, using k -nearest neighbor ( k NN) graphs to improve efficiency. We outline conditions for k NN-DBSCAN to match DBSCAN’s results and present a parallel implementation using OpenMP and MPI for shared and distributed memory systems. Testing on datasets up to 32 dimensions, we achieve remarkable scalability. Our implementation clusters one billion 3D points in under one second on 28K cores at TACC’s Frontera system. In a larger run, we cluster 65 billion points in 20 dimensions in under 40 seconds using 114,688 cores. Our method is up to 37× faster than state-of-the-art parallel DBSCAN on a 20-dimensional dataset with 4 million points. Code is available at https://github.com/ut-padas/knndbscan .

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.

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.

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.

Modern Python programs in high-performance computing call into compiled libraries and kernels for performance-critical tasks. However, effectively parallelizing these finer-grained, and often dynamic, kernels across modern heterogeneous platforms remains a challenge. First, we perform an experimental study to examine the impact of Python's Global Interpreter Lock (GIL), and potential speedups under a GIL-less PEP703 future, to guide runtime design. Using our optimized runtime, we explore scheduling tasks with constraints that require resources across multiple, potentially diverse, devices through the introduction of new programming abstractions and runtime mechanisms. We extend an existing Python tasking library, Parla, to augment its performance and add support for such multi-device tasks. Our experimental analysis, using tasks graphs from synthetic and real applications, shows at least a 3 ×(and up to 6 ×) performance improvement over its predecessor in scenarios with high GIL contention. When scheduling multi-GPU tasks, we observe an 8x reduction in per-task launching overhead compared to a multi-process system.

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 [1]." 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 righthand sides. Large-scale numerical experiments are presented to show the performance of our Julia implementation.