Dynamic Low-rank Approximation of Full-Matrix Preconditioner for Training Generalized Linear Models

TL;DR

AdaGram introduces an efficient full-matrix adaptive optimizer that captures parameter correlations with low-rank approximations, enabling faster convergence in large-scale models while reducing computational costs.

Contribution

The paper presents AdaGram, a novel optimizer that efficiently approximates full-matrix preconditioning using low-rank structures and fast symmetric factorization.

Findings

AdaGram converges faster than diagonal adaptive methods on standard tasks.

It matches the performance of full-matrix methods with low-rank approximations.

Demonstrates scalability for large models with reduced computational overhead.

Abstract

Adaptive gradient methods like Adagrad and its variants are widespread in large-scale optimization. However, their use of diagonal preconditioning matrices limits the ability to capture parameter correlations. Full-matrix adaptive methods, approximating the exact Hessian, can model these correlations and may enable faster convergence. At the same time, their computational and memory costs are often prohibitive for large-scale models. To address this limitation, we propose AdaGram, an optimizer that enables efficient full-matrix adaptive gradient updates. To reduce memory and computational overhead, we utilize fast symmetric factorization for computing the preconditioned update direction at each iteration. Additionally, we maintain the low-rank structure of a preconditioner along the optimization trajectory using matrix integrator methods. Numerical experiments on standard machine learning tasks show that AdaGram converges faster or matches the performance of diagonal adaptive optimizers when using rank five and smaller rank approximations. This demonstrates AdaGram's potential as a scalable solution for adaptive optimization in large models.

Tables

Table 1. Table 1: Statistics of the UCI Datasets Used in This Work

Dataset	Instances	Features	Classes
Australian Credit	690	14	2
Heart Disease	303	13	2
Splice Junctions	3,186	60	3

Equations

$$\theta \leftarrow \theta -\eta G_t^{-\text{nicefrac}12}{g}_t,$$

$$G_t=ϵI+\sum _{\tau =1}^t{g}_{\tau }{g}_{\tau }^{\top }.$$

$$G_{t+1}=G_t+g_{t+1}{g}_{t+1}^{\top },$$

$$G_t+g_{t+1}{g}_{t+1}^{\top }=G_t^{\text{nicefrac}12}\left(I+{\bar{g}}_{t+1}{\bar{g}}_{t+1}^{\top }\right)G_t^{\text{nicefrac}12}=G_t^{\text{nicefrac}12}\left(I+\alpha {\bar{g}}_{t+1}{\bar{g}}_{t+1}^{\top }\right)\left(I+\alpha {\bar{g}}_{t+1}{\bar{g}}_{t+1}^{\top }\right)G_t^{\text{nicefrac}12},$$

$$1+\alpha \Vert \bar{g}{\Vert }^2={\left(1+\Vert \bar{g}{\Vert }^2\right)}^{\text{nicefrac}12}.$$

$$G_t^{\text{nicefrac}12}(I+\alpha {\bar{g}}_{t+1}{\bar{g}}_{t+1}^{\top })\ne (I+\alpha {\bar{g}}_{t+1}{\bar{g}}_{t+1}^{\top })G^{\text{nicefrac}12}\ne (G_t+g_{t+1}{g}_{t+1}^{\top }{)}^{\text{nicefrac}12}.$$

$$L_{t+1}^{-1}{g}_{t+1}=\frac{1}{\sqrt{1+\Vert L_t^{-1}{g}_{t+1}{\Vert }^2}}{L}_t^{-1}{g}_{t+1}.$$

$$L_{t+1}=L_t\left(I+{\alpha }_{t+1}{\bar{g}}_{t+1}{\bar{g}}_{t+1}^{\top }\right),$$

$$L_{t+1}^{-1}=\left(I-{\beta }_{t+1}{\bar{g}}_{t+1}{\bar{g}}_{t+1}^{\top }\right)L_t^{-1},$$

$$P_{t+1}$$

$$Q_{t+1}$$

$$L_{t+1}^{-1}=\left(I-P_{t+1}{Q}_{t+1}^{\top }\right)L_0^{-1}.$$

$$L_{t+1}^{-1}$$

$$=(I-P_t{Q}_t^T-{\beta }_{t+1}{\bar{g}}_{t+1}{\bar{g}}_{t+1}^T+P_t{Q}_t^T{\beta }_{t+1}{\bar{g}}_{t+1}{\bar{g}}_{t+1}^T)L_0^{-1}$$

$$=(I-\underset{P_{t+1}{Q}_{t+1}^{\top }}{\underbrace{(P_t{Q}_t^{\top }+{\beta }_{t+1}{\bar{g}}_{t+1}{\bar{g}}_{t+1}^{\top }\left(I-P_t{Q}_t^{\top }\right))}}).$$

$$A_t=P_t{Q}_t^{\top }$$

$$\Delta A_{t+1}={\beta }_{t+1}{\bar{g}}_{t+1}{\bar{g}}_{t+1}^{\top }\left(I-P_t{Q}_t^{\top }\right).$$

$$P_{t+1}{Q}_{t+1}^{\top }=A_{t+1}=A_t+\Delta A_{t+1}.$$

$$P_{t+1}{Q}_{t+1}^{\top }=A_{t+1}=\mu A_t+(1-\mu )\Delta A_{t+1}.$$

$$\mathbf{H}=\frac{1}{B}{\mathbf{X}}^T\mathbf{PX}.$$

$${\nabla }_{\theta }{J}_B(\theta )=\frac{1}{b}\sum _{i\in B}\left(\sigma ({\theta }^T{x}^{(i)})-y^{(i)}\right)x^{(i)}.$$

$$L_{t+1}^{-1}g=\frac{1}{\sqrt{1+\Vert L_t^{-1}g{\Vert }^2}}\cdot L_t^{-1}g$$

$$G_0=ϵI⟹L_0=\frac{1}{\sqrt{ϵ}}I$$

$$L_1{L}_1^{\top }{g}_1=G_1{g}_1=\left(ϵI+g_1{g}_1^{\top }\right)g_1=\left(ϵ+\Vert g_1{\Vert }^2\right)g_1,$$

$$L_1^{-1}{g}_1=\frac{1}{\sqrt{\left(ϵ+\Vert g_1{\Vert }^2\right)}}{g}_1$$

$$G_{t+1}=G_t+g_{t+1}{g}_{t+1}^{\top }.$$

$$G+g{g}^{\top }=L_t\left(I+\bar{g}{\bar{g}}^{\top }\right)L_t^{\top },$$

$$\left(I+\bar{g}{\bar{g}}^{\top }\right)=P{P}^{\top }.$$

$$P=I+\alpha \bar{g}{\bar{g}}^{\top },$$

$$G+g{g}^{\top }=L_t\left(I+\alpha \bar{g}{\bar{g}}^{\top }\right){\left(I+\alpha \bar{g}{\bar{g}}^{\top }\right)}^{\top }{L}_t^{\top },$$

$$L_{t+1}=L_t\left(I+\alpha \bar{g}{\bar{g}}^{\top }\right).$$

$$L_{t+1}^{-1}=\left(I-\frac{\alpha }{1+\alpha \Vert \bar{g}{\Vert }^2}\bar{g}{\bar{g}}^{\top }\right)L_t^{-1}.$$

Full text

1 Introduction

The landscape of stochastic optimization has been fundamentally transformed by the introduction of adaptive gradient methods, which automatically adjust learning rates based on the geometry of the optimization problem. In classical gradient-based optimization, the selection of a learning rate $\eta$ is critical; an improper choice can lead to slow convergence or divergence of the optimization trajectory. A common desideratum is for the step size to decrease as the optimizer approaches a minimum. To enforce this, various learning rate decay schedulers are employed. These schedulers systematically reduce the learning rate $\eta$ over time. Common strategies include step decay, where the learning rate is reduced by a fixed factor at specific epochs; exponential decay, which applies a continuous reduction of the form $\eta_{t}=\eta_{0}\exp(-kt)$; and cosine annealing, which follows a sinusoidal trajectory from the initial rate to a minimum value. While effective, these schedulers are typically non-adaptive, as their evolution is pre-determined and does not respond to the observed gradient statistics during training.

A parallel line of inquiry for improving upon vanilla SGD is the incorporation of momentum. Rather than relying solely on the current gradient $g_{t}$ at step $t$, momentum methods accumulate an exponentially decaying moving average of past gradients, often termed velocity $v_{t}$. The parameter update is then a function of this velocity, i.e., $\theta_{t+1}=\theta_{t}-\eta v_{t}$. This has the effect of dampening high-frequency oscillations along directions of high curvature while accelerating progress along directions with persistent, low-variance gradients. This interpretation becomes especially evident if momentum-based methods are analyzed from the ODE discretization perspective [1, 2, 3]. The canonical example is SGD with Momentum, with a notable refinement being the Nesterov Accelerated Gradient (NAG) [4], which calculates the gradient after applying a preliminary momentum step, providing a more effective ”lookahead” mechanism.

Adaptive methods, as introduced, directly address the limitation of a single, static learning rate by maintaining per-parameter learning rates scaled by historical gradient information. These methods address a critical limitation of traditional stochastic gradient descent (SGD): the use of a single learning rate for all parameters, regardless of their individual optimization characteristics. The pioneering method, AdaGrad, introduced in [5], accumulates the sum of squared gradients for each parameter and normalizes the update by the square root of this sum. While highly effective for sparse problems, its monotonically decreasing learning rate can prematurely halt training. RMSProp refines this by using an exponentially decaying moving average of squared gradients, preventing the aggressive decay. The most prominent adaptive method, Adam (Adaptive Moment Estimation) [6], synergistically combines the momentum concept (a first-moment estimate) with the per-parameter scaling of RMSProp (a second-moment estimate [7]), establishing itself as a default optimizer for a vast range of deep learning applications.

In this work, we exclude methods enriched with momentum and exponential decay techniques from our scope and develop a method based solely on Adagrad that considers feature correlations. The method, introduced in this work, can be further developed with techniques used in more sophisticated adaptive algorithms (Adam [6], RMSProp [7], etc.)

Our contributions are:

• •

AdaGram, an optimization algorithm that enables efficient full-matrix adaptive gradient updates. We use the fast symmetric factorization technique, introduced in [8], for the computation of the preconditioned update direction at each iteration. Additionally, we link the iterative nature of the preconditioner with matrix integrator theory and decomposition update techniques to avoid costly recomputation of the matrix decomposition at each step.

• •

An empirical validation of AdaGram’s performance and efficiency. We conduct numerical experiments on Generalized Linear Models (GLMs), specifically a logistic regression task using a single-layer perceptron. The results, presented in Section 4, demonstrate the practical advantages of our method and motivate its generalization to more complex, arbitrary model architectures.

2 Related work

Classical stochastic gradient descent (SGD) faces several fundamental limitations that impact its optimization performance. To address these, the full-matrix AdaGrad algorithm, introduced in [5] and represented in Algorithm 3, marked a significant step forward. By capturing the complete second-order information from the gradient history, including correlations and cross-parameter dependencies, this approach enables sophisticated preconditioning that accounts for the geometry of the loss landscape. However, the same work also recognized the practical challenge of this method and proposed a diagonal preconditioning strategy as a computationally feasible alternative. In this variant, each parameter receives an individualized learning rate based only on its own gradient history. While efficient, this simplification comes at the cost of ignoring the rich information contained in parameter interactions.

The computational intractability of full-matrix methods in high-dimensional settings stems from their prohibitive complexity. These approaches necessitate matrix inversion and square root operations, resulting in an $O(n^{3})$ computational cost per iteration. This cubic scaling renders such methods impractical for large-scale applications, where diagonal approximations with linear $O(n)$ complexity remain the standard.

Consequently, much research has focused on making full-matrix methods practical. Initial explorations sought to bypass the direct computation of the inverse square root through scalable approximations. For instance, the study [9] proposed maintaining a low-rank factorization of the preconditioning matrix. While this avoids cubic complexity, it introduces challenges in maintaining an accurate factorization. Similarly, the study [10] explored online matrix estimation to approximate the inverse Hessian directly but found that the high variance of stochastic gradients led to unstable estimates, degrading optimizer performance.

A significant conceptual advance came with the Shampoo algorithm [11], which introduced a novel preconditioning strategy based on Kronecker factorization. Instead of approximating the full covariance matrix, Shampoo leverages the tensor structure of gradients in deep learning by maintaining separate, smaller preconditioners for the input and output dimensions. The core innovation is applying a $-\frac{1}{4}$ power to each factor, such that their combined effect approximates the inverse square root of the full covariance matrix. This dramatically reduces computational overhead while capturing crucial geometric information. More recent variants like SOAP [12] have built upon this foundation, achieving near-linear time complexity and making them competitive with Adam in wall-clock time.

Other promising research directions have emerged concurrently. The paper [13] revisited the full-matrix challenge by combining efficient regularization with low-rank approximations. However, despite strong theoretical convergence guarantees, the algorithm did not outperform established benchmarks in practical wall-clock time. Another approach, proposed in [14], involves removing the square root from the update rule entirely. This modification not only improves empirical performance but also endows the resulting algorithm with scale-invariance, a highly desirable property for stabilizing training.

To address the aforementioned limitations, we propose an efficient algorithm for the iterative computation of the inverse square root of the gradient second-moment matrix. Our approach leverages recent advances in numerical linear algebra, exploiting the low-rank structure inherent in gradient accumulation over finite time windows. The method circumvents the computational bottleneck of direct matrix inversion through specialized techniques for fast matrix function evaluation, thereby making a full-matrix approach feasible.

The foundational principles explored in this work are broadly applicable. These ideas can be integrated with or seen as extensions to a wide range of established first-order adaptive methods — including Adam [6], RMSprop [7], Adadelta [15], and Nadam [16]—offering a rich avenue for future hybrid algorithms.

3 Adagram: adaptive optimizer based on the low-rank integrator

This section presents the mathematical framework and algorithmic components for constructing an adaptive optimizer that enables efficient full-matrix adaptive gradient updates, leveraging the special structure of a preconditioner matrix.

3.1 Preliminaries

Following the seminal work [5], we define the general schedule for adaptive gradient step updates as:

[TABLE]

where $\theta\in\mathbb{R}^{n}$ is a parameter vector that we optimize with the preconditioned gradient descent, $g_{t}$ is a gradient at step $t$, and $\eta$ is a scalar denoting the learning rate. An updatable preconditioner matrix $G_{t}$ has the form

[TABLE]

It can be equivalently rewritten via the following rank-1 update series

[TABLE]

initialized by $G_{0}=\epsilon I$, where $\epsilon$ is a scalar and $I$ is an $n\times n$ identity matrix.

We denote an ordinary $L_{2}$ vector norm as $\|\cdot\|$ and use vec() notation to denote the flattening a matrix (a tensor in general case) into a column vector. By construction, the existence of an inverse of matrix $G_{t}$ is guaranteed by the addition of an $\epsilon$ to the matrix main diagonal on the step [math]. The square root exists as the matrix $G_{t}$ is positive semi-definite by design.

Commonly, only the diagonal part of $G_{t}$ is considered during optimization to ensure fast parameter updates at each step. On the other hand, the full-matrix adaptive optimization can potentially bring faster convergence by incorporating the richer second-order information into the convergence process. However, direct implementation of such an approach incurs a prohibitive $O(n^{3})$ computational cost of inverting the full matrix with additional $O(n^{2})$ memory overhead for storing it.

To address the aforementioned computational challenges, we propose a novel algorithm that utilizes the incremental nature of the matrix $G_{t}$ along with the special dynamical low-rank representation for efficiently updating the preconditioner. The method maintains the theoretical advantages of full-matrix preconditioning while achieving better computational efficiency. It enables additional control over the balance between the second-order accuracy and practical scalability in large-scale optimization. To achieve this, we revisit the problem under the isometrical equivalence condition, which helps lift strict constraints that prevent the derivation of an efficient computational form of the preconditioner.

3.2 Revisiting the square root computation

Consider the task of incremental computation of the square root of the matrix $G_{t}$ at two consecutive time steps. Having computed the matrix square root $G_{t}^{\nicefrac{{1}}{{2}}}$ at time $t$, one can derive the following form for the factors of $G_{t+1}$ at time $t+1$:

[TABLE]

where $\bar{g}_{t+1}=G_{t}^{-1}g_{t+1}$ $\alpha$ is a scalar constant that allows factorizing the term $\left(I+\bar{g}\bar{g}^{\top}\right)$ into the product of two matrices, it can be found from the identity:

[TABLE]

From Eq. 4 it follows that even if we knew the exact symmetric matrix square root of $G_{t}$ at some moment $t$, this knowledge could not be directly reused for obtaining the matrix square root at the next step as it would require matrices $G_{t}^{\nicefrac{{1}}{{2}}}$ and $(I+\bar{g}\bar{g}^{\top})$ to commute, which does not hold in the general case, i.e.:

[TABLE]

However, it also hints that computing the non-symmetric triangular factors $L_{t}$ of the Cholesky decomposition $G=L_{t}L_{t}^{\top}$ instead of symmetric $G^{\nicefrac{{1}}{{2}}}$ offers a viable solution. Although Cholesky factors exhibit different properties compared to symmetric matrix square roots, the substitution still maintains equivalence under the isometry condition, which is directly verified by the expression $\|G^{-\nicefrac{{1}}{{2}}}g\|^{2}=g^{\top}G^{-1}g=g^{\top}(LL^{\top})^{-1}g=\|L^{-1}g\|^{2}$.

Conveniently, as we will demonstrate below, both the Cholesky factor $L_{t}$ and its inverse $L_{t}^{-1}$ can be obtained analytically for any $t$ in the form of a matrix factorization with special structure. We will show how to improve the scalability of this new representation further using a dynamical low-rank approximation scheme.

3.3 Efficient Cholesky Factor Computation

Effective optimization requires computation of the preconditioned gradient, which involves a matrix-vector product. We demonstrate that this preconditioned vector can be computed recursively. Furthermore, through matrix decomposition techniques, computational complexity can be substantially reduced by operating exclusively on low-rank matrices and vectors.

Theorem 1** (Recursive computation of a preconditioned gradient).**

Let $g_{t+1}$ be the gradient on step $t+1$ and $L_{t}$ be the Cholesky factor of a matrix $G$ on step $t$. Then the preconditioned gradient on step $t+1$ can be calculated as follows:

[TABLE]

Proof Sketch.

Following the symmetric factorization approach similarly to Eq. 4, one gets an explicit recursive expression for the Cholesky factor of $G_{t+1}$:

[TABLE]

where $\bar{g}_{t+1}=L_{t}^{-1}g_{t+1}$ and $\alpha_{t+1}$ satisfies Eq. 5.

Applying the Sherman-Woodberry-Morrison identity to Eq. 8 yields:

[TABLE]

where $\beta_{t+1}=\frac{\alpha_{t+1}}{1+\alpha_{t+1}\|\bar{g}_{t+1}\|^{2}}$. The full proof is presented in Appendix A.2. ∎

The scheme provided above significantly simplifies the preconditioning process. However, in practice, keeping the matrix $L_{t}$ in memory requires $O(n^{2})$ space; moreover, multiplication by a vector still results in $2n^{2}+n$, leading to undesirable computational costs. This estimation can be improved by introducing low-rank factors, which we will update “on-the-fly” with the projector splitting algorithm.

Unraveling the recursive definition of Eq. 9 gives rise to the following lemma.

Lemma 1** (Recursive representation of the preconditioner inverse).**

Let $P_{t}$ and $Q_{t}$ be iteratively built matrices of the form:

[TABLE]

Then the preconditioner matrix $L^{-1}$ on step $t+1$ takes the form:

[TABLE]

Proof.

The proof is presented in Appendix A.3. ∎

Note that matrices $P_{t}$ and $Q_{t}$ grow with time, thus requiring an increasing amount of memory to support recursive computations.

3.4 Updating preconditioner with low-rank integrator

Following Lemma 1, we obtain a recursive procedure for computing the inverse of the Cholesky factor using the matrices $P_{t}$ and $Q_{t}$. The primary complication of this approach is the unbounded growth in the dimensionality of these factors over time.

One way to restrict the growth is to perform a low-rank approximation of the $P_{t}Q_{t}^{\top}$ term via SVD. There are schemes that improve the efficiency of this step by utilizing the special low-rank structure of the updates that allow avoiding recomputation of the full SVD at each step, for example, the algorithm given in [17]. Another way is to perform updates using a dynamic low-rank approximation, via matrix integrator approach introduced in [18].

Let us have a closer look on the Lemma 1 base step:

[TABLE]

Let

[TABLE]

and

[TABLE]

Then,

[TABLE]

The resulting form can be interpreted as a matrix integration process. The method allows for modeling matrix dynamics over steps by recalculating the rank-r SVD decomposition $P_{0}Q_{0}^{\top}=U_{0}S_{0}V_{0}^{\top}$ with respect to the update from Eq. 12.

For the implementation, we used the full practical algorithm from section 3.2 of the [18], which is given in Algorithm 1.

A key advantage of this approach lies in its adaptive update mechanism. The magnitude of decomposition modification is directly proportional to the magnitude of the gradient update. This property is particularly beneficial for optimization, as it enables the method to “forget” about the large initial updates that are typically irrelevant to later epochs.

Additionally, a hyperparameter $\mu$ can be incorporated in Eq. 13 to control the relative weighting between current gradient updates and historical information accumulated from previous epochs. The optimal value of $\mu$ should be determined before training based on the characteristics of the loss landscape.

[TABLE]

The resulting AdaGram optimizer is presented in Algorithm 2.

4 Numerical experiments

This section presents our experimental evaluation of the AdaGram algorithm in the following Section. All experiments were conducted on a cluster with Tesla V100-SXM3-32GB, Python 3.9, and PyTorch 2.5.1.

4.1 Training generalized linear models

For the consistent hypothesis confirmation, we consider only generalized linear models (GLM), specifically the logistic regression case. Below, we summarize some classical results for these models.

A framework of generalized linear models introduced in [19] generalizes traditional linear regression models by allowing the response variable to follow distributions from the exponential family. The linear combination of predictors is linked to the mean of the response through a monotonic link function.

In binary logistic regression, the target variable $y$ is modeled by a Bernoulli distribution with probability $p=\sigma(\theta^{\top}x)=\frac{\exp(\theta^{\top}x)}{1+\exp(\theta^{\top}x)}$. The Hessian matrix $\mathbf{H}$ of the cross-entropy loss function for the batch case looks as follows:

[TABLE]

This formula resembles the AdaGrad preconditioner matrix, given in Eq. 3 and motivates such preconditioner choice.

4.2 Experimental Setup

For empirical validation of the proposed algorithms, we adapt several UCI [20] datasets (Heart, Australian, and Splice) from the experimental framework of [14]. Their main characteristics are shown in the Table 1.

The baseline methods are the following: Vanilla SGD, Diagonal AdaGrad [5], Shampoo [11], KATE [14]. Hyperparameter selection for each dataset was performed via grid search over batch size, learning rate, and numerical stability parameter $\epsilon$. For AdaGram variants, additional grid search was conducted over rank $r$ and memory parameter $\mu$. The configuration that minimized the objective function value at the final training epoch was selected for all reported experiments. The source code for the experiments is available in the following GitHub repository https://github.com/tnmtvv/AdaGram.

4.3 Evaluation on synthetic datasets

These experiments aim to establish key algorithmic properties through controlled evaluation. The incorporation of the full gradient covariance matrix is designed to capture inter-parameter correlations, and we seek to validate this capability for our proposed methods empirically. In the case of logistic regression, gradients are linked to the samples as follows

[TABLE]

Thus, we can control the complexity of the loss landscape by changing the feature correlations. Following [21], we introduce the following structures of the synthetic datasets and compare the optimizer’s behavior for them:

• •

Isotropic – a dataset with independent features. We expect simple methods, such as vanilla SGD and Adagrad, to be the fastest to converge in this scenario.

• •

Data with simple correlations – dataset with a three-diagonal correlation matrix, which induces structured diagonal dominant correlations. In this case, we expect to see the compatible or supreme convergence speed of the full-matrix approximating methods over the simple ones.

• •

Data with complex correlation – we consider a borderline case with highly correlated features and expect to see the visible gain in full-matrix approximating algorithms’ performance.

The corresponding correlation matrices are visualized in Figure 1.

The performance of our proposed methods on synthetic data with varying correlation structures is presented in Figure 2. As expected, on isotropic data where feature correlations are absent, our methods are outperformed by simpler optimizers like AdaGrad and SGD. However, as correlations are introduced, the advantages of our approach become apparent. On data with simple correlations, our AdaGram variants are highly competitive. In the complex correlations case, our methods demonstrate superior final performance, converging to a lower objective value than baselines.

4.4 Evaluation on real datasets

To explore AdaGram performance on real-world data, we conduct experiments on three UCI benchmark datasets: Australian, Heart, and Splice. The principal characteristics of these datasets are summarized in Table 1. Our primary objective is to analyze convergence behavior with respect to wall-clock time, while also examining the time-to-optimum metric for each optimizer to achieve its best performance.

Figure 3 demonstrates the evolution of test loss as a function of wall-clock time, establishing that our SVD and projector splitting variants consistently achieve superior or equivalent performance compared to all baseline methods. Notably, grid search optimization reveals that the optimal rank for the AdaGram_PS variant remains consistently low (rank 1 or 2 across all datasets). This low-rank property is theoretically significant, as it demonstrates that the method effectively captures essential gradient covariance structure while maintaining minimal computational complexity. The observed efficiency indicates favorable scalability characteristics for high-dimensional optimization scenarios.

Plots, demonstrated in Figure 4, indicate a clear advantage of our proposed methods. Across the datasets, they consistently reach a higher or equal final accuracy than the baseline optimizers. Crucially, they achieve this performance with no significant time overhead; in the case of the Splice dataset, convergence is notably faster.

5 Conclusion

In this work, we introduced AdaGram, a novel optimization algorithm that implements full-matrix adaptive preconditioning with high efficiency. By leveraging a fast symmetric factorization technique from [8] and establishing a connection to matrix integrator theory, our method circumvents the prohibitive computational cost typically associated with full-preconditioner updates. Specifically, AdaGram avoids the need for costly matrix decomposition and inversion at each iteration, making it a practical approach for capturing gradient correlations.

Our empirical evaluation on Generalized Linear Models, particularly for a logistic regression task using a single-layer perceptron, demonstrates the efficiency and potential of AdaGram. The promising results validate our methodology and strongly motivate future work. Key directions include extending the AdaGram framework to more complex, arbitrary model architectures and investigating its performance on a broader range of machine learning tasks.

Appendix A Appendix

A.1 Algorithms

A.2 Proof of Theorem 1 (Recursive computation of a preconditioned gradient)

[TABLE]

Proof.

Let us consider a base case $t=0$. For $t=0$ we have:

[TABLE]

The identity immediately verifies the equality:

[TABLE]

That implies that:

[TABLE]

Let’s prove the proposition for some time step $t>0$. By the means of Eq. (3) we have

[TABLE]

Consider the following relation (indices are omitted for better comprehension):

[TABLE]

where $\bar{g}=L_{t}^{-1}g$. We can further symmetrically factorize the middle term:

[TABLE]

By [8]:

[TABLE]

where $\alpha$ is a constant satisfying the relation 5 which leads to

[TABLE]

or

[TABLE]

By applying the Sherman-Morrison-Woodbury identity, one gets

[TABLE]

From here,

[TABLE]

where we used (5) to simplify expression within brackets. Expanding $\bar{g}$ in the last equation finally gives:

[TABLE]

which proves the relation (7) for arbitrary timestep $t\geq 0$. ∎

A.3 Proof of Lemma (Recursive representation of the preconditioner inverse))

Proof.

We proceed by induction on $t$. Base case: For $t=0$, we have $L_{0}^{-1}=(I-P_{0}Q_{0}^{T})L_{0}^{-1}$ with $P_{0}=Q_{0}=\emptyset$, which holds trivially.

Inductive step: Assume the formula holds for time $t$. From the Sherman-Morrison formula, we have:

[TABLE]

By the inductive hypothesis, $L_{t}^{-1}=(I-P_{t}Q_{t}^{T})L_{0}^{-1}$. Substituting:

[TABLE]

Define $\tilde{g}_{t+1}=(I-Q_{t}P_{t}^{T})\bar{g}_{t+1}$. Then:

[TABLE]

Using matrix distributive property and the definitions of $P_{t+1}$ and $Q_{t+1}$, we can verify that:

[TABLE]

where the updated matrices are:

[TABLE]

∎