Matrix-free Method Research Articles

Simulations of X-ray focusing by zone plates in rotationally symmetric optical field utilizing the matrix-free Finite Difference Beam Propagation Method

We present the use of a finite difference method based on Crank-Nicholson scheme and recurrence scheme for computationally efficient simulation of the X-ray propagation through a zone plate. By introducing boundary and central conditions and by avoiding large matrix operations, the method achieves considerable speed, little memory occupation and low background noise. Accommodating refractive index profiles of arbitrary shape, it can be applied to assist optimizing X-ray zone plates and understanding focusing mechanism.

Self-similarity and recurrence in stability spectra of near-extreme Stokes waves

We consider steady surface waves in an infinitely deep two-dimensional ideal fluid with potential flow, focusing on high-amplitude waves near the steepest wave with a 120 $^{\circ }$ corner at the crest. The stability of these solutions with respect to coperiodic and subharmonic perturbations is studied, using new matrix-free numerical methods. We provide evidence for a plethora of conjectures on the nature of the instabilities as the steepest wave is approached, especially with regards to the self-similar recurrence of the stability spectrum near the origin of the spectral plane.

Sparsity-Independent Lyapunov Exponent in the Sachdev-Ye-Kitaev Model.

The saturation of a recently proposed universal bound on the Lyapunov exponent has been conjectured to signal the existence of a gravity dual. This saturation occurs in the low-temperature limit of the dense Sachdev-Ye-Kitaev (SYK) model, N Majorana fermions with q body (q>2) infinite-range interactions. We calculate certain out-of-time-order correlators (OTOCs) for N≤64 fermions for a highly sparse SYK model and find no significant dependence of the Lyapunov exponent on sparsity up to near the percolation limit where the Hamiltonian breaks up into blocks. This provides strong support to the saturation of the Lyapunov exponent in the low-temperature limit of the sparse SYK. A key ingredient to reaching N=64 is the development of a novel quantum spin model simulation library that implements highly optimized matrix-free Krylov subspace methods on graphical processing units. This leads to a significantly lower simulation time as well as vastly reduced memory usage over previous approaches, while using modest computational resources. Strong sparsity-driven statistical fluctuations require both the use of a much larger number of disorder realizations with respect to the dense limit and a careful finite size scaling analysis. The saturation of the bound in the sparse SYK points to the existence of a gravity analog that would enlarge substantially the number of field theories with this feature.

Mixed node's residual descent method for hyperelastic problem analysis

Geometric nonlinearities, material nonlinearities, and volume locking are the notable challenges faced in hyperelastic analysis. Traditional methods in this regard are complex and laborious for implementation as they require linearization and formulation of global matrix equations while simultaneously addressing volumetric locking. A mixed node's residual descent method (NRDM) proposed herein can effectively address the numerical challenges associated with geometric nonlinearities, material nonlinearities, force nonlinearities, and incompressibility. First, the implementation of the NRDM is considerably simplified as unlike traditional incremental–iterative methods, it's a matrix-free iterative method that does not require incremental linear equations. Second, the NRDM addresses the geometric and material nonlinearities with relative ease as it flexibly describes the deformation with the initial configuration or assumed deformed configuration as the reference frame. Third, the NRDM prevents the occurrence of volumetric locking by employing hydrostatic pressure as an independent variable. Furthermore, the NRDM can easily treat the force nonlinearities and boundary nonlinearities by controlling the relation between load and deformation during iteration. Moreover, a notable accuracy of the NRDM is confirmed through numerical verifications, and several critical matters are discussed, including the scheme for adjusting the basic independent variables, treatment of displacement boundaries, scheme for imposing loads, and computational parameter setting.

A matrix-free parallel two-level deflation preconditioner for two-dimensional heterogeneous Helmholtz problems

We propose a matrix-free parallel two-level deflation method combined with the Complex Shifted Laplacian Preconditioner (CSLP) for two-dimensional heterogeneous Helmholtz problems encountered in seismic exploration, antennas, and medical imaging. These problems pose challenges in terms of accuracy and convergence due to scalability issues with numerical solvers. Motivated by the limitations imposed by excessive computational time and memory constraints when employing a sequential solver with constructed matrices, we parallelize the two-level deflation method without constructing any matrices. Our approach utilizes preconditioned Krylov subspace methods and approximates the CSLP preconditioner with a parallel geometric multigrid V-cycle. For the two-level deflation, standard inter-grid deflation vectors and further high-order deflation vectors are considered. As another main contribution, the matrix-free Galerkin coarsening approach and a novel re-discretization scheme as well as high-order finite-difference schemes on the coarse grid are studied to obtain wavenumber-independent convergence. The optimal settings for an efficient coarse-grid problem solver are investigated. Numerical experiments of model problems show that the wavenumber independence has been obtained for medium wavenumbers. The matrix-free parallel framework shows satisfactory weak and strong parallel scalability.

SmashGP: Large-Scale Spatial Modeling via Matrix-Free Gaussian Processes

Gaussian processes are essential for spatial data analysis. Not only do they allow the prediction of unknown values, but they also allow for uncertainty quantification. However, in the era of big data, directly using Gaussian processes has become computationally infeasible as cubic run times are required for dense matrix decomposition and inversion. Various alternatives have been proposed to reduce the computational burden of directly fitting Gaussian processes. These alternatives rely on assumptions on the underlying structure of the covariance or precision matrices, such as sparsity or low-rank. In contrast, this article uses hierarchical matrices and matrix-free methods to enable the computation of Gaussian processes for large spatial datasets by exploiting the underlying kernel properties. The proposed framework, smashGP, represents the covariance matrix as an H 2 matrix in O ( n ) time and is able to estimate the unknown parameters of the model and predict the values of spatial observations at unobserved locations in O ( n log n ) time thanks to fast matrix-vector products. Additionally, it can be parallelized to take full advantage of shared-memory computing environments. With simulations and case studies, we illustrate the advantage of smashGP to model large-scale spatial datasets. Supplementary materials for this article are available online.

3D cell cultures exhibit increased drug resistance, particularly to taxanes.

e12599 Background: Cell culture platforms could be categorized into two main types: 2D and 3D. 3D cell cultures, known for mimicking the physiological cellular environment, are commonly preferred by researchers to assess drug responses as they accurately reproduce aspects of the in vivo cell behavior. We hypothesized that dose-response curves of 2D and 3D platforms may be drug dependent, with some being more sensitive or resistant depending on the platform. Methods: MCF-7 cancer cells were cultured and plated using a 2D flat bottom well plate platform or using a matrix-based and a matrix-free 3D platform. Cells were treated with chemotherapies and targeted therapies for four days and cellular sensitivity was assessed through a viability assay. Four-parameter logistic (4PL) dose response curves were generated and IC50 values (µM/Cmax) were quantified to understand the differences in dose response curves between the 2D and 3D cellular platforms. Results: While the dose-response curves generally correlate between 2D vs. 3D matrix-based (r2=0.82, p<0.001) and 2D vs. 3D matrix-free methods (r2=0.78, p<0.001), noteworthy variations were observed for certain drugs, particularly paclitaxel (2D IC50= 66; 3D matrix-based IC50 = ≥100; 3D matrix-free IC50 = ≥100) and docetaxel (2D IC50 = ≥34; 3D matrix-based IC50 = ≥100; 3D matrix-free IC50 = ≥100). Additionally, a general increase in resistance was evident when cells formed 3D structures compared to the 2D monolayer condition, as indicated by regression analysis (Y-intercept = 0.55, CI95% 0.32 to 0.78; slope = 0.47, CI95% 0.20 to 0.75). Conclusions: Our findings reveal that taxanes such as paclitaxel and docetaxel exhibit distinctive resistance patterns in 3D structures. In addition, a substantial increase in cellular resistance is observed in 3D cellular models when compared to the 2D assay for most drugs. Despite the overall correlation, these drug-specific differences should be considered. Thus, the next generation of an extreme drug resistance assay, aimed at enhancing the quality of life for patients by avoiding ineffective therapies, might preferentially employ a 3D platform to tailor patients’ cancer therapies.

Matrix-Free Monolithic Multigrid Methods for Stokes and Generalized Stokes Problems

Matrix-Free Monolithic Multigrid Methods for Stokes and Generalized Stokes Problems

Truncated Gauss-Newton full-waveform inversion of pure quasi-P waves in vertical transverse isotropic media

Full-waveform inversion (FWI) uses the full information of seismic data to obtain a quantitative estimation of subsurface physical parameters. Anisotropic FWI has the potential to recover high-resolution velocity and anisotropy parameter models, which are critical for imaging the long-offset and wide-azimuth data. We develop an acoustic anisotropic FWI method based on a simplified pure quasi P-wave (qP-wave) equation, which can be solved efficiently and is beneficial for the subsequent inversion. Using the inverse Hessian operator to precondition the functional gradients helps to reduce the parameter tradeoff in the multi-parameter inversion. To balance the accuracy and efficiency, we extend the truncated Gauss-Newton (TGN) method into FWI of pure qP-waves in vertical transverse isotropic (VTI) media. The inversion is performed in a nested way: a linear inner loop and a nonlinear outer loop. We derive the formulation of Hessian-vector products for pure qP-waves in VTI media based on the Lagrange multiplier method and compute the model update by solving a Gauss-Newton linear system via a matrix-free conjugate gradient method. A suitable preconditioner and the Eisenstat and Walker stopping criterion for the inner iterations are used to accelerate the convergence and avoid prohibitive computational cost. We test the proposed FWI method on several synthetic data sets. Inversion results reveal that the pure acoustic VTI FWI exhibits greater accuracy than the conventional pseudoacoustic VTI FWI. Additionally, the TGN method proves effective in mitigating the parameter crosstalk and increasing the accuracy of anisotropy parameters.

Protocol for high throughput 3D drug screening of patient derived melanoma and renal cell carcinoma

High Throughput Screening (HTS) with 3D cell models is possible thanks to the recent progress and development in 3D cell culture technologies. Results from multiple studies have demonstrated different drug responses between 2D and 3D cell culture. It is now widely accepted that 3D cell models more accurately represent the physiologic conditions of tumors over 2D cell models. However, there is still a need for more accurate tests that are scalable and better imitate the complex conditions in living tissues. Here, we describe ultrahigh throughput 3D methods of drug response profiling in patient derived primary tumors including melanoma as well as renal cell carcinoma that were tested against the NCI oncologic set of FDA approved drugs. We also tested their autologous patient derived cancer associated fibroblasts, varied the in-vitro conditions using matrix vs matrix free methods and completed this in both 3D vs 2D rendered cancer cells. The result indicates a heterologous response to the drugs based on their genetic background, but not on their maintenance condition. Here, we present the methods and supporting results of the HTS efforts using these 3D of organoids derived from patients. This demonstrated the possibility of using patient derived 3D cells for HTS and expands on our screening capabilities for testing other types of cancer using clinically approved anti-cancer agents to find drugs for potential off label use.

A matrix-free parallel solution method for the three-dimensional heterogeneous Helmholtz equation

A matrix-free parallel solution method for the three-dimensional heterogeneous Helmholtz equation

Fast Matrix-Free Methods for Model-Based Personalized Synthetic MR Imaging

Synthetic Magnetic Resonance (MR) imaging predicts images at new design parameter settings from a few observed MR scans. Model-based methods, that use both the physical and statistical properties underlying the MR signal and its acquisition, can predict images at any setting from as few as three scans, allowing it to be used in individualized patient- and anatomy-specific contexts. However, the estimation problem in model-based synthetic MR imaging is ill-posed and so regularization, in the form of correlated Gaussian markov random fields, is imposed on the voxel-wise spin-lattice relaxation time, spin-spin relaxation time and the proton density underlying the MR image. We develop theoretically sound but computationally practical matrix-free estimation methods for synthetic MR imaging. Our evaluations demonstrate superior performance of our methods in currently-used clinical settings when compared to existing model-based and deep learning methods. Moreover, unlike deep learning approaches, our fast methodology can synthesize needed images during patient visits, with good estimation and prediction accuracy and consistency. An added strength of our model-based approach, also developed and illustrated here, is the accurate estimation of standard errors of regional contrasts in the synthesized images. A R package symr implements our methodology. Supplementary materials for this article are available online.

A correction function-based kernel-free boundary integral method for elliptic PDEs with implicitly defined interfaces

This work addresses a novel version of the kernel-free boundary integral (KFBI) method for solving elliptic PDEs with implicitly defined irregular boundaries and interfaces. We focus on boundary value problems and interface problems, which are reformulated into boundary integral equations and solved with the matrix-free GMRES method. In the KFBI method, evaluating boundary and volume integrals only requires solving equivalent but much simpler interface problems in a bounding box, for which fast solvers such as FFTs and geometric multigrid methods are applicable. For the simple interface problem, a correction function is introduced for both the evaluation of right-hand side correction terms and the interpolation of a non-smooth potential function. A mesh-free collocation method is proposed to compute the correction function near the interface. The new method avoids complicated derivation for derivative jumps of the solution and is easy to implement, especially for the fourth-order method in three space dimensions. Various numerical examples are presented, including challenging cases such as high-contrast coefficients, arbitrarily close interfaces and heterogeneous interface problems. The reported numerical results verify that the proposed method is both accurate and efficient.

Fundamental data structures for matrix-free finite elements on hybrid tetrahedral grids

This paper presents efficient data structures for the implementation of matrix-free finite element methods on block-structured, hybrid tetrahedral grids. It provides a complete categorization of all geometric sub-objects that emerge from the regular refinement of the unstructured, tetrahedral coarse grid and describes efficient iteration patterns and analytical linearization functions for the mapping of coefficients to memory addresses. This foundation enables the implementation of fast, extreme-scalable, matrix-free, iterative solvers, and in particular geometric multigrid methods by design. Their application to the variable-coefficient Stokes system subject to an enriched Galerkin discretization and to the curl-curl problem discretized with Nédélec edge elements showcases the flexibility of the implementation. Finally, the solution of a curl-curl problem with 1.6 ⋅ 10 11 (more than one hundred billion) unknowns on more than 32,000 processes with a matrix-free full multigrid solver demonstrates its extreme-scalability.

Jacobian-free variational method for computing connecting orbits in nonlinear dynamical systems.

One approach for describing spatiotemporal chaos is to study the unstable invariant sets embedded in the chaotic attractor of the system. While equilibria, periodic orbits, and invariant tori can be computed using existing methods, the numerical identification of heteroclinic and homoclinic connections between them remains challenging. We propose a robust matrix-free variational method for computing connecting orbits between equilibrium solutions. Instead of a common shooting-based approach, we view the identification of a connecting orbit as a minimization problem in the space of smooth curves in the state space that connect the two equilibria. In this approach, the deviation of a connecting curve from an integral curve of the vector field is penalized by a non-negative cost function. Minimization of the cost function deforms a trial curve until, at a global minimum, a connecting orbit is obtained. The method has no limitation on the dimension of the unstable manifold at the origin equilibrium and does not suffer from exponential error amplification associated with time-marching a chaotic system. Owing to adjoint-based minimization techniques, no Jacobian matrices need to be constructed. Therefore, the memory requirement scales linearly with the size of the problem, allowing the method to be applied to high-dimensional dynamical systems. The robustness of the method is demonstrated for the one-dimensional Kuramoto-Sivashinsky equation.

Node's residual descent method for linear elastic boundary value problems

This study proposes a node's residual descent method (NRDM) for the linear elasticity boundary value problem. In this method, the boundary information and deformation states are attached to the nodes, and the partial differential equations (PDEs) are then transformed into a generalized difference form with the first-order generalized finite difference (GFD) scheme. Furthermore, this study thoroughly describes the NRDM, confirms the accuracy and convergence by numerical examples, performs an error analysis, and discusses the influence of computational parameters. Results demonstrate that NRDM has high computational accuracy and superliner convergence rate. The double derivation is implemented to solve the second-order PDE boundary value problem with the first-order GFD scheme, thus simplifying the requirement for star connectivity and reducing the computational complexity. NRDM is a matrix-free iterative method that eliminates the task of assembling large global sparse matrices, thus avoiding the numerous matrix–matrix product operations. Rather than solving large systems of linear equations, many extremely small matrix multiplications are present during iterations, which are extremely friendly to parallel computing and can be adapted to the computation of ultra-multi-degree-of-freedom problems. The computational accuracy and convergence rate are considerably improved in NRDM compared to the classical GFD method.

A 3D frequency-domain electromagnetic solver employing a high order compact finite-difference scheme

Three-dimensional (3D) electromagnetic (EM) modeling is an important tool for geophysical applications. In EM induction methods the modified Helmholtz equation is often used to describe scattered or residual electric fields in three dimensions. Throughout this paper, a high order compact finite difference scheme for the solution of that equation for vertical magnetic dipole source (VMD) is presented. The approximation of the residual electric field intensity using a fourth order compact finite difference (FD) discretizer is achieved with the solution of a block linear system, the coefficient matrices are large and sparse with a particular structure, implying the application of a matrix-free Krylov subspace method for an efficient numerical solution. The proposed solver is being examined using a number of test problems in uniform and non-uniform grid spacing where numerical and analytical solutions for the homogeneous half-space cases are being compared.

Fast and robust parameter estimation with uncertainty quantification for the cardiac function

Background and objectivesParameter estimation and uncertainty quantification are crucial in computational cardiology, as they enable the construction of digital twins that faithfully replicate the behavior of physical patients. Many model parameters regarding cardiac electromechanics and cardiovascular hemodynamics need to be robustly fitted by starting from a few, possibly non-invasive, noisy observations. Moreover, short execution times and a small amount of computational resources are required for the effective clinical translation. MethodsIn the framework of Bayesian statistics, we combine Maximum a Posteriori estimation and Hamiltonian Monte Carlo to find an approximation of model parameters and their posterior distributions. Fast simulations and minimal memory requirements are achieved by using an accurate and geometry-specific Artificial Neural Network surrogate model for the cardiac function, matrix–free methods, automatic differentiation and automatic vectorization. Furthermore, we account for the surrogate modeling error and measurement error. ResultsWe perform three different in silico test cases, ranging from the ventricular function to the entire cardiocirculatory system, involving whole-heart mechanics, arterial and venous hemodynamics. By employing a single central processing unit on a standard laptop, we attain highly accurate estimations for all model parameters in short computational times. Furthermore, we obtain posterior distributions that contain the true values inside the 90% credibility regions. ConclusionsMany model parameters regarding the entire cardiovascular system can be fastly and robustly identified with minimal hardware requirements. This can be achieved when a small amount of non-invasive data is available and when high levels of signal-to-noise ratio are present in the quantities of interest. With these features, our approach meets the requirements for clinical exploitation, while being compliant with Green Computing practices.

Accelerating high-order mesh optimization using finite element partial assembly on GPUs

In this paper we present a new GPU-oriented mesh optimization method based on high-order finite elements. Our approach relies on node movement with fixed topology, through the Target-Matrix Optimization Paradigm (TMOP) and uses a global nonlinear solve over the whole computational mesh, i.e., all mesh nodes are moved together. A key property of the method is that the mesh optimization process is recast in terms of finite element operations, which allows us to utilize recent advances in the field of GPU-accelerated high-order finite element algorithms. For example, we reduce data motion by using tensor factorization and matrix-free methods, which have superior performance characteristics compared to traditional full finite element matrix assembly and offer advantages for GPU-based HPC hardware. We describe the major mathematical components of the method along with their efficient GPU-oriented implementation. In addition, we propose an easily reproducible mesh optimization test that can serve as a performance benchmark for the mesh optimization community.

Accelerating boundary element methods in wideband frequency sweep analysis by matrix-free model order reduction

In frequency-domain acoustic simulations, a wideband frequency sweep analysis is often required. This becomes particularly challenging for boundary element methods (BEM), as the system matrix has a complex dependency on frequency. Previous attempts have been made to accelerate the frequency sweep analysis for the conventional BEM i.e., without fast multipole or ℋ-matrices accelerations, which restrains the applications to limited scale problems due to the dense system matrix. This paper investigates and incorporates a rational Krylov projection based model order reduction (MOR) technology i.e., matrix-free method into various BEM formulations to address the critical challenge of accelerating frequency sweep analysis. Operating on the transfer functions of the BEM system’s input and output, the matrix-free MOR avoids the direct interference with the BEM kernels, which makes it feasible and attractive to work with existing fast BEM techniques such as ℋ-matrices and fast multipole method. The combined holistic BEM solution retains the merits of the acceleration from both MOR and fast BEM techniques. Significant speed-ups have been achieved by the matrix-free MOR for wideband BEM simulations. Numerical examples are presented to demonstrate the capability of the presented holistic approach in simulating various acoustic problems, including interior–exterior models, complex geometries and frequency dependent boundary conditions.