Search

Randomly pivoted Cholesky (RPCholesky) is an algorithm for constructing a low-rank approximation of a positive-semidefinite matrix using a small number of columns. This paper develops an accelerated version of RPCholesky that employs block matrix computations and rejection sampling to efficiently simulate the execution of the original algorithm. For the task of approximating a kernel matrix, the accelerated algorithm can run over $40\times$ faster. The paper contains implementation details, theoretical guarantees, experiments on benchmark data sets, and an application to computational chemistry., Comment: 26 pages, 4 figures, 5 pages of supplementary material

Markov state modeling has gained popularity in various scientific fields due to its ability to reduce complex time series data into transitions between a few states. Yet, current frameworks are limited by assuming a single Markov chain describes the data, and they suffer an inability to discern heterogeneities. As a solution, this paper proposes a variational expectation-maximization algorithm that identifies a mixture of Markov chains in a time-series data set. The method is agnostic to the definition of the Markov states, whether data-driven (e.g. by spectral clustering) or based on domain knowledge. Variational EM efficiently and organically identifies the number of Markov chains and dynamics of each chain without expensive model comparisons or posterior sampling. The approach is supported by a theoretical analysis and numerical experiments, including simulated and observational data sets based on ${\tt Last.fm}$ music listening, ultramarathon running, and gene expression. The results show the new algorithm is competitive with contemporary mixture modeling approaches and powerful in identifying meaningful heterogeneities in time series data.

Rare event schemes require an approximation of the probability of the rare event as a function of system state. Finding an appropriate reaction coordinate is typically the most challenging aspect of applying a rare event scheme. Here we develop an artificial intelligence (AI) based reaction coordinate that effectively predicts which of a limited number of simulations of the Solar System will go unstable using a convolutional neural network classifier. The performance of the algorithm does not degrade significantly even 3.5 billion years before the instability. We overcome the class imbalance intrinsic to rare event problems using a combination of minority class oversampling, increased minority class weighting, and pulling multiple non-overlapping training sequences from simulations. Our success suggests that AI may provide a promising avenue for developing reaction coordinates without detailed theoretical knowledge of the system., Comment: 3 pages, 1 figure

Recently, a class of algorithms combining classical fixed point iterations with repeated random sparsification of approximate solution vectors has been successfully applied to eigenproblems with matrices as large as $10^{108} \times 10^{108}$. So far, a complete mathematical explanation for their success has proven elusive. Additionally, the methods have not been extended to linear system solves. In this paper we propose a new scheme based on repeated random sparsification that is capable of solving linear systems in extremely high dimensions. We provide a complete mathematical analysis of this new algorithm. Our analysis establishes a faster-than-Monte Carlo convergence rate and justifies use of the scheme even when the solution vector itself is too large to store., Comment: 27 pages, 2 figures

This survey explores modern approaches for computing low-rank approximations of high-dimensional matrices by means of the randomized SVD, randomized subspace iteration, and randomized block Krylov iteration. The paper compares the procedures via theoretical analyses and numerical studies to highlight how the best choice of algorithm depends on spectral properties of the matrix and the computational resources available. Despite superior performance for many problems, randomized block Krylov iteration has not been widely adopted in computational science. The paper strengthens the case for this method in three ways. First, it presents new pseudocode that can significantly reduce computational costs. Second, it provides a new analysis that yields simple, precise, and informative error bounds. Last, it showcases applications to challenging scientific problems, including principal component analysis for genetic data and spectral clustering for molecular dynamics data., Comment: 60 pages, 14 figures

Mercury's orbit can destabilize, generally resulting in a collision with either Venus or the Sun. Chaotic evolution can cause g1 to decrease to the approximately constant value of g5 and create a resonance. Previous work has approximated the variation in g1 as stochastic diffusion, which leads to a phenomological model that can reproduce the Mercury instability statistics of secular and N-body models on timescales longer than 10 Gyr. Here we show that the diffusive model underpredicts the Mercury instability probability by a factor of 3-10,000 on timescales less than 5 Gyr, the remaining lifespan of the Solar System. This is because g1 exhibits larger variations on short timescales than the diffusive model would suggest. To better model the variations on short timescales, we build a new subdiffusive phenomological model for g1. Subdiffusion is similar to diffusion but exhibits larger displacements on short timescales and smaller displacements on long timescales. We choose model parameters based on the behavior of the g1 trajectories in the N-body simulations, leading to a tuned model that can reproduce Mercury instability statistics from 1-40 Gyr. This work motivates fundamental questions in Solar System dynamics: Why does subdiffusion better approximate the variation in g1 than standard diffusion? Why is there an upper bound on g1, but not a lower bound that would prevent it from reaching g5?, Comment: accepted at ApJ

This paper investigates two randomized preconditioning techniques for solving kernel ridge regression (KRR) problems with a medium to large number of data points ($10^4 \leq N \leq 10^7$), and it introduces two new methods with state-of-the-art performance. The first method, RPCholesky preconditioning, accurately solves the full-data KRR problem in $O(N^2)$ arithmetic operations, assuming sufficiently rapid polynomial decay of the kernel matrix eigenvalues. The second method, KRILL preconditioning, offers an accurate solution to a restricted version of the KRR problem involving $k \ll N$ selected data centers at a cost of $O((N + k^2) k \log k)$ operations. The proposed methods solve a broad range of KRR problems, making them ideal for practical applications., Comment: 29 pages, 11 figures

The implicit trace estimation problem asks for an approximation of the trace of a square matrix, accessed via matrix-vector products (matvecs). This paper designs new randomized algorithms, XTrace and XNysTrace, for the trace estimation problem by exploiting both variance reduction and the exchangeability principle. For a fixed budget of matvecs, numerical experiments show that the new methods can achieve errors that are orders of magnitude smaller than existing algorithms, such as the Girard-Hutchinson estimator or the Hutch++ estimator. A theoretical analysis confirms the benefits by offering a precise description of the performance of these algorithms as a function of the spectrum of the input matrix. The paper also develops an exchangeable estimator, XDiag, for approximating the diagonal of a square matrix using matvecs., Comment: 31 pages, 8 figures

A major problem in numerical weather prediction (NWP) is the estimation of high-dimensional covariance matrices from a small number of samples. Maximum likelihood estimators cannot provide reliable estimates when the overall dimension is much larger than the number of samples. Fortunately, NWP practitioners have found ingenious ways to boost the accuracy of their covariance estimators by leveraging the assumption that the correlations decay with spatial distance. In this work, Bayesian statistics is used to provide a new justification and analysis of the practical NWP covariance estimators. The Bayesian framework involves manipulating distributions over symmetric positive definite matrices, and it leads to two main findings: (i) the commonly used "hybrid estimator" for the covariance matrix has a naturally Bayesian interpretation; (ii) the very commonly used "Schur product estimator" is not Bayesian, but it can be studied and understood within the Bayesian framework. As practical implications, the Bayesian framework shows how to reduce the amount of tuning required for covariance estimation, and it suggests that efficient covariance estimation should be rooted in understanding and penalizing conditional correlations, rather than correlations., Comment: 20 pages, 4 figures

The long-term stability of the Solar System is an issue of significant scientific and philosophical interest. The mechanism leading to instability is Mercury's eccentricity being pumped up so high that Mercury either collides with Venus or is scattered into the Sun. Previously, only three five-billion-year $N$-body ensembles of the Solar System with thousands of simulations have been run to assess long-term stability. We generate two additional ensembles, each with 2750 members, and make them publicly available at \texttt{https://archive.org/details/@dorianabbot}. We find that accurate Mercury instability statistics can be obtained by (1) including only the Sun and the 8 planets, (2) using a simple Wisdom-Holman scheme without correctors, (3) using a basic representation of general relativity, and (4) using a time step of 3.16 days. By combining our Solar System ensembles with previous ensembles we form a 9,601-member ensemble of ensembles. In this ensemble of ensembles, the logarithm of the frequency of a Mercury instability event increases linearly with time between 1.3 and 5 Gyr, suggesting that a single mechanism is responsible for Mercury instabilities in this time range and that this mechanism becomes more active as time progresses. Our work provides a robust estimate of Mercury instability statistics over the next five billion years, outlines methodologies that may be useful for exoplanet system investigations, and provides two large ensembles of publicly available Solar System integrations that can serve as testbeds for theoretical ideas as well as training sets for artificial intelligence schemes., Comment: accepted at ApJ

The use of neural network parametrizations to represent the ground state in variational Monte Carlo (VMC) calculations has generated intense interest in recent years. However, as we demonstrate in the context of the periodic Heisenberg spin chain, this approach can produce unreliable wave function approximations. One of the most obvious signs of failure is the occurrence of random, persistent spikes in the energy estimate during training. These energy spikes are caused by regions of configuration space that are over-represented by the wave function density, which are called ``spurious modes'' in the machine learning literature. After exploring these spurious modes in detail, we demonstrate that a collective-variable-based penalization yields a substantially more robust training procedure, preventing the formation of spurious modes and improving the accuracy of energy estimates. Because the penalization scheme is cheap to implement and is not specific to the particular model studied here, it can be extended to other applications of VMC where a reasonable choice of collective variable is available., Comment: 12 pages, 13 figures

The randomly pivoted partial Cholesky algorithm (RPCholesky) computes a factorized rank-k approximation of an N x N positive-semidefinite (psd) matrix. RPCholesky requires only (k + 1) N entry evaluations and O(k^2 N) additional arithmetic operations, and it can be implemented with just a few lines of code. The method is particularly useful for approximating a kernel matrix. This paper offers a thorough new investigation of the empirical and theoretical behavior of this fundamental algorithm. For matrix approximation problems that arise in scientific machine learning, experiments show that RPCholesky matches or beats the performance of alternative algorithms. Moreover, RPCholesky provably returns low-rank approximations that are nearly optimal. The simplicity, effectiveness, and robustness of RPCholesky strongly support its use in scientific computing and machine learning applications., Comment: 40 pages, 4 figures

We present a stable and systematically improvable quantum Monte Carlo (QMC) approach to calculating excited-state energies, which we implement using our fast randomized iteration method for the full configuration interaction problem (FCI-FRI). Unlike previous excited-state quantum Monte Carlo methods, our approach, which is an asymmetric variant of subspace iteration, avoids the use of dot products of random vectors and instead relies upon trial vectors to maintain orthogonality and estimate eigenvalues. By leveraging recent advances, we apply our method to calculate ground- and excited-state energies of strongly correlated molecular systems in large active spaces, including the carbon dimer with 8 electrons in 108 orbitals (8e,108o), an oxo-Mn(salen) transition metal complex (28e,28o), ozone (18e,87o), and butadiene (22e,82o). In the majority of these test cases, our approach yields total excited-state energies that agree with those from state-of-the-art methods -- including heat-bath CI, the density matrix renormalization group approach, and FCIQMC -- to within sub-milliHartree accuracy. In all cases, estimated excitation energies agree to within about 0.1 eV., Comment: 16 pages, 5 figures, 3 tables

Atmospheric regime transitions are highly impactful as drivers of extreme weather events, but pose two formidable modeling challenges: predicting the next event (weather forecasting), and characterizing the statistics of events of a given severity (the risk climatology). Each event has a different duration and spatial structure, making it hard to define an objective "average event." We argue here that transition path theory (TPT), a stochastic process framework, is an appropriate tool for the task. We demonstrate TPT's capacities on a wave-mean flow model of sudden stratospheric warmings (SSWs) developed by Holton and Mass (1976), which is idealized enough for transparent TPT analysis but complex enough to demonstrate computational scalability. Whereas a recent article (Finkel et al. 2021) studied near-term SSW predictability, the present article uses TPT to link predictability to long-term SSW frequency. This requires not only forecasting forward in time from an initial condition, but also \emph{backward in time} to assess the probability of the initial conditions themselves. TPT enables one to condition the dynamics on the regime transition occurring, and thus visualize its physical drivers with a vector field called the \emph{reactive current}. The reactive current shows that before an SSW, dissipation and stochastic forcing drive a slow decay of vortex strength at lower altitudes. The response of upper-level winds is late and sudden, occurring only after the transition is almost complete from a probabilistic point of view. This case study demonstrates that TPT quantities, visualized in a space of physically meaningful variables, can help one understand the dynamics of regime transitions., Comment: 18 pages, 7 figures (main text), 19 pages, 1 figure (supplement). Accepted for publication in the Journal of the Atmospheric Sciences

Variational Monte Carlo (VMC) is an approach for computing ground-state wavefunctions that has recently become more powerful due to the introduction of neural network-based wavefunction parametrizations. However, efficiently training neural wavefunctions to converge to an energy minimum remains a difficult problem. In this work, we analyze optimization and sampling methods used in VMC and introduce alterations to improve their performance. First, based on theoretical convergence analysis in a noiseless setting, we motivate a new optimizer that we call the Rayleigh-Gauss-Newton method, which can improve upon gradient descent and natural gradient descent to achieve superlinear convergence at no more than twice the computational cost. Second, in order to realize this favorable comparison in the presence of stochastic noise, we analyze the effect of sampling error on VMC parameter updates and experimentally demonstrate that it can be reduced by the parallel tempering method. In particular, we demonstrate that RGN can be made robust to energy spikes that occur when the sampler moves between metastable regions of configuration space. Finally, putting theory into practice, we apply our enhanced optimization and sampling methods to the transverse-field Ising and XXZ models on large lattices, yielding ground-state energy estimates with remarkably high accuracy after just 200 parameter updates., Comment: 12 pages, 8 figures

Due to the chaotic nature of planetary dynamics, there is a non-zero probability that Mercury's orbit will become unstable in the future. Previous efforts have estimated the probability of this happening between 3 and 5 billion years in the future using a large number of direct numerical simulations with an N-body code, but were not able to obtain accurate estimates before 3 billion years in the future because Mercury instability events are too rare. In this paper we use a new rare event sampling technique, Quantile Diffusion Monte Carlo (QDMC), to estimate that the probability of a Mercury instability event in the next 2 billion years is approximately $10^{-4}$ in the REBOUND N-body code. We show that QDMC provides unbiased probability estimates at a computational cost of up to 100 times less than direct numerical simulation. QDMC is easy to implement and could be applied to many problems in planetary dynamics in which it is necessary to estimate the probability of a rare event., Comment: Abbot et al 2021 ApJ 923 236

Traditional numerical methods for calculating matrix eigenvalues are prohibitively expensive for high-dimensional problems. Iterative random sparsification methods allow for the estimation of a single dominant eigenvalue at reduced cost by leveraging repeated random sampling and averaging. We present a general approach to extending such methods for the estimation of multiple eigenvalues and demonstrate its performance for several benchmark problems in quantum chemistry., Comment: 31 pages, 8 figures

Rare events arising in nonlinear atmospheric dynamics remain hard to predict and attribute. We address the problem of forecasting rare events in a prototypical example, Sudden Stratospheric Warmings (SSWs). Approximately once every other winter, the boreal stratospheric polar vortex rapidly breaks down, shifting midlatitude surface weather patterns for months. We focus on two key quantities of interest: the probability of an SSW occurring, and the expected lead time if it does occur, as functions of initial condition. These \emph{optimal forecasts} concretely measure the event's progress. Direct numerical simulation can estimate them in principle, but is prohibitively expensive in practice: each rare event requires a long integration to observe, and the cost of each integration grows with model complexity. We describe an alternative approach using integrations that are \emph{short} compared to the timescale of the warming event. We compute the probability and lead time efficiently by solving equations involving the transition operator, which encodes all information about the dynamics. We relate these optimal forecasts to a small number of interpretable physical variables, suggesting optimal measurements for forecasting. We illustrate the methodology on a prototype SSW model developed by Holton and Mass (1976) and modified by stochastic forcing. While highly idealized, this model captures the essential nonlinear dynamics of SSWs and exhibits the key forecasting challenge: the dramatic separation in timescales between a single event and the return time between successive events. Our methodology is designed to fully exploit high-dimensional data from models and observations, and has the potential to identify detailed predictors of many complex rare events in meteorology., Comment: 26 pages, 7 figures, major revision after original. Accepted to Monthly Weather Review, American Meteorological Society

We explore whether splitting and killing methods can improve the accuracy of Markov chain Monte Carlo (MCMC) estimates of rare event probabilities, and we make three contributions. First, we prove that "weighted ensemble" is the only splitting and killing method that provides asymptotically consistent estimates when combined with MCMC. Second, we prove a lower bound on the asymptotic variance of weighted ensemble's estimates. Third, we give a constructive proof and numerical examples to show that weighted ensemble can approach this optimal variance bound, in many cases reducing the variance of MCMC estimates by multiple orders of magnitude., Comment: 30 pages, 9 figures

We introduce a new Markov chain Monte Carlo (MCMC) sampler for infinite-dimensional inverse problems. Our new sampler is based on the affine invariant ensemble sampler, which uses interacting walkers to adapt to the covariance structure of the target distribution. We extend this ensemble sampler for the first time to infinite-dimensional function spaces, yielding a highly efficient gradient-free MCMC algorithm. Because our new ensemble sampler does not require gradients or posterior covariance estimates, it is simple to implement and broadly applicable.

One approach to analyzing the dynamics of a physical system is to search for long-lived patterns in its motions. This approach has been particularly successful for molecular dynamics data, where slowly decorrelating patterns can indicate large-scale conformational changes. Detecting such patterns is the central objective of the variational approach to conformational dynamics (VAC), as well as the related methods of time-lagged independent component analysis and Markov state modeling. In VAC, the search for slowly decorrelating patterns is formalized as a variational problem solved by the eigenfunctions of the system's transition operator. VAC computes solutions to this variational problem by optimizing a linear or nonlinear model of the eigenfunctions using time series data. Here, we build on VAC's success by addressing two practical limitations. First, VAC can give poor eigenfunction estimates when the lag time parameter is chosen poorly. Second, VAC can overfit when using flexible parameterizations such as artificial neural networks with insufficient regularization. To address these issues, we propose an extension that we call integrated VAC (IVAC). IVAC integrates over multiple lag times before solving the variational problem, making its results more robust and reproducible than VAC's., Comment: 16 pages, 7 figures

Dynamical spectral estimation is a well-established numerical approach for estimating eigenvalues and eigenfunctions of the Markov transition operator from trajectory data. Although the approach has been widely applied in biomolecular simulations, its error properties remain poorly understood. Here we analyze the error of a dynamical spectral estimation method called "the variational approach to conformational dynamics" (VAC). We bound the approximation error and estimation error for VAC estimates. Our analysis establishes VAC's convergence properties and suggests new strategies for tuning VAC to improve accuracy., Comment: 34 pages, 7 figures

We present three modifications to our recently introduced fast randomized iteration method for full configuration interaction (FCI-FRI) and investigate their effects on the method's performance for Ne, H$_2$O, and N$_2$. The initiator approximation, originally developed for full configuration interaction quantum Monte Carlo, significantly reduces statistical error in FCI-FRI when few samples are used in compression operations, enabling its application to larger chemical systems. The semi-stochastic extension, which involves exactly preserving a fixed subset of elements in each compression, improves statistical efficiency in some cases but reduces it in others. We also developed a new approach to sampling excitations that yields consistent improvements in statistical efficiency and reductions in computational cost. We discuss possible strategies based on our findings for improving the performance of stochastic quantum chemistry methods more generally., Comment: 13 pages, 5 figures

We introduce a family of methods for the full configuration interaction problem in quantum chemistry, based on the fast randomized iteration (FRI) framework [L.-H. Lim and J. Weare, SIAM Rev. 59, 547 (2017)]. These methods, which we term "FCI-FRI," stochastically impose sparsity during iterations of the power method and can be viewed as a generalization of full configuration interaction quantum Monte Carlo (FCIQMC) without walkers. In addition to the multinomial scheme commonly used to sample excitations in FCIQMC, we present a systematic scheme where excitations are not sampled independently. Performing ground-state calculations on five small molecules at fixed cost, we find that the systematic FCI-FRI scheme is 11 to 45 times more statistically efficient than the multinomial FCI-FRI scheme, which is in turn 1.4 to 178 times more statistically efficient than the original FCIQMC algorithm., Comment: 19 pages, 7 figures

Direct computer simulation of intense tropical cyclones (TCs) in weather models is limited by computational expense. Intense TCs are rare and have small-scale structures, making it difficult to produce large ensembles of storms at high resolution. Further, models often fail to capture the process of rapid intensification, which is a distinguishing feature of many intense TCs. Understanding rapid intensification is especially important in the context of global warming, which may increase the frequency of intense TCs. To better leverage computational resources for the study of rapid intensification, we introduce an action minimization algorithm applied to the WRF and WRFPLUS models. Action minimization nudges the model into forming more intense TCs than it otherwise would; it does so via the maximum likelihood path in a stochastic formulation of the model, thereby allowing targeted study of intensification mechanisms. We apply action minimization to simulations of Hurricanes Danny (2015) and Fred (2009) at 6 km resolution to demonstrate that the algorithm consistently intensifies TCs via physically plausible pathways. We show an approximately ten-fold computational savings using action minimization to study the tail of the TC intensification distribution. Further, for Hurricanes Danny and Fred, action minimization produces perturbations that preferentially reduce low-level shear as compared to upper-level shear, at least above a threshold of approximately $4 \mathrm{\ m \ s^{-1}}$. We also demonstrate that asymmetric, time-dependent patterns of heating can cause significant TC intensification beyond symmetric, azimuthally-averaged heating and find a regime of non-linear response to asymmetric heating that has not been extensively studied in previous work.

Extreme mesoscale weather, including tropical cyclones, squall lines, and floods, can be enormously damaging and yet challenging to simulate; hence, there is a pressing need for more efficient simulation strategies. Here we present a new rare event sampling algorithm called Quantile Diffusion Monte Carlo (Quantile DMC). Quantile DMC is a simple-to-use algorithm that can sample extreme tail behavior for a wide class of processes. We demonstrate the advantages of Quantile DMC compared to other sampling methods and discuss practical aspects of implementing Quantile DMC. To test the feasibility of Quantile DMC for extreme mesoscale weather, we sample extremely intense realizations of two historical tropical cyclones, 2010 Hurricane Earl and 2015 Hurricane Joaquin. Our results demonstrate Quantile DMC's potential to provide low-variance extreme weather statistics while highlighting the work that is necessary for Quantile DMC to attain greater efficiency in future applications., Comment: 18 pages, 9 figures

Sequential Monte Carlo (SMC) is a class of algorithms that approximate high-dimensional expectations of a Markov chain. SMC algorithms typically include a resampling step. There are many possible ways to resample, but the relative advantages of different resampling schemes remains poorly understood. Here, a theoretical framework for comparing resampling schemes is presented. The framework uses resampling matrices to provide a simple description for the SMC resampling step. The framework identifies the matrix resampling scheme that gives the lowest possible error. The framework leads to new asymptotic error formulas that can be used to compare different resampling schemes., Comment: 36 pages, 3 figures

Provided that sufficient resources are deployed, we can look forward to an extraordinary future in which we will characterize potentially habitable planets. Until now, we have had to base interpretations of observations on habitability hypotheses that have remained untested. To test these theories observationally, we propose a statistical comparative planetology approach to questions of planetary habitability. The key objective of this approach will be to make quick and cheap measurements of critical planetary characteristics on a large sample of exoplanets, exploiting statistical marginalization to answer broad habitability questions. This relaxes the requirement of obtaining multiple types of data for a given planet, as it allows us to test a given hypothesis from only one type of measurement using the power of an ensemble. This approach contrasts with a "systems science" approach, where a few planets would be extensively studied with many types of measurements. A systems science approach is associated with a number of difficulties which may limit overall scientific return, including: the limited spectral coverage and noise of instruments, the diversity of exoplanets, and the extensive list of potential false negatives and false positives. A statistical approach could also be complementary to a systems science framework by providing context to interpret extensive measurements on planets of particular interest. We strongly recommend future missions with a focus on exoplanet characterization, and with the capability to study large numbers of planets in a homogenous way, rather than exclusively small, intense studies directed at a small sample of planets., Comment: White paper submitted in response to the solicitation of feedback for the "2020 Decadal Survey" by the National Academy of Sciences