Fundamental Bias in Inverting Random Sampling Matrices with Application to Sub-sampled Newton

Fundamental Bias in Inverting Random Sampling Matrices with Application to Sub-sampled Newton

Speakers: Chengmei Niu, Zhenyu Liao, Zenan Ling, Michael Mahoney (presented by Liam Hutchinson, University of Melbourne)

Conference: ICML 2025

YouTube: https://slideslive.com/39044017

Overview

This talk, presented by Liam Hutchinson on behalf of authors Chengmei Niu, Zhenyu Liao, Zenan Ling, and Michael Mahoney, addresses a critical theoretical and practical challenge in randomized numerical linear algebra (RNLA): the inherent bias when inverting matrices formed via random sampling. Specifically, the presentation delves into the fundamental bias in inverting random sampling matrices and demonstrates its significant application to sub-sampled Newton methods in large-scale optimization. The core problem arises because while sketching matrices can provide unbiased estimates of the original Gram matrix, their inverses—crucial for solving normal equations and powering second-order optimization algorithms—often exhibit a significant, unaddressed bias.

The research presents a novel and precise derivation of this inversion bias for general random sampling procedures. This breakthrough enables the development of debiased sampling procedures that yield an epsilon-delta unbiased estimator for the inverse Gram matrix. The practical import of this theoretical advancement is showcased through its integration into sub-sampled Newton methods, where it leads to demonstrably faster and, notably, problem-independent convergence rates—a significant achievement in the field of stochastic optimization. This work is pivotal for researchers and practitioners grappling with the computational burdens of large-scale machine learning, offering a pathway to more robust and efficient second-order optimization techniques.

The significance of this work extends beyond theoretical elegance. In an era dominated by massive datasets, the computational cost of traditional linear algebra operations, particularly those involving full Hessian matrices, becomes prohibitive. RNLA offers a powerful paradigm to mitigate these costs by introducing probabilistic elements, trading some precision for substantial computational savings. However, the statistical guarantees of these randomized methods are paramount. By precisely characterizing and correcting the inversion bias, this research strengthens the theoretical foundations of RNLA, making it a more reliable tool for high-stakes applications in machine learning, scientific computing, and beyond.

Background

▶ Watch: Introduction to Randomized Numerical Linear Algebra (0:00)

The explosion of data in modern machine learning has pushed traditional numerical linear algebra methods to their computational limits. Operations like matrix inversion or solving large linear systems, which scale polynomially with the number of data points or features, quickly become intractable. This challenge has fueled the growth of Randomized Numerical Linear Algebra (RNLA), a field dedicated to developing algorithms that leverage randomness to approximate linear algebra operations more efficiently, often with strong statistical guarantees. The underlying premise is that in many machine learning contexts, high numerical precision is not strictly necessary, allowing for a trade-off where computational cost is reduced by introducing probabilistic elements.

A canonical problem in this domain is solving overdetermined linear systems of the form $Ax=b$, where $A$ is an $n \times d$ matrix (typically $n \gg d$, with $n$ being the number of data points and $d$ the number of features) and $b$ is a vector of length $n$. When $n$ becomes very large, computing the exact solution via the normal equations ($A^T A x = A^T b$) requiring the inversion of the $d \times d$ Gram matrix $A^T A$, becomes computationally expensive. To circumvent this, RNLA often employs sketching matrices, denoted as $S$. A sketching matrix $S \in \mathbb{R}^{m \times n}$ (where $m \ll n$) is applied to the original system, transforming it into a smaller, sketched system: $(SA)x_{tilde} = SB$. Solving this reduced system is significantly cheaper, and the goal is for $x_{tilde}$ to be a good approximation of the true solution $x$.

Two primary methods for constructing sketching matrices are commonly used:

1. Random Sampling: This involves randomly selecting rows from the original matrix $A$. This can be done via uniform sampling (each row has an equal probability of being chosen) or importance sampling (rows are sampled with probabilities proportional to their "importance," often quantified by leverage scores). Random sampling typically results in a very sparse sketching matrix, making it highly computationally efficient.
2. Random Projections: This involves multiplying $A$ by a dense random matrix $S$ whose entries are drawn from a specific distribution. Common choices include i.i.d. Gaussian or sub-Gaussian elements. Other methods falling under this category include sub-sampled randomized Walsh-Hadamard Fourier transforms, which leverage fast transforms for computational efficiency.

For these sketching methods to be useful, they need to come with statistical guarantees. Two crucial properties are typically sought:

• Unbiasedness: The expectation of the sketched Gram matrix, $E[(SA)^T (SA)]$, should be equal to the exact Gram matrix $A^T A$. This ensures that, on average, the sketched system reflects the properties of the original system.
• Epsilon-Delta Subspace Embedding Property: This stronger guarantee states that not only is the sketched Gram matrix unbiased, but it is also "epsilon-close" to the exact Gram matrix with high probability. Formally, for some small $\epsilon > 0$ and small $\delta > 0$, the spectral norm of the difference between the sketched and exact Gram matrices (possibly scaled) is bounded by $\epsilon$ with probability at least $1-\delta$. This property ensures that the geometry of the data is largely preserved in the lower-dimensional sketch.

However, a critical challenge arises when we move from the Gram matrix itself to its inverse. While many sketching methods can provide an unbiased or approximately unbiased estimate of $A^T A$, the inverse sketched Gram matrix, $((SA)^T (SA))^{-1}$, is generally not an unbiased estimator of $(A^T A)^{-1}$. This inversion bias is particularly problematic because the inverse Gram matrix is fundamental to solving the normal equations and is a core component in many second-order optimization methods, such as Newton's method, where it approximates the inverse Hessian. The presence of this bias often leads to slow convergence rates when these randomized methods are incorporated into stochastic optimization procedures.

Existing work in this area has made progress, particularly for projection-based methods. For instance, it is known that if $S$ consists of i.i.d. Gaussian random variables, the expectation of the inverse of the sketched Gram matrix can be precisely characterized using techniques from Random Matrix Theory (RMT), specifically by analyzing the resolvent object. This often leads to a simple constant-factor debiasing that can correct the inverse. Similar results have been found for sub-Gaussian projections and leverage score sparsified embedding procedures.

The significant gap, which this research addresses, lies with random sampling methods. While these methods are highly computationally efficient due to the sparsity of their sketching matrices, prior art established that, in general, no simple constant-factor debiasing can achieve an epsilon-delta unbiased estimator for the inverse Gram matrix when using random sampling. This indicated a need for more sophisticated and refined analyses to achieve the desired "win-win" of both computational efficiency and strong statistical guarantees for random sampling.

Key Findings

▶ Watch: Two Main Sketching Methods: Random Sampling vs. Projections (0:40)

The central contribution of this work is the groundbreaking derivation of a precise analytical expression for the inversion bias inherent in inverse sketched Gram matrices formed by general random sampling procedures. This is a fundamental theoretical breakthrough that resolves a long-standing challenge in randomized numerical linear algebra.

The key findings can be summarized as follows:

• Precise Bias Characterization: The research provides a definitive mathematical expression for the expectation of the inverse sketched Gram matrix, $((SA)^T (SA))^{-1}$, demonstrating that it is an epsilon-delta unbiased estimator of the exact inverse Gram matrix, $(A^T A)^{-1}$, but only after applying a specific diagonal correction factor. This factor, referred to as "equation one" in the talk, represents the exact nature of the bias.
• Matrix-Level Debiasing Procedure: Leveraging the derived bias expression, the authors propose a novel matrix-level debiasing procedure. Unlike previous methods that sought simple scalar corrections, this approach directly utilizes the structure of the bias to correct the sampling matrix, thereby achieving the desired epsilon-delta unbiased estimator for the inverse Gram matrix. This debiasing is more intricate than constant-factor approaches but is crucial for general random sampling.
• Discovery of Constant-Factor Debiasing for Specific Sampling Schemes: While general random sampling requires matrix-level debiasing, the work identifies specific instances where constant-factor debiasing is indeed achievable. This includes:
• Approximate Leverage Sampling: Here, a debiasing factor similar to $m/(m-d)$ (where $m$ is sketch size, $d$ is feature dimension) is applicable, akin to projection methods. However, this comes with an additional error term stemming from the approximation inherent in leverage score estimation.
• Sub-sampled Randomized Walsh-Hadamard Transforms: Similar constant-factor debiasing results are shown to hold for these specific types of sketching matrices.
• Application to Sub-sampled Newton Methods: The debiased inverse Gram matrix is directly applied to enhance sub-sampled Newton optimization procedures. By providing a more accurate and unbiased approximation of the inverse Hessian (or its Gauss-Newton form), the debiased iterations exhibit superior performance.
• Problem-Independent Convergence Rates: A highly significant result is that the debiased sub-sampled Newton procedure achieves better convergence rates than all other existing sub-sampled Newton methods in the literature. Crucially, this is reported as the first instance of a problem-independent convergence rate for sub-sampled Newton methods, implying a robust and universally applicable performance guarantee, regardless of the specific problem structure within a broad class of functions.
• Empirical Validation: Numerical experiments on standard datasets (MNIST, CIFAR-10) confirm the efficacy of the proposed algorithm. The debiased sub-sampled Newton (particularly with approximate leverage sampling, "ALev") demonstrates the best convergence complexity tradeoff compared to first-order and other second-order baselines.

In essence, this research provides both the theoretical underpinnings and practical methods to overcome a critical statistical hurdle in RNLA, particularly for the computationally attractive random sampling schemes, thereby paving the way for more efficient and robust large-scale optimization.

Technical Deep Dive

▶ Watch: Challenge: No Constant Factor Debiasing for Random Sampling (1:30)

The technical core of this work revolves around a precise characterization and correction of the bias in the inverse of sketched Gram matrices formed by random sampling. Consider an optimization problem involving a smooth function $f: \mathbb{R}^d \to \mathbb{R}$ that we wish to minimize over a convex set $C$. Many such problems, particularly in machine learning (e.g., least squares, logistic regression), involve a Hessian matrix that can be decomposed into a Gauss-Newton form, often involving the Gram matrix $A^T A$, where $A$ represents the data matrix.

The normal equations for a linear system $Ax=b$ are given by $A^T A x = A^T b$. The solution involves $(A^T A)^{-1}$. In a large-scale setting, computing $A^T A$ and its inverse is prohibitive. Randomized sketching aims to approximate this by forming $S A$, where $S \in \mathbb{R}^{m \times n}$ is a sketching matrix with $m \ll n$. The sketched Gram matrix is then $G_S = (SA)^T (SA)$. While $E[G_S]$ might be proportional to $A^T A$ (potentially after scaling), the fundamental problem addressed is that $E[G_S^{-1}] \neq (A^T A)^{-1}$. This expectation mismatch is the inversion bias.

The authors define an epsilon-delta unbiased estimator as a random matrix $\tilde{C}$ that is an unbiased estimator of $C$ (the true inverse Gram matrix in this context) if there exists an event $\zeta$ such that, with high probability (i.e., with probability at least $1-\delta$), the conditional expectation $E[\tilde{C} | \zeta]$ is close to $C$ (within $\epsilon$ error). Previous work showed that for projection methods (e.g., i.i.d. Gaussian $S$) and leverage score sparsified embeddings, such unbiasedness (or approximate unbiasedness) could be achieved with a simple constant-factor debiasing, often of the form $m/(m-d)$.

The challenge for random sampling methods is that the sketching matrix $S$ is typically very sparse, constructed by selecting rows of an identity matrix according to a sampling distribution. This sparsity, while computationally beneficial, complicates the statistical analysis of $G_S^{-1}$. The key technical achievement is the derivation of a precise expression for the inversion bias for these general random sampling procedures. This expression reveals that $G_S^{-1}$ is an epsilon-delta unbiased estimator of $(A^T A)^{-1}$ if a specific diagonal correction factor is applied. While the exact form of "equation one" was not provided in the talk, its existence implies a structured bias that can be analytically characterized. This correction is not a simple scalar multiplier but rather a matrix-level debiasing procedure, suggesting that the bias depends on the specific diagonal entries of the underlying matrices or their properties.

The methodology for deriving this bias expression heavily relies on Random Matrix Theory (RMT) techniques, particularly those involving the resolvent. The resolvent of a matrix $M$ is defined as $(M - zI)^{-1}$ for a complex scalar $z$. RMT provides powerful tools and concentration inequalities for analyzing such objects, and by setting $z=0$, one can recover the inverse of $M$. This sophisticated machinery allows for a rigorous characterization of the expectation of the inverse of random matrices, even when they arise from sparse sampling.

Beyond the general case, the research also identifies specific random sampling schemes where a simpler constant-factor debiasing can be achieved. For approximate leverage sampling schemes, where rows are sampled based on approximate leverage scores, a factor similar to $m/(m-d)$ is found. However, this introduces an additional error term into the inversion bias, directly attributable to the approximation error in estimating leverage scores. Similar results are reported for sub-sampled randomized Walsh-Hadamard transforms, which combine random sampling with fast orthogonal transforms.

The most compelling application of this debiasing procedure is to sub-sampled Newton methods. Newton's method is a powerful second-order optimization algorithm that uses the inverse Hessian (or its approximation) to determine search directions. In large-scale settings, computing the full Hessian and its inverse is infeasible. Sub-sampled Newton methods approximate the Hessian (and its inverse) using a randomly sampled subset of data. By integrating the newly derived debiased inverse Gram matrix (which approximates the inverse Hessian in Gauss-Newton settings), the proposed debiased sub-sampled Newton procedure can construct a more accurate and statistically sound approximation of the inverse Hessian. This leads to a significant improvement in the local convergence rate. Crucially, the authors claim this is the first time a problem-independent convergence rate has been achieved for sub-sampled Newton methods, indicating a robust theoretical guarantee that holds broadly across different problem instances.

Experimental Setup & Results

▶ Watch: Main Result: Precise Inversion Bias for Random Sampling (1:40)

To demonstrate the practical efficacy of their debiased sub-sampled Newton algorithm, the authors conducted a series of numerical experiments, comparing its performance against established baselines on standard machine learning datasets.

Datasets: The experiments utilized two widely recognized image classification datasets:

• MNIST: A dataset of handwritten digits, commonly used for benchmarking classification algorithms.
• CIFAR-10: A dataset of 32x32 color images across 10 classes, posing a more complex classification challenge than MNIST.

Baselines: The performance of the proposed debiased sub-sampled Newton procedure was evaluated against several types of optimization algorithms:

• First-order baselines: These typically include methods like Stochastic Gradient Descent (SGD) or Adam, which rely solely on gradient information.
• Newton-less procedures: This category likely refers to other optimization methods that do not explicitly compute or approximate the Hessian, perhaps quasi-Newton methods that build up Hessian approximations over iterations, or other forms of stochastic optimization.
• Other second-order methods: This would encompass existing sub-sampled Newton variants or other approaches that attempt to leverage curvature information.

Algorithm Tested: The primary algorithm under investigation was the debiased sub-sampled Newton ALev procedure. This refers to the debiased sub-sampled Newton method specifically employing approximate leverage sampling for sketching. As highlighted in the technical deep dive, approximate leverage sampling is one of the specific cases where a constant-factor debiasing can be applied, making it computationally attractive.

Metrics: The primary metric for evaluation was the convergence complexity tradeoff. This typically involves plotting the objective function value (or loss) against the computational effort (e.g., number of iterations, number of gradient evaluations, or wall-clock time). A better tradeoff implies achieving a lower objective value with less computational cost.

Headline Results: The experimental findings consistently showed that:

• Second-order methods outperform first-order methods: As expected for many problems, algorithms that leverage curvature information (second-order methods) generally converged faster or to a better solution than first-order or Newton-less baselines. The talk explicitly states that "first-order baselines and Newton-less perform much worse in these cases than second-order methods."
• Debiased sub-sampled Newton ALev achieves superior performance: Among the second-order methods, the introduced debiased sub-sampled Newton ALev procedure "seems to offer the best convergence complexity tradeoff." This indicates that by precisely correcting the inversion bias, the algorithm can achieve a desired level of accuracy with fewer computational steps or in less time compared to other existing second-order stochastic optimization methods.

Specifics Not Provided: Due to the concise nature of the conference talk, specific numerical details such as exact objective function values, iteration counts, speedup factors, exact model sizes, GPU types used, or precise hyperparameter settings were not elaborated upon. Similarly, detailed ablation studies breaking down the impact of different components of the debiasing procedure were not presented in this brief overview. However, the clear statement about achieving the "best convergence complexity tradeoff" on standard datasets provides strong evidence of its practical utility. The claim of achieving the "first problem-independent convergence rates" further underscores the theoretical robustness validated by these experiments.

Practical Implications

▶ Watch: Numerical Experiments Demonstrating Algorithm Efficacy (2:05)

The findings presented in "Fundamental Bias in Inverting Random Sampling Matrices with Application to Sub-sampled Newton" carry significant practical implications across various roles in the machine learning ecosystem. By addressing a foundational statistical issue, this work enhances the reliability and efficiency of large-scale optimization.

For Practitioners (e.g., ML Engineers, Data Scientists):

• More Robust Second-Order Optimization: Practitioners can now leverage sub-sampled Newton methods with greater confidence. The debiasing procedure means these powerful second-order optimizers are less prone to the convergence slowdowns caused by statistical bias, potentially leading to faster training times and better model performance for complex models.
• Improved Efficiency for Large Datasets: When working with massive datasets where full Hessian computations are infeasible, randomized sketching is a necessity. This work makes random sampling—a computationally very efficient sketching method due to its sparsity—a more viable option for second-order methods, offering a better balance between speed and statistical guarantees.
• Stronger Theoretical Guarantees: The achievement of problem-independent convergence rates is a major benefit. It suggests that the debiased sub-sampled Newton method will perform consistently well across a wider range of optimization problems, reducing the need for extensive hyperparameter tuning or method switching based on problem specifics.

For Infrastructure Teams (e.g., MLOps Engineers, Cloud Architects):

• Optimized Resource Utilization: Faster and more stable convergence directly translates to more efficient use of computational resources (GPUs, TPUs, CPU clusters). If models train faster, infrastructure costs can be reduced, or more experiments can be run within the same budget.
• Enabling Advanced Algorithms: By making sub-sampled Newton methods more reliable and predictable, this research could encourage wider adoption of these sophisticated algorithms. Infrastructure teams might need to support frameworks or libraries that implement these debiased methods, potentially requiring specific hardware or software stack optimizations.

For Model Builders (e.g., Researchers, Algorithm Developers):

• Foundation for New Algorithms: This work provides a rigorous statistical foundation for building new randomized algorithms that rely on inverse matrix approximations. Future research could explore applying similar debiasing principles to other areas of RNLA or other stochastic algorithms.
• Enhanced Algorithm Design: The precise characterization of bias offers insights into the specific conditions under which different debiasing strategies (matrix-level vs. constant-factor) are necessary or sufficient. This informs the design of more tailored and effective randomized algorithms.

Tradeoffs and Limitations:

• Computational Cost of Debiasing: While random sampling itself is efficient, the "matrix-level debiasing procedure" can be more computationally expensive than simple constant-factor debiasing. The practical overhead of computing and applying this diagonal correction factor needs to be carefully weighed against the benefits of improved convergence. For very high-dimensional problems, computing this correction might still be a bottleneck if it requires operations scaling poorly with $d$.
• Approximation Errors: For methods like approximate leverage sampling, while constant-factor debiasing is possible, it introduces additional error. Practitioners must understand the interplay between the approximation quality (e.g., of leverage scores) and the overall convergence behavior.
• Generalizability of "Problem-Independent": While a significant achievement, "problem-independent" typically refers to independence from specific problem parameters within a class of problems (e.g., smooth convex functions with Lipschitz Hessian). It does not imply universal applicability to all possible optimization problems. Understanding the precise scope of this independence is crucial.
• Implementation Complexity: Implementing these debiased methods might require more sophisticated numerical routines and a deeper understanding of RNLA compared to standard first-order optimizers. This could pose a barrier to adoption for teams without specialized expertise.

In summary, this research provides a crucial tool for enhancing the performance and reliability of large-scale machine learning optimization. It moves the field closer to a reality where the computational efficiency of randomized methods is consistently matched by robust theoretical and practical guarantees.

Key Takeaways

• Fundamental Bias Identified: There is an inherent statistical bias when inverting matrices formed by random sampling procedures in randomized numerical linear algebra.
• Precise Bias Expression Derived: The research provides a novel and exact analytical expression for this inversion bias for general random sampling methods, a significant theoretical advancement.
• Effective Debiasing Procedures: A matrix-level debiasing procedure is introduced to correct this bias, enabling the inverse sketched Gram matrix to serve as an epsilon-delta unbiased estimator of the exact inverse Gram matrix. For specific cases like approximate leverage sampling, a simpler constant-factor debiasing is also shown to be applicable.
• Improved Sub-sampled Newton Convergence: Applying this debiasing to sub-sampled Newton methods leads to demonstrably better convergence rates compared to existing methods.
• First Problem-Independent Convergence Rates: A major breakthrough is the achievement of the first problem-independent convergence rates for sub-sampled Newton methods, offering robust and generalizable performance guarantees.
• Leverages Random Matrix Theory: The analytical derivation of the bias expression relies on sophisticated techniques from Random Matrix Theory (RMT), particularly involving the resolvent.

About the Speaker(s)

The work presented, "Fundamental Bias in Inverting Random Sampling Matrices with Application to Sub-sampled Newton," is authored by Chengmei Niu, Zhenyu Liao, Zenan Ling, and Michael Mahoney. The presentation at ICML 2025 was delivered by Liam Hutchinson from the University of Melbourne on behalf of the authors.

Michael Mahoney is a well-known and highly respected figure in the field of Randomized Numerical Linear Algebra (RNLA) and large-scale data analysis. His extensive work has focused on developing and analyzing randomized algorithms for fundamental problems in linear algebra, machine learning, and scientific computing, often addressing challenges related to massive datasets. While specific affiliations for Chengmei Niu, Zhenyu Liao, and Zenan Ling were not detailed in the brief talk, their collaboration with Mahoney suggests expertise in related areas of theoretical machine learning, optimization, and numerical analysis. Liam Hutchinson's role as presenter from the University of Melbourne indicates his involvement or expertise in the domain, likely as a researcher contributing to or familiar with the project.