The randomized Kaczmarz (RK) method is a well-known approach for solving linear least-squares problems with a large number of rows. RK accesses and processes just one row at a time, leading to exponentially fast convergence for consistent linear systems. However, RK fails to converge to the least-squares solution for inconsistent systems. This work presents a simple fix: average the RK iterates produced in the tail part of the algorithm. The proposed tail-averaged randomized Kaczmarz (TARK) converges for both consistent and inconsistent least-squares problems at a polynomial rate, which is known to be optimal for any row-access method. An extension of TARK also leads to efficient solutions for ridge-regularized least-squares problems.

Two widely used randomized algorithms are the sketch-and-solve method for least-squares regression and the randomized SVD for low-rank approximation. These algorithms apply a random embedding to compress a target matrix, and they perform computations on the compressed matrix to save computational cost. This paper asks, what is the optimal random embedding in these algorithms? Also, what is the sharpest possible error bound for the optimal embedding? The paper proves that a random orthonormal matrix is minimax optimal for the sketch-and-solve algorithm while any rotation-invariant embedding is minimax optimal for the randomized SVD. Following these results, the paper obtains the best possible error bounds for sketched least-squares and the randomized SVD. Last, empirical experiments provide evidence of universality phenomena, in which several random embeddings lead to similar accuracy to the optimal embeddings in practice.

A significant challenge in molecular dynamics (MD) simulations is ensuring that sampled configurations converge to the equilibrium or nonequilibrium stationary distribution of interest. Lack of convergence constrains the estimation of free energies and of rates and mechanisms for molecular transitions. Here, we introduce the "Randomized ITErative trajectory reWeighting" (RiteWeight) algorithm to estimate a stationary distribution from unconverged simulation data. This method iteratively reweights trajectory segments in a self-consistent way by solving for the stationary distribution of a Markov state model (MSM), updating segment weights, and employing a new random clustering in each iteration. The repeated clustering mitigates the configuration-space discretization error inherent in existing trajectory reweighting techniques and yields quasi-continuous configuration-space distributions. RiteWeight accurately recovers the stationary distribution even without requiring the Markov property at the cluster level. We present mathematical analysis of the RiteWeight fixed point. We empirically validate the method using both synthetic MD Trp-cage trajectories, for which the stationary solution is exactly calculable, and standard atomistic MD Trp-cage trajectories, which are extracted from a long reference simulation. In both test systems, RiteWeight corrects flawed distributions and generates accurate observables for equilibrium and nonequilibrium steady states. The results highlight the value of correcting the underlying trajectory distribution rather than using a standard MSM.

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 . So far, a complete mathematical explanation for this success has proven elusive.The family of methods has not yet been extended to the important case of linear system solves. In this paper, we propose a new scheme based on repeated random sparsification that is capable of solving sparse linear systems in arbitrarily 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 is too large to store as a dense vector.

Markov state models (MSMs) have become ubiquitous tools for analyzing molecular dynamics (MD) simulations because of their simple, powerful premise: although complete MD sampling may be impossible, the MSM can "stitch together" transition probabilities derived from local sampling to provide a global picture of kinetics and mechanisms. In the standard MSM framework, the available MD data is organized into a single transition matrix, which is then used to estimate all observables at a lag time chosen so the coarse-grained dynamics are approximately Markovian. This approach leads to avoidable model bias and motivates long lag times that obscure short-timescale processes of interest. In contrast, this paper shows how to obtain unbiased coarse-grained observables at any fixed lag time and for any fixed coarse-graining in the limit of infinite, properly weighted data. The central idea is to replace the single-matrix framework with two transition matrices -- one representing equilibrium dynamics and another representing source-sink recycling dynamics -- and use the correct matrix or matrices to estimate the matched dynamical observables.

Markov state modeling has gained popularity in various scientific fields since it reduces complex time-series data sets into transitions between a few states. Yet common Markov state modeling frameworks assume a single Markov chain describes the data, so they suffer from an inability to discern heterogeneities. As an alternative, this paper models time-series data using a mixture of Markov chains, and it automatically determines the number of mixture components using the variational expectation-maximization algorithm. Variational EM simultaneously identifies the number of Markov chains and the dynamics of each chain without expensive model comparisons or posterior sampling. As a theoretical contribution, this paper identifies the natural limits of Markov chain mixture modeling by proving a lower bound on the classification error. This paper then presents numerical experiments where variational EM achieves performance consistent with the theoretically optimal error scaling. The experiments are based on synthetic and observational data sets including Last.fm music listening, ultramarathon running, and gene expression. In each of the three data sets, variational EM leads to the identification of meaningful heterogeneities.

Weighted ensemble (WE) is an enhanced path-sampling method that is conceptually simple, widely applicable, and statistically exact. In a WE simulation, an ensemble of trajectories is periodically pruned or replicated to enhance sampling of rare transitions and improve estimation of mean first passage times (MFPTs). However, poor choices of the parameters governing pruning and replication can lead to high-variance MFPT estimates. Our previous work [J. Chem. Phys. 158, 014108 (2023)] presented an optimal WE parameterization strategy and applied it in low-dimensional example systems. The strategy harnesses estimated local MFPTs from different initial configurations to a single target state. In the present work, we apply the optimal parameterization strategy to more challenging, high-dimensional molecular models, namely, synthetic molecular dynamics (MD) models of Trp-cage folding and unfolding, as well as atomistic MD models of NTL9 folding in high-friction and low-friction continuum solvents. In each system we use WE to estimate the MFPT for folding or unfolding events. We show that the optimal parameterization reduces the variance of MFPT estimates in three of four systems, with dramatic improvement in the most challenging atomistic system. Overall, the parameterization strategy improves the accuracy and reliability of WE estimates for the kinetics of biophysical processes.

Neural network parametrizations have increasingly been used to represent the ground and excited states in variational Monte Carlo (VMC) with promising results. However, traditional VMC methods only optimize the wave function in regions of peak probability. The wave function is uncontrolled in the tails of the probability distribution, which can limit the accuracy of the trained wave function. To improve the approximation accuracy in the probability tails, this paper interprets VMC as a gradient flow in the space of wave functions, followed by a projection step. From this perspective, arbitrary probability distributions can be used in the projection step, allowing the user to prioritize accuracy in different regions of state space. Motivated by this theoretical perspective, the paper tests a weighted VMC method on the antiferromagnetic Heisenberg model for a periodic spin chain. Compared to traditional VMC, weighted VMC reduces the error in the ground state energy by a factor of 2, and it reduces the errors in the local energies away from the mode by large factors of 10^{2}–10^{4}.

The partial pivoted Cholesky approximation accurately represents matrices that are close to being low-rank. Meanwhile, the Vecchia approximation accurately represents matrices with inverse Cholesky factors that are close to being sparse. What happens if a partial Cholesky approximation is combined with a Vecchia approximation of the residual? This paper shows how the sum is exactly a Vecchia approximation of the original matrix with an augmented sparsity pattern. Thus, Vecchia approximations subsume a class of existing matrix approximations and have broad applicability.

Neural network parametrizations have increasingly been used to represent the ground and excited states in variational Monte Carlo (VMC) with promising results. However, traditional VMC methods only optimize the wave function in regions of peak probability. The wave function is uncontrolled in the tails of the probability distribution, which can limit the accuracy of the trained wavefunction approximation. To improve the approximation accuracy in the probability tails, this paper interprets VMC as a gradient flow in the space of wave functions, followed by a projection step. From this perspective, arbitrary probability distributions can be used in the projection step, allowing the user to prioritize accuracy in different regions of state space. Motivated by this theoretical perspective, the paper tests a new weighted VMC method on the antiferromagnetic Heisenberg model for a periodic spin chain. Compared to traditional VMC, weighted VMC reduces the error in the ground state energy by a factor of 2 and it reduces the errors in the local energies away from the mode by large factors of 10^2–10^4.

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 https://archive.org/details/@dorianabbot . We find that accurate Mercury instability statistics can be obtained by (1) including only the Sun and the eight 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 9601-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 test beds for theoretical ideas as well as training sets for artificial intelligence schemes.

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.

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.

Markov state modeling has gained popularity in various scientific fields since it reduces complex time-series data sets into transitions between a few states. Yet common Markov state modeling frameworks assume a single Markov chain describes the data, so they suffer from an inability to discern heterogeneities. As an alternative, this paper models time-series data using a mixture of Markov chains, and it automatically determines the number of mixture components using the variational expectation-maximization algorithm.Variational EM simultaneously identifies the number of Markov chains and the dynamics of each chain without expensive model comparisons or posterior sampling. As a theoretical contribution, this paper identifies the natural limits of Markov state mixture modeling by proving a lower bound on the classification error. It then presents numerical experiments where variational EM achieves performance consistent with the theoretically optimal error scaling. The experiments are based on synthetic and observational data sets including Last.fm music listening, ultramarathon running, and gene expression. In each of the three data sets, variational EM leads to the identification of meaningful heterogeneities.

In the study of stochastic systems, the committor function describes the probability that a system starting from an initial configuration x will reach a set B before a set A. This paper introduces an efficient and interpretable algorithm for approximating the committor, called the "fast committor machine" (FCM). The FCM uses simulated trajectory data to build a kernel-based model of the committor. The kernel function is constructed to emphasize low-dimensional subspaces that optimally describe the A to B transitions. The coefficients in the kernel model are determined using randomized linear algebra, leading to a runtime that scales linearly with the number of data points. In numerical experiments involving a triple-well potential and alanine dipeptide, the FCM yields higher accuracy and trains more quickly than a neural network with the same number of parameters. The FCM is also more interpretable than the neural net.

Mercury's orbit can destabilize, resulting in a collision with either Venus or the Sun. Chaotic evolution can cause g_1 to decrease to the approximately constant value of g_5 and create a resonance. Previous work has approximated the variation in g_1 as stochastic diffusion, which leads to a 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 g_1 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 model for g_1 including a quadratic spring potential above a certain value of g_1, which we refer to as a soft upper boundary. 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 short-time behavior of the g_1 trajectories in the N-body simulations, leading to a tuned model that can reproduce Mercury instability statistics from 1-40 Gyr. This work motivates several questions in planetary dynamics: Why does subdiffusion better approximate the variation in g_1 than standard diffusion? Why is a soft upper boundary condition on g_1 an appropriate approximation? Why is there an upper bound on g_1, but not a lower bound that would prevent it from reaching g_5?

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.

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.

Weighted ensemble (WE) simulation runs multiple weighted trajectories that are periodically replicated or pruned to spread computing effort across conformation space and exploit parallelization. The method can provide enhanced sampling power as well as statistically unbiased estimates for the rate constant of rare events, independent of hyperparameter choice. However, the run-to-run variance of WE simulations can be greatly affected by the choice of hyperparameters: the bins, which partition the configurational space, and the allocation, which determines how many trajectories are to be maintained in each bin.