VAMO: Efficient Zeroth-Order Variance Reduction for SGD with Faster Convergence

## Abstract

Abstract Optimizing large-scale nonconvex problems, common in deep learning, demands balancing rapid convergence with computational efficiency. First-order (FO) optimizers, which serve as today’s baselines, provide fast convergence and good generalization but often incur high computation and memory costs due to the large size of modern models. Conversely, zeroth-order (ZO) algorithms reduce this burden using estimated gradients, yet their slow convergence in high-dimensional settings limits practicality. We introduce VAMO (VAriance-reduced Mixed-gradient Optimizer), a stochastic variance-reduced method that extends mini-batch SGD with full-batch ZO gradients under an SVRG-style framework. VAMO’s hybrid design utilizes a two-point ZO estimator to achieve a dimension-agnostic convergence rate of 𝒪 ​ ( 1 / T + 1 / b ) \mathcal{O}(1/T+1/b) , where T T is the number of iterations and b b is the batch-size, surpassing the dimension-dependent slowdown of purely ZO methods and significantly improving over SGD’s 𝒪 ​ ( 1 / T ) \mathcal{O}(1/\sqrt{T}) rate. Additionally, we propose a multi-point variant that mitigates the O ​ ( 1 / b ) O(1/b) error by adjusting the number of estimation points to balance convergence and cost. Importantly, VAMO achieves these gains with smaller dynamic memory requirements than many FO baselines, making it particularly attractive for edge deployment. Experiments including traditional neural network training and LLM finetuning confirm that VAMO not only outperforms established FO and ZO methods, but also does so with a light memory footprint.

## 1 Introduction

First-order (FO) optimization methods, particularly Stochastic Gradient Descent (SGD), have been applied in training a wide range of machine learning models. For large-scale problems, variance-reduced (VR) techniques, such as the Stochastic Variance Reduced Gradient (SVRG) algorithm (Johnson & Zhang, 2013; Allen-Zhu & Yuan, 2016; Reddi et al., 2016), offer significant improvements, achieving faster convergence rates, $O(1/T)$, compared to the rate of SGD $O(1/\sqrt{T})$ (Reddi et al., 2016). However, in recent years, extremely large models—such as Large Language Models (LLMs) with billions of parameters—have become increasingly prevalent in machine learning. When training these models, traditional variance reduction methods like SVRG face a major challenge: they require periodically computing the full gradient over the entire dataset, which requires high cost for such large-scale models (Reddi et al., 2016). For LLMs, this step results in prohibitive computational and memory overhead, severely hindering efficient training of large models.

Zeroth-order (ZO) optimization methods present an appealing alternative in this context, as they completely bypass the need for explicit gradient calculations, estimating gradients using only function value queries (Nesterov & Spokoiny, 2017; Liu et al., 2018b). This gradient-free characteristic drastically reduces per-iteration computational cost and memory footprint, making ZO methods attractive for resource-constrained training of LLMs (Malladi et al., 2023; Gautam et al., 2024). Despite these advantages, ZO methods typically exhibit slower theoretical convergence rates than FO methods and, critically, often suffer from a strong dependence on the parameter dimension $d$ (Duchi et al., 2015; Nesterov & Spokoiny, 2017). Given the vast dimensionality of modern LLMs, this dependence can render pure ZO approaches impractically slow. This creates a clear dilemma for large model training: FO methods offer desirable convergence, but suffer from high gradient costs, while ZO methods are cheaper per step but often too slow and scale poorly with the parameter dimension. This naturally raises the question: can we devise a hybrid strategy that combines the strengths of both FO and ZO methods, thereby overcoming their limitations for efficient training of large models?

In this work, we propose VAMO (VAriance-reduced Mixed-gradient Optimizer), a new adaptive hybrid algorithm specifically designed to navigate this dilemma in large-scale non-convex optimization. Our approach integrates FO and ZO techniques within the SVRG framework, aiming to maintain SVRG’s fast convergence while substantially mitigating its computational burden. A major breakthrough here is the replacement of the prohibitively expensive full FO gradient $\nabla f(\hat{\mathbf{x}})$ at SVRG snapshots with an efficient ZO gradient estimate $\hat{\nabla}f(\hat{\mathbf{x}})$, which significantly reduces the computation. This leads to a convergence rate of $O(1/T+1/b)$, significantly outperforming FO-SGD, matching the rate of FO-SVRG only with an additional error of $O(1/b)$. The adaptability of VAMO is then enhanced through several key innovations proposed in this work. First, we extend VAMO with a multi-point ZO gradient estimator for the snapshot gradient $\hat{\nabla}f(\hat{\mathbf{x}})$. Compared to the standard two-point variant, this extension reduces the additional $O(1/b)$ error term and allows for balancing performance and cost flexibly by varying the number of query directions. Second, we introduce a mixing coefficient $\alpha$ into the update rule. While incorporating the ZO full gradient can effectively reduce variance, it inevitably introduces estimation error; the coefficient $\alpha$ provides fine-grained control over this balance by adaptively weighting the FO stochastic gradient and the ZO correction term. Together, these mechanisms make VAMO highly adaptable across diverse optimization scenarios. Importantly, despite this hybrid design, our gradient estimator preserves the unbiasedness of FO-SVRG, distinguishing VAMO from many biased ZO methods and enabling a rigorous convergence analysis with stronger theoretical guarantees.

## 2 Related Work

#### First-order optimization.

While Stochastic Gradient Descent (SGD) (Robbins & Monro, 1951) remains a foundational algorithm in machine learning, its convergence can be slow in large-scale settings due to high gradient variance (Johnson & Zhang, 2013). Addressing this, variance-reduced (VR) methods (Gower et al., 2020), notably SVRG (Johnson & Zhang, 2013; Allen-Zhu & Yuan, 2016; Reddi et al., 2016) and SAGA (Defazio et al., 2014), represent a significant theoretical advancement. These algorithms reduce variance by leveraging gradients from past iterates or periodically computing full-batch gradients and achieve faster convergence rates (e.g., linear convergence under certain assumptions) compared to SGD (Johnson & Zhang, 2013; Reddi et al., 2016). Despite these theoretical benefits, a primary practical limitation is the substantial computational and memory cost associated with full-batch gradients. This overhead can become prohibitive as model sizes and datasets scale. Another common extension of SGD is the use of adaptive step-sizes, as in Adagrad (Duchi et al., 2011) and ADAM (Kingma & Ba, 2014). However, the convergence properties of these adaptive methods remain debated and can be highly sensitive to hyper-parameter choices (Reddi et al., 2019; Défossez et al., 2020; Zhang et al., 2022). To provide a clearer and more interpretable comparison, we focus on standard baselines such as SGD and SVRG on our theoretical analysis, which better isolate the effects of our proposed modifications.

#### Zeroth-order optimization.

Zeroth-order (ZO) optimization approximates gradients using function evaluations instead of explicit gradient computation, offering reduced computational and memory overhead (Zhang et al., 2024). This advantage makes them attractive for large-scale problems like fine-tuning of large language models (Malladi et al., 2023; Gautam et al., 2024; Zhang et al., 2024; Tang et al., 2024; Ling et al., 2024), and their convergence properties are theoretically studied (Duchi et al., 2015; Jamieson et al., 2012; Nesterov & Spokoiny, 2017; Liu et al., 2018b). However, the convergence rates of ZO methods often degrade with increasing parameter dimension $d$, which makes them slow for large models compared to FO methods (Liu et al., 2018b; a; Wang et al., 2018). Even recent applications like MeZO (Malladi et al., 2023) and MeZO-SVRG (Gautam et al., 2024) are constrained by large parameter dimension. While ZO optimization is crucial for black-box scenarios (Chen et al., 2017; Tu et al., 2019), tasks like fine-tuning often have accessible gradients whose computation is merely expensive. This motivates our hybrid approach, which aims to combine ZO’s efficiency with FO’s faster convergence by strategically incorporating both types of gradient information, thereby avoiding the high computational cost associated with FO methods, while also mitigating the performance degradation that ZO methods often suffer in high-dimensional settings.

#### Hybrid Zeroth-Order and First-Order Algorithms.

Combining the strengths of FO and ZO optimization is a relatively recent and underexplored direction, with limited established theoretical analysis. The goal is to enjoy faster convergence while using less computational resource via ZO techniques (Zhang et al., 2024; Li et al., 2024). Early explorations include schemes like applying ZO to shallower model layers and FO to deeper ones (Zhang et al., 2024), or concurrently computing and combining FO-SGD and ZO-SGD updates at each step, as in Addax (Li et al., 2024). However, these initial hybrid approaches may have limitations; for instance, the theoretical convergence rate of Addax still exhibits dependence on the problem dimension $d$ (Li et al., 2024), hindering its effectiveness for large-scale models. Furthermore, many early hybrid strategies often lack explicit mechanisms to adaptively tune the balance between FO accuracy and ZO query efficiency in response to varying computational resources or specific problem demands. The scarcity of hybrid strategies that offer both theoretical convergence and controlled adaptability underscores the novelty of our work.
To contextualize our contributions, Table [1](https://arxiv.org/html/2505.13954v2#S2.T1) summarizes the convergence rates and computational complexities of our proposed methods, referred to as VAMO and VAMO (multi-point) in the table alongside several FO and ZO algorithms. We provide a detailed explanation of Table [1](https://arxiv.org/html/2505.13954v2#S2.T1) in Appendix [B](https://arxiv.org/html/2505.13954v2#A2).

**Table 1: Summary of convergence rate and computational complexity of our proposals given $T$ total iterations. $n$ represents the total number of samples or component functions, $d$ is the parameter dimension, $b$ denotes the mini-batch size, $S$ is the number of epochs or outer loops (for SVRG-type methods, $T\approx Sm$ where $m$ is the number of inner iterations per epoch), and $q$ signifies the number of query directions used for ZO estimation.**
| Method | Grad. estimator | Stepsize | Convergence rate (worst case as $b<n$) | Computational complexity |
| --- | --- | --- | --- | --- |
| ZO-SVRG | Gradient Estimate | $O(1/d)$ | $O(d/T+1/b)$ | $O(nS+bT)$ |
| FO-SGD | Explicit Gradient | $O(1/\sqrt{T})$ | $O(1/\sqrt{T})$ | $O(bdT)$ |
| FO-SVRG | Explicit Gradient | $O(1)$ | $O(1/T)$ | $O(dnS+bdT)$ |
| VAMO | Mixed Gradient | $O(1)$ | $O(1/T+1/b)$ | $O(nS+bdT)$ |
| VAMO(multi-point) | Mixed Gradient | $O(1)$ | $O(1/T+(1-q/d)^{2}/b)$ | $O(qnS+bdT)$ |

## 3 Preliminaries

We consider the following nonconvex finite-sum optimization problem:

$$ $\min_{\mathbf{x}\in\mathbb{R}^{d}}\quad f(\mathbf{x}):=\frac{1}{n}\sum_{i=1}^{n}f_{i}(\mathbf{x}),$ (1) $$

where $\{f_{i}(\mathbf{x})\}_{i=1}^{n}$ are $n$ individual nonconvex cost functions. Note that  equation [1](https://arxiv.org/html/2505.13954v2#S3.E1) is the generic form of many machine learning problems such as training neural networks, since this is the natural form arising from empirical risk minimization (ERM). Next we introduce assumptions we will make throughout the paper and provide the background of ZO gradient estimate.

### 3.1 Assumptions

Throughout this paper, we make the following standard assumptions on the objective function components $f_{i}(\mathbf{x})$. Let $d$ be the dimension of the optimization variable $\mathbf{x}$.

###### Assumption 1 (L-smooth) .

Each function $f_{i}:\mathbb{R}^{d}\to\mathbb{R}$ is $L$-smooth for $i\in[n]:=\{1,2,\ldots,n\}$. That is, for any $\mathbf{x},\mathbf{y}\in\mathbb{R}^{d}$, there exists a constant $L>0$ such that:

$$ $\|\nabla f_{i}(\mathbf{x})-\nabla f_{i}(\mathbf{y})\|_{2}\leq L\|\mathbf{x}-\mathbf{y}\|_{2}$ $$

This also implies that the full objective function $f(\mathbf{x})=\frac{1}{n}\sum_{i=1}^{n}f_{i}(\mathbf{x})$ is $L$-smooth.

###### Assumption 2 (Bounded Variance) .

The variance of the stochastic gradients is bounded. Specifically, for any $\mathbf{x}\in\mathbb{R}^{d}$, there exists a constant $\sigma^{2}\geq 0$ such that:

$$ $\frac{1}{n}\sum_{i=1}^{n}\|\nabla f_{i}(\mathbf{x})-\nabla f(\mathbf{x})\|_{2}^{2}\leq\sigma^{2}$ $$

Here, $\nabla f_{i}(\mathbf{x})$ is the gradient of a single component function, which can be viewed as a stochastic gradient of $f(\mathbf{x})$ if $i$ is chosen uniformly at random from $[n]$.

These assumptions are standard in the analysis of stochastic optimization algorithms for nonconvex problems (Ghadimi & Lan, 2013; Reddi et al., 2016; Nesterov & Spokoiny, 2017).

###### Assumption 1 (L-smooth) .

###### Assumption 2 (Bounded Variance) .

### 3.2 Convergence Notion

We study the nonconvex optimization problem in equation [1](https://arxiv.org/html/2505.13954v2#S3.E1), which is common in modern machine learning.
For convex problems, convergence is typically measured by the expected suboptimality $\mathbb{E}[f(\mathbf{x}_{T})-f(\mathbf{x}^{*})]$.
In contrast, nonconvex problems often contain multiple local minima and saddle points (Dauphin et al., 2014; Balasubramanian & Ghadimi, 2022), making global optimality intractable. Consequently, convergence is evaluated by the first-order stationarity condition via the expected squared gradient norm $\mathbb{E}[\|\nabla f(\mathbf{x})\|_{2}^{2}]$.
An algorithm is deemed convergent when this metric approaches zero or falls below a tolerance $\epsilon$ (Liu et al., 2020; Mu et al., 2024), which is the criterion used in our theoretical guarantees.

### 3.3 ZO Gradient Estimation

Consider an individual cost function $f_{i}:\mathbb{R}^{d}\to\mathbb{R}$ that satisfies the conditions in Assumption [1](https://arxiv.org/html/2505.13954v2#Thmassumption1). The ZO approach estimates gradients using only function evaluations. A commonly used two-point ZO gradient estimator for $f_{i}(\mathbf{x})$ is defined as (Spall, 1992; Nesterov & Spokoiny, 2017):

$$ $\hat{\nabla}f_{i}(\mathbf{x})=\frac{d}{\mu}\left[f_{i}(\mathbf{x}+\mu\mathbf{u}_{i})-f_{i}(\mathbf{x})\right]\mathbf{u}_{i},\quad\text{for }i\in[n],$ (2) $$

where $d$ is the dimensionality of the parameter vector $\mathbf{x}$, $\mu>0$ is a small smoothing parameter, and $\{\mathbf{u}_{i}\}_{i=1}^{n}$ are i.i.d. random direction vectors drawn uniformly from the unit Euclidean sphere in $\mathbb{R}^{d}$ (i.e., $\mathbf{u}_{i}\sim U(\mathbb{S}^{d-1})$) (Flaxman et al., 2004; Shamir, 2017; Gao et al., 2018).

In general, for $\mu>0$, the ZO gradient estimator $\hat{\nabla}f_{i}(\mathbf{x})$ is a biased approximation of the true gradient $\nabla f_{i}(\mathbf{x})$. The bias tends to decrease as $\mu\to 0$. However, in practical implementations, choosing an excessively small $\mu$ can render the function difference $f_{i}(\mathbf{x}+\mu\mathbf{u}_{i})-f_{i}(\mathbf{x})$ highly sensitive to numerical errors or system noise or numerical precision issues, potentially failing to accurately represent the local change in the function (Lian et al., 2016). A key property of the ZO estimator is that for $\mu>0$, it provides an unbiased estimate of the gradient of a smoothed version of $f_{i}$, often denoted $f_{i,\mu}(\mathbf{x})=\mathbb{E}_{\mathbf{v}}[f_{i}(\mathbf{x}+\mu\mathbf{v})]$ (where $\mathbf{v}$ is a random vector from a unit ball or sphere), i.e., $\mathbb{E}_{\mathbf{u}_{i}}[\hat{\nabla}f_{i}(\mathbf{x})]=\nabla f_{i,\mu}(\mathbf{x})$ (Liu et al., 2020).

To reduce the error of the ZO gradient estimate, a multi-point version can be employed. Instead of using a single random direction $\mathbf{u}_{i}$, $q\geq 1$ i.i.d. random directions $\{\mathbf{u}_{i,j}\}_{j=1}^{q}$ are sampled for each $f_{i}$. Since estimating along each direction requires two function queries, the multi-point ZO gradient estimator involves a total of $2q$ function queries and is defined as (Duchi et al., 2015; Liu et al., 2020):

$$ $\hat{\nabla}f_{i}(\mathbf{x})=\frac{d}{\mu q}\sum_{j=1}^{q}\left[f_{i}(\mathbf{x}+\mu\mathbf{u}_{i,j})-f_{i}(\mathbf{x})\right]\mathbf{u}_{i,j},\quad\text{for }i\in[n].$ (3) $$

We refer to this as the multi-point ZO gradient estimate throughout the paper.

### 3.4 Notations

In this paper, we denote $\nabla f(x),\nabla f_{i}(x)$ as FO gradients of $f(x)$ and $f_{i}(x)$, respectively and $\hat{\nabla}f(x),\hat{\nabla}f_{i}(x)$ as their ZO variants. $\mathbb{E}[\cdot]$ operates as the usual mathematical expectation, and $\mathcal{I}$ is a mini-batch of indices sampled from $[n]:={1,\dots,n}$, with size $b=|\mathcal{I}|$. $\|\cdot\|_{2}$ denotes the Euclidean (i.e., $\ell_{2}$) norm.

## 4 Hybrid FO and ZO Stochastic Variance Reduction (VAMO)

### 4.1 From SVRG and ZO-SVRG to Hybrid SVRG

The principles of FO-SVRG and ZO-SVRG have been extensively explored in optimization literature (Johnson & Zhang, 2013; Reddi et al., 2016; Liu et al., 2018b; Ji et al., 2019). FO-SVRG, in particular, is known to achieve a linear convergence rate $O(1/T)$ for non-convex problems under certain conditions, significantly outperforming the convergence rate of FO-SGD (Reddi et al., 2016). The key step of FO-SVRG involves leveraging a full gradient $\nabla f(\hat{\mathbf{x}})$, computed periodically at the snapshot $\hat{\mathbf{x}}$, to construct a variance-reduced stochastic gradient estimate (Johnson & Zhang, 2013):

$$ $\hat{\mathbf{g}}_{\text{FO-SVRG}}=\nabla f_{\mathcal{I}}(\mathbf{x})-\nabla f_{\mathcal{I}}(\hat{\mathbf{x}})+\nabla f(\hat{\mathbf{x}}),$ (4) $$

where $\nabla f_{\mathcal{I}}(\mathbf{x})=\frac{1}{b}\sum_{i\in\mathcal{I}}\nabla f_{i}(\mathbf{x})$ is the mini-batch stochastic gradient from a subset $\mathcal{I}\subseteq[n]$ of size $b$. A crucial property of $\hat{\mathbf{g}}_{\text{FO-SVRG}}$ is that $\hat{\mathbf{g}}_{\text{FO-SVRG}}$ is an unbiased gradient estimate of $\nabla f(\mathbf{x})$, $\mathbb{E}[\hat{\mathbf{g}}_{\text{FO-SVRG}}]=\nabla f(\mathbf{x})$ (Johnson & Zhang, 2013).

In the ZO setting, ZO-SVRG adapts the SVRG structure by replacing all explicit gradient computations with ZO estimates derived from function evaluations:

$$ $\hat{\mathbf{g}}_{\text{ZO-SVRG}}=\hat{\nabla}f_{\mathcal{I}}(\mathbf{x})-\hat{\nabla}f_{\mathcal{I}}(\hat{\mathbf{x}})+\hat{\nabla}f(\hat{\mathbf{x}}),$ (5) $$

where $\hat{\nabla}f_{\mathcal{I}}(\mathbf{x})=(1/b)\sum_{i\in\mathcal{I}}\hat{\nabla}f_{i}(\mathbf{x})$, $\hat{\nabla}f(\mathbf{x})=\hat{\nabla}f_{[n]}(\mathbf{x})$, and $\hat{\nabla}f_{i}(\mathbf{x})$ is a ZO gradient estimate, as defined in Section [3.3](https://arxiv.org/html/2505.13954v2#S3.SS3). While structurally similar, a key distinction is that $\hat{\mathbf{g}}_{\text{ZO-SVRG}}$ is generally a biased estimate of $\nabla f(\mathbf{x})$ due to the inherent bias of $\hat{\nabla}f_{i}(\mathbf{x})$ relative to $\nabla f_{i}(\mathbf{x})$ (Liu et al., 2020). This bias significantly complicates its convergence analysis compared to FO-SVRG.

VAMO (Algorithm [1](https://arxiv.org/html/2505.13954v2#alg1)) is motivated by the high cost of computing the full gradient $\nabla f(\hat{\mathbf{x}})$ in SVRG for large-scale models, and introduces a hybrid gradient estimator that combines FO and ZO components to reduce this overhead:

$$ $\hat{\mathbf{g}}=\nabla f_{\mathcal{I}}(\mathbf{x})-\alpha\left(\hat{\nabla}f_{\mathcal{I}}(\hat{\mathbf{x}})-\hat{\nabla}f(\hat{\mathbf{x}})\right).$ (6) $$

Here, $\nabla f_{\mathcal{I}}(\mathbf{x})$ is the standard FO mini-batch stochastic gradient at the current iterate $\mathbf{x}$, while $\hat{\nabla}f(\hat{\mathbf{x}})$ is the ZO estimate of the full gradient at the snapshot $\hat{\mathbf{x}}$.

A critical design choice in  equation [6](https://arxiv.org/html/2505.13954v2#S4.E6) is the construction of the variance-reduction term $\alpha\left(\hat{\nabla}f_{\mathcal{I}}(\hat{\mathbf{x}})-\hat{\nabla}f(\hat{\mathbf{x}})\right)$. To preserve the desirable unbiased property of FO-SVRG, we ensure that the expectation of the ZO variance correction term, $\mathbb{E}\left[\hat{\nabla}f_{\mathcal{I}}(\hat{\mathbf{x}})-\hat{\nabla}f(\hat{\mathbf{x}})\right]$, is zero. Using ZO estimates for both terms within the parentheses, $\hat{\nabla}f_{\mathcal{I}}(\hat{\mathbf{x}})$ and $\hat{\nabla}f(\hat{\mathbf{x}})$, is key to this property and also contributes to computational savings at the snapshot. Maintaining this unbiasedness is pivotal, as it allows for a more tractable convergence analysis similar to FO-SVRG, distinguishing our approach from many ZO algorithms that contend with biased estimators. Furthermore, VAMO introduces a novel mixing coefficient $\alpha>0$. This parameter allows for explicit control over the influence of the ZO variance correction term. We will provide a detailed theoretical and empirical analysis of $\alpha$ in Appendix [J.2](https://arxiv.org/html/2505.13954v2#A10.SS2).

*The introduction of this hybrid structure, particularly the ZO estimation at snapshot and the mixing coefficient $\alpha$, means that the convergence analysis of VAMO cannot be trivially inherited from existing FO-SVRG or ZO-SVRG analyses.* It requires a dedicated theoretical investigation to characterize its behavior and prove its convergence guarantees, which constitutes a core part of our contribution in Section [4.2](https://arxiv.org/html/2505.13954v2#S4.SS2). This distinct analytical challenge underscores the theoretical novelty of our work.

Figure: Algorithm 1 VAMO $(T,m,\{\eta_{k}\},b,\bar{\mathbf{x}}_{0},\mu,\alpha)$

### 4.2 Convergence analysis

In this section, we present the convergence analysis for VAMO using the two-point ZO gradient estimate in  equation [2](https://arxiv.org/html/2505.13954v2#S3.E2). Our analysis is based on an upper bound on the expected squared gradient norm $\mathbb{E}\bigl[\|\nabla f(\bar{\mathbf{x}})\|_{2}^{2}\bigr]$, as shown in Theorem [1](https://arxiv.org/html/2505.13954v2#Thmtheorem1). As discussed in Section [3.2](https://arxiv.org/html/2505.13954v2#S3.SS2), for non-convex objectives, a small value of $\mathbb{E}\bigl[\|\nabla f(\bar{\mathbf{x}})\|_{2}^{2}\bigr]$ implies convergence to a stationary point.

###### Theorem 1 .

Under the assumptions in Section [3.1](https://arxiv.org/html/2505.13954v2#S3.SS1), and the two-point ZO gradient estimate is used. The output $\bar{\mathbf{x}}$ of Algorithm [1](https://arxiv.org/html/2505.13954v2#alg1) satisfies:

$$ $\displaystyle\mathbb{E}\bigl[\|\nabla f(\bar{\mathbf{x}})\|_{2}^{2}\bigr]$ $\displaystyle\leq\frac{\mathbb{E}[f(\bar{\mathbf{x}}_{0})-f^{*}]}{T\bar{\gamma}}+\frac{S\chi_{m}}{T\bar{\gamma}},$ (7) $$

where $T=Sm$, $f^{*}=\min_{\mathbf{x}}f(\mathbf{x})$, $\bar{\gamma}=\min_{k\in[m]}\gamma_{k}$, and $\chi_{m}=\sum_{k=0}^{m-1}\chi_{k}$.
$\gamma_{k}$ and $\chi_{k}$ are coefficients which depend on $\{\eta_{k},\mu,b,d,\alpha\}$. The proof is provided in Appendix [E](https://arxiv.org/html/2505.13954v2#A5).

Compared to FO-SVRG, Theorem [1](https://arxiv.org/html/2505.13954v2#Thmtheorem1) has an additional error $(S\chi_{m}/(T\bar{\gamma}))$ due to the use of the ZO gradient estimator. To obtain a clear dependence on these parameters and explore deeper insights into convergence, we simplify equation [7](https://arxiv.org/html/2505.13954v2#S4.E7) to suit the specific parameter settings, as shown below.

###### Corollary 1 .

Suppose parameters are set as

$$ $\mu=\frac{1}{\sqrt{T}},\quad\eta_{k}=\eta=\frac{\rho}{L},\quad\alpha=\frac{1}{d},$ (8) $$

with $\beta_{k}=\beta=L$, where $0<\rho\leq 1$ is a universal constant independent of $b$, $d$, $\alpha$, $L$ and $T$. Then Theorem [1](https://arxiv.org/html/2505.13954v2#Thmtheorem1) implies

$$ $\frac{\mathbb{E}[f(\bar{\mathbf{x}}_{0})-f^{*}]}{T\bar{\gamma}}\leq O\biggl(\frac{1}{T}\biggr),\quad\frac{S\chi_{m}}{T\bar{\gamma}}\leq O\biggl(\frac{1}{bT}+\frac{1}{b}\biggr),$ (9) $$

yielding the convergence rate:

$$ $\mathbb{E}\bigl[\|\nabla f(\bar{\mathbf{x}})\|_{2}^{2}\bigr]\leq O\biggl(\frac{1}{T}+\frac{1}{bT}+\frac{1}{b}\biggr).$ (10) $$

The proof is provided in Appendix [F](https://arxiv.org/html/2505.13954v2#A6).

From Corollary [1](https://arxiv.org/html/2505.13954v2#Thmcorollary1), we can observe that one advantage of VAMO is that, compared to previous ZO algorithms, the value of smoothing parameters $\mu$ is less restrictive. For example, ZO-SVRG required $\mu\leq O(1/\sqrt{dT})$, and ZO-SGD required $\mu\leq O(1/d\sqrt{T})$ (Liu et al., 2018b). Compared to FO-SGD, the algorithm achieves an improved rate of $O(1/T)$ rather than the rate of $O(1/\sqrt{T})$. Compared to ZO algorithms, the convergence rate is independent of the parameter dimension $d$. Compared to FO-SVRG, VAMO suffers an additional error of $O(1/b)$ inherited from $(S\chi_{m}/(T\bar{\gamma}))$ in equation [1](https://arxiv.org/html/2505.13954v2#Thmtheorem1).

###### Theorem 1 .

###### Corollary 1 .

### 4.3 Memory Efficiency of VAMO

Building on the convergence analysis of VAMO in the Section [4.2](https://arxiv.org/html/2505.13954v2#S4.SS2), where we showed that VAMO achieves an $O(1/T)$ convergence rate, we now analyze its memory efficiency in large-scale training.

In large-scale training, the total memory footprint of an algorithm is dominated by three components: model parameters $|\mathbf{x}|$, optimizer states, and dynamic memory for gradient computation (Zhang et al., 2024). Here $|\mathbf{x}|$ denotes the memory required to store the model parameters in full precision, and $|\mathbf{x}_{l}|$ is that of the parameters in layer $l$. The symbol $|\mathbf{a}_{l}|$ denotes the memory of activations or intermediate results of one single sample at layer $l$.

For FO methods, dynamic memory is dominated by intermediate results, scaling roughly as $\sum_{l}b\cdot|\mathbf{a}_{l}|$ in Table [1](https://arxiv.org/html/2505.13954v2#S2.T1) for a mini-batch of size $b$. FO-SVRG (Johnson & Zhang, 2013) further requires computing the full gradient at the snapshot point; while this can be accumulated over smaller mini-batches to reduce peak memory of intermediate results to $\sum_{l}b\cdot|\mathbf{a}_{l}|$, it increases computational cost per epoch due to multiple passes over all $n$ samples. If the full FO gradient is computed at once, the dynamic memory of intermediate results will increase to $\sum_{l}n\cdot|\mathbf{a}_{l}|$, which poses a substantial memory burden for large-scale settings. Therefore, FO-SVRG must trade off compute and memory when computing full FO gradients. In contrast, VAMO leverages full ZO gradient estimates, which do not require backpropagation, and thus do not need to store intermediate results. Even for estimation of full batch gradients at once, dynamic memory only reaches $\max_{l}|\mathbf{x}_{l}|$ , independent of $b$ , the same as ZO-SGD, which is much smaller than $\sum_{l}b\cdot|\mathbf{a}_{l}|$. Therefore, VAMO do not need to introduce a compute–memory trade-off like FO-SVRG in large-scale settings. Regarding optimizer states, VAMO only needs $|\mathbf{x}|$ for storing full ZO gradient at the snapshot, substantially smaller than Adagrad and Adam which need around $2|\mathbf{x}|$ or $3|\mathbf{x}|$ to store first- or second-momentum. In training large models, the optimizer states of Adam (Kingma & Ba, 2014) and Adagrad (Duchi et al., 2011) can occupy a substantial portion of memory, posing a significant challenge in resource-constrained environments. Table [3](https://arxiv.org/html/2505.13954v2#A1.T3) summarizes the three memory components for different optimizers and a more detailed theoretical and empirical analysis of memory efficiency is provided in Appendix [A](https://arxiv.org/html/2505.13954v2#A1).

## 5 VAMO with Multi-Point ZO Gradient Estimation

Building upon the VAMO algorithm previously introduced with a two-point ZO gradient estimator, this section presents its multi-point ZO estimation variant. This extension is a key component of VAMO’s adaptive design, as adjusting the number of random directions $q$ in ZO estimate allows for explicit tuning of the trade-off between computational cost and the precision of the ZO-based variance reduction, thereby directly influencing convergence performance.

###### Theorem 2 .

Suppose assumptions A1 and A2 hold, and the multi-point ZO gradient estimate is used in Algorithm [1](https://arxiv.org/html/2505.13954v2#alg1).
The gradient norm bound in equation [7](https://arxiv.org/html/2505.13954v2#S4.E7) yields the simplified convergence rate:

$$ $\mathbb{E}\bigl[\|\nabla f(\bar{\mathbf{x}})\|_{2}^{2}\bigr]\leq O\biggl(\frac{1}{T}+\frac{1}{bT}+\frac{1}{b}\Bigl(1-\frac{q}{d}\Bigr)^{2}\biggr).$ (11) $$

With parameter choices $\mu=\frac{1}{q\sqrt{T}}$, $\eta=\frac{\rho}{L}$, $\alpha=\frac{q}{d}$. The proof is provided in Appendix [G](https://arxiv.org/html/2505.13954v2#A7).

By contrast with Corollary [1](https://arxiv.org/html/2505.13954v2#Thmcorollary1), it can be seen from  equation [11](https://arxiv.org/html/2505.13954v2#S5.E11) that the use of multi-point version of VAMO reduces the error $O(1/b)$ in  equation [1](https://arxiv.org/html/2505.13954v2#Thmcorollary1) by leveraging multiple $q$ direction sampling, while increasing the computational cost accordingly. If $q=d$, the algorithm’s convergence rate become comparable to FO-SVRG. A comprehensive summary and comparison of the computational complexities and convergence rates of our proposed VAMO methods against various FO and ZO algorithms can be found in Table [1](https://arxiv.org/html/2505.13954v2#S2.T1) presented in Section [2](https://arxiv.org/html/2505.13954v2#S2). Briefly, the mixing coefficient $\alpha$ governs the trade-off between the FO stochastic gradient and the ZO variance-correction: larger $\alpha$ increases reliance on ZO information while smaller $\alpha$ favors FO updates to mitigate the error of ZO estimation. A detailed theoretical and empirical analysis of $\alpha$ is provided in Appendix [J.2](https://arxiv.org/html/2505.13954v2#A10.SS2).

###### Theorem 2 .

## 6 Applications and Experiments

In this section, we present empirical results to validate the effectiveness and adaptability of the proposed VAMO through three experiments: (i) an adaptability study on a synthetic task varying the number of ZO query directions $q$, (ii) a benchmark on MNIST classification (LeCun et al., 1998), and (iii) large-scale fine-tuning of GPT-2, GPT-2 (Radford et al., 2019) and RoBERTa (Liu et al., 2019) on SST-2 (Socher et al., 2013) and MNLI (Williams et al., 2017). We benchmark against popular FO methods (FO-SGD (Robbins & Monro, 1951), FO-Adagrad (Duchi et al., 2011), FO-Adam (Kingma & Ba, 2014)) and ZO methods (ZO-SGD (Ghadimi & Lan, 2013), ZO-SVRG (Liu et al., 2018b)).
Full setups and hyperparameters are deferred to Appendix [I](https://arxiv.org/html/2505.13954v2#A9).

Figure: (a)
Refer to caption: x1.png

#### Adaptability Experiment.

To highlight VAMO’s adaptive nature, we tested its multi-point ZO estimation strategy on a synthetic non-convex least-squares task. We compared variants using $q\in\{1,3,5\}$ query directions against FO-SGD.
Fig. [1(a)](https://arxiv.org/html/2505.13954v2#S6.F1.sf1) presents the training loss convergence. Consistent with our theoretical analysis in Table [1](https://arxiv.org/html/2505.13954v2#S2.T1), all VAMO variants achieve an $O(1/T)$ convergence rate, outperforming FO-SGD’s $O(1/\sqrt{T})$ rate. The figure clearly illustrates VAMO’s adaptability: increasing $q$ improves convergence performance, effectively mitigating the additional $O(1/b)$ error term associated with the two-point ($q=1$) variant. This aligns with the theoretical prediction that this error term diminishes for larger $q$ (scaling towards $O((1-q/d)^{2}/b$). These results empirically validate our theory and highlight VAMO’s practical ability to adaptively trade computational cost for enhanced convergence by managing ZO estimation error, a key aspect of its adaptive design.

#### Multiclass Classification.

On the MNIST benchmark, we trained a Multi-Layer Perceptron (MLP) and compared VAMO (two-point version, $q=1$) against FO optimizer (FO-SGD) and ZO optimizers (ZO-SGD and ZO-SVRG).
As shown in Fig. [1(b)](https://arxiv.org/html/2505.13954v2#S6.F1.sf2), VAMO significantly outperforms purely ZO methods, converging faster and reaching a lower final loss.
Its performance is also highly competitive with FO-SGD, underscoring the practical effectiveness of our hybrid strategy.

**Table 2: We measure peak GPU memory consumption (in GB) when fine-tuning RoBERTa-Large and GPT-2 Medium models under varying batch sizes (bs) with a fixed context length (cl=128).**
|  | RoBERTa-Large | GPT-2-Medium |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Method | bs = 16 | bs = 32 | bs = 64 | bs = 16 | bs = 32 | bs = 64 |
| FO-SGD | 4.63 | 6.72 | 10.83 | 7.08 | 11.50 | 20.33 |
| FO-Adagrad | 7.21 | 10.46 | 16.54 | 10.25 | 15.74 | 27.57 |
| FO-Adam | 8.62 | 11.59 | 17.66 | 11.59 | 17.05 | 29.07 |
| VAMO | 5.95 | 7.83 | 11.96 | 8.40 | 12.62 | 21.46 |

#### Fine-Tuning Experiments.

To further assess VAMO’s practical utility and its advantages in complex, large-scale settings, we further evaluate VAMO on fine-tuning pre-trained large models, including GPT-2, GPT-2 Medium, and RoBERTa-large, on the SST-2 and MNLI datasets, comparing against standard FO optimizers (FO-SGD, FO-Adagrad, FO-Adam) and ZO optimizers (ZO-SGD, ZO-SVRG).
Table [2](https://arxiv.org/html/2505.13954v2#S6.T2) presents the peak GPU memory usage of different algorithms for batch sizes of 16, 32, and 64, with the context length fixed at 128. Fig. [2](https://arxiv.org/html/2505.13954v2#S6.F2) presents the training loss of different algorithms over both iteration steps and wall-clock time.
VAMO achieves markedly faster and more stable convergence than FO-SGD, while maintaining a similar memory footprint.
This improvement is theoretically supported by Table [1](https://arxiv.org/html/2505.13954v2#S2.T1): VAMO enjoys a convergence rate of $\mathcal{O}(1/T)$, superior to the $\mathcal{O}(1/\sqrt{T})$ rate of FO-SGD. Compared with ZO methods, VAMO not only yields better empirical performance but also avoids dependence on the parameter dimension $d$, demonstrating stronger scalability to large models. When compared to Adagrad, VAMO achieves a similar level of performance but with far lower memory overhead, further validating its practicality.
Although VAMO underperforms Adam in terms of empirical convergence speed, it remains much more memory-efficient, as shown in Table [2](https://arxiv.org/html/2505.13954v2#S6.T2). Overall, these results confirm that VAMO provides a favorable balance between convergence efficiency and memory cost, making it a competitive optimizer for fine-tuning large models.

Figure: Figure 2: Convergence comparison on fine-tuning GPT-2, GPT-2 medium and RoBERTa-large on SST-2 and MNLI datasets against pure FO and ZO methods.
Refer to caption: x3.png

## 7 Conclusion

In this paper, we propose a hybrid FO and ZO variance-reduced algorithm, VAMO, for nonconvex optimization. We demonstrate that compared to FO-SGD, our algorithm improves the convergence rate from $O(1/\sqrt{T})$ to a linear rate of $O(1/T)$, achieving convergence performance similar to FO-SVRG. Compared to ZO algorithms, our method maintains convergence performance independent of the parameter dimension $d$, making it effective for optimizing high-dimensional problems. However, due to the use of two-point ZO gradient estimation, our convergence result includes an additional error term $O(1/b)$. To mitigate this, we introduce a multi-point ZO gradient estimation variant, which reduces this error. Unlike previous purely FO or ZO methods, our hybrid approach leverages the advantages of both, enabling a more flexible trade-off between efficiency and convergence performance. This makes it more adaptable to real-world applications with complex constraints. Our theoretical analysis and empirical evaluations, including comparisons with state-of-the-art methods, demonstrate the effectiveness of our approach.

## References

- Allen-Zhu & Yuan (2016)
Zeyuan Allen-Zhu and Yang Yuan.
Improved svrg for non-strongly-convex or sum-of-non-convex objectives.
In *International conference on machine learning*, pp. 1080–1089. PMLR, 2016.
- Balasubramanian & Ghadimi (2022)
Krishnakumar Balasubramanian and Saeed Ghadimi.
Zeroth-order nonconvex stochastic optimization: Handling constraints, high dimensionality, and saddle points.
*Foundations of Computational Mathematics*, 22(1):35–76, 2022.
- Chen et al. (2017)
Pin-Yu Chen, Huan Zhang, Yash Sharma, Jinfeng Yi, and Cho-Jui Hsieh.
Zoo: Zeroth order optimization based black-box attacks to deep neural networks without training substitute models.
In *Proceedings of the 10th ACM workshop on artificial intelligence and security*, pp. 15–26, 2017.
- Dauphin et al. (2014)
Yann N Dauphin, Razvan Pascanu, Caglar Gulcehre, Kyunghyun Cho, Surya Ganguli, and Yoshua Bengio.
Identifying and attacking the saddle point problem in high-dimensional non-convex optimization.
*Advances in neural information processing systems*, 27, 2014.
- Defazio et al. (2014)
Aaron Defazio, Francis Bach, and Simon Lacoste-Julien.
Saga: A fast incremental gradient method with support for non-strongly convex composite objectives.
*Advances in neural information processing systems*, 27, 2014.
- Défossez et al. (2020)
Alexandre Défossez, Léon Bottou, Francis Bach, and Nicolas Usunier.
A simple convergence proof of adam and adagrad.
*arXiv preprint arXiv:2003.02395*, 2020.
- Duchi et al. (2011)
John Duchi, Elad Hazan, and Yoram Singer.
Adaptive subgradient methods for online learning and stochastic optimization.
*Journal of machine learning research*, 12(7), 2011.
- Duchi et al. (2015)
John C Duchi, Michael I Jordan, Martin J Wainwright, and Andre Wibisono.
Optimal rates for zero-order convex optimization: The power of two function evaluations.
*IEEE Transactions on Information Theory*, 61(5):2788–2806, 2015.
- Flaxman et al. (2004)
Abraham D Flaxman, Adam Tauman Kalai, and H Brendan McMahan.
Online convex optimization in the bandit setting: gradient descent without a gradient.
*arXiv preprint cs/0408007*, 2004.
- Gao et al. (2018)
Xiang Gao, Bo Jiang, and Shuzhong Zhang.
On the information-adaptive variants of the admm: an iteration complexity perspective.
*Journal of Scientific Computing*, 76:327–363, 2018.
- Gautam et al. (2024)
Tanmay Gautam, Youngsuk Park, Hao Zhou, Parameswaran Raman, and Wooseok Ha.
Variance-reduced zeroth-order methods for fine-tuning language models.
*arXiv preprint arXiv:2404.08080*, 2024.
- Ghadimi & Lan (2013)
Saeed Ghadimi and Guanghui Lan.
Stochastic first-and zeroth-order methods for nonconvex stochastic programming.
*SIAM journal on optimization*, 23(4):2341–2368, 2013.
- Gower et al. (2020)
Robert M Gower, Mark Schmidt, Francis Bach, and Peter Richtárik.
Variance-reduced methods for machine learning.
*Proceedings of the IEEE*, 108(11):1968–1983, 2020.
- Jamieson et al. (2012)
Kevin G Jamieson, Robert Nowak, and Ben Recht.
Query complexity of derivative-free optimization.
*Advances in Neural Information Processing Systems*, 25, 2012.
- Ji et al. (2019)
Kaiyi Ji, Zhe Wang, Yi Zhou, and Yingbin Liang.
Improved zeroth-order variance reduced algorithms and analysis for nonconvex optimization.
In *International conference on machine learning*, pp. 3100–3109. PMLR, 2019.
- Johnson & Zhang (2013)
Rie Johnson and Tong Zhang.
Accelerating stochastic gradient descent using predictive variance reduction.
*Advances in neural information processing systems*, 26, 2013.
- Kingma & Ba (2014)
Diederik P Kingma and Jimmy Ba.
Adam: A method for stochastic optimization.
*arXiv preprint arXiv:1412.6980*, 2014.
- LeCun et al. (1998)
Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.
Gradient-based learning applied to document recognition.
*Proceedings of the IEEE*, 86(11):2278–2324, 1998.
- Lei et al. (2017)
Lihua Lei, Cheng Ju, Jianbo Chen, and Michael I Jordan.
Non-convex finite-sum optimization via scsg methods.
*Advances in neural information processing systems*, 30, 2017.
- Li et al. (2024)
Zeman Li, Xinwei Zhang, Peilin Zhong, Yuan Deng, Meisam Razaviyayn, and Vahab Mirrokni.
Addax: Utilizing zeroth-order gradients to improve memory efficiency and performance of sgd for fine-tuning language models.
*arXiv preprint arXiv:2410.06441*, 2024.
- Lian et al. (2016)
Xiangru Lian, Huan Zhang, Cho-Jui Hsieh, Yijun Huang, and Ji Liu.
A comprehensive linear speedup analysis for asynchronous stochastic parallel optimization from zeroth-order to first-order.
*Advances in neural information processing systems*, 29, 2016.
- Ling et al. (2024)
Zhenqing Ling, Daoyuan Chen, Liuyi Yao, Yaliang Li, and Ying Shen.
On the convergence of zeroth-order federated tuning for large language models.
In *Proceedings of the 30th ACM SIGKDD Conference on Knowledge Discovery and Data Mining*, pp. 1827–1838, 2024.
- Liu et al. (2018a)
Liu Liu, Minhao Cheng, Cho-Jui Hsieh, and Dacheng Tao.
Stochastic zeroth-order optimization via variance reduction method.
*arXiv preprint arXiv:1805.11811*, 2018a.
- Liu et al. (2018b)
Sijia Liu, Bhavya Kailkhura, Pin-Yu Chen, Paishun Ting, Shiyu Chang, and Lisa Amini.
Zeroth-order stochastic variance reduction for nonconvex optimization.
*Advances in neural information processing systems*, 31, 2018b.
- Liu et al. (2020)
Sijia Liu, Pin-Yu Chen, Bhavya Kailkhura, Gaoyuan Zhang, Alfred O Hero III, and Pramod K Varshney.
A primer on zeroth-order optimization in signal processing and machine learning: Principals, recent advances, and applications.
*IEEE Signal Processing Magazine*, 37(5):43–54, 2020.
- Liu et al. (2019)
Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov.
Roberta: A robustly optimized bert pretraining approach.
*arXiv preprint arXiv:1907.11692*, 2019.
- Malladi et al. (2023)
Sadhika Malladi, Tianyu Gao, Eshaan Nichani, Alex Damian, Jason D Lee, Danqi Chen, and Sanjeev Arora.
Fine-tuning language models with just forward passes.
*Advances in Neural Information Processing Systems*, 36:53038–53075, 2023.
- Mu et al. (2024)
Huaiyi Mu, Yujie Tang, and Zhongkui Li.
Variance-reduced gradient estimator for nonconvex zeroth-order distributed optimization.
*arXiv preprint arXiv:2409.19567*, 2024.
- Nesterov & Spokoiny (2017)
Yurii Nesterov and Vladimir Spokoiny.
Random gradient-free minimization of convex functions.
*Foundations of Computational Mathematics*, 17(2):527–566, 2017.
- Radford et al. (2019)
Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, Ilya Sutskever, et al.
Language models are unsupervised multitask learners.
*OpenAI blog*, 1(8):9, 2019.
- Reddi et al. (2016)
Sashank J Reddi, Ahmed Hefny, Suvrit Sra, Barnabas Poczos, and Alex Smola.
Stochastic variance reduction for nonconvex optimization.
In *International conference on machine learning*, pp. 314–323. PMLR, 2016.
- Reddi et al. (2019)
Sashank J Reddi, Satyen Kale, and Sanjiv Kumar.
On the convergence of adam and beyond.
*arXiv preprint arXiv:1904.09237*, 2019.
- Robbins & Monro (1951)
Herbert Robbins and Sutton Monro.
A stochastic approximation method.
*The annals of mathematical statistics*, pp. 400–407, 1951.
- Shamir (2017)
Ohad Shamir.
An optimal algorithm for bandit and zero-order convex optimization with two-point feedback.
*Journal of Machine Learning Research*, 18(52):1–11, 2017.
- Socher et al. (2013)
Richard Socher, Alex Perelygin, Jean Wu, Jason Chuang, Christopher D Manning, Andrew Y Ng, and Christopher Potts.
Recursive deep models for semantic compositionality over a sentiment treebank.
In *Proceedings of the 2013 conference on empirical methods in natural language processing*, pp. 1631–1642, 2013.
- Spall (1992)
James C Spall.
Multivariate stochastic approximation using a simultaneous perturbation gradient approximation.
*IEEE transactions on automatic control*, 37(3):332–341, 1992.
- Tang et al. (2024)
Xinyu Tang, Ashwinee Panda, Milad Nasr, Saeed Mahloujifar, and Prateek Mittal.
Private fine-tuning of large language models with zeroth-order optimization.
*arXiv preprint arXiv:2401.04343*, 2024.
- Tu et al. (2019)
Chun-Chen Tu, Paishun Ting, Pin-Yu Chen, Sijia Liu, Huan Zhang, Jinfeng Yi, Cho-Jui Hsieh, and Shin-Ming Cheng.
Autozoom: Autoencoder-based zeroth order optimization method for attacking black-box neural networks.
In *Proceedings of the AAAI conference on artificial intelligence*, volume 33, pp. 742–749, 2019.
- Wang et al. (2018)
Yining Wang, Simon Du, Sivaraman Balakrishnan, and Aarti Singh.
Stochastic zeroth-order optimization in high dimensions.
In *International conference on artificial intelligence and statistics*, pp. 1356–1365. PMLR, 2018.
- Williams et al. (2017)
Adina Williams, Nikita Nangia, and Samuel R Bowman.
A broad-coverage challenge corpus for sentence understanding through inference.
*arXiv preprint arXiv:1704.05426*, 2017.
- Zhang et al. (2024)
Yihua Zhang, Pingzhi Li, Junyuan Hong, Jiaxiang Li, Yimeng Zhang, Wenqing Zheng, Pin-Yu Chen, Jason D Lee, Wotao Yin, Mingyi Hong, et al.
Revisiting zeroth-order optimization for memory-efficient llm fine-tuning: A benchmark.
*arXiv preprint arXiv:2402.11592*, 2024.
- Zhang et al. (2022)
Yushun Zhang, Congliang Chen, Naichen Shi, Ruoyu Sun, and Zhi-Quan Luo.
Adam can converge without any modification on update rules.
*Advances in neural information processing systems*, 35:28386–28399, 2022.

## Appendix A Detailed Discussion on Memory Efficiency

We have made a brief analysis of VAMO’s memory efficiency in Section [4.3](https://arxiv.org/html/2505.13954v2#S4.SS3). In this section, we analyse VAMO’s memory profile in detail and compare it against standard FO and ZO optimizers, theoretically and empirically.

### A.1 Theoretical Analysis

We adapt the memory analysis framework of Zhang et al. (2024) to compare the peak memory consumption of different optimizers. Table [3](https://arxiv.org/html/2505.13954v2#A1.T3) summarizes the decomposition into three components: the model (Weight Mem.), the optimizer states (Opt. State Mem.), and dynamic allocation for computing gradients and optimization (Dynamic Mem.).

**Table 3: Comparison of the instant peak memory consumption of different optimizers when fine-tuning the full model. Here $|\mathbf{x}|$ denotes the memory required to store the model parameters in full precision, and $|\mathbf{x}_{l}|$ is that of the parameters in layer $l$. The symbol $|\mathbf{a}_{l}|$ denotes the memory of intermediate activations of one single sample at layer $l$, where $b$ and $n$ correspond to the mini-batch size and the total dataset size, respectively. In Dynamic Mem., $b\cdot|\mathbf{a}_{l}|$ represents the memory of saving mini-batch activations at layer $l$, while $|\mathbf{x}_{l}|$ also accounts for temporarily saved gradients, as they have the same size as the corresponding parameters (Zhang et al., 2024).**
| Optimizer | Weight Mem. | Dynamic Mem. (Grad.&Opt.) | Opt. State Mem. |
| --- | --- | --- | --- |
| FO-SGD | $|\mathbf{x}|$ | $\sum_{l}\max\{b\cdot|\mathbf{a}_{l}|,|\mathbf{x}_{l}|\}$ | 0 |
| FO-SVRG (accumulation) | $|\mathbf{x}|$ | $\sum_{l}\max\{b\cdot|\mathbf{a}_{l}|,|\mathbf{x}_{l}|\}$ | $2|\mathbf{x}|$ |
| FO-SVRG (full batch at once) | $|\mathbf{x}|$ | $\sum_{l}\max\{n\cdot|\mathbf{a}_{l}|,|\mathbf{x}_{l}|\}$ | $2|\mathbf{x}|$ |
| FO-Adam | $|\mathbf{x}|$ | $\sum_{l}\max\{b\cdot|\mathbf{a}_{l}|,|\mathbf{x}_{l}|\}$ | $3|\mathbf{x}|$ |
| FO-Adagrad | $|\mathbf{x}|$ | $\sum_{l}\max\{b\cdot|\mathbf{a}_{l}|,|\mathbf{x}_{l}|\}$ | $2|\mathbf{x}|$ |
| VAMO | $|\mathbf{x}|$ | $\sum_{l}\max\{b\cdot|\mathbf{a}_{l}|,|\mathbf{x}_{l}|\}$ | $|\mathbf{x}|$ |
| ZO-SGD | $|\mathbf{x}|$ | $\max_{l}|\mathbf{x}_{l}|$ | 0 |

The advantage of ZO methods is the removal of backpropagation, which eliminates the need to store intermediate activations ${|\mathbf{a}_{l}|}$ for each input example during the forward pass. In FO optimizers, the total activation memory grows as $b\cdot\sum_{l}|\mathbf{a}_{l}|$ for a mini-batch of size $b$, while ZO methods completely avoid this dependence on the batch size. This reduction in dynamic memory cost becomes especially significant when $b$ is large in large-scale settings. Moreover, ZO gradients can be estimated in a layer-wise manner, so only the gradients of the currently processed layer need to be stored, further lowering the peak memory footprint compared to FO methods (Zhang et al., 2024). For ZO-SGD, dynamic memory is dominated by a temporary copy of the parameters of a single layer for gradient estimation, resulting in a peak memory of $\max_{l}|\mathbf{x}_{l}|$, independent of batch size, and no extra optimizer state is required.

In contrast, all FO methods require storing intermediate activations or other temporary variables to compute gradients,
which results in dynamic memory that grows with the mini-batch size.
FO (Johnson & Zhang, 2013) additionally needs to compute the full gradient at the snapshot point.
Although this full gradient can be obtained by accumulating gradients over smaller mini-batches,
reducing dynamic memory of intermediate results per step to $b\cdot|\mathbf{a}_{l}|$, it comes at the expense of higher computational cost per epoch due to multiple passes over all mini-batches.
Computing the full gradient in a single full batch would substantially increase dynamic memory of intermediate results to $n\cdot|\mathbf{a}_{l}|$, requiring a compute-memory trade-off as commonly done in FO-SVRG.
In contrast, VAMO leverages ZO estimation to compute the full gradient without storing intermediate activations,
so even when computing over the full batch at once, the peak dynamic memory is only $\max_{l}|\mathbf{x}_{l}|$,
the same as ZO-SGD, and no compute–memory trade-off is required.
This makes VAMO highly memory-efficient for large mini-batches or large models.

Regarding optimizer state, adaptive methods such as Adam (Kingma & Ba, 2014) or Adagrad (Duchi et al., 2011) require additional memory for first- or second-order momentum,
resulting in $2|\mathbf{x}|$–$3|\mathbf{x}|$ extra memory per step. FO-SVRG also increases optimizer state by $2|\mathbf{x}|$ for storing the snapshot and full gradient (Johnson & Zhang, 2013). VAMO, on the other hand, only needs $|\mathbf{x}|$ for the full ZO gradient used in variance reduction, which is substantially smaller than Adam, Adagrad, or FO-SVRG. In memory-constrained scenarios, this makes VAMO a more favorable choice for fine-tuning large models.

### A.2 Empirical Results Analysis

To empirically validate our theoretical analysis, we conducted full-parameter fine-tuning of a pre-trained RoBERTa-large model (Liu et al., 2019) and GPT-2 medium model (Radford et al., 2019) on the MNLI dataset (Williams et al., 2017) with FP32 (see Section [6](https://arxiv.org/html/2505.13954v2#S6)). We measured the peak GPU memory consumption of VAMO, FO-SGD, FO-Adagrad and FO-Adam by varying the batch size (with sequence length fixed at 128). The detailed results are presented in Table [2](https://arxiv.org/html/2505.13954v2#S6.T2) in Section [6](https://arxiv.org/html/2505.13954v2#S6).

We observe that VAMO’s peak memory consumption is consistently only marginally higher than FO-SGD, with the difference stemming from the overhead of storing the ZO full gradient for variance reduction. Importantly, as the batch size increases, the scaling trend of VAMO almost parallels that of FO-SGD, confirming that the intermediate results for computing mini-batch FO gradients dominates the memory cost in both methods. Notably, the ZO full gradient in VAMO does not require storing additional intermediate activations, further preventing the memory blow-up typically associated with FO-SVRG. In contrast, FO-Adam requires significantly more memory due to the additional storage of first- and second-moment estimates, resulting in an optimizer state of size $3|\mathbf{x}|$. This overhead makes FO-Adam substantially less memory efficient, especially for larger models, whereas VAMO maintains a footprint close to FO-SGD while avoiding the prohibitive costs of FO-SVRG. Taken together, these results highlight VAMO’s practical advantage: for a negligible increase in memory compared to FO-SGD, it achieves markedly better convergence, as shown in Section [6](https://arxiv.org/html/2505.13954v2#S6), while remaining far more memory efficient than FO-Adam. This balance makes VAMO particularly well suited for memory-constrained large-scale fine-tuning tasks.

## Appendix B Detailed Discussing of Convergence and Complexity

In this section, we provide additional discussion of the convergence rates and computational complexities summarized in Table [1](https://arxiv.org/html/2505.13954v2#S2.T1) (see Section [2](https://arxiv.org/html/2505.13954v2#S2)).
Table [1](https://arxiv.org/html/2505.13954v2#S2.T1) summarizes the convergence rates and computational complexities of our proposed methods, referred to as VAMO and VAMO (multi-point) in the table alongside several FO and ZO algorithms.

For ZO methods, ZO-SVRG is listed with a complexity of $O(nS+bT)$ function queries. Among FO methods, FO-SGD has the lowest computational cost $O(bdT)$ but also exhibits the slowest convergence rate of $O(1/\sqrt{T})$. FO-SVRG improves convergence to $O(1/T)$ but increases the cost to $O(dnS+bdT)$ due to full gradient computations. Our proposed VAMO maintains a complexity of $O(nS+bdT)$, similar to ZO-SVRG in terms of $nS$ but replacing the $ndS$ full gradient cost of FO-SVRG with a cheaper $nS$ ZO estimation cost for snapshots, while achieving a fast $O(1/T+1/b)$ convergence rate. This makes its computational complexity significantly slower than FO-SVRG, especially when $d$ is large. The VAMO (multi-point) variant has a complexity of $O(qnS+bdT)$. Here, increasing $q$ (the number of ZO sampling directions) leads to higher complexity but also improves the convergence rate to $O(1/T+(1-q/d)^{2}/b)$, reducing the $O(1/b)$ error term and making its performance more comparable to FO-SVRG, particularly if $q\ll d$. This demonstrates that our proposed methods provide a flexible and often more efficient trade-off between computational cost and convergence performance compared to existing pure FO or ZO approaches. Our work further develops such an adaptive hybrid approach by specifically integrating ZO estimation within the SVRG structure, aiming to reduce the full-gradient cost while preserving strong convergence guarantees independent of dimensionality.

## Appendix C ZO gradient estimator

###### Lemma 1 .

Under the assumptions in Section [3.1](https://arxiv.org/html/2505.13954v2#S3.SS1), and define $f_{\mu}=\mathbb{E}_{\mathbf{u}\sim U_{b}}[f(\mathbf{x}+\mu\mathbf{u})]$ where $U_{b}$ is the uniform distribution over the unit Euclidean ball. Then:

- (i)
$f_{\mu}$ is $L$-smooth with
$\nabla f_{\mu}(\mathbf{x})=\mathbb{E}_{\mathbf{u}}\bigl[\hat{\nabla}f(\mathbf{x})\bigr].$
(12)
- (ii)
For any $\mathbf{x}\in\mathbb{R}^{d}$:
$|f_{\mu}(\mathbf{x})-f(\mathbf{x})|\leq\frac{L\mu^{2}}{2},$
(13)
$\|\nabla f_{\mu}(\mathbf{x})-\nabla f(\mathbf{x})\|_{2}^{2}\leq\frac{\mu^{2}L^{2}d^{2}}{4},$
(14)
$\frac{1}{2}\|\nabla f(\mathbf{x})\|_{2}^{2}-\frac{\mu^{2}L^{2}d^{2}}{4}\leq\|\nabla f_{\mu}(\mathbf{x})\|_{2}^{2}\leq 2\|\nabla f(\mathbf{x})\|_{2}^{2}+\frac{\mu^{2}L^{2}d^{2}}{2}.$
(15)
- (iii)
For any $\mathbf{x}\in\mathbb{R}^{d}$:
$\mathbb{E}_{\mathbf{u}}\bigl[\|\hat{\nabla}f(\mathbf{x})-\nabla f_{\mu}(\mathbf{x})\|_{2}^{2}\bigr]\leq 2d\|\nabla f(\mathbf{x})\|_{2}^{2}+\frac{\mu^{2}L^{2}d^{2}}{2}.$
(16)

###### Proof.

See the proof of Lemma 1 in (Liu et al., 2018b)
∎

###### Lemma 2 .

Under the conditions of Lemma [1](https://arxiv.org/html/2505.13954v2#Thmlemma1):

- (i)
For any $\mathbf{x}\in\mathbb{R}^{d}$:
$\nabla f_{\mu}(\mathbf{x})=\mathbb{E}_{\mathbf{u}}\bigl[\hat{\nabla}f(\mathbf{x})\bigr].$
(17)
where $\hat{\nabla}f(\mathbf{x})$ is the multi-point gradient estimate.
- (ii)
For any $\mathbf{x}\in\mathbb{R}^{d}$:
$\mathbb{E}\bigl[\|\hat{\nabla}f(\mathbf{x})-\nabla f_{\mu}(\mathbf{x})\|_{2}^{2}\bigr]\leq\frac{2d}{q}\|\nabla f(\mathbf{x})\|_{2}^{2}+\frac{\mu^{2}L^{2}d^{2}}{2q}.$
(18)

###### Proof.

See the proof of Lemma 2 in (Liu et al., 2018b)
∎

###### Lemma 1 .

###### Proof.

###### Lemma 2 .

###### Proof.

## Appendix D Second-Order Moment of the Hybrid Gradient Estimator

The primary goal of our convergence analysis is to establish theoretical guarantees for VAMO in solving non-convex optimization problems. Specifically, we aim to bound the expected squared norm of the gradient, $\mathbb{E}[\|\nabla f(\bar{\mathbf{x}})\|_{2}^{2}]$, as shown in Theorem [1](https://arxiv.org/html/2505.13954v2#Thmtheorem1). Due to the hybrid structure of the gradient estimator $\mathbf{v}_{k}^{s}$ used in VAMO, directly analyzing the final convergence metric is challenging. As a key intermediate step, we first derive an upper bound on the second-order moment $\mathbb{E}[\|\mathbf{v}_{k}^{s}\|_{2}^{2}]$.

###### Proposition 1 .

Under the assumptions in Section [3.1](https://arxiv.org/html/2505.13954v2#S3.SS1), and two-point ZO gradient estimate is used in Algorithm [1](https://arxiv.org/html/2505.13954v2#alg1). The blended gradient $\mathbf{v}_{k}^{s}$ in Step 7 of Algorithm [1](https://arxiv.org/html/2505.13954v2#alg1) satisfies,

$$ $\begin{split}\mathbb{E}\bigl[\|{\mathbf{v}}_{k}^{s}\|_{2}^{2}\bigr]&\leq 4\biggl(2\alpha^{2}-2\alpha+1+\frac{24d\delta_{n}}{b}\alpha^{2}\biggr)\mathbb{E}\bigl[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\bigr]\\ &\quad+\frac{12\delta_{n}(4d+1)L^{2}}{b}\alpha^{2}\mathbb{E}\bigl[\|\mathbf{x}_{0}^{s}-\mathbf{x}_{k}^{s}\|_{2}^{2}\bigr]\\ &\quad+\frac{9\delta_{n}}{b}d^{2}L^{2}\mu^{2}\alpha^{2}+\frac{4\sigma^{2}}{b}\biggl(24d\delta_{n}\alpha^{2}+(1-\alpha)^{2}\biggr),\end{split}$ (19) $$

where $\delta_{n}=1$ if the mini-batch contains $i.i.d$. samples from $[n]$ with replacement, and $\delta_{n}=I(b<n)$ if samples are randomly selected without replacement. Here $I(b<n)$ is 1 if $b<n$, and 0 if $b=n$.

###### Proof.

In Algorithm [1](https://arxiv.org/html/2505.13954v2#alg1), we recall that the mini-batch $\mathcal{I}$ is chosen uniformly randomly (with replacement). It is known from Lemma [1](https://arxiv.org/html/2505.13954v2#Thmlemma1) and Lemma [3](https://arxiv.org/html/2505.13954v2#Thmlemma3) that

$$ $\displaystyle\mathbb{E}_{\mathcal{I}_{k}}\left[\nabla{f}_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})-\hat{\nabla}{f}_{\mathcal{I}_{k}}(\mathbf{x}_{0}^{s})\right]=\nabla f(\mathbf{x}_{k}^{s})-\hat{\nabla}f(\mathbf{x}_{0}^{s}).$ (20) $$

We then rewrite $\mathbf{v}_{k}^{s}$ as

$$ $\displaystyle\mathbf{v}_{k}^{s}$ $\displaystyle=\left(1-\alpha\right)\nabla{f}_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})$ (21) $\displaystyle+\alpha\left(\nabla{f}_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})-\hat{\nabla}{f}_{\mathcal{I}_{k}}(\mathbf{x}_{0}^{s})-\mathbb{E}_{\mathcal{I}_{k}}\left[\nabla{f}_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})-\hat{\nabla}{f}_{\mathcal{I}_{k}}(\mathbf{x}_{0}^{s})\right]+\nabla f(\mathbf{x}_{k}^{s})\right)$ $$

Taking the expectation of $\left\|\mathbf{v}_{k}^{s}\right\|_{2}^{2}$ with respect to all the random variables, we have

$$ $\displaystyle\mathbb{E}\left[\left\|\mathbf{v}_{k}^{s}\right\|_{2}^{2}\right]$ $\displaystyle\leq 2\left(1-\alpha\right)^{2}\mathbb{E}\left[\left\|\nabla{f}_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})\right\|_{2}^{2}\right]$ (22) $\displaystyle+2\alpha^{2}\mathbb{E}\left[\left\|\nabla{f}_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})-\hat{\nabla}{f}_{\mathcal{I}_{k}}(\mathbf{x}_{0}^{s})-\mathbb{E}_{\mathcal{I}_{k}}\left[\nabla{f}_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})-\hat{\nabla}{f}_{\mathcal{I}_{k}}(\mathbf{x}_{0}^{s})\right]+\nabla f(\mathbf{x}_{k}^{s})\right\|_{2}^{2}\right]$ $\displaystyle\leq 4\alpha^{2}\mathbb{E}\left[\left\|\nabla{f}_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})-\hat{\nabla}{f}_{\mathcal{I}_{k}}(\mathbf{x}_{0}^{s})-\mathbb{E}_{\mathcal{I}_{k}}\left[\nabla{f}_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})-\hat{\nabla}{f}_{\mathcal{I}_{k}}(\mathbf{x}_{0}^{s})\right]\right\|_{2}^{2}\right]$ $\displaystyle+4\alpha^{2}\mathbb{E}\left[\left\|\nabla f(\mathbf{x}_{k}^{s})\right\|_{2}^{2}\right]+2\left(1-\alpha\right)^{2}\mathbb{E}\left[\left\|\nabla{f}_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})\right\|_{2}^{2}\right]$ $$

where the first inequality holds due to Lemma [4](https://arxiv.org/html/2505.13954v2#Thmlemma4). Based on equation [20](https://arxiv.org/html/2505.13954v2#A4.E20), we note that the following holds

$$ $\displaystyle\sum_{i=1}^{n}\left\{\nabla f_{i}(\mathbf{x}_{k}^{s})-\hat{\nabla}f_{i}(\mathbf{x}_{0}^{s})-\mathbb{E}_{\mathcal{I}_{k}}\left[\nabla f_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})-\hat{\nabla}f_{\mathcal{I}_{k}}(\mathbf{x}_{0}^{s})\right]\right\}$ (23) $\displaystyle=n(\nabla f(\mathbf{x}_{k}^{s})-\hat{\nabla}f(\mathbf{x}_{0}^{s}))-n(\nabla f(\mathbf{x}_{k}^{s})-\hat{\nabla}f(\mathbf{x}_{0}^{s}))=0.$ $$

Based on equation [23](https://arxiv.org/html/2505.13954v2#A4.E23) and applying Lemma [1](https://arxiv.org/html/2505.13954v2#Thmlemma1) and Lemma [3](https://arxiv.org/html/2505.13954v2#Thmlemma3), the first term at the right hand side (RHS) of equation [22](https://arxiv.org/html/2505.13954v2#A4.E22) yields

$$ $\displaystyle\mathbb{E}\left[\left\|\nabla f_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})-\hat{\nabla}f_{\mathcal{I}_{k}}(\mathbf{x}_{0}^{s})-\mathbb{E}_{\mathcal{I}_{k}}\left[\nabla f_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})-\hat{\nabla}f_{\mathcal{I}_{k}}(\mathbf{x}_{0}^{s})\right]\right\|_{2}^{2}\right]$ (24) $\displaystyle\leq\frac{\delta_{n}}{bn}\sum_{i=1}^{n}\mathbb{E}\left[\|{\nabla}f_{i}(\mathbf{x}_{k}^{s})-\hat{\nabla}f_{i}(\mathbf{x}_{0}^{s})-({\nabla}f(\mathbf{x}_{k}^{s})-\hat{\nabla}f(\mathbf{x}_{0}^{s}))\|_{2}^{2}\right]$ $\displaystyle=\mathbb{E}\left[\frac{\delta_{n}}{b}\left(\frac{1}{n}\sum_{i=1}^{n}\|{\nabla}f_{i}(\mathbf{x}_{k}^{s})-\hat{\nabla}f_{i}(\mathbf{x}_{0}^{s})\|_{2}^{2}-\|{\nabla}f(\mathbf{x}_{k}^{s})-\hat{\nabla}f(\mathbf{x}_{0}^{s})\|_{2}^{2}\right)\right]$ $\displaystyle\leq\frac{\delta_{n}}{bn}\sum_{i=1}^{n}\mathbb{E}\left[\left\|\nabla f_{i}(\mathbf{x}_{k}^{s})-\hat{\nabla}f_{i}(\mathbf{x}_{0}^{s})\right\|_{2}^{2}\right].$ $$

where the first inequality holds due to Lemma [1](https://arxiv.org/html/2505.13954v2#Thmlemma1) and Lemma [3](https://arxiv.org/html/2505.13954v2#Thmlemma3) (taking the expectation with respect to mini-batch $\mathcal{I}$), we define $\delta_{n}$ as

$$ $\displaystyle\delta_{n}=\left\{\begin{array}[]{ll}1&\text{if }\mathcal{I}\text{ contains i.i.d. samples with replacement (Lemma~\ref{lem:lemma3})}\\ I(b<n)&\text{if }\mathcal{I}\text{ contains samples without replacement (Lemma~\ref{lem:lemma4})}.\end{array}\right.$ (25) $$

Substituting equation [24](https://arxiv.org/html/2505.13954v2#A4.E24) into equation [22](https://arxiv.org/html/2505.13954v2#A4.E22), we obtain

$$ $\displaystyle\mathbb{E}\left[\left\|\mathbf{v}_{k}^{s}\right\|_{2}^{2}\right]$ $\displaystyle\leq 2\left(1-\alpha\right)^{2}\mathbb{E}\left[\left\|\nabla{f}_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})\right\|_{2}^{2}\right]$ (26) $\displaystyle+\frac{4\alpha^{2}\delta_{n}}{bn}\sum_{i=1}^{n}\mathbb{E}\left[\left\|\nabla f_{i}(\mathbf{x}_{k}^{s})-\hat{\nabla}f_{i}(\mathbf{x}_{0}^{s})\right\|_{2}^{2}\right]+4\alpha^{2}\mathbb{E}\left[\left\|\nabla f(\mathbf{x}_{k}^{s})\right\|_{2}^{2}\right].$ $$

Similar to Lemma [1](https://arxiv.org/html/2505.13954v2#Thmlemma1), we introduce a smoothing function $f_{i,\mu}$ of $f_{i}$, and continue to bound the second term at the right hand side (RHS) of equation [26](https://arxiv.org/html/2505.13954v2#A4.E26). This yields

$$ $\displaystyle\mathbb{E}\left[\|{\nabla}f_{i}(\mathbf{x}_{k}^{s})-\hat{\nabla}f_{i}(\mathbf{x}_{0}^{s})\|_{2}^{2}\right]$ (27) $\displaystyle\leq 3\mathbb{E}\left[\|{\nabla}f_{i}(\mathbf{x}_{k}^{s})-\nabla f_{i,\mu}(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]+3\mathbb{E}\left[\|{\nabla}f_{i,\mu}(\mathbf{x}_{0}^{s})-\hat{\nabla}f_{i}(\mathbf{x}_{0}^{s})\|_{2}^{2}\right]$ $\displaystyle+3\mathbb{E}\left[\|\nabla f_{i,\mu}(\mathbf{x}_{k}^{s})-\nabla f_{i,\mu}(\mathbf{x}_{0}^{s})\|_{2}^{2}\right]$ $\displaystyle\leq 6d\mathbb{E}[\|\nabla f_{i}(\mathbf{x}_{0}^{s})\|_{2}^{2}]+\frac{9}{4}L^{2}d^{2}\mu^{2}+3\mathbb{E}\left[\|\nabla f_{i,\mu}(\mathbf{x}_{k}^{s})-\nabla f_{i,\mu}(\mathbf{x}_{0}^{s})\|_{2}^{2}\right]$ $$

Since both $f_{i}$ and $f_{i,\mu}$ are L-smooth (Lemma [1](https://arxiv.org/html/2505.13954v2#Thmlemma1)), we have

$$ $\displaystyle\mathbb{E}\left[\|\nabla f_{i,\mu}(\mathbf{x}_{k}^{s})-\nabla f_{i,\mu}(\mathbf{x}_{0}^{s})\|_{2}^{2}\right]\leq L^{2}\mathbb{E}\left[\|\mathbf{x}_{k}^{s}-\mathbf{x}_{0}^{s}\|_{2}^{2}\right],$ (28) $\displaystyle\mathbb{E}\left[\|\nabla f_{i}(\mathbf{x}_{0}^{s})\|_{2}^{2}\right]\leq 2\mathbb{E}\left[\|\nabla f_{i}(\mathbf{x}_{0}^{s})-\nabla f_{i}(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]+2\mathbb{E}\left[\|\nabla f_{i}(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]$ $\displaystyle\leq 2L^{2}\mathbb{E}\left[\|\mathbf{x}_{0}^{s}-\mathbf{x}_{k}^{s}\|_{2}^{2}\right]+2\mathbb{E}\left[\|\nabla f_{i}(\mathbf{x}_{k}^{s})\|_{2}^{2}\right].$ $$

We obtain

$$ $\displaystyle\mathbb{E}\left[\|{\nabla}f_{i}(\mathbf{x}_{k}^{s})-\hat{\nabla}f_{i}(\mathbf{x}_{0}^{s})\|_{2}^{2}\right]$ (29) $\displaystyle\leq 2d\mathbb{E}[\|\nabla f_{i}(\mathbf{x}_{k}^{s})\|_{2}^{2}]+(2d+3)L^{2}\mathbb{E}\left[\|\mathbf{x}_{0}^{s}-\mathbf{x}_{k}^{s}\|_{2}^{2}\right]+\frac{9}{4}L^{2}d^{2}\mu^{2}$ $\displaystyle\leq 4d\mathbb{E}\left[\|\nabla f_{i}(\mathbf{x}_{k}^{s})-\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]+4d\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]$ $\displaystyle+(2d+3)L^{2}\mathbb{E}\left[\|\mathbf{x}_{0}^{s}-\mathbf{x}_{k}^{s}\|_{2}^{2}\right]+\frac{9}{4}L^{2}d^{2}\mu^{2}$ $\displaystyle\leq 4d\sigma^{2}+4d\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]+(2d+3)L^{2}\mathbb{E}\left[\|\mathbf{x}_{0}^{s}-\mathbf{x}_{k}^{s}\|_{2}^{2}\right]+\frac{9}{4}L^{2}d^{2}\mu^{2},$ $$

where the last inequality holds due to Assumption in Section [3.1](https://arxiv.org/html/2505.13954v2#S3.SS1).
We bound the first term at the right hand side (RHS) of equation [26](https://arxiv.org/html/2505.13954v2#A4.E26). This yields

$$ $\displaystyle\mathbb{E}\left[\left\|\nabla{f}_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})\right\|_{2}^{2}\right]$ $\displaystyle\leq 2\mathbb{E}\left[\left\|\nabla{f}_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})-\nabla{f}(\mathbf{x}_{k}^{s})\right\|_{2}^{2}\right]+2\mathbb{E}\left[\left\|\nabla{f}(\mathbf{x}_{k}^{s})\right\|_{2}^{2}\right]$ (30) $\displaystyle\leq\frac{2}{b}\sigma^{2}+2\mathbb{E}\left[\left\|\nabla{f}(\mathbf{x}_{k}^{s})\right\|_{2}^{2}\right]$ $$

Therefore, we have

$$ $\displaystyle\mathbb{E}\left[\|{\mathbf{v}}_{k}^{s}\|_{2}^{2}\right]$ $\displaystyle\leq\frac{4(1-\alpha)^{2}}{b}\sigma^{2}+4(1-\alpha)^{2}\mathbb{E}\left[\left\|\nabla{f}(\mathbf{x}_{k}^{s})\right\|_{2}^{2}\right]$ (31) $\displaystyle+\frac{12\delta_{n}(4d+1)L^{2}}{b}\alpha^{2}\mathbb{E}\left[\|\mathbf{x}_{0}^{s}-\mathbf{x}_{k}^{s}\|_{2}^{2}\right]+\left(4+\frac{96d\delta_{n}}{b}\right)\alpha^{2}\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]$ $\displaystyle+\frac{9\delta_{n}}{b}d^{2}L^{2}\mu^{2}\alpha^{2}+\frac{96d\sigma^{2}\delta_{n}}{b}\alpha^{2}.$ $\displaystyle=4\left(2\alpha^{2}-2\alpha+1+\frac{24d\delta_{n}}{b}\alpha^{2}\right)\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]$ $\displaystyle+\frac{12\delta_{n}(4d+1)L^{2}}{b}\alpha^{2}\mathbb{E}\left[\|\mathbf{x}_{0}^{s}-\mathbf{x}_{k}^{s}\|_{2}^{2}\right]$ $\displaystyle+\frac{9\delta_{n}}{b}d^{2}L^{2}\mu^{2}\alpha^{2}+\frac{4\sigma^{2}}{b}\left(24d\delta_{n}\alpha^{2}+(1-\alpha)^{2}\right).$ $$

∎

The bound on $\mathbb{E}[\|\mathbf{v}_{k}^{s}\|_{2}^{2}]$, detailed in Proposition [1](https://arxiv.org/html/2505.13954v2#Thmproposition1), plays a central role in our analysis. It enables us to control the error accumulation during the optimization process and ultimately leads to the convergence rate stated in Theorem [1](https://arxiv.org/html/2505.13954v2#Thmtheorem1). Based on Proposition [1](https://arxiv.org/html/2505.13954v2#Thmproposition1), Theorem [1](https://arxiv.org/html/2505.13954v2#Thmtheorem1) provides the convergence rate of VAMO in terms of an upper bound on $\mathbb{E}\bigl[\|\nabla{f}(\bar{\mathbf{x}})\|_{2}^{2}\bigr]$ at the solution $\bar{\mathbf{x}}$.

###### Proposition 1 .

###### Proof.

## Appendix E Proof of Theorem 1

###### Proof.

Since $f$ is L-smooth (Lemma [1](https://arxiv.org/html/2505.13954v2#Thmlemma1)), from Lemma [5](https://arxiv.org/html/2505.13954v2#Thmlemma5) we have

$$ $\displaystyle f(\mathbf{x}_{k}^{s+1})$ $\displaystyle\leq f(\mathbf{x}_{k}^{s})+\langle\nabla f(\mathbf{x}_{k}^{s}),\mathbf{x}_{k}^{s+1}-\mathbf{x}_{k}^{s}\rangle+\frac{L}{2}\|\mathbf{x}_{k}^{s+1}-\mathbf{x}_{k}^{s}\|_{2}^{2}$ (32) $\displaystyle=f(\mathbf{x}_{k}^{s})-\eta_{k}\langle\nabla f(\mathbf{x}_{k}^{s}),\mathbf{v}_{k}^{s}\rangle+\frac{L}{2}\eta_{k}^{2}\|\mathbf{v}_{k}^{s}\|_{2}^{2}$ $$

where the last equality holds due to $x_{k+1}^{s}=x_{k}^{s}-\eta_{k}v_{k}^{s}$. Since $x_{k}^{s}$ and $x_{0}^{s}$ are independent of $\mathcal{I}$ and random directions $u$ used for ZO gradient estimates, from equation [12](https://arxiv.org/html/2505.13954v2#A3.E12) we obtain

$$ $\displaystyle\mathbb{E}_{\mathbf{u},\mathcal{I}_{k}}\left[{\mathbf{v}}_{k}^{s}\right]=$ $\displaystyle\mathbb{E}_{\mathbf{u},\mathcal{I}_{k}}\left[{\nabla}f_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})-\alpha\left(\hat{\nabla}f_{\mathcal{I}_{k}}(\mathbf{x}_{0}^{s})-\hat{\nabla}f(\mathbf{x}_{0}^{s})\right)\right]$ (33) $\displaystyle=$ $\displaystyle\nabla f(\mathbf{x}_{k}^{s})-\alpha\left(\nabla f_{\mu}(\mathbf{x}_{0}^{s})-\nabla f_{\mu}(\mathbf{x}_{0}^{s})\right)=\nabla f(\mathbf{x}_{k}^{s}).$ $$

Combining equation [32](https://arxiv.org/html/2505.13954v2#A5.E32) and equation [33](https://arxiv.org/html/2505.13954v2#A5.E33), we have

$$ $\displaystyle\mathbb{E}\left[f(\mathbf{x}_{k+1}^{s})\right]\leq\mathbb{E}\left[f(\mathbf{x}_{k}^{s})\right]-\eta_{k}\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]+\frac{L}{2}\eta_{k}^{2}\mathbb{E}\left[\|{\mathbf{v}}_{k}^{s}\|_{2}^{2}\right],$ (34) $$

where the expectation is taken with respect to all random variables.
At RHS of equation [34](https://arxiv.org/html/2505.13954v2#A5.E34), the upper bound on $\mathbb{E}\left[\|{\mathbf{v}}_{k}^{s}\|_{2}^{2}\right]$ is given by Proposition [1](https://arxiv.org/html/2505.13954v2#Thmproposition1),

$$ $\displaystyle\mathbb{E}\left[\|{\mathbf{v}}_{k}^{s}\|_{2}^{2}\right]$ $\displaystyle\leq 4\left(2\alpha^{2}-2\alpha+1+\frac{24d\delta_{n}}{b}\alpha^{2}\right)\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]$ (35) $\displaystyle+\frac{12\delta_{n}(4d+1)L^{2}}{b}\alpha^{2}\mathbb{E}\left[\|\mathbf{x}_{0}^{s}-\mathbf{x}_{k}^{s}\|_{2}^{2}\right]$ $\displaystyle+\frac{9\delta_{n}}{b}d^{2}L^{2}\mu^{2}\alpha^{2}+\frac{4\sigma^{2}}{b}\left(24d\delta_{n}\alpha^{2}+(1-\alpha)^{2}\right).$ $$

In equation [34](https://arxiv.org/html/2505.13954v2#A5.E34), we further bound $\mathbb{E}\left[\|\mathbf{x}_{k+1}^{s}-\mathbf{x}_{0}^{s}\|_{2}^{2}\right]$ as,

$$ $\displaystyle\mathbb{E}\left[\|\mathbf{x}_{k+1}^{s}-\mathbf{x}_{0}^{s}\|_{2}^{2}\right]=\mathbb{E}\left[\|\mathbf{x}_{k+1}^{s}-\mathbf{x}_{k}^{s}+\mathbf{x}_{k}^{s}-\mathbf{x}_{0}^{s}\|_{2}^{2}\right]$ (36) $\displaystyle=\eta_{k}^{2}\mathbb{E}\left[\|{\mathbf{v}}_{k}^{s}\|_{2}^{2}\right]+\mathbb{E}\left[\|\mathbf{x}_{k}^{s}-\mathbf{x}_{0}^{s}\|_{2}^{2}\right]-2\eta_{k}\mathbb{E}\left[\langle{\mathbf{v}}_{k}^{s},\mathbf{x}_{k}^{s}-\mathbf{x}_{0}^{s}\rangle\right]$ $\displaystyle=\eta_{k}^{2}\mathbb{E}\left[\|{\mathbf{v}}_{k}^{s}\|_{2}^{2}\right]+\mathbb{E}\left[\|\mathbf{x}_{k}^{s}-\mathbf{x}_{0}^{s}\|_{2}^{2}\right]-2\eta_{k}\mathbb{E}\left[\langle\nabla f(\mathbf{x}_{k}^{s}),\mathbf{x}_{k}^{s}-\mathbf{x}_{0}^{s}\rangle\right]$ $\displaystyle\leq\eta_{k}^{2}\mathbb{E}\left[\|{\mathbf{v}}_{k}^{s}\|_{2}^{2}\right]+\mathbb{E}\left[\|\mathbf{x}_{k}^{s}-\mathbf{x}_{0}^{s}\|_{2}^{2}\right]+2\eta_{k}\mathbb{E}\left[\frac{1}{2\beta_{k}}\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}+\frac{\beta_{k}}{2}\|\mathbf{x}_{k}^{s}-\mathbf{x}_{0}^{s}\|_{2}^{2}\right],$ $$

We introduce a Lyapunov function with respect to $f_{\mu}$,

$$ $\displaystyle R_{k}^{s}=\mathbb{E}\left[f(\mathbf{x}_{k}^{s})+c_{k}\|\mathbf{x}_{k}^{s}-\mathbf{x}_{0}^{s}\|_{2}^{2}\right],$ (37) $$

for some $c_{k}>0$, Substituting equation [34](https://arxiv.org/html/2505.13954v2#A5.E34) and equation [36](https://arxiv.org/html/2505.13954v2#A5.E36) into $R_{k+1}^{s}$, we obtain

$$ $\displaystyle R_{k+1}^{s}$ $\displaystyle=\mathbb{E}\left[f(\mathbf{x}_{k+1}^{s})+c_{k+1}\|\mathbf{x}_{k+1}^{s}-\mathbf{x}_{0}^{s}\|_{2}^{2}\right]$ (38) $\displaystyle\leq\mathbb{E}\left[f(\mathbf{x}_{k}^{s})-\eta_{k}\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}+\frac{L}{2}\eta_{k}^{2}\|{\mathbf{v}}_{k}^{s}\|_{2}^{2}\right]+\mathbb{E}\left[c_{k+1}\eta_{k}^{2}\|{\mathbf{v}}_{k}^{s}\|_{2}^{2}+c_{k+1}\|\mathbf{x}_{k}^{s}-\mathbf{x}_{0}^{s}\|_{2}^{s}\right]$ $\displaystyle+\mathbb{E}\left[\frac{c_{k+1}\eta_{k}}{\beta_{k}}\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}+c_{k+1}\beta_{k}\eta_{k}\|\mathbf{x}_{k}^{s}-\mathbf{x}_{0}^{s}\|_{2}^{2}\right]$ $\displaystyle=\mathbb{E}\left[f(\mathbf{x}_{k}^{s})\right]-\left(\eta_{k}-\frac{c_{k+1}\eta_{k}}{\beta_{k}}\right)\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]$ $\displaystyle+\left(c_{k+1}+c_{k+1}\beta_{k}\eta_{k}\right)\mathbb{E}\left[\|\mathbf{x}_{k}^{s}-\mathbf{x}_{0}^{s}\|_{2}^{2}\right]+\left(\frac{L}{2}\eta_{k}^{2}+c_{k+1}\eta_{k}^{2}\right)\mathbb{E}\left[\|{\mathbf{v}}_{k}^{s}\|_{2}^{2}\right].$ $$

Moreover, substituting equation [35](https://arxiv.org/html/2505.13954v2#A5.E35) into equation [38](https://arxiv.org/html/2505.13954v2#A5.E38), we have

$$ $\displaystyle R_{k+1}^{s}\leq$ $\displaystyle\mathbb{E}\left[f(\mathbf{x}_{k}^{s})\right]\bm{-}\left(\eta_{k}-\frac{c_{k\bm{+}1}\eta_{k}}{\beta_{k}}\right)\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]+\left(c_{k\bm{+}1}+c_{k\bm{+}1}\beta_{k}\eta_{k}\right)\mathbb{E}\left[\|\mathbf{x}_{k}^{s}-\mathbf{x}_{0}^{s}\|_{2}^{2}\right]$ (39) $\displaystyle+\left(\frac{L}{2}\eta_{k}^{2}+c_{k+1}\eta_{k}^{2}\right)\frac{12(4d+1)L^{2}\delta_{n}}{b}\alpha^{2}\mathbb{E}\left[\|\mathbf{x}_{k}^{s}-\mathbf{x}_{0}^{s}\|_{2}^{2}\right]$ $\displaystyle+4\left(\frac{L}{2}\eta_{k}^{2}+c_{k+1}\eta_{k}^{2}\right)\left(2\alpha^{2}-2\alpha+1+\frac{24d\delta_{n}}{b}\alpha^{2}\right)\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]$ $\displaystyle+\left(\frac{L}{2}\eta_{k}^{2}+c_{k+1}\eta_{k}^{2}\right)\left(\frac{9\delta_{n}}{b}d^{2}L^{2}\mu^{2}\alpha^{2}+\frac{4\sigma^{2}}{b}\left(24d\delta_{n}\alpha^{2}+(1-\alpha)^{2}\right)\right).$ $$

The definition of $c_{k}$ is given by

$$ $c_{k}=c_{k+1}+\beta_{k}\eta_{k}c_{k+1}+\left(\frac{L}{2}\eta_{k}^{2}+c_{k+1}\eta_{k}^{2}\right)\frac{12(4d+1)L^{2}\delta_{n}}{b}\alpha^{2}$ (40) $$

Based on the definition of $c_{k}$ and the definition of $R_{k}^{s}$ in equation [37](https://arxiv.org/html/2505.13954v2#A5.E37), we can simplify the inequality equation [39](https://arxiv.org/html/2505.13954v2#A5.E39) as

$$ $\displaystyle R_{k+1}^{s}$ $\displaystyle\leq R_{k}^{s}-\left(\eta_{k}-\frac{c_{k+1}\eta_{k}}{\beta_{k}}\right)\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]$ (41) $\displaystyle+4\left(\frac{L}{2}\eta_{k}^{2}+c_{k+1}\eta_{k}^{2}\right)\left(2\alpha^{2}-2\alpha+1+\frac{24d\delta_{n}}{b}\alpha^{2}\right)\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]$ $\displaystyle+\left(\frac{L}{2}\eta_{k}^{2}+c_{k+1}\eta_{k}^{2}\right)\left(\frac{9\delta_{n}}{b}d^{2}L^{2}\mu^{2}\alpha^{2}+\frac{4\sigma^{2}}{b}\left(24d\delta_{n}\alpha^{2}+(1-\alpha)^{2}\right)\right)$ $\displaystyle=R_{k}^{s}-\gamma_{k}\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]+\chi_{k},$ $$

where $\gamma_{k}$ and $\chi_{k}$ are coefficients given by

$$ $\displaystyle\gamma_{k}=\left(1-\frac{c_{k+1}}{\beta_{k}}\right)\eta_{k}-4\left(\frac{L}{2}+c_{k+1}\right)\left(2\alpha^{2}-2\alpha+1+\frac{24d\delta_{n}}{b}\alpha^{2}\right)\eta_{k}^{2},$ (42) $\displaystyle\chi_{k}=\left(\frac{L}{2}+c_{k+1}\right)\left(\frac{9\delta_{n}}{b}d^{2}L^{2}\mu^{2}\alpha^{2}+\frac{4\sigma^{2}}{b}\left(24d\delta_{n}\alpha^{2}+(1-\alpha)^{2}\right)\right)\eta_{k}^{2}$ $$

Taking a telescopic sum for equation [42](https://arxiv.org/html/2505.13954v2#A5.E42), we obtain

$$ $\displaystyle R_{m}^{s}\leq R_{0}^{s}-\sum_{k=0}^{m-1}\gamma_{k}\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]+\chi_{m},$ (43) $$

where $\chi_{m}=\sum_{k=0}^{m-1}\chi_{k}$. It is known from equation [37](https://arxiv.org/html/2505.13954v2#A5.E37) that,

$$ $\displaystyle R_{0}^{s}=\mathbb{E}\left[f(\mathbf{x}_{0}^{s})\right],\quad R_{m}^{s}=\mathbb{E}\left[f(\mathbf{x}_{m}^{s})\right],$ (44) $$

where the last equality used the fact that $c_{m}=0$, since $\bar{\mathbf{x}}_{s\bm{-}1}=\mathbf{x}_{0}^{s}$ and $\bar{\mathbf{x}}_{s}=\mathbf{x}_{m}^{s}$, we obtain

$$ $\displaystyle R_{0}^{s}-R_{m}^{s}=\mathbb{E}\left[f(\bar{\mathbf{x}}_{s-1})-f(\bar{\mathbf{x}}_{s})\right].$ (45) $$

Telescoping the sum for $s=1,2,\ldots,S$, we obtain,

$$ $\displaystyle\sum_{s=1}^{S}\sum_{k=0}^{m-1}\gamma_{k}\mathbb{E}[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}]\leq\mathbb{E}[f(\bar{\mathbf{x}}_{0})-f(\bar{\mathbf{x}}_{S})]+S\chi_{m}.$ (46) $$

let $\bar{\gamma}=\min_{k}\gamma_{k}$ and we choose $\bar{\mathbf{x}}$ uniformly random from $\{\{\mathbf{x}_{k}^{s}\}_{k=0}^{m-1}\}_{s=1}^{S}$, then we obtain

$$ $\displaystyle\mathbb{E}[\|\nabla f(\bar{\mathbf{x}})\|_{2}^{2}]\leq\frac{\mathbb{E}[f(\bar{\mathbf{x}}_{0})-f^{*}]}{T\bar{\gamma}}+\frac{S\chi_{m}}{T\bar{\gamma}}.$ (47) $$

∎

###### Proof.

## Appendix F Proof of Corollary 1

###### Proof.

We start by rewriting $c_{k}$ in equation [40](https://arxiv.org/html/2505.13954v2#A5.E40) as

$$ $\displaystyle c_{k}=(1+\theta)c_{k+1}+\frac{6(1+4d)L^{3}\delta_{n}\eta^{2}}{b}\alpha^{2}$ (48) $$

where $\theta=\beta\eta+\frac{12(1+4d)L^{2}\delta_{n}\eta^{2}}{b}\alpha^{2}$. The recursive formula equation [48](https://arxiv.org/html/2505.13954v2#A6.E48) implies that $c_{k}\leq c_{0}$ for any $k$, and

$$ $\displaystyle c_{0}=\frac{6(1+4d)L^{3}\delta_{n}\eta^{2}\alpha^{2}}{b}\frac{(1+\theta)^{m}-1}{\theta}.$ (49) $$

Based on the choice of $\eta=\frac{\rho}{L}$, $\alpha=\frac{1}{d}$, and $\beta=L$, we have

$$ $\displaystyle\theta=\rho+\frac{12(4d+1)\delta_{n}\rho^{2}}{bd^{2}}$ (50) $$

where we have used the fact that $\delta_{n}\leq 1$, Substituting equation [50](https://arxiv.org/html/2505.13954v2#A6.E50) into equation [49](https://arxiv.org/html/2505.13954v2#A6.E49), we have

$$ $\displaystyle c_{k}\leq c_{0}$ $\displaystyle=\frac{6(1+4d)L^{3}\delta_{n}\alpha^{2}}{b}\frac{\eta^{2}}{\theta}[(1+\theta)^{m}-1]=\frac{6(1+4d)L\rho\delta_{n}}{bd^{2}+12(4d+1)\delta_{n}\rho}[(1+\theta)^{m}-1]$ (51) $\displaystyle\leq\frac{30L\rho\delta_{n}}{bd}[(1+\theta)^{m}-1]\leq\frac{30L\rho\delta_{n}}{bd}(e-1)\leq\frac{60L\rho\delta_{n}}{bd},$ $$

where the third inequality holds since $(1+\theta)^{m}\leq(1+\frac{31\rho}{d})^{m},(1+1/a)^{a}\leq\lim_{a\to\infty}(1+\frac{1}{a})^{a}=e\mathrm{~for~}a>0$, and the last inequality loosely uses the notion ’$\leq$’ since $e<3$.
We recall from equation [41](https://arxiv.org/html/2505.13954v2#A5.E41) that

$$ $\displaystyle\bar{\gamma}=\min_{0\leq k\leq m-1}\left\{\left(1-\frac{c_{k+1}}{\beta_{k}}\right)\eta_{k}-4\left(\frac{L}{2}+c_{k+1}\right)\left(2\alpha^{2}-2\alpha+1+\frac{24d\delta_{n}}{b}\alpha^{2}\right)\eta_{k}^{2}\right\}.$ (52) $$

Since $\eta_{k}=\eta$, $\beta_{k}=\beta$ and $\eta_{k}=\eta,\beta_{k}=\beta$, we have

$$ $\displaystyle\bar{\gamma}\geq\left(1-\frac{c_{0}}{\beta}\right)\eta-4\left(\frac{L}{2}+c_{0}\right)\left(2\alpha^{2}-2\alpha+1+\frac{24d\delta_{n}}{b}\alpha^{2}\right)\eta^{2}.$ (53) $$

From equation [51](https://arxiv.org/html/2505.13954v2#A6.E51) and the definition of $\beta$, we have

$$ $\displaystyle\frac{c_{0}}{\beta}\leq\frac{60\rho}{bd},$ (54) $$

and

$$ $\displaystyle\left(\frac{L}{2}+c_{0}\right)\left(2\alpha^{2}-2\alpha+1+\frac{24d\delta_{n}}{b}\alpha^{2}\right)\eta$ (55) $\displaystyle\leq\left(\frac{L}{2}+\frac{60L\rho}{bd}\right)\left(\frac{2}{d^{2}}-\frac{2}{d}+1+\frac{24\delta_{n}}{bd}\right)\frac{\rho}{L}$ $\displaystyle\leq\rho\left(1+\frac{24}{bd}\right)$ $$

Substituting equation LABEL:eq:gamma_first_bound and equation [55](https://arxiv.org/html/2505.13954v2#A6.E55) into equation [53](https://arxiv.org/html/2505.13954v2#A6.E53), we obtain

$$ $\displaystyle\bar{\gamma}\geq\eta\left(1-\frac{60\rho}{bd}-4\rho-\frac{96\rho}{bd}\right)\geq\eta\left(1-\frac{156\rho}{bd}-4\rho\right),$ (56) $$

where we have used the fact that $b<d$. Moreover, if we set $\rho\leq\frac{1}{160}$ , then $\bar{\gamma}>0$. In other words, the current parameter setting is valid for Theorem [1](https://arxiv.org/html/2505.13954v2#Thmtheorem1). Upon defining a universal constant $z_{0}=1-\frac{156\rho}{bd}-4\rho$, we have

$$ $\displaystyle\bar{\gamma}\geq\eta z_{0}$ (57) $$

Next, we find the upper bound on $\chi_{m}$ in equation [41](https://arxiv.org/html/2505.13954v2#A5.E41) given the current parameter setting and $c_{k}\leq c_{0}$,

$$ $\displaystyle\chi_{m}\leq m\left(\frac{L}{2}+c_{0}\right)\left(\frac{9\delta_{n}}{b}d^{2}L^{2}\mu^{2}\alpha^{2}+\frac{4\sigma^{2}}{b}\left(24d\delta_{n}\alpha^{2}+(1-\alpha)^{2}\right)\right)\eta^{2}$ (58) $$

Based on $\bar{\gamma}\geq\eta z_{0}$ and $c_{0}\leq 60L\rho\delta_{n}\leq\frac{L}{2}$, we have

$$ $\displaystyle\frac{\chi_{m}}{\bar{\gamma}}\leq m\rho\left(\frac{9\delta_{n}}{bz_{0}}d^{2}L^{2}\mu^{2}\alpha^{2}+\frac{4\sigma^{2}}{bz_{0}}\left(24d\delta_{n}\alpha^{2}+(1-\alpha)^{2}\right)\right)$ (59) $$

since $T=Sm$, and $\mu=\frac{1}{\sqrt{T}}$,the above inequality yields

$$ $\displaystyle\frac{S\chi_{m}}{T\bar{\gamma}}\leq\frac{9\rho L^{2}\delta_{n}}{z_{0}bT}+\frac{4\sigma^{2}}{bz_{0}}\left(\frac{24\delta_{n}}{d}+(1-\frac{1}{d})^{2}\right)=O\left(\frac{1}{Tb}+\frac{1}{b}\right),$ (60) $$

where in the big $O$ notation, we only keep the dominant terms and ignore the constant numbers that are independent of $d$, $b$, and $T$ .
Substituting equation [57](https://arxiv.org/html/2505.13954v2#A6.E57) and equation [60](https://arxiv.org/html/2505.13954v2#A6.E60) into equation [7](https://arxiv.org/html/2505.13954v2#S4.E7), we have

$$ $\displaystyle\mathbb{E}[\|\nabla f(\overline{\mathbf{x}})\|_{2}^{2}]\leq\frac{[f(\bar{\mathbf{x}}_{0})-f^{*}]}{Tz_{0}}\frac{L}{\rho}+\frac{S\chi_{m}}{T\bar{\gamma}}=O\left(\frac{1}{T}+\frac{1}{bT}+\frac{1}{b}\right).$ (61) $$

∎

###### Proof.

## Appendix G Proof of Theorem 2

###### Proof.

Motivated by Proposition [1](https://arxiv.org/html/2505.13954v2#Thmproposition1), we first bound $\left\|\mathbf{v}_{k}^{s}\right\|_{2}^{2}$, Following, we have

$$ $\displaystyle\mathbb{E}\left[\left\|\mathbf{v}_{k}^{s}\right\|_{2}^{2}\right]$ $\displaystyle\leq 2\left(1-\alpha\right)^{2}\mathbb{E}\left[\left\|\nabla{f}_{\mathcal{I}_{k}}(\mathbf{x}_{k}^{s})\right\|_{2}^{2}\right]$ (62) $\displaystyle+\frac{4\alpha^{2}\delta_{n}}{bn}\sum_{i=1}^{n}\mathbb{E}\left[\left\|\nabla f_{i}(\mathbf{x}_{k}^{s})-\hat{\nabla}f_{i}(\mathbf{x}_{0}^{s})\right\|_{2}^{2}\right]+4\alpha^{2}\mathbb{E}\left[\left\|\nabla f(\mathbf{x}_{k}^{s})\right\|_{2}^{2}\right].$ $$

Following together with equation [18](https://arxiv.org/html/2505.13954v2#A3.E18), we can obtain that

$$ $\displaystyle\mathbb{E}\left[\|{\nabla}f_{i}(\mathbf{x}_{k}^{s})-\hat{\nabla}f_{i}(\mathbf{x}_{0}^{s})\|_{2}^{2}\right]$ (63) $\displaystyle\leq\frac{24d}{q}\sigma^{2}+\frac{24d}{q}\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]$ $\displaystyle+\left(3+\frac{12d}{q}\right)L^{2}\mathbb{E}\left[\|\mathbf{x}_{0}^{s}-\mathbf{x}_{k}^{s}\|_{2}^{2}\right]+\left(\frac{3}{4}+\frac{3}{2q}\right)L^{2}d^{2}\mu^{2},$ $$

Substituting equation [63](https://arxiv.org/html/2505.13954v2#A7.E63) and equation [30](https://arxiv.org/html/2505.13954v2#A4.E30) into equation [62](https://arxiv.org/html/2505.13954v2#A7.E62), we have:

$$ $\displaystyle\mathbb{E}\left[\|{\mathbf{v}}_{k}^{s}\|_{2}^{2}\right]$ $\displaystyle\leq 4\left(\left(2+\frac{24d\delta_{n}}{qb}\right)\alpha^{2}-2\alpha+1\right)\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]$ (64) $\displaystyle+\frac{12L^{2}\delta_{n}}{b}\left(1+\frac{4d}{q}\right)\alpha^{2}\mathbb{E}\left[\|\mathbf{x}_{0}^{s}-\mathbf{x}_{k}^{s}\|_{2}^{2}\right]+\frac{3\delta_{n}}{b}\left(1+\frac{2}{q}\right)L^{2}d^{2}\mu^{2}\alpha^{2}$ $\displaystyle+\frac{4\sigma^{2}}{b}\left(\frac{24d\delta_{n}\alpha^{2}}{q}+(1-\alpha)^{2}\right).$ $$

Substituting equation [64](https://arxiv.org/html/2505.13954v2#A7.E64) into equation [38](https://arxiv.org/html/2505.13954v2#A5.E38), we have:

$$ $\displaystyle R_{k+1}^{s}\leq$ $\displaystyle\mathbb{E}\left[f(\mathbf{x}_{k}^{s})\right]\bm{-}\left(\eta_{k}-\frac{c_{k\bm{+}1}\eta_{k}}{\beta_{k}}\right)\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]+\left(c_{k\bm{+}1}+c_{k\bm{+}1}\beta_{k}\eta_{k}\right)\mathbb{E}\left[\|\mathbf{x}_{k}^{s}-\mathbf{x}_{0}^{s}\|_{2}^{2}\right]$ (65) $\displaystyle+\left(\frac{L}{2}\eta_{k}^{2}+c_{k+1}\eta_{k}^{2}\right)\frac{12L^{2}\delta_{n}}{b}\left(1+\frac{4d}{q}\right)\alpha^{2}\mathbb{E}\left[\|\mathbf{x}_{k}^{s}-\mathbf{x}_{0}^{s}\|_{2}^{2}\right]$ $\displaystyle+4\left(\frac{L}{2}\eta_{k}^{2}+c_{k+1}\eta_{k}^{2}\right)\left(\left(2+\frac{24d\delta_{n}}{qb}\right)\alpha^{2}-2\alpha+1\right)\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]$ $\displaystyle+\left(\frac{L}{2}\eta_{k}^{2}+c_{k+1}\eta_{k}^{2}\right)\frac{3\delta_{n}}{b}\left(1+\frac{2}{q}\right)L^{2}d^{2}\mu^{2}\alpha^{2}$ $\displaystyle+\left(\frac{L}{2}\eta_{k}^{2}+c_{k+1}\eta_{k}^{2}\right)\frac{4\sigma^{2}}{b}\left(\frac{24d\delta_{n}\alpha^{2}}{q}+(1-\alpha)^{2}\right)$ $$

Based on the definition of $c_{k}=\left(1+\beta_{k}\eta_{k}\right)c_{k\bm{+}1}+\left(\frac{L}{2}+c_{k+1}\right)\left(1+\frac{4d}{q}\right)\frac{12L^{2}\delta_{n}\eta_{k}^{2}\alpha^{2}}{b}$ and $R_{k}^{s}$ given by equation [37](https://arxiv.org/html/2505.13954v2#A5.E37), we can simplify equation [65](https://arxiv.org/html/2505.13954v2#A7.E65) to

$$ $\displaystyle R_{k+1}^{s}$ $\displaystyle\leq R_{k}^{s}\bm{-}\left(\eta_{k}-\frac{c_{k\bm{+}1}\eta_{k}}{\beta_{k}}\right)\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]$ (66) $\displaystyle+4\left(\frac{L}{2}\eta_{k}^{2}+c_{k+1}\eta_{k}^{2}\right)\left(\left(2+\frac{24d\delta_{n}}{qb}\right)\alpha^{2}-2\alpha+1\right)\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]$ $\displaystyle+\left(\frac{L}{2}\eta_{k}^{2}+c_{k+1}\eta_{k}^{2}\right)\frac{3\delta_{n}}{b}\left(1+\frac{2}{q}\right)L^{2}d^{2}\mu^{2}\alpha^{2}$ $\displaystyle+\left(\frac{L}{2}\eta_{k}^{2}+c_{k+1}\eta_{k}^{2}\right)\frac{4\sigma^{2}}{b}\left(\frac{24d\delta_{n}\alpha^{2}}{q}+(1-\alpha)^{2}\right)$ $\displaystyle\leq R_{k}^{s}-\gamma_{k}\mathbb{E}\left[\|\nabla f(\mathbf{x}_{k}^{s})\|_{2}^{2}\right]+\chi_{k},$ $$

where $\gamma_{k}$ and $\chi_{k}$ are defined coefficients in Theorem [2](https://arxiv.org/html/2505.13954v2#Thmtheorem2).
Based on equation [66](https://arxiv.org/html/2505.13954v2#A7.E66) and the following argument in, we can achieve

$$ $\displaystyle\mathbb{E}[\|\nabla f(\bar{\mathbf{x}})\|_{2}^{2}]\leq\frac{\mathbb{E}[f(\bar{\mathbf{x}}_{0})-f^{*}]}{T\bar{\gamma}}+\frac{S\chi_{m}}{T\bar{\gamma}}.$ (67) $$

The rest of the proof is similar to the proof of Corollary [1](https://arxiv.org/html/2505.13954v2#Thmcorollary1) with the added complexity of the parameter $q$.
Let $\theta=\beta\eta_{k}+\left(1+\frac{4d}{q}\right)\frac{12L^{2}\delta_{n}\alpha^{2}}{b}\eta_{k}^{2}$, and $c_{k}=c_{k+1}(1+\theta)+\left(1+\frac{4d}{q}\right)\frac{6L^{3}\delta_{n}\eta_{k}^{2}\alpha^{2}}{b}$. This leads to:

$$ $\displaystyle c_{0}=\left(1+\frac{4d}{q}\right)\frac{6L^{3}\delta_{n}\eta^{2}\alpha^{2}}{b}\frac{(1+\theta)^{m}-1}{\theta}$ (68) $$

Let $\eta=\frac{\rho}{L}$, $\alpha=\frac{q}{d}$, $\beta=L$, and $q\leq d$ we have:

$$ $\displaystyle\theta=\rho+\left(q+4d\right)\frac{12\delta_{n}q\rho^{2}}{bd^{2}}\leq\rho+2\rho\left(\frac{q^{2}}{d^{2}}+4\frac{q}{d}\right)\leq\rho+\frac{60\rho q}{d}$ (69) $$

Substituting equation [69](https://arxiv.org/html/2505.13954v2#A7.E69) into equation [68](https://arxiv.org/html/2505.13954v2#A7.E68), we have:

$$ $\displaystyle c_{k}\leq c_{0}$ $\displaystyle=\left(1+\frac{4d}{q}\right)\frac{6L^{3}\delta_{n}\eta^{2}\alpha^{2}}{b}\frac{(1+\theta)^{m}-1}{\theta}$ (70) $\displaystyle=\frac{6(q+4d)L\delta_{n}\rho q}{bd^{2}+12(q+4d)\delta_{n}\rho q}\left[(1+\theta)^{m}-1\right]$ $\displaystyle\leq\frac{6(q+4d)L\delta_{n}\rho q}{bd^{2}}\left[(1+\theta)^{m}-1\right]$ $\displaystyle\leq\frac{30L\delta_{n}\rho q}{bd}\left(e-1\right)=\frac{60L\delta_{n}\rho q}{bd},$ $$

where the second inequality holds since $q\leq d$, and the first inequality holds if $m=\lceil\frac{1}{\rho+\frac{108\rho q}{d}}\rceil$
Because we define $\bar{\gamma}=\min_{k}\gamma_{k}$, we have

$$ $\displaystyle\bar{\gamma}\geq\eta-\frac{c_{0}\eta}{\beta}-4\left(\frac{L}{2}\eta^{2}+c_{0}\eta^{2}\right)\left(\left(2+\frac{24d\delta_{n}}{qb}\right)\alpha^{2}-2\alpha+1\right)$ (71) $$

From equation [70](https://arxiv.org/html/2505.13954v2#A7.E70), we have,

$$ $\displaystyle\frac{c_{0}}{\beta}\leq\frac{60\rho q}{bd}$ (72) $$

Because $\eta=\frac{\rho}{L}$, $\alpha=\frac{q}{d}$, and $q\leq d$ we have

$$ $\displaystyle\left(\frac{L}{2}\eta+c_{0}\eta\right)\left(\left(2+\frac{24d\delta_{n}}{qb}\right)\alpha^{2}-2\alpha+1\right)$ (73) $\displaystyle\leq\left(\frac{\rho}{2}+\frac{60\rho^{2}q}{bd}\right)\left(\frac{2q^{2}}{d^{2}}+\frac{24q}{bd}-\frac{2q}{d}+1\right)$ $\displaystyle\leq\rho\left(\frac{24q}{bd}+1\right)\leq\rho\left(\frac{24}{b}+1\right)$ $$

The second inequality holds if we let $\rho\leq\frac{1}{120}$
Substituting equation [73](https://arxiv.org/html/2505.13954v2#A7.E73) and equation [72](https://arxiv.org/html/2505.13954v2#A7.E72) into equation [71](https://arxiv.org/html/2505.13954v2#A7.E71), we can get

$$ $\displaystyle\bar{\gamma}\geq\eta\left(1-\frac{60\rho q}{bd}-4\rho\left(\frac{24}{b}+1\right)\right)=\eta z_{0},$ (74) $$

where $z_{0}>0$, and $\bar{\gamma}$ is a universal constant that is independent of $T$, $b$ and $d$.
Then, we bound $\chi_{m}=\sum_{k}\chi_{k}$

$$ $\displaystyle\chi_{m}$ $\displaystyle\leq m\left(\frac{L}{2}\eta^{2}+c_{0}\eta^{2}\right)\frac{3\delta_{n}}{b}\left(1+\frac{2}{q}\right)L^{2}d^{2}\mu^{2}\alpha^{2}$ (75) $\displaystyle+m\left(\frac{L}{2}\eta^{2}+c_{0}\eta^{2}\right)\frac{4\sigma^{2}}{b}\left(\frac{24d\delta_{n}\alpha^{2}}{q}+(1-\alpha)^{2}\right)$ $$

Because $c_{0}\leq\frac{60L\rho q}{bd}\leq\frac{L}{2}$ if $\rho\leq\frac{1}{120}$, this yields

$$ $\displaystyle\frac{\chi_{m}}{\bar{\gamma}}$ $\displaystyle\leq\frac{\rho}{z_{0}}\frac{3\delta_{n}}{b}\left(1+\frac{2}{q}\right)L^{2}\mu^{2}q^{2}$ (76) $\displaystyle+\frac{\rho}{z_{0}}\frac{4\sigma^{2}}{b}\left(\frac{24q\delta_{n}}{d}+(1-\frac{q}{d})^{2}\right)$ $$

Since $T=Sm$ and $\mu=\frac{1}{q\sqrt{T}}$, we have

$$ $\displaystyle\frac{S\chi_{m}}{T\bar{\gamma}}$ $\displaystyle\leq\frac{\rho}{z_{0}}\frac{3\delta_{n}}{b}\left(1+\frac{2}{q}\right)\frac{L^{2}}{T}$ (77) $\displaystyle+\frac{\rho}{z_{0}}\frac{4\sigma^{2}}{b}\left(\frac{24q\delta_{n}}{d}+(1-\frac{q}{d})^{2}\right)$ $\displaystyle\leq O\left(\frac{1}{bT}+\frac{1}{b}\left(1-\frac{q}{d}\right)^{2}\right)$ $$

Substituting equation [74](https://arxiv.org/html/2505.13954v2#A7.E74) and equation [77](https://arxiv.org/html/2505.13954v2#A7.E77) into equation [7](https://arxiv.org/html/2505.13954v2#S4.E7), we have

$$ $\displaystyle\mathbb{E}[\|\nabla f(\bar{\mathbf{x}})\|_{2}^{2}]\leq\frac{\mathbb{E}[f(\bar{\mathbf{x}}_{0})-f^{*}]}{Tz_{0}}\frac{L}{\rho}+\frac{S\chi_{m}}{T\bar{\gamma}}=O\left(\frac{1}{T}+\frac{1}{bT}+\frac{1}{b}\left(1-\frac{q}{d}\right)^{2}\right)$ (78) $$

∎

###### Proof.

### G.1 Auxiliary Lemmas

###### Lemma 3 .

Let $\{\mathbf{z}_{i}\}^{n}_{i=1}$ be a sequence of n vectors. Let $\mathcal{I}$ be a mini-batch of size b, which contains i.i.d. samples selected uniformly randomly (with replacement) from $[n]$.

$$ $\displaystyle\mathbb{E}_{\mathcal{I}}\left[\frac{1}{b}\sum_{i\in\mathcal{I}}\mathbf{z}_{i}\right]=\frac{1}{n}\sum_{j=1}^{n}\mathbf{z}_{j}.$ (79) $$

When $\sum_{i=1}^{n}\mathbf{z}_{i}=\mathbf{0}$, then

$$ $\displaystyle\mathbb{E}_{\mathcal{I}}\left[\left\|\frac{1}{b}\sum_{i\in\mathcal{I}}\mathbf{z}_{i}\right\|_{2}^{2}\right]=\frac{1}{bn}\sum_{i=1}^{n}\|\mathbf{z}_{i}\|_{2}^{2}.$ (80) $$

###### Proof.

See the proof of Lemma 4 in (Liu et al., 2018b).
∎

###### Lemma 4 .

Let $\{\mathbf{z}_{i}\}^{n}_{i=1}$ be a sequence of n vectors. Let $\mathcal{I}$ be a uniform random mini-batch of $[n]$ with size b (no replacement in samples). Then

$$ $\displaystyle\mathbb{E}_{\mathcal{I}}\left[\frac{1}{b}\sum_{i\in\mathcal{I}}\mathbf{z}_{i}\right]=\frac{1}{n}\sum_{j=1}^{n}\mathbf{z}_{j}.$ (81) $$

When $\sum_{i=1}^{n}\mathbf{z}_{i}=\mathbf{0}$, then

$$ $\displaystyle\mathbb{E}_{\mathcal{I}}\left[\left\|\frac{1}{b}\sum_{i\in\mathcal{I}}\mathbf{z}_{i}\right\|_{2}^{2}\right]=\frac{\mathcal{I}(b<n)}{bn}\sum_{i=1}^{n}\|\mathbf{z}_{i}\|_{2}^{2}.$ (82) $$

where $I$ is an indicator function, which is equal to 1 if $b<n$ and 0 if $b=n$.

###### Proof.

See the proof of Lemma A.1 in (Lei et al., 2017).
∎

###### Lemma 5 .

For variables $\{\mathbf{z}_{i}\}^{n}_{i=1}$, we have

$$ $\displaystyle\left\|\sum_{i=1}^{n}\mathbf{z}_{i}\right\|_{2}^{2}\leq n\sum_{i=1}^{n}\|\mathbf{z}_{i}\|_{2}^{2}.$ (83) $$

###### Proof.

See the proof of Lemma 6 in (Liu et al., 2018b).
∎

###### Lemma 6 .

if f is L-smooth, then for any $\mathbf{x},\mathbf{y}\in\mathbb{R}^{d}$

$$ $\displaystyle|f(\mathbf{x})-f(\mathbf{y})-\langle\nabla f_{i}(\mathbf{y}),\mathbf{x}-\mathbf{y}\rangle|\leq\frac{L}{2}\|\mathbf{x}-\mathbf{y}\|_{2}^{2}.$ (84) $$

###### Proof.

This is a direct consequence of Lemma A.2 in (Lei et al., 2017).
∎

###### Lemma 3 .

###### Proof.

###### Lemma 4 .

###### Proof.

###### Lemma 5 .

###### Proof.

###### Lemma 6 .

###### Proof.

## Appendix H Analysis of Zeroth-Order Gradient Estimation Error

This section details bounds on the expected squared error of the ZO gradient estimators used in our work. We consider a ZO gradient estimator $\hat{\nabla}f_{i}(\mathbf{x})$ for a component function $f_{i}(\mathbf{x})$, which approximates the true gradient $\nabla f_{i}(\mathbf{x})$ with an estimation error $\omega_{i}(\mathbf{x})$, such that $\hat{\nabla}f_{i}(\mathbf{x})=\nabla f_{i}(\mathbf{x})+\omega_{i}(\mathbf{x})$. The characteristics of the expected squared error, $\mathbb{E}[\|\omega_{i}(\mathbf{x})\|_{2}^{2}]$, are presented below.

For the two-point ZO gradient estimator of $f_{i}(\mathbf{x})$, as defined in Equation equation [2](https://arxiv.org/html/2505.13954v2#S3.E2) in the main text, the expected squared error is bounded by:

$$ $\mathbb{E}[\|\omega_{i}(\mathbf{x})\|_{2}^{2}]\leq O(d)\|\nabla f_{i}(\mathbf{x})\|_{2}^{2}+O(\mu^{2}L^{2}d^{2}).$ (85) $$

Here, $d$ is the problem dimension, $\mu$ is the smoothing parameter, and $L$ is the smoothness constant associated with $f_{i}$.

Subsequently, for the multi-point ZO gradient estimator of $f_{i}(\mathbf{x})$ using $2q$ query points, as defined in Equation equation [3](https://arxiv.org/html/2505.13954v2#S3.E3) in the main text, the expected squared error is bounded by:

$$ $\mathbb{E}[\|\omega_{i}(\mathbf{x})\|_{2}^{2}]\leq O(d/q)\|\nabla f_{i}(\mathbf{x})\|_{2}^{2}+O(\mu^{2}L^{2}d^{2}).$ (86) $$

The detailed proofs for these bounds can be found in Proposition 2 of (Liu et al., 2018b).

## Appendix I Experiment Setup

### I.1 Adaptability Experiment

The primary objective of this adaptability experiment was to empirically investigate the impact of the number of ZO query points ($q$) on the performance of VAMO and to validate the theoretical benefits of its multi-point ZO estimation strategy (see Section [5](https://arxiv.org/html/2505.13954v2#S5)). The optimization problem was a finite-sum non-convex least-squares objective: $f(\mathbf{x})=\frac{1}{n}\sum_{i=1}^{n}(h(\mathbf{x};\mathbf{z}_{i})-y_{i})^{2}$. We configured this synthetic task with $n=1000$ individual component functions and a parameter dimension of $d=100$. The function $h(\mathbf{x};\cdot)$ was parameterized using a simple neural network with a non-convex activation function to ensure the overall non-convexity of the loss landscape.

In this setup, VAMO variants utilizing $q\in\{1,3,5\}$ query directions for the multi-point ZO gradient estimator were compared against the classical FO-SGD algorithm. A mini-batch size of $b=8$ was consistently applied across all methods. Learning rates for both VAMO (for each $q$ setting) and FO-SGD were individually tuned by selecting the best performing value from the range $[10^{-2},10^{-1}]$. For all VAMO variants, the ZO smoothing parameter was fixed at $\mu=10^{-3}$. The mixing coefficient $\alpha$ for VAMO was also tuned for each value of $q$, guided by the theoretical insights on balancing FO and ZO information discussed in Appendix [J.2](https://arxiv.org/html/2505.13954v2#A10.SS2).

### I.2 MNIST Classification Task

For the MNIST multi-class image classification task (LeCun et al., 1998), we trained a Multi-Layer Perceptron (MLP) to evaluate VAMO against established baselines. The MLP architecture consisted of an input layer receiving flattened $28\times 28$ pixel images (784 dimensions), followed by two hidden layers with 32 and 16 units respectively, both employing ReLU activation functions. The final output layer comprised 10 units corresponding to the digit classes, and the network was trained using a standard cross-entropy loss function. Images were normalized to the range $[0,1]$.

Our proposed VAMO algorithm, configured with a single ZO query direction ($q=1$), was benchmarked against pure first-order FO-SGD (Robbins & Monro, 1951) and pure zeroth-order methods, ZO-SGD and ZO-SVRG  (Liu et al., 2018b; Ghadimi & Lan, 2013). For all methods, the mini-batch size was set to $b=4$. The learning rates were independently tuned for each method, selected from the range $[10^{-4},10^{-3}]$ for the ZO methods, and $[10^{-3},10^{-2}]$ for FO methods and VAMO. For VAMO with $q=1$, we fixed the mixing coefficient at $\alpha=0.1$ and used a ZO smoothing parameter of $\mu=10^{-3}$. All models were trained for 10 epochs.

### I.3 Fine-Tuning Experiments

To further examine VAMO in realistic large-scale scenarios, we conducted fine-tuning experiments on three language models under the MultiNLI (MNLI) dataset (Williams et al., 2017) for a three-way natural language inference task. In all settings, the training and validation sets were subsampled to 256 and 128 examples, respectively, with a maximum input sequence length of 128 tokens. All models were fine-tuned for 1500 epochs with full precision (FP32) on a single NVIDIA RTX5880 49GB GPU. All experiments involved full-parameter fine-tuning.

For fair comparisons, VAMO was compared against representative FO methods including FO-SGD (Robbins & Monro, 1951), FO-Adagrad (Duchi et al., 2011) and FO-Adam (Kingma & Ba, 2014), and ZO methods including ZO-SGD (Ghadimi & Lan, 2013) and ZO-SVRG (Liu et al., 2018b). We set the batch size as 32. The learning rates of FO methods were tuned within $[10^{-4},10^{-3}]$, while those of ZO methods were tuned within $[10^{-6},10^{-5}]$. For VAMO, the learning rate was also tuned in the range $[10^{-4},10^{-3}]$. We adopted the same learning rate schedule as illustrated in Algorithm [2](https://arxiv.org/html/2505.13954v2#alg2). For ZO methods, we adopted a smoothing parameter $\mu=10^{-3}$ and $q=1$ (two-point version). VAMO was configured with the the smoothing parameter $\mu=10^{-3}$, $q=1$, and inner loop length $m=10$.

The main difference is in setting the mixing coefficient $\alpha$. $\alpha$ was set to $0.01$ for the GPT-2 experiments, while a lower value $\alpha=0.005$ was used for RoBERTa-Large and GPT-2 Medium. We used a smaller $\alpha$ for larger models because the ZO gradient estimates exhibit higher inherent error in higher-dimensional models, so reducing $\alpha$ limits the contribution of the noisier ZO correction and helps stabilize the training. A detailed analysis of $\alpha$ selection and the error of the ZO estimation is provided in Appendix [J.2](https://arxiv.org/html/2505.13954v2#A10.SS2) and Appendix [H](https://arxiv.org/html/2505.13954v2#A8).

Figure: Algorithm 2 Learning Rate Scheduling for VAMO

## Appendix J Ablation on Key Hyperparameters

In this section, we investigate the effect of two critical hyperparameters in VAMO: the number of inner loop iterations $m$ and the mixing coefficient $\alpha$. Both parameters play a central role in balancing the trade-offs of the algorithm—$m$ controls the frequency of full-batch ZO gradient corrections, while $\alpha$ adjusts the relative contribution of ZO-based variance reduction. To better understand their influence, we conduct ablation studies on fine-tuning GPT-2 Medium (Radford et al., 2019) under the MNLI dataset (Williams et al., 2017), examining both convergence speed and stability.

### J.1 Role of m m (Inner Loop Iterations)

The parameter $m$ is a significant role in VAMO as it governs the frequency of the fullbatch ZO gradient updates. Specifically, a smaller $m$ means that the full-batch ZO gradient estimate is computed more frequently, which can lead to a more effective reduction in gradient variance. To better understand the trade-off between the effectiveness of these updates and their computational overhead, we perform an ablation study on $m$.

We consider the task of fine-tuning a GPT-2 model on the MNLI dataset . For this study, we fixed the other key hyperparameters ($\alpha=0.1$, $q=1$, and $\mu=10^{-3}$) while varying $m$. The results, illustrated in Figure [3](https://arxiv.org/html/2505.13954v2#A10.F3), show that increasing the frequency of fullbatch ZO gradient update ($2\leq m\leq 10$) generally enhances the performance of VAMO. However, this trend does not continue indefinitely; when a full batch update becomes too infrequent (e.g., $m\geq 20$), the variance correction becomes stale and less effective, leading to a performance degradation. This suggests an optimal range for $m$ that balances the benefits of frequent variance correction against its computational cost.

### J.2 Role of α \alpha (Mixing Coefficient)

#### Theoretical Analysis.

The mixing coefficient $\alpha$ in the VAMO update in  equation [6](https://arxiv.org/html/2505.13954v2#S4.E6) critically balances the FO stochastic gradient $\nabla f_{\mathcal{I}}(\mathbf{x})$ against the ZO variance correction term $\hat{\nabla}f_{\mathcal{I}}(\hat{\mathbf{x}})-\hat{\nabla}f(\hat{\mathbf{x}})$. The optimal choice for $\alpha$ directly depends on the estimation error $\omega_{i}$ inherent in the ZO gradient components ($\hat{\nabla}f_{i}(\mathbf{x})=\nabla f_{i}(\mathbf{x})+\omega_{i}$). As established in the literature (Liu et al., 2018b; 2020) and detailed in Appendix [H](https://arxiv.org/html/2505.13954v2#A8), the expected squared ZO error $\mathbb{E}[\|\omega_{i}\|_{2}^{2}]$ typically scales as $\mathcal{O}(d)$ for two-point estimates and $\mathcal{O}(d/q)$ for multi-point estimates using $q$ random directions. This relationship dictates that $\alpha$ should reflect the trustworthiness of the ZO estimates: when the ZO error is substantial (e.g., large $d$, small $q$), a smaller $\alpha$ is warranted to prevent amplifying this error. Conversely, when ZO estimates are more reliable (e.g., larger $q$ reducing error), a larger $\alpha$ can more aggressively leverage the variance reduction. This principled inverse relationship between ZO error magnitude (influenced by $d$ and $q$) and the appropriate scale of $\alpha$ is key. While specific forms like $\alpha\propto 1/d$ or $\alpha\propto q/d$ analyzed in our theoretical sections (e.g., Corollary [1](https://arxiv.org/html/2505.13954v2#Thmcorollary1) and Theorem [2](https://arxiv.org/html/2505.13954v2#Thmtheorem2)) illustrate this adaptive trend, the core insight is that $\alpha$ must be adjusted to counterbalance the ZO estimator’s error profile. Such adaptability enables VAMO to effectively navigate the trade-off between computational cost and convergence performance, a central aspect of its practical utility.

#### Empirical Results Analysis.

The mixing coefficient $\alpha$ controls how much weight VAMO assigns to the ZO-based variance reduction term relative to the FO gradient. A larger $\alpha$ strengthens the effect of variance correction but also amplifies the noise of the ZO estimator. To study this trade-off, we fine-tuned GPT-2 on the MNLI dataset, fixing other hyperparameters ($m=10$, $q=1$, $\mu=10^{-3}$) while varying $\alpha$. As shown in Figure [3](https://arxiv.org/html/2505.13954v2#A10.F3), smaller values (e.g. $\alpha=0.01$) provide a good balance, leading to fast convergence without being overly sensitive to noisy ZO gradients. However, if $\alpha$ is set too small (e.g., $\alpha=0.005$), the variance reduction term becomes underutilized, diminishing its benefit. In contrast, setting $\alpha$ too high (e.g., $\alpha=0.2$) magnifies estimation error, destabilizes training, and slows convergence. These results align with our theoretical analysis, which highlights the need for careful tuning of $\alpha$ in high-dimensional settings.

Figure: (a)
Refer to caption: x9.png

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 4
presentation: 3
questions: see weakenesses.
soundness: 3
strengths: 1. I like the discussion of memory: the paper argues VAMO’s snapshot uses forward-only ZO passes, so peak dynamic memory resembles ZO-SGD rather than FO-SVRG, and optimizer state is lighter than Adam/Adagrad. Included tables and a clear decomposition (weights / states / dynamics) are helpful.
2. The main bound removes the typical d dependence of ZO methods and improves over FO-SGD’s O(1/\sqrt{T}) in theory.
summary: The paper proposes VAMO, a hybrid variance-reduced optimizer that blends a mini-batch first-order (FO) gradient with a zeroth-order (ZO) SVRG-style correction computed at snapshot points. The key idea is to replace the expensive full-batch FO snapshot gradient with a full-batch ZO estimate and to weight the correction with a mixing coefficient \alpha. The authors prove a  convergence rate of O(1/T+1/b) independent of dimension d. Experiments span a synthetic task, MNIST MLP, and fine-tuning GPT-2 / GPT-2-Medium / RoBERTa-Large on SST-2 and MNLI, with GPU memory comparisons versus FO-SGD/Adagrad/Adam.
weaknesses: 1.ZO methods inherently trade off performance, memory, and wall-clock. The paper treats ZO snapshots as “cheap,” but a full-batch snapshot still requires multiple forward passes per direction over the entire dataset (or many mini-batches). In practice, if memory is the bottleneck then ZO can help. However, runtime can balloon unless the number of directions q is tiny—yet shrinking q raises estimator variance and hurts convergence. Therefore, without an accounting of function evaluations (FEs) and throughput (e.g., tokens/sec), it’s unclear when VAMO is actually preferable to tuned FO baselines (Adam/Lion/Adafactor) or to ZO-SVRG variants with different q.

2. ZO can also be performed when the random direction is Gaussian, i.e., let u(x;\theta) = E_{\delta\sim N(0,\sigma^2 I) f(x+\delta;\theta), then \nabla_x u(x;\theta) = E_{\delta\sim N(0,\sigma^2 I) [\delta/\sigma^2  f(x+\delta;\theta). In high dimension, Gaussian vs. coordinate-wise directions can yield different estimator variance and smoothing bias. Yet this is not compared in the paper.

3. Fine-tuning results present training loss (showing convergence) and memory costs. However, I would expect also the task metrics (accuracy/F1 for SST-2/MNLI) and stability stats (divergence/NaNs) to show the effectiveness of VAMO.

4. SVRG-style methods hinge on the outer-loop frequency S and inner length m. The paper would benefit from an empirical study of how often to recompute ZO snapshots (and with what q) under a fixed compute budget; otherwise, it’s hard to see when VAMO is preferable to carefully tuned FO-SGD/Adam or existing ZO-SVRG variants.

---

confidence: 3
contribution: 2
initial_rating: 6
presentation: 3
questions: See weakness.
soundness: 2
strengths: 1. VAMO achieves a fast, linear convergence rate of $\mathcal{O}(1/T)$, which is an improvement over the $\mathcal{O}(1/\sqrt{T})$ rate of standard SGD.
2. A key advantage is that its convergence rate is independent of the model's parameter dimension $d$. This allows it to overcome the "curse of dimensionality" that makes purely Zeroth-Order (ZO) methods impractical for large models.
3. This paper provides a strong theoretical guarantee. VAMO's gradient estimator is designed to be unbiased.
summary: The research presents VAMO as a hybrid optimizer that works with big models under memory restrictions. The optimizer uses First-Order (FO) speed together with Zeroth-Order (ZO) efficiency by substituting the expensive SVRG variance-reduction algorithm step with a low-cost ZO gradient estimate. The algorithm provides fast dimension-independent convergence speed that outperforms SGD yet requires memory levels similar to SGD and less than Adam. The empirical results show VAMO runs at a slower pace than Adam but provides an attractive trade-off between performance and memory usage for restricted resource scenarios.
weaknesses: 1. The experimental results demonstrate that VAMO achieves better performance than SGD, but it fails to match the convergence speed and training efficiency of the Adam optimizer during large-scale fine-tuning.
2. The theoretical analysis shows that VAMO's convergence rate, while faster than SGD's, includes an additional error term of $\mathcal{O}(1/b)$ that is not present in the purely First-Order FO-SVRG algorithm.
3. The paper introduces a multi-point variant to minimize the additional error term. The improvement requires additional computational resources, which increase with the number of query points $q$, thus users need to decide between faster convergence and higher processing expenses.

---

confidence: 4
contribution: 2
initial_rating: 2
presentation: 3
questions: 1. The text in Section 4.3 (Lines 325-328) claims that VAMO does not need to store intermediate
activations. However, Step 7 of Algorithm 1 computes an FO gradient ∇fIk (xsk), which necessitates storing activations. Your own Table 3 and Table 2 confirm that VAMO has the same dynamic memory as FO-SGD. Can you please clarify this fundamental contradiction? Is the premise of the paper’s memory-saving benefits not incorrect?

2. The O(1/T + 1/b) rate converges to a noise ball, not a stationary point. Why is this
misleadingly presented as "similar to FO-SVRG" (which converges to 0) and an improvement
over FO-SGD (which also converges to 0)? Figure 1b seems to confirm that VAMO(q = 1) converges
to a worse solution.

3. The inner loop (Step 7) computes both an FO-grad (cost O(bd)) and a ZO-grad (cost O(bq)).
Why is the O(bqT) term missing from the complexity analysis in Table 1?

4. Given that VAMO(q = d) has the same rate and computation as FO-SVRG, why was this
Comparison (which is the only fair comparison of two O(1/T ) methods) is not included in the
experiments?

5. How sensitive are convergence and stability to the choice of α, μ, and q? Could you provide
empirical ablations or adaptive scheduling strategies, validating these settings across different
model dimensions?
soundness: 2
strengths: Novel Hybrid Estimator: The proposed gradient estimator (Eq. 6) is novel. The insight to
use a ZO-based correction term α(ˆ∇fIk (ˆx) − ˆ∇f(ˆx)), which has zero expectation, is technically
sound. This ensures the full estimator vsk is an unbiased estimator of the true current gradient
∇f(xsk), which is an elegant property for the analysis.

2. Dimension-Independent ZO-Hybrid Rate: The theoretical analysis successfully breaks the
dimension-dependency curse of pure ZO methods. Achieving a rate of O(1/T + 1/b) (Corollary
1) that is independent of d is a strong theoretical contribution for a method that incorporates
ZO components.

3. Strong Empirical Performance vs. FO-SGD: The experiments, particularly in Figure 2
(LLM finetuning), show that VAMO (with q = 1) converges significantly faster in terms of both
steps and wall-clock time than FO-SGD. This empirically validates the theoretical advantage of
the O(1/T + ...) rate over FO-SGD’s O(1/√T) rate for reaching a certain precision.
summary: This paper proposes VAMO, a new stochastic optimization algorithm for large-scale non-convex optimization.
The algorithm aims to bridge the gap between first-order (FO) methods, which have fast
convergence but high computational and memory costs, and zeroth-order (ZO) methods, which are
memory-light but suffer from dimension-dependent and slow convergence.
VAMO is a two-loop "SVRG-style" algorithm. Its central idea is to replace the computationally
expensive full-batch FO gradient (∇f(ˆx)) used in FO-SVRG with a full-batch ZO gradient estimate
(ˆ∇f(ˆx)). The inner-loop update is a novel hybrid estimator: vsk
= ∇fIk (xsk) − α(ˆ∇fIk (ˆx) − ˆ∇f(ˆx)).
This construction cleverly uses a ZO-based variance correction term that has zero expectation, making
vsk an unbiased estimator of the true gradient ∇f(xsk).
The authors provide a theoretical analysis showing the two-point (q = 1) version of VAMO achieves a
convergence rate of O(1/T +1/b), which is dimension-independent, a significant improvement over typical
ZO methods. They also propose a multi-point (q > 1) variant with a rate of O(1/T +(1−q/d)2/b),
which can converge to a stationary point if q = d. The primary claims are that VAMO achieves FOSVRG-
like convergence speed with significantly lower computational and, critically, dynamic memory
costs. Experiments on neural networks and LLM finetuning are presented to support these claims.
weaknesses: Despite its theoretical novelty, the paper’s core claims about its practical advantages are based on a
series of critical, and in some cases contradictory, flaws in the analysis of its computational and memory
costs.

1. Fundamentally Misleading Convergence Claims: The paper repeatedly claims an O(1/T )
rate, equating it with FO-SVRG (e.g., Abstract: "significantly improving over SGD’s... rate";
Conclusion: "achieving convergence performance similar to FO-SVRG"). This is a misrepresentation.
The actual derived rate is O(1/T + 1/b). This is not convergence to a stationary point.
It is linear convergence to a noise ball of size O(1/b). This non-vanishing error term, which
dominates at large T, means VAMO (q = 1) cannot achieve high-precision solutions and may, in
fact, converge to a worse solution than FO-SGD (which does converge to ϵ = 0, albeit slower).
This is a crucial distinction that is glossed over.

2. Contradictory and Factually Incorrect Memory Analysis: This is the paper’s most severe
flaw. The paper is motivated by reducing the high dynamic memory of FO methods (e.g., storing
activations for backpropagation). The inner-loop update (Algorithm 1, Step 7) explicitly requires
a mini-batch FO gradient, ∇fIk (xsk). Computing this FO gradient requires backpropagation
and storing intermediate activations, leading to a dynamic memory cost of O(b ·P|al|). This is
the exact same dynamic memory cost as FO-SGD. The paper’s text (e.g., Section 4.3, Lines 325-
328: "thus do not need to store intermediate results", "dynamic memory only reaches maxl |xl|")
is factually incorrect and fundamentally misunderstands the cost of its own algorithm. The
data in the paper confirm this. Table 3 correctly lists VAMO’s dynamic memory as
Pl max{b ·|al|, |xl|}, identical to FO-SGD. Table 2 empirically shows VAMO’s memory (e.g., 21.46 GB) is
almost identical to FO-SGD’s (20.33 GB). The small increase is expected, as VAMO = FO-SGD
(dynamic) + |x| (snapshot state). The central claim that VAMO is a low-memory algorithm
(relative to FO-SGD) is false.

3. No Asymptotic Computational Advantage for O(1/T ) Convergence: To achieve the true
O(1/T ) convergence (i.e., remove the O(1/b) error), one must use the multi-point variant and set
q = d (Theorem 2).The computational complexity of VAMO (q = d) is O(qnS + (bd + bq)T) =
O(dnS + (bd + bd)T) = O(dnS + bdT ). Therefore, to achieve the same convergence rate as FOSVRG,
VAMO requires the identical asymptotic computational complexity. The paper’s claim
of computational efficiency is only valid when comparing the non-converging q = 1 version to the
converging FO-SVRG, which is an apples-to-oranges comparison.

4. Incorrect Computational Complexity in Table 1: The complexity for VAMO (q = 1) is
listed as O(nS +bdT ). This is incorrect. The inner loop (Step 7) computes two gradients: ∇fIk
(cost O(bd)) and ˆ∇fIk (cost O(bq), or O(b) for q = 1). The correct complexity is O(qnS +(bd+
bq)T). For q = 1, this is O(nS + (bd + b)T), which is O(nS + bdT ) only if d ≫ 1. This O(bqT)
term is missing from the table.

5. Empirical Evidence of Non-Convergence: In Figure 1b, VAMO(q = 1) clearly converges to
a final training loss that is visibly *higher* than that of FO-SGD. This plot empirically confirms
the O(1/b) noise ball limitation. This is a poor result, as it fails to match the solution quality of
the baseline FO-SGD, yet this is not discussed.

---

confidence: 3
contribution: 2
initial_rating: 4
presentation: 2
questions: - See weaknesses.
soundness: 2
strengths: - Combining FO and ZO in SVRG is an interesting idea and the FO order component improves the performance of the algorithm compared with other ZO methods.
- The presentation of the paper is easy to follow with detailed comparison with prior works.
- Proofs seem all good.
summary: The paper proposes VAMO - a variance reduction optimization method that combines both zero-order gradient estimates and gradient to produce the update direction. Concretely, the algorithm uses the full batch and batch gradient estimates at the checkpoint (obtained from the standard zero-order method) and combines with the batch gradient at the current iterate ($g = \nabla f(x, B) - \alpha (\hat{\nabla} f (x_{\text{cpt}}, B) - \hat{\nabla}f(x_{\text{cpt}}))$). This gives a trade off between convergence and computational complexity compared with ZO-SVRG. In terms of convergence, the rate of VAMO is $1/T + 1/b$ where $b$ is the batch size. This has the extra  $1/b$ term is inherited from ZO-SVRG due to the bias in the estimator, but the first term is improved by a factor $d$. The computational complexity on the other hand is increased from $nS+bT$ to $nS+dbT$ (due to the gradient computation). An extension to this method is to use multi-point query, which can trade off between the improvement in the additional error term and the computational complexity. Experiments show that VAMO has faster convergence than other ZO methods and better memory footprint than FO methods.
weaknesses: I have several concerns below.
- I'm not sure I understand the memory analysis in Section 4.3, B.1 and Table 3. 
  - For VAMO, why is the memory for the optimizer states only $|x|$? The algorithm needs to store $\hat{\nabla}(x_{\text{cpt}})$ but also $x_{\text{cpt}}$ to compute the estimate $\hat{\nabla}(x_{\text{cpt}}, B)$ (see also line 3-4 in Alg. 1). I don't understand how to reduce the memory to only $|x|$? If my understanding is true, I don't see an improvement in the memory for the proposed method.
  - For Adagrad, why is the optimizer states $2|x|$? Doesn't the algorithm only store the accumulation of the gradients per coordinate so the memory needed is only $|x|$? Same for Adam, shouldn't it be $2|x|$?
- The experiment appears quite weak. First of all, the paper only reports the training loss. While optimization algorithms only optimize the training loss, we also care about the test loss and accuracy. Second, one main motivation the paper mentioned is to overcome limitations of FO and ZO methods for efficient training of large models. However, the set of experiments is quite limited. MNIST is an outdated dataset as an optimization benchmark. Usually, for optimization papers, starting with CIFAR-10/100 will show clearer impacts. I also suggest that the authors use a similar experiment setup in the paper MeZO and add more experiments.