Publications

Optimization in deep learning remains poorly understood. A key difficulty is that optimizers exhibit complex oscillatory dynamics, referred to as "edge of stability," which cannot be captured by traditional optimization theory. In this paper, we show that the path taken by an oscillatory optimizer can often be captured by a central flow: a differential equation which directly models the time-averaged (i.e. smoothed) optimization trajectory. We empirically show that these central flows can predict long-term optimization trajectories for generic neural networks with a high degree of numerical accuracy. By interpreting these flows, we are able to understand how gradient descent makes progress even as the loss sometimes goes up; how adaptive optimizers ``adapt'' to the local loss landscape; and how adaptive optimizers implicitly seek out regions of weight space where they can take larger steps. These insights (and others) are not apparent from the optimizers' update rules, but are revealed by the central flows. Therefore, we believe that central flows constitute a promising tool for reasoning about optimization in deep learning.

We present a generalization of the discrete Lehmann representation (DLR) to three-point correlation and vertex functions in imaginary time and Matsubara frequency. The representation takes the form of a linear combination of judiciously chosen exponentials in imaginary time, and products of simple poles in Matsubara frequency, which are universal for a given temperature and energy cutoff. We present a systematic algorithm to generate compact sampling grids, from which the coefficients of such an expansion can be obtained by solving a linear system. We show that the explicit form of the representation can be used to evaluate diagrammatic expressions involving infinite Matsubara sums, such as polarization functions or self-energies, with controllable, high-order accuracy. This collection of techniques establishes a framework through which methods involving three-point objects can be implemented robustly, with a substantially reduced computational cost and memory footprint.

Notwithstanding the evidence against them, classical variational phase-field models continue to be used and pursued in an attempt to describe fracture nucleation in elastic brittle materials. In this context, the main objective of this paper is to provide a comprehensive review of the existing evidence against such a class of models as descriptors of fracture nucleation. To that end, a review is first given of the plethora of experimental observations of fracture nucleation in nominally elastic brittle materials under quasi-static loading conditions, as well as of classical variational phase-field models, without and with energy splits. These models are then confronted with the experimental observations. The conclusion is that they cannot possibly describe fracture nucleation in general. This because classical variational phase-field models cannot account for material strength as an independent macroscopic material property. The last part of the paper includes a brief summary of a class of phase-field models that can describe fracture nucleation. It also provides a discussion of how pervasively material strength has been overlooked in the analysis of fracture at large, as well as an outlook into the modeling of fracture nucleation beyond the basic setting of elastic brittle materials.

A new scheme is proposed to construct an (n)-times differentiable function extension of an (n)-times differentiable function defined on a smooth domain, (D,) in (d)-dimensions. The extension scheme relies on an explicit formula consisting of a linear combination of (n+1) function values in (D,) which extends the function along directions normal to the boundary. Smoothness tangent to the boundary is automatic. The performance of the scheme is illustrated by using function extension as part of a numerical solver for the Poisson equation on domains with complex geometry in both two and three dimensions. Although the cost of extending the function increases mildly with the extension order, it remains a small fraction of the overall algorithm. Moreover, the modest additional work required for high order function extensions leads to considerably more accurate solutions of the partial differential equation.

Mamba has recently emerged as a promising alternative to Transformers, demonstrating competitive performance in many language modeling tasks with linear-time computational complexity. Theoretical characterization of Mamba has largely focused on its approximation power for solving certain tasks through specific constructions. However, it remains unclear whether Mamba trained with gradient descent can learn such constructions. As a first step to address this gap, we perform a mechanistic study of simplified Mamba models on associative recall tasks. By analyzing the learned model weights and the hidden state evolution, we uncover the mechanisms used by simplified Mamba models to perform associative recall. We complement our study with theoretical analysis on the optimization dynamics of simplified Mamba models that give rise to such mechanisms.

The randomized singular value decomposition (SVD) has become a popular approach to computing cheap, yet accurate, low-rank approximations to matrices due to its efficiency and strong theoretical guarantees. Recent work by Boullé and Townsend [Found. Comput. Math., 23 (2023), pp. 709–739] presents an infinite-dimensional analogue of the randomized SVD to approximate Hilbert–Schmidt operators. However, many applications involve computing low-rank approximations to symmetric positive semi-definite matrices. In this setting, it is well established that the randomized Nyström approximation is usually preferred over the randomized SVD. This paper explores an infinite-dimensional analogue of the Nyström approximation to compute low-rank approximations to non-negative self-adjoint trace-class operators. We present an analysis of the method and, along the way, improve the existing infinite-dimensional bounds for the randomized SVD. Our analysis yields bounds on the expected value and tail bounds for the Nyström approximation error in the operator, trace, and Hilbert–Schmidt norms. Numerical experiments on integral operators arising from Gaussian process sampling and Bayesian inverse problems are used to validate the proposed infinite-dimensional Nyström algorithm.

Extracting conformational heterogeneity from cryo-electron microscopy (cryo-EM) images is particularly challenging for flexible biomolecules, where traditional 3D classification approaches often fail. Over the past few decades, advancements in experimental and computational techniques have been made to tackle this challenge, especially Bayesian-based approaches that provide physically interpretable insights into cryo-EM heterogeneity. To reduce the computational cost for Bayesian approaches, and building upon previously developed Fourier–Bessel image-representation methods, we created CryoLike, computationally efficient software for evaluating image-to-structure (or image-to-volume) likelihoods across large image data sets, packaged in a user-friendly Python workflow.

We introduce a Fourier-based fast algorithm for Gaussian process regression in low dimensions. It approximates a translationally invariant covariance kernel by complex exponentials on an equispaced Cartesian frequency grid of \(M\) nodes. This results in a weight-space \(M\times M\) system matrix with Toeplitz structure, which can thus be applied to a vector in \({\mathcal O}(M \log{M})\) operations via the fast Fourier transform (FFT), independent of the number of data points \(N\). The linear system can be set up in \({\mathcal O}(N+M \log{M})\) operations using nonuniform FFTs. This enables efficient massive-scale regression via an iterative solver, even for kernels with fat-tailed spectral densities (large \(M\)). We provide bounds on both kernel approximation and posterior mean errors. Numerical experiments for squared-exponential and Matérn kernels in one, two, and three dimensions often show 1–2 orders of magnitude acceleration over state-of-the-art rank-structured solvers at comparable accuracy. Our method allows two-dimensional Matérn-\(\frac{3}{2}\) regression from \(N=10^9\) data points to be performed in two minutes on a standard desktop, with posterior mean accuracy \(10^{-3}\). This opens up spatial statistics applications 100 times larger than previously possible.

The Gibbs sampler, also known as the coordinate hit-and-run algorithm, is a Markov chain that is widely used to draw samples from probability distributions in arbitrary dimensions. At each iteration of the algorithm, a randomly selected coordinate is resampled from the distribution that results from conditioning on all the other coordinates. We study the behavior of the Gibbs sampler on the class of log-smooth and strongly log-concave target distributions supported on ℝ

We study the problem of gradient descent learning of a single-index target function f∗(x) = σ∗(⟨x,θ⟩) under isotropic Gaussian data in Rd, where the unknown link function σ∗ : R → R has information exponent p (defined as the lowest degree in the Hermite expansion). Prior works showed that gradientbased training of neural networks can learn this target with n ≳ dΘ(p) samples, and such complexity is predicted to be necessary by the correlational statistical query lower bound. Surprisingly, we prove that a two-layer neural network optimized by an SGD-based algorithm (on the squared loss) learns f∗ with a complexity that is not governed by the information exponent. Specifically, for arbitrary polynomial single-index models, we establish a sample and runtime complexity of n ≃ T = Θ(d·polylogd), where Θ(·) hides a constant only depending on the degree of σ∗; this dimension dependence matches the information theoretic limit up to polylogarithmic factors. More generally, we show that n ≳ d(p∗−1)∨1 samples are sufficient to achieve low generalization error, where p∗ ≤ p is the generative exponent of the link function. Core to our analysis is the reuse of minibatch in the gradient computation, which gives rise to higher-order information beyond correlational queries.