Major Computational Bottleneck Research Articles

Accelerating large-scale multi-scalar multiplication in Zk-SNARK through exploiting its multilevel parallelism

In the context of zk-SNARK, MSM emerges as a major computational bottleneck, particularly due to its high computational and memory overhead. In this work, we exploit multiple levels of parallelism to effectively accelerate largescale MSM in zk-SNARKs. Firstly, a distributed parameter generation method is proposed in this paper to replace that of centralized method. Based on this methodology, the paper realizes a system with extraordinary scalability, capable of computing large-scale MSMs. Subsequently, the approach presented in this paper elevates the computation of Bellperson, the most prominent zero-knowledge proof system in real-world applications, from a single-node computation to a clustering mode, significantly enhancing its computational performance – a crucial advancement for practical applications. Finally, we implement a multi-level, fully parallelised MSM computing system by leveraging hierarchical sub-task partitioning and cross-node communication optimization, thereby thoroughly exploiting parallelism at diverse granularities. Experimental results show that in the cluster scenario, the proposed approach achieves acceleration ratios of approximately 3.60 and 6.50 times compared to the cutting-edge heterogeneous version Bellperson in dual-node and quad-node settings, respectively. On a single node, the proposed optimization approach achieves an acceleration ratio of 1.38 times compared to the current State-of-the-Art MSM calculation module of cuZK, outperforms the industry-popular Bellman by 186 times and the Bellperson by 1.96 times.

A structure translation model for crystal compounds

High-throughput virtual screening for crystals aims to discover new materials by evaluating the property of every virtual candidate in the database exhaustively. During this process, the major computational bottleneck is the costly structural relaxation of each hypothetical material on the large-scale dataset using density functional theory (DFT) calculations. Here, we present a generative domain translation framework that maps the unrelaxed structural domains to the relaxed domains, enabling data-driven structural translations. The model predicts the materials formation energy with a small mean absolute error without DFT relaxations, and furthermore can produce the atomic coordinates consistent with the DFT relaxed structures. The utility of the proposed concept is not restricted to the structural domains, and we expect that it can be extended to translate the domain of easy-to-compute properties into the domain of more difficult properties.

Faster Convergence in Deep-Predictive-Coding Networks to Learn Deeper Representations.

Deep-predictive-coding networks (DPCNs) are hierarchical, generative models. They rely on feed-forward and feedback connections to modulate latent feature representations of stimuli in a dynamic and context-sensitive manner. A crucial element of DPCNs is a forward-backward inference procedure to uncover sparse, invariant features. However, this inference is a major computational bottleneck. It severely limits the network depth due to learning stagnation. Here, we prove why this bottleneck occurs. We then propose a new forward-inference strategy based on accelerated proximal gradients. This strategy has faster theoretical convergence guarantees than the one used for DPCNs. It overcomes learning stagnation. We also demonstrate that it permits constructing deep and wide predictive-coding networks. Such convolutional networks implement receptive fields that capture well the entire classes of objects on which the networks are trained. This improves the feature representations compared with our lab's previous nonconvolutional and convolutional DPCNs. It yields unsupervised object recognition that surpass convolutional autoencoders and is on par with convolutional networks trained in a supervised manner.

Digital in-memory stochastic computing architecture for vector-matrix multiplication

The applications of the Artificial Intelligence are currently dominating the technology landscape. Meanwhile, the conventional Von Neumann architectures are struggling with the data-movement bottleneck to meet the ever-increasing performance demands of these data-centric applications. Moreover, The vector-matrix multiplication cost, in the binary domain, is a major computational bottleneck for these applications. This paper introduces a novel digital in-memory stochastic computing architecture that leverages the simplicity of the stochastic computing for in-memory vector-matrix multiplication. The proposed architecture incorporates several new approaches including a new stochastic number generator with ideal binary-to-stochastic mapping, a best seeding approach for accurate-enough low stochastic bit-precisions, a hybrid stochastic-binary accumulation approach for vector-matrix multiplication, and the conversion of conventional memory read operations into on-the-fly stochastic multiplication operations with negligible overhead. Thanks to the combination of these approaches, the accuracy analysis of the vector-matrix multiplication benchmark shows that scaling down the stochastic bit-precision from 16-bit to 4-bit achieves nearly the same average error (less than 3%). The derived analytical model of the proposed in-memory stochastic computing architecture demonstrates that the 4-bit stochastic architecture achieves the highest throughput per sub-array (122 Ops/Cycle), which is better than the 16-bit stochastic precision by 4.36x, while still maintaining a small average error of 2.25%.

A linear programming approach to difference-of-convex piecewise linear approximation

We address the problem of finding continuous piecewise linear (CPWL) approximations of deterministic functions of any dimension that satisfy any predefined error-tolerance, while keeping the number of polytopes that partition the approximation domain low. Specifically, we focus on overcoming the major computational bottleneck of the CPWL Approximation Algorithm (CPWL-AA) that has been proposed in the recent literature. CPWL-AA uses the difference-of-convex CPWL representation to search CPWL approximations which can partition the approximation domain to have polytopes of any shape. A computational bottleneck of the method is to solve a mixed-integer linear program (MILP) in which the number of binary variables is large for many problems of practical interest. In this paper, we overcome this by introducing a method that obtains a high quality solution of the MILP by iteratively solving a linear program (LP). We further reduce the computational expense by developing a method that treats some constraints in the LP problem as lazy constraints. Through a computational study we demonstrate that the proposed methods substantially reduce the computation time of CPWL-AA, while maintaining high quality CPWL approximations. With this, we demonstrate that we can generate CPWL approximations that satisfy predefined error-tolerances on functions of up to five dimensions within reasonable solution times.

Accelerating the density-functional tight-binding method using graphical processing units.

Acceleration of the density-functional tight-binding (DFTB) method on single and multiple graphical processing units (GPUs) was accomplished using the MAGMA linear algebra library. Two major computational bottlenecks of DFTB ground-state calculations were addressed in our implementation: the Hamiltonian matrix diagonalization and the density matrix construction. The code was implemented and benchmarked on two different computer systems: (1) the SUMMIT IBM Power9 supercomputer at the Oak Ridge National Laboratory Leadership Computing Facility with 1-6 NVIDIA Volta V100 GPUs per computer node and (2) an in-house Intel Xeon computer with 1-2 NVIDIA Tesla P100 GPUs. The performance and parallel scalability were measured for three molecular models of 1-, 2-, and 3-dimensional chemical systems, represented by carbon nanotubes, covalent organic frameworks, and water clusters.

Resta-like preconditioning for self-consistent field iterations in the linearized augmented planewave method

Convergence in self-consistent-field cycles can be a major computational bottleneck of density-functional theory calculations. We propose a Resta-like preconditioning method for full-potential all-electron calculations in the linearized augmented planewave method to smoothly converge to self-consistency. We implemented this preconditioner in the exciting code and apply it to the two semiconducting systems of MoS2 slabs and P-rich GaP(100) surfaces as well as the metallic system Au(111), containing a sufficiently large amount of vacuum. Two magnetic systems of 19-atoms Fe and Co are also considered. Our calculations demonstrate that the implemented scheme performs reliably as well as more efficiently regardless of system size, suppressing long-range charge sloshing. While the suitability of this preconditioning higher for semiconducting systems, the convergence for metals is, depending on the system, only slightly de- or increased and thus still trustworthy to apply. Furthermore, a mixing algorithm with the preconditioner shows an improvement over that with the Kerker preconditioner for the investigated semiconducting systems.

Finite-Size Error Cancellation in Diffusion Monte Carlo Calculations of Surface Chemistry.

The accurate prediction of reaction mechanisms in heterogeneous (surface) catalysis is one of the central challenges in computational chemistry. Quantum Monte Carlo methods─Diffusion Monte Carlo (DMC) in particular─are being recognized as higher-accuracy, albeit more computationally expensive, alternatives to Density Functional Theory (DFT) for energy predictions of catalytic systems. A major computational bottleneck in the broader adoption of DMC for catalysis is the need to perform finite-size extrapolations by simulating increasingly large periodic cells (supercells) to eliminate many-body finite-size effects and obtain energies in the thermodynamic limit. Here, we show that it is possible to significantly reduce this computational cost by leveraging the cancellation of many-body finite-size errors that accompanies the evaluation of energy differences when calculating quantities like adsorption (binding) energies and mapping potential energy surfaces. We analyze the cancellation and convergence of many-body finite-size errors in two well-known adsorbate/slab systems, H2O/LiH(001) and CO/Pt(111). Based on this analysis, we identify strategies for obtaining binding energies in the thermodynamic limit that optimally utilize error cancellation to balance accuracy and computational efficiency. Using one such strategy, we then predict the correct order of adsorption site preference on CO/Pt(111), a challenging problem for a wide range of density functionals. Our accurate and inexpensive DMC calculations are found to unambiguously recover the top > bridge > hollow site order, in agreement with experimental observations. We proceed to use this DMC method to map the potential energy surface of CO hopping between Pt(111) adsorption sites. This reveals the existence of an L-shaped top-bridge-hollow diffusion trajectory characterized by energy barriers that provide an additional kinetic justification for experimental observations of CO/Pt(111) adsorption. Overall, this work demonstrates that it is routinely possible to achieve order-of-magnitude speedups and memory savings in DMC calculations by taking advantage of error cancellation in the calculation of energy differences that are ubiquitous in heterogeneous catalysis and surface chemistry more broadly.

Programmable Quantum Annealers as Noisy Gibbs Samplers

Drawing independent samples from high-dimensional probability distributions represents the major computational bottleneck for modern algorithms, including powerful machine learning frameworks such as deep learning. The quest for discovering larger families of distributions for which sampling can be efficiently realized has inspired an exploration beyond established computing methods and turning to novel physical devices that leverage the principles of quantum computation. Quantum annealing embodies a promising computational paradigm that is intimately related to the complexity of energy landscapes in Gibbs distributions, which relate the probabilities of system states to the energies of these states. Here, we study the sampling properties of physical realizations of quantum annealers which are implemented through programmable lattices of superconducting flux qubits. Comprehensive statistical analysis of the data produced by these quantum machines shows that quantum annealers behave as samplers that generate independent configurations from low-temperature noisy Gibbs distributions. We show that the structure of the output distribution probes the intrinsic physical properties of the quantum device such as effective temperature of individual qubits and magnitude of local qubit noise, which result in a non-linear response function and spurious interactions that are absent in the hardware implementation. We anticipate that our methodology will find widespread use in characterization of future generations of quantum annealers and other emerging analog computing devices.

End-to-end orientation estimation from 2D cryo-EM images.

Cryo-electron microscopy (cryo-EM) is a Nobel Prize-winning technique for determining high-resolution 3D structures of biological macromolecules. A 3D structure is reconstructed from hundreds of thousands of noisy 2D projection images. However, existing 3D reconstruction methods are still time-consuming, and one of the major computational bottlenecks is recovering the unknown orientation of the particle in each 2D image. The dominant methods typically exploit an expensive global search on each image to estimate the missing orientations. Here, a novel end-to-end supervised learning method is introduced to directly recover the missing orientations from 2D cryo-EM images. A neural network is used to approximate the mapping from images to orientations. A robust loss function is proposed for optimizing the parameters of the network, which can handle both asymmetric and symmetric 3D structures. Experiments on synthetic data sets with various symmetry types confirm that the neural network is capable of recovering orientations from 2D cryo-EM images, and the results on a real cryo-EM data set further demonstrate its potential under more challenging imaging conditions.

Multigrid reduction preconditioning framework for coupled processes in porous and fractured media

Many subsurface engineering applications involve tight-coupling between fluid flow, solid deformation, fracturing, and similar processes. To better understand the complex interplay of different governing equations, and therefore design efficient and safe operations, numerical simulations are widely used. Given the relatively long time-scales of interest, fully-implicit time-stepping schemes are often necessary to avoid time-step stability restrictions. A major computational bottleneck for these methods, however, is the linear solver. These systems are extremely large and ill-conditioned. Because of the wide range of processes and couplings that may be involved – e.g. formation and propagation of fractures, deformation of the solid porous medium, viscous flow of one or more fluids in the pores and fractures, complicated well sources and sinks, etc. – it is difficult to develop general-purpose but scalable linear solver frameworks. This challenge is further aggravated by the range of different discretization schemes that may be adopted, which have a direct impact on the linear system structure. To address this obstacle, we describe a flexible strategy based on multigrid reduction (MGR) that can produce purely algebraic preconditioners for a wide spectrum of relevant physics and discretizations. We demonstrate that MGR, guided by physics and theory in block preconditioning, can tackle several distinct and challenging problems, notably: a hybrid discretization of single-phase flow, compositional multiphase flow with complex wells, and hydraulic fracturing simulations. Extension to other systems can be handled quite naturally. We demonstrate the efficiency and scalability of the resulting solvers through numerical examples of difficult, field-scale problems.

High performance implementation of the hierarchical likelihood for generalized linear mixed models: an application to estimate the potassium reference range in massive electronic health records datasets

BackgroundConverting electronic health record (EHR) entries to useful clinical inferences requires one to address the poor scalability of existing implementations of Generalized Linear Mixed Models (GLMM) for repeated measures. The major computational bottleneck concerns the numerical evaluation of multivariable integrals, which even for the simplest EHR analyses may involve millions of dimensions (one for each patient). The hierarchical likelihood (h-lik) approach to GLMMs is a methodologically rigorous framework for the estimation of GLMMs that is based on the Laplace Approximation (LA), which replaces integration with numerical optimization, and thus scales very well with dimensionality.MethodsWe present a high-performance, direct implementation of the h-lik for GLMMs in the R package TMB. Using this approach, we examined the relation of repeated serum potassium measurements and survival in the Cerner Real World Data (CRWD) EHR database. Analyzing this data requires the evaluation of an integral in over 3 million dimensions, putting this problem beyond the reach of conventional approaches. We also assessed the scalability and accuracy of LA in smaller samples of 1 and 10% size of the full dataset that were analyzed via the a) original, interconnected Generalized Linear Models (iGLM), approach to h-lik, b) Adaptive Gaussian Hermite (AGH) and c) the gold standard for multivariate integration Markov Chain Monte Carlo (MCMC).ResultsRandom effects estimates generated by the LA were within 10% of the values obtained by the iGLMs, AGH and MCMC techniques. The H-lik approach was 4–30 times faster than AGH and nearly 800 times faster than MCMC. The major clinical inferences in this problem are the establishment of the non-linear relationship between the potassium level and the risk of mortality, as well as estimates of the individual and health care facility sources of variations for mortality risk in CRWD.ConclusionsWe found that the direct implementation of the h-lik offers a computationally efficient, numerically accurate approach for the analysis of extremely large, real world repeated measures data via the h-lik approach to GLMMs. The clinical inference from our analysis may guide choices of treatment thresholds for treating potassium disorders in the clinic.

Modeling the α- and β-resorcinol phase boundary via combination of density functional theory and density functional tight-binding.

The ability to predict not only what organic crystal structures might occur but also the thermodynamic conditions under which they are the most stable would be extremely useful for discovering and designing new organic materials. The present study takes a step in that direction by predicting the temperature- and pressure-dependent phase boundary between the α and β polymorphs of resorcinol using density functional theory (DFT) and the quasi-harmonic approximation. To circumvent the major computational bottleneck associated with computing a well-converged phonon density of states via the supercell approach, a recently developed approximation is employed, which combines a supercell phonon density of states from dispersion-corrected third-order density functional tight binding [DFTB3-D3(BJ)] with frequency corrections derived from a smaller B86bPBE-XDM functional DFT phonon calculation on the crystallographic unit cell. This mixed DFT/DFTB quasi-harmonic approach predicts the lattice constants and unit cell volumes to within 1%-2% at lower pressures. It predicts the thermodynamic phase boundary in almost perfect agreement with the experiment, although this excellent agreement does reflect fortuitous cancellation of errors between the enthalpy and entropy of transition.

Parareal with a learned coarse model for robotic manipulation

A key component of many robotics model-based planning and control algorithms is physics predictions, that is, forecasting a sequence of states given an initial state and a sequence of controls. This process is slow and a major computational bottleneck for robotics planning algorithms. Parallel-in-time integration methods can help to leverage parallel computing to accelerate physics predictions and thus planning. The Parareal algorithm iterates between a coarse serial integrator and a fine parallel integrator. A key challenge is to devise a coarse model that is computationally cheap but accurate enough for Parareal to converge quickly. Here, we investigate the use of a deep neural network physics model as a coarse model for Parareal in the context of robotic manipulation. In simulated experiments using the physics engine Mujoco as fine propagator we show that the learned coarse model leads to faster Parareal convergence than a coarse physics-based model. We further show that the learned coarse model allows to apply Parareal to scenarios with multiple objects, where the physics-based coarse model is not applicable. Finally, we conduct experiments on a real robot and show that Parareal predictions are close to real-world physics predictions for robotic pushing of multiple objects. Code (https://doi.org/10.5281/zenodo.3779085) and videos (https://youtu.be/wCh2o1rf-gA) are publicly available.

An adaptive method for speeding up the numerical integration of chemical mechanisms in atmospheric chemistry models: application to GEOS-Chem version 12.0.0

Abstract. The major computational bottleneck in atmospheric chemistry models is the numerical integration of the stiff coupled system of kinetic equations describing the chemical evolution of the system as defined by the model chemical mechanism (typically over 100 coupled species). We present an adaptive method to greatly reduce the computational cost of that numerical integration in global 3-D models while maintaining high accuracy. Most of the atmosphere does not in fact require solving for the full chemical complexity of the mechanism, so considerable simplification is possible if one can recognize the dynamic continuum of chemical complexity required across the atmospheric domain. We do this by constructing a limited set of reduced chemical mechanisms (chemical regimes) to cover the range of atmospheric conditions and then pick locally and on the fly which mechanism to use for a given grid box and time step on the basis of computed production and loss rates for individual species. Application to the GEOS-Chem global 3-D model for oxidant–aerosol chemistry in the troposphere and stratosphere (full mechanism of 228 species) is presented. We show that 20 chemical regimes can largely encompass the range of conditions encountered in the model. Results from a 2-year GEOS-Chem simulation shows that our method can reduce the computational cost of chemical integration by 30 %–40 % while maintaining accuracy better than 1 % and with no error growth. Our method retains the full complexity of the original chemical mechanism where it is needed, provides the same model output diagnostics (species production and loss rates, reaction rates) as the full mechanism, and can accommodate changes in the chemical mechanism or in model resolution without having to reconstruct the chemical regimes.

Mstar – a fast parallelized algorithmically regularized integrator with minimum spanning tree coordinates

ABSTRACTWe present the novel algorithmically regularized integration method mstar for high-accuracy (|ΔE/E| ≳ 10−14) integrations of N-body systems using minimum spanning tree coordinates. The twofold parallelization of the $\mathcal {O}(N_\mathrm{part}^2)$ force loops and the substep divisions of the extrapolation method allow for a parallel scaling up to NCPU = 0.2 × Npart. The efficient parallel scaling of mstar makes the accurate integration of much larger particle numbers possible compared to the traditional algorithmic regularization chain (ar-chain) methods, e.g. Npart = 5000 particles on 400 CPUs for 1 Gyr in a few weeks of wall-clock time. We present applications of mstar on few particle systems, studying the Kozai mechanism and N-body systems like star clusters with up to Npart = 104 particles. Combined with a tree or fast multipole-based integrator, the high performance of mstar removes a major computational bottleneck in simulations with regularized subsystems. It will enable the next-generation galactic-scale simulations with up to 109 stellar particles (e.g. $m_\star = 100 \, \mathrm{M}_\odot$ for an $M_\star = 10^{11} \, \mathrm{M}_\odot$ galaxy), including accurate collisional dynamics in the vicinity of nuclear supermassive black holes.

Exact two-component equation-of-motion coupled-cluster singles and doubles method using atomic mean-field spin-orbit integrals.

A new semi-atomic-orbital- based algorithm for a two-component spin-orbit (SO) equation-of-motion coupled-cluster singles and doubles (EOM-CCSD) method using mean-field SO integrals is reported. The new algorithm removes the major computational bottlenecks of a SO-EOM-CCSD calculation associated with the evaluation, storage, and processing of the H¯ab,cd elements in the similarity-transformed Hamiltonian involving four virtual orbital labels. The partial recovery of spin symmetry in the present algorithm reduces the storage requirement by an order of magnitude and the floating point operation count for the evaluation of the ladder-like term by a factor of three to four. EOM-CCSD calculations of excited states in the triiodide ion (I3 -) using the exact two-component Hamiltonian in combination with atomic mean-field SO integrals (X2CAMF) are reported as a validation of the implementation and also as a demonstration of the capability of the new algorithm to correlate extended virtual spaces. X2CAMF-EOM-CCSD calculations of the ground and excited states in As2, Sb2, and Bi2 are also presented and compared with the available experimental studies. An analysis based on the computed spectroscopic constants as well as the compositions of the excited-state wavefunctions strongly supports a new assignment for the lowest 2u and 0u - levels in the photoelectron spectrum of Bi2.

All-atom computations with irreversible Markov chains.

We apply the irreversible event-chain Monte Carlo (ECMC) algorithm to the simulation of dense all-atom systems with long-range Coulomb interactions. ECMC is event-driven and exactly samples the Boltzmann distribution. It neither uses time-step approximations nor spatial cutoffs on the range of the interaction potentials. Most importantly, it need not evaluate the total Coulomb potential and thus circumvents the major computational bottleneck of traditional approaches. It only requires the derivatives of the two-particle Coulomb potential, for which we discuss mutually consistent choices. ECMC breaks up the total interaction potential into factors. For particle systems made up of neutral dipolar molecules, we demonstrate the superior performance of dipole-dipole factors that do not decompose the Coulomb potential beyond the two-molecule level. We demonstrate that these long-range factors can nevertheless lead to local lifting schemes, where subsequently moved particles are mostly close to each other. For the simple point-charge water model with flexible molecules (SPC/Fw), which combines the long-ranged intermolecular Coulomb potential with hydrogen-oxygen bond-length vibrations, a flexible hydrogen-oxygen-hydrogen bond angle, and Lennard-Jones oxygen-oxygen potentials, we break up the potential into factors containing between two and six particles. For this all-atom liquid-water model, we demonstrate that the computational complexity of ECMC scales very well with the system size. This is achieved in a pure particle-particle framework, without the interpolating mesh required for the efficient implementation of other modern Coulomb algorithms. Finally, we discuss prospects and challenges for ECMC and outline several future applications.

Propagating Uncertainty in Power Flow With the Alternating Direction Method of Multipliers

In this paper, we develop an optimization-based method to propagate input and parametric uncertainty to the power flow solution. The approach is based on maximizing and minimizing quadratic approximations of the power-flow states as a function of the uncertainties subject to inequality constraints that capture all possible values the uncertain elements can take. A major computational bottleneck in such an approach is that the formulation of the quadratic approximations requires the solution of sensitivities (up to second order) from algebraic equations that are derived from the power flow equations. We demonstrate how decoupling assumptions based on the form and function of power networks can be applied to facilitate computations in this regard. The formulated quadratic programs are non-convex in general, and we adopt the Alternating Direction Method of Multipliers to solve them. Conditions for convergence in this non-convex setting are established leveraging recent advances in optimization theory. Numerical simulations for the matpower 1354-bus test system are provided to validate the accuracy and demonstrate the scalability of the approach.

Inference of species phylogenies from bi-allelic markers using pseudo-likelihood.

MotivationPhylogenetic networks represent reticulate evolutionary histories. Statistical methods for their inference under the multispecies coalescent have recently been developed. A particularly powerful approach uses data that consist of bi-allelic markers (e.g. single nucleotide polymorphism data) and allows for exact likelihood computations of phylogenetic networks while numerically integrating over all possible gene trees per marker. While the approach has good accuracy in terms of estimating the network and its parameters, likelihood computations remain a major computational bottleneck and limit the method’s applicability.ResultsIn this article, we first demonstrate why likelihood computations of networks take orders of magnitude more time when compared to trees. We then propose an approach for inference of phylogenetic networks based on pseudo-likelihood using bi-allelic markers. We demonstrate the scalability and accuracy of phylogenetic network inference via pseudo-likelihood computations on simulated data. Furthermore, we demonstrate aspects of robustness of the method to violations in the underlying assumptions of the employed statistical model. Finally, we demonstrate the application of the method to biological data. The proposed method allows for analyzing larger datasets in terms of the numbers of taxa and reticulation events. While pseudo-likelihood had been proposed before for data consisting of gene trees, the work here uses sequence data directly, offering several advantages as we discuss.Availability and implementationThe methods have been implemented in PhyloNet (http://bioinfocs.rice.edu/phylonet).