Christopher Musco

NeurIPS Conference 2025 Conference Paper

Distance Adaptive Beam Search for Provably Accurate Graph-Based Nearest Neighbor Search

• Yousef Al-Jazzazi
• Haya Diwan
• Jinrui Gou
• Cameron Musco
• Christopher Musco
• Torsten Suel

Nearest neighbor search is central in machine learning, information retrieval, and databases. For high-dimensional datasets, graph-based methods such as HNSW, DiskANN, and NSG have become popular thanks to their empirical accuracy and efficiency. These methods construct a directed graph over the dataset and perform beam search on the graph to find nodes close to a given query. While significant work has focused on practical refinements and theoretical understanding of graph-based methods, many questions remain. We propose a new distance-based termination condition for beam search to replace the commonly used condition based on beam width. We prove that, as long as the search graph is navigable, our resulting Adaptive Beam Search method is guaranteed to approximately solve the nearest-neighbor problem, establishing a connection between navigability and the performance of graph-based search. We also provide extensive experiments on our new termination condition for both navigable graphs and approximately navigable graphs used in practice, such as HNSW and Vamana graphs. We find that Adaptive Beam Search outperforms standard beam search over a range of recall values, data sets, graph constructions, and target number of nearest neighbors. It thus provides a simple and practical way to improve the performance of popular methods.

SODA Conference 2025 Conference Paper

Faster Linear Systems and Matrix Norm Approximation via Multi-level Sketched Preconditioning

• Michal Derezinski
• Christopher Musco
• Jiaming Yang

SODA Conference 2025 Conference Paper

Improved Spectral Density Estimation via Explicit and Implicit Deflation

• Rajarshi Bhattacharjee
• Rajesh Jayaram
• Cameron Musco
• Christopher Musco
• Archan Ray

ICLR Conference 2025 Conference Paper

Matrix Product Sketching via Coordinated Sampling

• Majid Daliri
• Juliana Freire
• Danrong Li
• Christopher Musco

We revisit the well-studied problem of approximating a matrix product, $\bv{A}^T\bv{B}$, based on small space sketches $\mathcal{S}(\bv{A})$ and $\mathcal{S}(\bv{B})$ of $\bv{A} \in \R^{n \times d}$ and $\bv{B}\in \R^{n \times m}$. We are interested in the setting where the sketches must be computed independently of each other, except for the use of a shared random seed. We prove that, when $\bv{A}$ and $\bv{B}$ are sparse, methods based on \emph{coordinated random sampling} can outperform classical linear sketching approaches, like Johnson-Lindenstrauss Projection or CountSketch. For example, to obtain Frobenius norm error $\epsilon\|\bv{A}\|_F\|\bv{B}\|_F$, coordinated sampling requires sketches of size $O(s/\epsilon^2)$ when $\bv{A}$ and $\bv{B}$ have at most $s \leq d,m$ non-zeros per row. In contrast, linear sketching leads to sketches of size $O(d/\epsilon^2)$ and $O(m/\epsilon^2)$ for $\bv{A}$ and $\bv{B}$. We empirically evaluate our approach on two applications: 1) distributed linear regression in databases, a problem motivated by tasks like dataset discovery and augmentation, and 2) approximating attention matrices in transformer-based language models. In both cases, our sampling algorithms yield an order of magnitude improvement over linear sketching.

SODA Conference 2025 Conference Paper

Near-optimal hierarchical matrix approximation from matrix-vector products

• Tyler Chen
• Feyza Duman Keles
• Diana Halikias
• Cameron Musco
• Christopher Musco
• David Persson

ICLR Conference 2025 Conference Paper

Provably Accurate Shapley Value Estimation via Leverage Score Sampling

• Christopher Musco
• R. Teal Witter

Originally introduced in game theory, Shapley values have emerged as a central tool in explainable machine learning, where they are used to attribute model predictions to specific input features. However, computing Shapley values exactly is expensive: for a model with $n$ features, $O(2^n)$ model evaluations are necessary. To address this issue, approximation algorithms are widely used. One of the most popular is the Kernel SHAP algorithm, which is model agnostic and remarkably effective in practice. However, to the best of our knowledge, Kernel SHAP has no strong non-asymptotic complexity guarantees. We address this issue by introducing *Leverage SHAP*, a light-weight modification of Kernel SHAP that provides provably accurate Shapley value estimates with just $O(n\log n)$ model evaluations. Our approach takes advantage of a connection between Shapley value estimation and agnostic active learning by employing *leverage score sampling*, a powerful regression tool. Beyond theoretical guarantees, we find that Leverage SHAP achieves an approximately 50% reduction in error compared to the highly optimized implementation of Kernel SHAP in the widely used SHAP library [Lundberg & Lee, 2017].

NeurIPS Conference 2025 Conference Paper

Regression-adjusted Monte Carlo Estimators for Shapley Values and Probabilistic Values

• R. Teal Witter
• Yurong Liu
• Christopher Musco

With origins in game-theory, probabilistic values like Shapley values, Banzhaf values, and semi-values have emerged as a central tool in explainable AI. They are used for feature attribution, data attribution, data valuation, and more. Since all of these values require exponential time to compute exactly, research has focused on efficient approximation methods using two techniques: Monte Carlo sampling and linear regression formulations. In this work, we present a new way of combining both of these techniques. Our approach is more flexible than prior algorithms, allowing for linear regression to be replaced with any function family whose probabilistic values can be computed efficiently. This allows us to harness the accuracy of tree-based models like XGBoost, while still producing unbiased estimates. From experiments across eight datasets, we find that our methods give state-of-the-art performance for estimating probabilistic values. For Shapley values, the error of our methods is up to $6\times$ lower than Permutation SHAP (the most popular Monte Carlo method), $2. 75\times$ lower than Kernel SHAP (the most popular linear regression method), and $1. 75\times$ lower than Leverage SHAP (the prior state-of-the-art Shapley value estimator). For more general probabilistic values, we can obtain error up to $60\times$ lower than prior work.

AAAI Conference 2024 Conference Paper

A Simple and Practical Method for Reducing the Disparate Impact of Differential Privacy

• Lucas Rosenblatt
• Julia Stoyanovich
• Christopher Musco

Differentially private (DP) mechanisms have been deployed in a variety of high-impact social settings (perhaps most notably by the U.S. Census). Since all DP mechanisms involve adding noise to results of statistical queries, they are expected to impact our ability to accurately analyze and learn from data, in effect trading off privacy with utility. Alarmingly, the impact of DP on utility can vary significantly among different sub-populations. A simple way to reduce this disparity is with stratification. First compute an independent private estimate for each group in the data set (which may be the intersection of several protected classes), then, to compute estimates of global statistics, appropriately recombine these group estimates. Our main observation is that naive stratification often yields high-accuracy estimates of population-level statistics, without the need for additional privacy budget. We support this observation theoretically and empirically. Our theoretical results center on the private mean estimation problem, while our empirical results center on extensive experiments on private data synthesis to demonstrate the effectiveness of stratification on a variety of private mechanisms. Overall, we argue that this straightforward approach provides a strong baseline against which future work on reducing utility disparities of DP mechanisms should be compared.

NeurIPS Conference 2024 Conference Paper

Benchmarking Estimators for Natural Experiments: A Novel Dataset and a Doubly Robust Algorithm

• R. Teal Witter
• Christopher Musco

Estimating the effect of treatments from natural experiments, where treatments are pre-assigned, is an important and well-studied problem. We introduce a novel natural experiment dataset obtained from an early childhood literacy nonprofit. Surprisingly, applying over 20 established estimators to the dataset produces inconsistent results in evaluating the nonprofits efficacy. To address this, we create a benchmark to evaluate estimator accuracy using synthetic outcomes, whose design was guided by domain experts. The benchmark extensively explores performance as real world conditions like sample size, treatment correlation, and propensity score accuracy vary. Based on our benchmark, we observe that the class of doubly robust treatment effect estimators, which are based on simple and intuitive regression adjustment, generally outperform other more complicated estimators by orders of magnitude. To better support our theoretical understanding of doubly robust estimators, we derive a closed form expression for the variance of any such estimator that uses dataset splitting to obtain an unbiased estimate. This expression motivates the design of a new doubly robust estimator that uses a novel loss function when fitting functions for regression adjustment. We release the dataset and benchmark in a Python package; the package is built in a modular way to facilitate new datasets and estimators. https: //github. com/rtealwitter/naturalexperiments

ICLR Conference 2024 Conference Paper

Improved Active Learning via Dependent Leverage Score Sampling

• Atsushi Shimizu
• Xiaoou Cheng 0001
• Christopher Musco
• Jonathan Weare

We show how to obtain improved active learning methods in the agnostic (adversarial noise) setting by combining marginal leverage score sampling with non-independent sampling strategies that promote spatial coverage. In particular, we propose an easily implemented method based on the \emph{pivotal sampling algorithm}, which we test on problems motivated by learning-based methods for parametric PDEs and uncertainty quantification. In comparison to independent sampling, our method reduces the number of samples needed to reach a given target accuracy by up to $50\%$. We support our findings with two theoretical results. First, we show that any non-independent leverage score sampling method that obeys a weak \emph{one-sided $\ell_{\infty}$ independence condition} (which includes pivotal sampling) can actively learn $d$ dimensional linear functions with $O(d\log d)$ samples, matching independent sampling. This result extends recent work on matrix Chernoff bounds under $\ell_{\infty}$ independence, and may be of interest for analyzing other sampling strategies beyond pivotal sampling. Second, we show that, for the important case of polynomial regression, our pivotal method obtains an improved bound of $O(d)$ samples.

NeurIPS Conference 2024 Conference Paper

Navigable Graphs for High-Dimensional Nearest Neighbor Search: Constructions and Limits

• Haya Diwan
• Jinrui Gou
• Cameron Musco
• Christopher Musco
• Torsten Suel

There has been significant recent interest in graph-based nearest neighbor search methods, many of which are centered on the construction of (approximately) "navigable" graphs over high-dimensional point sets. A graph is navigable if we can successfully move from any starting node to any target node using a greedy routing strategy where we always move to the neighbor that is closest to the destination according to the given distance function. The complete graph is obviously navigable for any point set, but the important question for applications is if sparser graphs can be constructed. While this question is fairly well understood in low-dimensions, we establish some of the first upper and lower bounds for high-dimensional point sets. First, we give a simple and efficient way to construct a navigable graph with average degree $O(\sqrt{n \log n })$ for any set of $n$ points, in any dimension, for any distance function. We compliment this result with a nearly matching lower bound: even under the Euclidean metric in $O(\log n)$ dimensions, a random point set has no navigable graph with average degree $O(n^{\alpha})$ for any $\alpha < 1/2$. Our lower bound relies on sharp anti-concentration bounds for binomial random variables, which we use to show that the {near-neighborhoods} of a set of random points do not overlap significantly, forcing any navigable graph to have many edges.

NeurIPS Conference 2024 Conference Paper

Nearly Optimal Approximation of Matrix Functions by the Lanczos Method

• Noah Amsel
• Tyler Chen
• Anne Greenbaum
• Cameron Musco
• Christopher Musco

Approximating the action of a matrix function $f(\vec{A})$ on a vector $\vec{b}$ is an increasingly important primitive in machine learning, data science, and statistics, with applications such as sampling high dimensional Gaussians, Gaussian process regression and Bayesian inference, principle component analysis, and approximating Hessian spectral densities. Over the past decade, a number of algorithms enjoying strong theoretical guarantees have been proposed for this task. Many of the most successful belong to a family of algorithms called Krylov subspace methods. Remarkably, a classic Krylov subspace method, called the Lanczos method for matrix functions (Lanczos-FA), frequently outperforms newer methods in practice. Our main result is a theoretical justification for this finding: we show that, for a natural class of rational functions, Lanczos-FA matches the error of the best possible Krylov subspace method up to a multiplicative approximation factor. The approximation factor depends on the degree of $f(x)$'s denominator and the condition number of $\vec{A}$, but not on the number of iterations $k$. Our result provides a strong justification for the excellent performance of Lanczos-FA, especially on functions that are well approximated by rationals, such as the matrix square root.

SODA Conference 2024 Conference Paper

On the Unreasonable Effectiveness of Single Vector Krylov Methods for Low-Rank Approximation

• Raphael A. Meyer
• Cameron Musco
• Christopher Musco

ICML Conference 2023 Conference Paper

Dimensionality Reduction for General KDE Mode Finding

• Xinyu Luo
• Christopher Musco
• Cas Widdershoven

Finding the mode of a high dimensional probability distribution $\mathcal{D}$ is a fundamental algorithmic problem in statistics and data analysis. There has been particular interest in efficient methods for solving the problem when $\mathcal{D}$ is represented as a mixture model or kernel density estimate, although few algorithmic results with worst-case approximation and runtime guarantees are known. In this work, we significantly generalize a result of (LeeLiMusco: 2021) on mode approximation for Gaussian mixture models. We develop randomized dimensionality reduction methods for mixtures involving a broader class of kernels, including the popular logistic, sigmoid, and generalized Gaussian kernels. As in Lee et al. ’s work, our dimensionality reduction results yield quasi-polynomial algorithms for mode finding with multiplicative accuracy $(1-\epsilon)$ for any $\epsilon > 0$. Moreover, when combined with gradient descent, they yield efficient practical heuristics for the problem. In addition to our positive results, we prove a hardness result for box kernels, showing that there is no polynomial time algorithm for finding the mode of a kernel density estimate, unless $\mathit{P} = \mathit{NP}$. Obtaining similar hardness results for kernels used in practice (like Gaussian or logistic kernels) is an interesting future direction.

SODA Conference 2023 Conference Paper

Near-Linear Sample Complexity for L p Polynomial Regression

• Raphael A. Meyer
• Cameron Musco
• Christopher Musco
• David P. Woodruff
• Samson Zhou

We study L p polynomial regression. Given query access to a function f: [−1, 1]→ℝ, the goal is to find a degree d polynomial q̂ such that, for a given parameter ε > 0 Here || · || p is the L p norm, ‖ g ‖ p = (∫ 1 −1 | g(t) | p dt ) 1/ p. We show that querying f at points randomly drawn from the Chebyshev measure on [-1, 1] is a near-optimal strategy for polynomial regression in all L p norms. In particular, to find q̂, it suffices to sample points from [-1, 1] with probabilities proportional to this measure. While the optimal sample complexity for polynomial regression was well understood for L 2 and L ∞, our result is the first that achieves sample complexity linear in d and error (1 + ε) for other values of p without any assumptions. Our result requires two main technical contributions. The first concerns p ≤ 2, for which we provide explicit bounds on the L p Lewis weight function of the infinite linear operator underlying polynomial regression. Using tools from the orthogonal polynomial literature, we show that this function is bounded by the Chebyshev density. Our second key contribution is to take advantage of the structure of polynomials to reduce the p > 2 case to the p ≤ 2 case. By doing so, we obtain a better sample complexity than what is possible for general p -norm linear regression problems, for which Ω (d p/2 ) samples are required.

NeurIPS Conference 2023 Conference Paper

Structured Semidefinite Programming for Recovering Structured Preconditioners

• Arun Jambulapati
• Jerry Li
• Christopher Musco
• Kirankumar Shiragur
• Aaron Sidford
• Kevin Tian

We develop a general framework for finding approximately-optimal preconditioners for solving linear systems. Leveraging this framework we obtain improved runtimes for fundamental preconditioning and linear system solving problems including: Diagonal preconditioning. We give an algorithm which, given positive definite $\mathbf{K} \in \mathbb{R}^{d \times d}$ with $\mathrm{nnz}(\mathbf{K})$ nonzero entries, computes an $\epsilon$-optimal diagonal preconditioner in time $\widetilde{O}(\mathrm{nnz}(\mathbf{K}) \cdot \mathrm{poly}(\kappa^\star, \epsilon^{-1}))$, where $\kappa^\star$ is the optimal condition number of the rescaled matrix. Structured linear systems. We give an algorithm which, given $\mathbf{M} \in \mathbb{R}^{d \times d}$ that is either the pseudoinverse of a graph Laplacian matrix or a constant spectral approximation of one, solves linear systems in $\mathbf{M}$ in $\widetilde{O}(d^2)$ time. Our diagonal preconditioning results improve state-of-the-art runtimes of $\Omega(d^{3. 5})$ attained by general-purpose semidefinite programming, and our solvers improve state-of-the-art runtimes of $\Omega(d^{\omega})$ where $\omega > 2. 3$ is the current matrix multiplication constant. We attain our results via new algorithms for a class of semidefinite programs (SDPs) we call matrix-dictionary approximation SDPs, which we leverage to solve an associated problem we call matrix-dictionary recovery.

FOCS Conference 2022 Conference Paper

Active Linear Regression for ℓp Norms and Beyond

• Cameron Musco
• Christopher Musco
• David P. Woodruff
• Taisuke Yasuda 0002

We study active sampling algorithms for linear regression, which aim to query only a small number of entries of a target vector and output a near minimizer to the objective function. For $\ell_{p}$ norm regression for any $ 0\lt p \lt \infty$, we give an algorithm based on Lewis weight sampling which outputs $\mathrm{a}(1+\epsilon)$-approximate solution using just $\tilde{O}(d/\epsilon^{2})$ queries to b for $p\in(0, 1)$, $\tilde{O}$ $(d/\epsilon)$ queries for $p\in(1, 2)$, and $\tilde{O}$ $(d^{p/2}/\epsilon^{p})$ queries for $p\in(2, \ \infty)$. For $p\in(0, 2)$, our bounds are optimal up to logarithmic factors, thus settling the query complexity for this range of p. For $p\in(2, \ \infty)$, our dependence on d is optimal, while our dependence on $\epsilon$ is off by at most a single $\epsilon$ factor, up to logarithmic factors. Our result resolves an open question of Chen and Dereziński, who gave near optimal bounds for the $\ell_{1}$ norm, but required at least $d^{2}/\epsilon^{2}$ samples for $\ell_{p}$ regression with $p\in(1, 2)$, and gave no bounds for $p\in(2, \ \infty)$ or $p\in(0, 1)$. We also provide the first total sensitivity upper bound for loss functions with at most degree p polynomial growth. This improves a recent result of Tukan, Maalouf, and Feldman. By combining this with our techniques for $\ell_{p}$ regression, we obtain the first active regression algorithms for such loss functions, including the important cases of the Tukey and Huber losses. This answers another question of Chen and Dereziński. Our sensitivity bounds also give improvements to a variety of previous results using sensitivity sampling, including Orlicz norm subspace embeddings, robust subspace approximation, and dimension reduction for smoothed p-norms. Finally, our active sampling results give the first sublinear time algorithms for Kronecker product regression under every $\ell_{p}$ norm. Previous results required reading the entire b vector in the kernel feature space. 1 1 Extended abstract; full version available at https: //arxiv. org/abs/2111. 04888.

ICLR Conference 2022 Conference Paper

Fast Regression for Structured Inputs

• Raphael A. Meyer
• Cameron Musco
• Christopher Musco
• David P. Woodruff
• Samson Zhou

We study the $\ell_p$ regression problem, which requires finding $\mathbf{x}\in\mathbb R^{d}$ that minimizes $\|\mathbf{A}\mathbf{x}-\mathbf{b}\|_p$ for a matrix $\mathbf{A}\in\mathbb R^{n \times d}$ and response vector $\mathbf{b}\in\mathbb R^{n}$. There has been recent interest in developing subsampling methods for this problem that can outperform standard techniques when $n$ is very large. However, all known subsampling approaches have run time that depends exponentially on $p$, typically, $d^{\mathcal{O}(p)}$, which can be prohibitively expensive. We improve on this work by showing that for a large class of common \emph{structured matrices}, such as combinations of low-rank matrices, sparse matrices, and Vandermonde matrices, there are subsampling based methods for $\ell_p$ regression that depend polynomially on $p$. For example, we give an algorithm for $\ell_p$ regression on Vandermonde matrices that runs in time $\mathcal{O}(n\log^3 n+(dp^2)^{0.5+\omega}\cdot\text{polylog}\,n)$, where $\omega$ is the exponent of matrix multiplication. The polynomial dependence on $p$ crucially allows our algorithms to extend naturally to efficient algorithms for $\ell_\infty$ regression, via approximation of $\ell_\infty$ by $\ell_{\mathcal{O}(\log n)}$. Of practical interest, we also develop a new subsampling algorithm for $\ell_p$ regression for arbitrary matrices, which is simpler than previous approaches for $p \ge 4$.

STOC Conference 2022 Conference Paper

Sublinear time spectral density estimation

• Vladimir Braverman
• Aditya Krishnan 0001
• Christopher Musco

We present a new sublinear time algorithm for approximating the spectral density (eigenvalue distribution) of an n × n normalized graph adjacency or Laplacian matrix. The algorithm recovers the spectrum up to є accuracy in the Wasserstein-1 distance in O ( n · (1/є)) time given sample access to the graph. This result compliments recent work by David Cohen-Steiner, Weihao Kong, Christian Sohler, and Gregory Valiant (2018), which obtains a solution with runtime independent of n , but exponential in 1/є. We conjecture that the trade-off between dimension dependence and accuracy is inherent. Our method is simple and works well experimentally. It is based on a Chebyshev polynomial moment matching method that employees randomized estimators for the matrix trace. We prove that, for any Hermitian A , this moment matching method returns an є approximation to the spectral density using just O (1/є) matrix-vector products with A . By leveraging stability properties of the Chebyshev polynomial three-term recurrence, we then prove that the method is amenable to the use of coarse approximate matrix-vector products. Our sublinear time algorithm follows from combining this result with a novel sampling algorithm for approximating matrix-vector products with a normalized graph adjacency matrix. Of independent interest, we show a similar result for the widely used kernel polynomial method (KPM), proving that this practical algorithm nearly matches the theoretical guarantees of our moment matching method. Our analysis uses tools from Jackson’s seminal work on approximation with positive polynomial kernels.

NeurIPS Conference 2021 Conference Paper

Dynamic Trace Estimation

• Prathamesh Dharangutte
• Christopher Musco

We study a dynamic version of the implicit trace estimation problem. Given access to an oracle for computing matrix-vector multiplications with a dynamically changing matrix A, our goal is to maintain an accurate approximation to A's trace using as few multiplications as possible. We present a practical algorithm for solving this problem and prove that, in a natural setting, its complexity is quadratically better than the standard solution of repeatedly applying Hutchinson's stochastic trace estimator. We also provide an improved algorithm assuming additional common assumptions on A's dynamic updates. We support our theory with empirical results, showing significant computational improvements on three applications in machine learning and network science: tracking moments of the Hessian spectral density during neural network optimization, counting triangles and estimating natural connectivity in a dynamically changing graph.

AAAI Conference 2021 Conference Paper

Graph Learning for Inverse Landscape Genetics

• Prathamesh Dharangutte
• Christopher Musco

The problem of inferring unknown graph edges from numerical data at a graph’s nodes appears in many forms across machine learning. We study a version of this problem that arises in the field of landscape genetics, where genetic similarity between organisms living in a heterogeneous landscape is explained by a weighted graph that encodes the ease of dispersal through that landscape. Our main contribution is an efficient algorithm for inverse landscape genetics, which is the task of inferring this graph from measurements of genetic similarity at different locations (graph nodes). Inverse landscape genetics is important in discovering impediments to species dispersal that threaten biodiversity and long-term species survival. In particular, it is widely used to study the effects of climate change and human development. Drawing on influential work that models organism dispersal using graph effective resistances (McRae 2006), we reduce the inverse landscape genetics problem to that of inferring graph edges from noisy measurements of these resistances, which can be obtained from genetic similarity data. Building on the NeurIPS 2018 work of Hoskins et al. (2018) on learning edges in social networks, we develop an efficient firstorder optimization method for solving this problem. Despite its non-convex nature, experiments on synthetic and real genetic data establish that our method provides fast and reliable convergence, significantly outperforming existing heuristics used in the field. By providing researchers with a powerful, general purpose algorithmic tool, we hope our work will have a positive impact on accelerating work on landscape genetics.

SODA Conference 2020 Conference Paper

Fast and Space Efficient Spectral Sparsification in Dynamic Streams

• Michael Kapralov
• Aida Mousavifar
• Cameron Musco
• Christopher Musco
• Navid Nouri
• Aaron Sidford
• Jakab Tardos

In this paper, we resolve the complexity problem of spectral graph sparcification in dynamic streams up to polylogarithmic factors. Using a linear sketch we design a streaming algorithm that uses Õ ( n ) space, and with high probability, recovers a spectral sparsifier from the sketch in Õ ( n ) time. 1 Prior results either achieved near optimal Õ ( n ) space, but Ω( n 2 ) recovery time [Kapralov et al. ‘14], or ran in o ( n 2 ) time, but used polynomially suboptimal space [Ahn et al ‘13]. Our main technical contribution is a novel method for recovering graph edges with high effective resistance from a linear sketch. We show how to do so in nearly linear time by ‘bucketing’ vertices of the input graph into clusters using a coarse approximation to the graph's effective resistance metric. A second main contribution is a new pseudorandom generator (PRG) for linear sketching algorithms. Constructed from a locally computable randomness extractor, our PRG stretches a seed of Õ (n) random bits polynomially in length with just log O (1) n run-time cost per evaluation. This improves on Nisan's commonly used PRG, which in our setting would require Õ ( n ) time per evaluation. Our faster PRG is essential to simultaneously achieving near optimal space and time complexity.

NeurIPS Conference 2020 Conference Paper

Fourier Sparse Leverage Scores and Approximate Kernel Learning

• Tamas Erdelyi
• Cameron Musco
• Christopher Musco

We prove new explicit upper bounds on the leverage scores of Fourier sparse functions under both the Gaussian and Laplace measures. In particular, we study s-sparse functions of the form $f(x) = \sum_{j=1}^s a_j e^{i \lambda_j x}$ for coefficients $a_j \in C$ and frequencies $\lambda_j \in R$. Bounding Fourier sparse leverage scores under various measures is of pure mathematical interest in approximation theory, and our work extends existing results for the uniform measure [Erd17, CP19a]. Practically, our bounds are motivated by two important applications in machine learning: 1. Kernel Approximation. They yield a new random Fourier features algorithm for approximating Gaussian and Cauchy (rational quadratic) kernel matrices. For low-dimensional data, our method uses a near optimal number of features, and its runtime is polynomial in the *statistical dimension* of the approximated kernel matrix. It is the first ``oblivious sketching method'' with this property for any kernel besides the polynomial kernel, resolving an open question of [AKM+17, AKK+20b]. 2. Active Learning. They can be used as non-uniform sampling distributions for robust active learning when data follows a Gaussian or Laplace distribution. Using the framework of [AKM+19], we provide essentially optimal results for bandlimited and multiband interpolation, and Gaussian process regression. These results generalize existing work that only applies to uniformly distributed data.

FOCS Conference 2020 Conference Paper

Near Optimal Linear Algebra in the Online and Sliding Window Models

• Vladimir Braverman
• Petros Drineas
• Cameron Musco
• Christopher Musco
• Jalaj Upadhyay
• David P. Woodruff
• Samson Zhou

We initiate the study of numerical linear algebra in the sliding window model, where only the most recent W updates in a stream form the underlying data set. Although many existing algorithms in the sliding window model use or borrow elements from the smooth histogram framework (Braverman and Ostrovsky, FOCS 2007), we show that many interesting linear-algebraic problems, including spectral and vector induced matrix norms, generalized regression, and low-rank approximation, are not amenable to this approach in the row-arrival model. To overcome this challenge, we first introduce a unified row-sampling based framework that gives randomized algorithms for spectral approximation, low-rank approximation/projection-cost preservation, and ℓ 1 -subspace embeddings in the sliding window model, which often use nearly optimal space and achieve nearly input sparsity runtime. Our algorithms are based on “reverse online” versions of offline sampling distributions such as (ridge) leverage scores, ℓ 1 sensitivities, and Lewis weights to quantify both the importance and the recency of a row; our structural results on these distributions may be of independent interest for future algorithmic design. Although our techniques initially address numerical linear algebra in the sliding window model, our row-sampling framework rather surprisingly implies connections to the well-studied online model; our structural results also give the first sample optimal (up to lower order terms) online algorithm for low-rank approximation/projection-cost preservation. Using this powerful primitive, we give online algorithms for column/row subset selection and principal component analysis that resolves the main open question of Bhaskara et al. (FOCS 2019). We also give the first online algorithm for ℓ 1 -subspace embeddings. We further formalize the connection between the online model and the sliding window model by introducing an additional unified framework for deterministic algorithms using a merge and reduce paradigm and the concept of online coresets, which we define as a weighted subset of rows of the input matrix that can be used to compute a good approximation to some given function on all of its prefixes. Our sampling based algorithms in the row-arrival online model yield online coresets, giving deterministic algorithms for spectral approximation, low-rank approximation/projection-cost preservation, and ℓ 1 - subspace embeddings in the sliding window model that use nearly optimal space.

SODA Conference 2020 Conference Paper

Sample Efficient Toeplitz Covariance Estimation

• Yonina C. Eldar
• Jerry Li 0001
• Cameron Musco
• Christopher Musco

NeurIPS Conference 2020 Conference Paper

The Statistical Cost of Robust Kernel Hyperparameter Turning

• Raphael Meyer
• Christopher Musco

This paper studies the statistical complexity of kernel hyperparameter tuning in the setting of active regression under adversarial noise. We consider the problem of finding the best interpolant from a class of kernels with unknown hyperparameters, assuming only that the noise is square-integrable. We provide finite-sample guarantees for the problem, characterizing how increasing the complexity of the kernel class increases the complexity of learning kernel hyperparameters. For common kernel classes (e. g. squared-exponential kernels with unknown lengthscale), our results show that hyperparameter optimization increases sample complexity by just a logarithmic factor, in comparison to the setting where optimal parameters are known in advance. Our result is based on a subsampling guarantee for linear regression under multiple design matrices which may be of independent interest.

STOC Conference 2019 Conference Paper

A universal sampling method for reconstructing signals with simple Fourier transforms

• Haim Avron
• Michael Kapralov
• Cameron Musco
• Christopher Musco
• Ameya Velingker
• Amir Zandieh

Reconstructing continuous signals based on a small number of discrete samples is a fundamental problem across science and engineering. We are often interested in signals with "simple'' Fourier structure -- e.g., those involving frequencies within a bounded range, a small number of frequencies, or a few blocks of frequencies -- i.e., bandlimited, sparse, and multiband signals, respectively. More broadly, any prior knowledge on a signal's Fourier power spectrum can constrain its complexity. Intuitively, signals with more highly constrained Fourier structure require fewer samples to reconstruct.

NeurIPS Conference 2018 Conference Paper

Inferring Networks From Random Walk-Based Node Similarities

• Jeremy Hoskins
• Cameron Musco
• Christopher Musco
• Babis Tsourakakis

Digital presence in the world of online social media entails significant privacy risks. In this work we consider a privacy threat to a social network in which an attacker has access to a subset of random walk-based node similarities, such as effective resistances (i. e. , commute times) or personalized PageRank scores. Using these similarities, the attacker seeks to infer as much information as possible about the network, including unknown pairwise node similarities and edges. For the effective resistance metric, we show that with just a small subset of measurements, one can learn a large fraction of edges in a social network. We also show that it is possible to learn a graph which accurately matches the underlying network on all other effective resistances. This second observation is interesting from a data mining perspective, since it can be expensive to compute all effective resistances or other random walk-based similarities. As an alternative, our graphs learned from just a subset of effective resistances can be used as surrogates in a range of applications that use effective resistances to probe graph structure, including for graph clustering, node centrality evaluation, and anomaly detection. We obtain our results by formalizing the graph learning objective mathematically, using two optimization problems. One formulation is convex and can be solved provably in polynomial time. The other is not, but we solve it efficiently with projected gradient and coordinate descent. We demonstrate the effectiveness of these methods on a number of social networks obtained from Facebook. We also discuss how our methods can be generalized to other random walk-based similarities, such as personalized PageRank scores. Our code is available at https: //github. com/cnmusco/graph-similarity-learning.

SODA Conference 2018 Conference Paper

Stability of the Lanczos Method for Matrix Function Approximation

• Cameron Musco
• Christopher Musco
• Aaron Sidford

Theoretically elegant and ubiquitous in practice, the Lanczos method can approximate f ( A ) x for any symmetric matrix A ∊ ℝ n × n, vector x ∊ ℝ n, and function f. In exact arithmetic, the method's error after k iterations is bounded by the error of the best degree- k polynomial uniformly approximating the scalar function f ( x ) on the range [λ min ( A ), λ max ( A )]. However, despite decades of work, it has been unclear if this powerful guarantee holds in finite precision. We resolve this problem, proving that when, Lanczos essentially matches the exact arithmetic guarantee if computations use roughly log( nC ║ A ║) bits of precision. Our proof extends work of Druskin and Knizhnerman [11], leveraging the stability of the classic Chebyshev recurrence to bound the stability of any polynomial approximating f ( x ). We also study the special case of f ( A ) = A –1 for positive definite A, where stronger guarantees hold for Lanczos. In exact arithmetic the algorithm performs as well as the best polynomial approximating 1/ x at each of A 's eigenvalues, rather than on the full range [λ min ( A ), λ max ( A )]. In seminal work, Greenbaum gives a natural approach to extending this bound to finite precision: she proves that finite precision Lanczos and the related conjugate gradient method match any polynomial approximating 1/ x in a tiny range around each eigenvalue [17]. For A –1, Greenbaum's bound appears stronger than our result. However, we exhibit matrices with condition number κ where exact arithmetic Lanczos converges in polylog( k ) iterations, but Greenbaum's bound predicts at best Ω( κ 1/5 ) iterations in finite precision. It thus cannot offer more than a polynomial improvement over the O ( κ 1/2 ) bound achievable via our result for general f ( A ). Our analysis bounds the power of stable approximating polynomials and raises the question of if they fully characterize the behavior of finite precision Lanczos in solving linear systems. If they do, convergence in less than poly( K ) iterations cannot be expected, even for matrices with clustered, skewed, or otherwise favorable eigenvalue distributions.

SODA Conference 2017 Conference Paper

Input Sparsity Time Low-rank Approximation via Ridge Leverage Score Sampling

• Michael B. Cohen
• Cameron Musco
• Christopher Musco

We present a new algorithm for finding a near optimal low-rank approximation of a matrix A in O ( n nz (A)) time. Our method is based on a recursive sampling scheme for computing a representative subset of A ‘s columns, which is then used to find a low-rank approximation. This approach differs substantially from prior O ( n nz (A)) time algorithms, which are all based on fast Johnson-Lindenstrauss random projections. Our algorithm matches the guarantees of the random projection methods while offering a number of advantages. In addition to better performance on sparse and structured data, sampling algorithms can be applied in settings where random projections cannot. For example, we give new streaming algorithms for the column subset selection and projection-cost preserving sample problems. Our method has also been used in the fastest algorithms for provably accurate Nyström approximation of kernel matrices [56].

ICML Conference 2017 Conference Paper

Random Fourier Features for Kernel Ridge Regression: Approximation Bounds and Statistical Guarantees

• Haim Avron
• Michael Kapralov
• Cameron Musco
• Christopher Musco
• Ameya Velingker
• Amir Zandieh

Random Fourier features is one of the most popular techniques for scaling up kernel methods, such as kernel ridge regression. However, despite impressive empirical results, the statistical properties of random Fourier features are still not well understood. In this paper we take steps toward filling this gap. Specifically, we approach random Fourier features from a spectral matrix approximation point of view, give tight bounds on the number of Fourier features required to achieve a spectral approximation, and show how spectral matrix approximation bounds imply statistical guarantees for kernel ridge regression.

NeurIPS Conference 2017 Conference Paper

Recursive Sampling for the Nystrom Method

• Cameron Musco
• Christopher Musco

We give the first algorithm for kernel Nystrom approximation that runs in linear time in the number of training points and is provably accurate for all kernel matrices, without dependence on regularity or incoherence conditions. The algorithm projects the kernel onto a set of s landmark points sampled by their ridge leverage scores, requiring just O(ns) kernel evaluations and O(ns^2) additional runtime. While leverage score sampling has long been known to give strong theoretical guarantees for Nystrom approximation, by employing a fast recursive sampling scheme, our algorithm is the first to make the approach scalable. Empirically we show that it finds more accurate kernel approximations in less time than popular techniques such as classic Nystrom approximation and the random Fourier features method.

ICML Conference 2016 Conference Paper

Principal Component Projection Without Principal Component Analysis

• Roy Frostig
• Cameron Musco
• Christopher Musco
• Aaron Sidford

We show how to efficiently project a vector onto the top principal components of a matrix, *without explicitly computing these components*. Specifically, we introduce an iterative algorithm that provably computes the projection using few calls to any black-box routine for ridge regression. By avoiding explicit principal component analysis (PCA), our algorithm is the first with no runtime dependence on the number of top principal components. We show that it can be used to give a fast iterative method for the popular principal component regression problem, giving the first major runtime improvement over the naive method of combining PCA with regression. To achieve our results, we first observe that ridge regression can be used to obtain a "smooth projection" onto the top principal components. We then sharpen this approximation to true projection using a low-degree polynomial approximation to the matrix step function. Step function approximation is a topic of long-term interest in scientific computing. We extend prior theory by constructing polynomials with simple iterative structure and rigorously analyzing their behavior under limited precision.

STOC Conference 2015 Conference Paper

Dimensionality Reduction for k-Means Clustering and Low Rank Approximation

• Michael B. Cohen
• Sam Elder
• Cameron Musco
• Christopher Musco
• Madalina Persu

We show how to approximate a data matrix A with a much smaller sketch ~A that can be used to solve a general class of constrained k-rank approximation problems to within (1+ε) error. Importantly, this class includes k-means clustering and unconstrained low rank approximation (i.e. principal component analysis). By reducing data points to just O(k) dimensions, we generically accelerate any exact, approximate, or heuristic algorithm for these ubiquitous problems. For k-means dimensionality reduction, we provide (1+ε) relative error results for many common sketching techniques, including random row projection, column selection, and approximate SVD. For approximate principal component analysis, we give a simple alternative to known algorithms that has applications in the streaming setting. Additionally, we extend recent work on column-based matrix reconstruction, giving column subsets that not only 'cover' a good subspace for A}, but can be used directly to compute this subspace. Finally, for k-means clustering, we show how to achieve a (9+ε) approximation by Johnson-Lindenstrauss projecting data to just O(log k/ε 2 ) dimensions. This is the first result that leverages the specific structure of k-means to achieve dimension independent of input size and sublinear in k.

NeurIPS Conference 2015 Conference Paper

Randomized Block Krylov Methods for Stronger and Faster Approximate Singular Value Decomposition

• Cameron Musco
• Christopher Musco

Since being analyzed by Rokhlin, Szlam, and Tygert and popularized by Halko, Martinsson, and Tropp, randomized Simultaneous Power Iteration has become the method of choice for approximate singular value decomposition. It is more accurate than simpler sketching algorithms, yet still converges quickly for any matrix, independently of singular value gaps. After ~O(1/epsilon) iterations, it gives a low-rank approximation within (1+epsilon) of optimal for spectral norm error. We give the first provable runtime improvement on Simultaneous Iteration: a randomized block Krylov method, closely related to the classic Block Lanczos algorithm, gives the same guarantees in just ~O(1/sqrt(epsilon)) iterations and performs substantially better experimentally. Our analysis is the first of a Krylov subspace method that does not depend on singular value gaps, which are unreliable in practice. Furthermore, while it is a simple accuracy benchmark, even (1+epsilon) error for spectral norm low-rank approximation does not imply that an algorithm returns high quality principal components, a major issue for data applications. We address this problem for the first time by showing that both Block Krylov Iteration and Simultaneous Iteration give nearly optimal PCA for any matrix. This result further justifies their strength over non-iterative sketching methods.

FOCS Conference 2014 Conference Paper

Single Pass Spectral Sparsification in Dynamic Streams

• Michael Kapralov
• Yin Tat Lee
• Cameron Musco
• Christopher Musco
• Aaron Sidford

We present the first single pass algorithm for computing spectral sparsifiers of graphs in the dynamic semi-streaming model. Given a single pass over a stream containing insertions and deletions of edges to a graph, G, our algorithm maintains a randomized linear sketch of the incidence matrix into dimension O(1/∈ 2n polylog(n)). Using this sketch, the algorithm can output a (1±∈) spectral sparsifier for G with high probability. While O(1/∈ 2n polylog(n)) space algorithms are known for computing cut sparsifiers in dynamic streams [1], [2] and spectral sparsifiers in insertion-only streams [3], prior to our work, the best known single pass algorithm for maintaining spectral sparsifiers in dynamic streams required sketches of dimension Ω(1/∈ 2n5/3 ). To achieve our result, we show that, using a coarse sparsifier of G and a linear sketch of G's incidence matrix, it is possible to sample edges by effective resistance, obtaining a spectral sparsifier of arbitrary precision. Sampling from the sketch requires a novel application of ℓ 2 /ℓ 2 sparse recovery, a natural extension of the ℓ 0 methods used for cut sparsifiers in [1]. Recent work of [2] on row sampling for matrix approximation gives a recursive approach for obtaining the required coarse sparsifiers. Under certain restrictions, our approach also extends to the problem of maintaining a spectral approximation for a general matrix A T A given a stream of updates to rows in A.