We thank the reviewer qodG for their thoughtful and constructive evaluation of our work. We address the concerns raised below.

W1 ($\tau$ may grow with $d$)

Thanks for this remark on the interplay between $\tau$ and $d$. First, we would like to highlight that the acceleration $\tau d$ $\to$ $\tau + d$ we achieve is significant regardless of the dependence between growth of $\tau$ and $d$. Second, we think these parameters are of independent nature. For example, mixing time of the chain may represent the complexity of the environment from which the data is sampled, while the problem dimension may be the number of parameters in the chosen model.

Another, more concrete example of an application, where the chain is unrelated to the problem dimension, is a token algorithm (see Section 1.1 from [1]). In short, it is a decentralized optimization algorithm, where a number of users each have a function $f_i(x)$: $R^d \to R$ stored locally, and the users can only communicate with their neighbors in some graph $G.$ Our goal is to minimize the average of these functions. To do so, we start a random walk over $G$, where the current vertex is said to receive the "token", representing the most up-to-date version of $x$. Each time the token moves, the receiving vertex updates it according to its $f_i$. This maps to our framework with random walk over $G$ becoming the Markov chain, and the expectation becoming averaging over all users. Thus $\tau$, the mixing time of the chain, depends solely on the topology of the network between users, while $d$ is the dimensionality of the local models.

[1] Even, M. (2023, July). Stochastic gradient descent under Markovian sampling schemes. In International Conference on Machine Learning (pp. 9412-9439). PMLR.

W2 (notations $J, M, l, B, t$)

That’s true, the presentation of MLMC technique was rather abrupt. We kept short definitions in Table 2, and the detailed description of these parameters will be added together with some intuition to the Section 2.1 in the final version, thanks! Here is the explanation draft.

To ease the understanding, let’s start with the case $l = 1, M = \infty, B = 1$ and enumerate baseline estimates as $g_1, g_2, \dots$.

Then, the MLMC estimate will be $g_1$ with prob. $1/2$; $g_1 + (g_3 - g_2)$ with prob. $1/4$; $g_1 + (g_4 + g_5 - g_2 - g_3)$ with prob. $1/8$ and so on. The parameter $M$ is the upper bound on the number of estimates used. The parameter $l$ transforms the base estimator into a sequence of $l$ base estimators, effectively stretching everything $l$ times. Finally, $B$ serves as a hyperparameter that can multiplicatively increase $l$.

The parameter $t$ is the step of finite difference scheme (see (1)-(3)). It is not tied to MLMC or other parts of the algorithm and the notation is classical.

We hope that after the revision this part is not going to confuse the reader.

W2 & Q1 (partitioned sequences usage)

You are right, the partitioned sequences are not used by the algorithm. They only appear in Section 2.1 as a motivation for $\tau + d$ behavior. The algorithm then uses a simpler procedure of choosing a random direction per sample. It also achieves a $\tau + d$ rate as we show in appendix (Lemma 5), but the proof is less intuitive. We are going to explicitly mention this in the final version, thanks.

W3 & Q3 (empirical comparison with first-order methods)

Zero-order (ZO) and first-order (FO) optimization methods are not directly comparable to each other, as they access different oracles. We conducted additional experiments assuming the following information-theoretic equivalence: one FO oracle call is equivalent to $d$ ZO oracle calls.

Specifically we compared our algorithm against Randomized Accelerated GD of [1] on the exact same quadratic problem as in our original experiments (see Section 3) with $\sigma_2^2 =$1e-3. Here is the average optimization error (initial error 10, all errors presented were multiplied by 1e5 for readability), over 1e4 runs, the number of ZO oracle calls is fixed at 2e3 (thus FO gets fewer iterations).

Zero-order case:

First-order case:

At $d=1$, there are no random directions to choose from, and the ZO algorithm is effectively equivalent to the FO one, except the true derivative is replaced with finite difference approximation, which explains why our results are worse.

At $d=5$, approximation errors no longer play a significant role, as FO algorithm gets 5 times less iterations, and our algorithm convincingly beats it.

At $d=50$, FO does not escape the linear regime within the allowed budget, while ZO performs dramatically better.

[1] Aleksandr Beznosikov, Sergey Samsonov, Marina Sheshukova, Alexander Gasnikov, Alexey Naumov, and Eric Moulines. First order methods with markovian noise: from acceleration to variational inequalities. Advances in Neural Information Processing Systems, 36, 2024.

S&W (lower bounds almost not discussed)

It was infeasible to keep everything in the main part, so we decided to move the discussion of lower bounds to Appendix F as additional content. In Section 3.2 we present main lower bounds and refer the reader to the appendix for details. If space permits, we will move the discussion into the main text.

Q2 (vanilla GD performance)

We are not aware of any existing analysis of vanilla SGD in a Markovian setting, but it is expected to have the same stochastic term if used with our gradient estimator, while the linear convergence would be worse. In two-point setting, for example, the expected rate is

However, the accelerated rate we achieve is important both theoretically to match deterministic lower bounds (see remark on line 295) and in practice, where accelerated methods are dominant (Momentum SGD, Nesterov SGD, Adam).

Many thanks to Reviewer Ryzo for their time, detailed review, and suggestions. We now turn to their specific concerns.

W (several results in the main text are just approximations based on intuition, and in particular I did not understand the rationale leading to the (d + tau) factors in the convergence rates)

While the main text emphasizes intuition to aid accessibility, all results — including those involving the $(d + \\tau)$ factors — are rigorously proved in the appendix with appropriate references in the main text. We intentionally placed technical details there to keep the main exposition focused. The $(d + \\tau)$ term is in fact central to our theoretical contribution, and no claims rely on unproven approximations. To address this particular misunderstanding, we note that the partitioned sequences we used to motivate the $(d + \tau)$ behavior are not utilized by the algorithm. The algorithm uses $g\_{rd}$, which is also introduced in the main part with some motivation, but a rigorous proof of its properties is available in Lemma 5 in Appendix (also see the Appendix C Notations and Definitions, as the notation is shortened).

Below, we clarify the notation and expand on the intuition behind this result to make the rationale more transparent. We hope this addresses the reviewer’s concern and helps in evaluating the contribution more fully.

Q0 (subsequence estimators dependence)

That is true that the resulting d estimators derived from $d$ intertwined subsequences are not independent of each other, thus trivially averaging them would not reduce the variance $d$-fold. However, as these $d$ subsequences each estimate a partial derivative [line 231] along its separate direction, we claim the reduction by restoring all $d$ coordinates of the gradient [lines 232-233] with the same accuracy as the naive estimator restored just one.

Let us present this explanation more formally:

The baseline estimator $g$ was able to produce, given a point $x$ and direction $e$, an estimate $g(x, Z, e) = p(x, Z, e) \cdot e$ with variance $\\mathbb{E}\_{e, Z} \\|g - \\nabla f\\|^2$ scaling as $\\frac{D\\tau}{n}$, where $n$ is the number of samples used and $D$ is the baseline variance factor (we omit the explicit value of $D$ here, as we focus on the acceleration only).

Denote by $e\_1, \\dots e\_d$ a random basis - the columns of a random rotation matrix $U$ sampled from the uniform (Haar) measure. First, we construct $d$ new estimators $g\_1, \\dots, g\_d$, each working on its own subchain with mixing time $\lceil \\frac{\\tau}{d}\rceil$, along its own direction ($g\_i = p\_i \cdot e\_i$), and using $\\frac{1}{d}$-fraction of all $n$ input samples. As you mentioned in the question, the variance of each of these estimators will scale as $\\frac{D \\lceil \\frac{\\tau}{d} \\rceil}{n/d} \sim \\frac{D\\max(d, \\tau)}{n}$.

Now, keeping in mind that these estimators may be dependent, we analyze the variance of the final estimate $g\_{final} = \\frac{g\_1 + \\dots + g\_d}{d}:$

As we see, while averaging we didn’t require the independence of $\\{g\_i\\}$, as orthogonal directions always have zero correlation. The final expression is very similar to the $\\mathbb{E}\_{e\_1, Z} \\|g\_1(x, Z) - \\nabla f\\|^2$ that we know scales as $\\frac{D\\max(d, \\tau)}{n}$, but the main term $\\mathbb{E}\_{e\_1, Z} \\|g\_1\\|^2$ experiences additional $\\frac{1}{d}$ factor. We show that even after $d$-fold reduction, it is still a dominant term:

Therefore the final reduction is $\\frac{\\max(\\tau, d)}{d\\tau} = \\min(d, \\tau)$.

Q1 ($f\_t$)

We are not sure about the requirement on samples you mention, as there are no estimators presented near the line 84. The function $f\_t$ is a purely theoretic concept with a couple of well known properties: ball smoothing of a Lipschitz function produces a smooth function; finite difference is an unbiased estimator of smoothed function’s gradient [22]. The function $f\_t$ is called the ”smoothed” version of the function $f$. The function $f\_t$ is only a virtual/theoretical object with good properties. Our algorithm does not use it in any way, but it will be used in the theoretical analysis.

[22] Alexander Gasnikov, Anton Novitskii, Vasilii Novitskii, Farshed Abdukhakimov, Dmitry Kamzolov, Aleksandr Beznosikov, Martin Takac, Pavel Dvurechensky, and Bin Gu. The power of first-order smooth optimization for black-box non-smooth problems. Proceedings of the 39th International Conference on Machine Learning, 2022.

Q2 (\pi in problem definition)

It is indeed the chain's invariant distribution [Line 170]. We will add a reference to Assumption 3 to the problem definition (4).

Q4 (total variation)

Yes, it is a total variation distance. We will add it to the Notations & Definitions section, thanks.

Q5 (squares in noise assumptions)

We kept squares to stay in line with other similar assumptions, which require squares either because there is an expectation on LHS, or there is an additional $\\delta^2$ on RHS.

Q6 (definition of $\\lesssim$)

$\\lesssim$ is defined in Table 1 as problem independent, thus the hidden factor is the same. We agree that “problem independent” is not ideal phrasing and will replace it with “independent of all other variables” / “universal”.

Q9 (minibatch variance lower bound)

Thanks for the detailed reading of this part, we fixed the missing $e$, as well as added one more transitional step in (8). Here is how the section now looks like:

In this section, we describe the main tools used to establish the $(d + \\tau)$-type scaling of the error rate. We will focus on reducing the variance and bias of gradient estimators using a specialized batching approach.

We begin by fixing a common building block of our gradient estimators at a point $x$ for both one-point and two-point feedback, as introduced in Section1.1:

These estimators exhibit a twofold randomness that affects how rapidly they concentrate around the true gradient, as we will discuss below.

For clarity, we focus our discussion on the one-point case, although our conclusions extend to the two-point case as well.

A widely used variance reduction technique is mini-batching, where one computes $F(x, Z\_i)$ over a batch of noise variables $\\{Z\_i\\}\_{i=1}^{n}$. The mini-batch gradient estimator is given by:

Let us estimate the scaling of its variance $\\mathbb{E}\_e \\mathbb{E}\_Z \\|\\hat{g}\_{mb} - \\nabla f\\|^2$ with the noise level $\\sigma\_1^2$.

As $\\mathbb{E}\_{Z}\\hat{p}\_{mb} \\approx \\frac{f(x + te) - f(x - te)}{2t} \\approx \\langle\\nabla f, e\\rangle$ we would like to estimate the following for any fixed direction $e$:

An important note here is that the last transition assumes a “worst case” scenario, when inequality in Lemma 1 is tight. For example, it is tight for the lazy markov chain discussed below. With that, we bound the variance:

As for the high-level idea if it is possible to control $\\|\\hat{g}\_{mb} - \\nabla f\\|$, we stress that $\\nabla f$ vanishes after the first step in (8), as we provide lower bound for variance. We also added the following remark:

We provide bounds in a “worst case” scenario, when inequality in Lemma 1 is tight. For example, it is tight for the lazy markov chain discussed below.

Q 3, 7, 8, 10, 11, 12

We fixed it, thank you!

Thanks to Reviewer xojL for their evaluation. Our responses follow.

W1,2 (parameters must be a-priori)

While we agree that this fact limits applying our algorithm in practice, and adaptive algorithms are an interesting direction of research, we feel that this limitation is ubiquitous in optimization and our theoretical contribution is significant enough on its own. We will propose resolving this issue as future work.

Q1 (weakened convexity assumptions)

While extending our result to other settings is obviously beneficial, we do not aim to cover all of them in this paper. This is because our main contribution is not the $(\\tau + d)$ oracle complexity in one regime, but rather the gradient estimate scheme with reduced variance that allows us to do that. This estimate is not tied to the strong convexity and we think it can be used with almost no changes in numerous other Markov ZOO problems, including nkXR’s suggestions (convex, non-convex) and Polyak-Łojasiewicz or weak convexity assumptions.

Q2 (weakened noise assumptions)

As long as Lemma 1 is true the proof should work the same. Investigating the sharp theoretical limits for the proposed proof scheme to work is an interesting direction and “Ergodic mirror descent” paper [14] suggests both Markovianity and uniform boundness of the oracle can be relaxed. Nevertheless, we argue that presenting the first result of this type in our limiting but clear setting is beneficial, as the main insights become easier to understand.

We thank Reviewer Exow for their time, work and positive feedback. Next, we answer the questions raised.

Q4 (upper bound on oracle calls)

Thanks for this question, it turned out to be insightful. We will add the discussion below to the final version of the paper.

The number of iterations $N$ to achieve the accuracy $\\varepsilon$ in expectation is deterministic. The number of oracle calls, however, is random, as the batch size $b_i$ on step $i$ comes from the truncated “log-geometric” distribution

We note here that $l \\sim \\log M$ and $\\mathbb{E} b_i \\sim l\\log M$ are (poly)logarithmic in $M$, but the variance is large, $\\mathbb{E} b_i^2 \\simeq l^2M$. Luckily, the total number of oracle calls $S_N = \\sum\\limits_{i=1}^{N} b_i$ is a finite sum of a large number of such i.i.d. variables and we may hope that the mean value of the sum will be typical (in the sense of high probability, similar to LLN or CLT). The main complication is that $M$ may depend on $N$, thus blocking the use of CLT. Any further discussion requires knowing the dependence between $N$ and $M$, and it can be recovered from the definition of $M$ and stepsize tuning lemma (Lemma 10). For big enough $N$, we can bound $M \\lesssim \\frac{N}{\\log N}$. With this, we apply the finite-time concentration inequality (Bernstein inequality) to our sum:

It shows the subexponential behaviour of the normalized deviation from the mean, thus confirming that the expectation is typical in the high-probability sense.

An important note here is that the normalization of error is of order $N$, and not $\\sqrt{N}$ as in CLT. Our attempt to use the Berry-Esseen type of results failed, as $\\mathbb{E}b_i \\gtrsim M^2$ is much larger than required $\\sqrt{N}$.

Q1 (choice of $p$ is obscured)

Agreed, thank you. In the final version we will clarify that $p$ either does not depend on the problem at all or $\\simeq B / (B + d)$, depending on the setting (one-point or two-point).

Q2 (parameter B is useless)

Fair remark, B overcomplicates understanding the presented results. Our goal was to highlight that the algorithm benefits (in terms of iteration complexity) from having bigger than necessary batches, which was condensed into remark on optimal iteration complexity (line 295). We will add a separate, easier to read corollary in the final version, with $B = 1$.

Q3 (formula for batch size)

True, and it is misrepresented by line 264 and the table of parameters. Batch size is $2^J l$, if $2^J \\leq M$, else $l$. Will be fixed, thanks!

We also plan to add more explanations to the MLMC procedure, here is the draft:

Let’s start with the case $l = 1, M = \\infty, B = 1$ and enumerate baseline estimates as $g_1, g_2, \\dots$.

Then, the MLMC estimate will be $g_1$ with prob. $1/2$, $g_1 + (g_3 - g_2)$ with prob. $1/4$, $g_1 + (g_4 + g_5 - g_2 - g_3)$ with prob. $1/8$ and so on. The parameter $M$ is the upper bound on the number of estimates used. The parameter $l$ transforms the base estimator into a sequence of $l$ base estimators, effectively stretching everything $l$ times. Finally, $B$ serves as a hyperparameter that can multiplicatively increase $l$.

Q3 (notation $\\hat{g}[l]$)

This is a slight abuse of notation of $\\hat{g}_{rd}$, with the rationale that definition would be too bulky if written completely straight. We think the situation could be improved with the explanation above.

Q5 (definition of $\\eta$)

The two terms in this row are equal and both can be used as a definition of $\\eta$. We will clarify this, thanks.

We appreciate Reviewer nkXR’s time, effort, and supportive feedback. Below, we address their questions.

Q1 (Markovianity to general mixing)

You are right, we do not use Markovianity at full potential. As long as Lemma 1 is true the proof should work the same. Investigating the sharp theoretical limits for the proposed proof scheme to work is an interesting direction and “Ergodic mirror descent” paper suggests both Markovianity and uniform boundness of the oracle can be relaxed. Nevertheless, we argue that presenting the first result of this type in our limiting but clear setting is beneficial, as the main insights become easier to understand. We are going to add this as a future work, thanks!

Q2 (strongly convex to general convex)

There is a classical regularization trick that transforms results in strongly convex setting into general convex by adding a quadratic part. For example, given a convex function $g$, we can run the optimization procedure querying the function $f(x) = g(x) + \\frac{\\mu}{2}\\|x\\|^2$ that is strongly convex. The solution produced by the algorithm may be only slightly worse for $g$ than for $f$:

Choosing $\\mu = \\frac{L}{N^2}$ and plugging in our convergence guarantee, we get the desired result:

There were no special reasons for choosing strong convexity except for personal choice.

Q3 (non-smooth)

We believe the question was “is it possible to relax strong convexity and smoothness at the same time?”. The same regularization trick should work in the non-smooth setting. If applied to the Theorem 3’ from Appendix B.2, the resulting rate is

While presenting new corollaries of the result is obviously beneficial, we do not aim to cover all of them in this paper. This is because our main contribution is not the $(\\tau + d)$ oracle complexity in one setting, but rather the gradient estimate scheme with reduced variance that allows us to do that. This estimate is not tied to the strong convexity and we think it can be used with almost no changes in numerous other Markov ZOO problems, including Q4 (non-convex) and xojL’s suggestions (Polyak-Łojasiewicz or weak convexity assumptions).

Q4 (strongly convex to non-convex)

There is no straightforward technique to transfer results from strongly convex setting to non-convex that we are aware of, but we believe the $\\tau + d$ behaviour will remain the same (see answer to the Q3). Particular recent works [1, 2] proposed a method for both non-smooth and non-convex stochastic optimization and may serve ground for the development of our technique.

[1] Guy Kornowski and Ohad Shamir. 2024. An algorithm with optimal dimension-dependence for zero-order nonsmooth nonconvex stochastic optimization. J. Mach. Learn. Res. 25, 1, Article 122 (January 2024).

[2] Yuyang Qiu, Uday Shanbhag, and Farzad Yousefian. Zeroth-Order Methods for Nondifferentiable, Nonconvex, and Hierarchical Federated Optimization, Advances in Neural Information Processing Systems, 36, 2023.

S&W (related insights in distributed optimization)

This connection to distributed optimization is new, thanks! We agree that delayed gradients behave similar to the Markovian case. Nevertheless, the speedup from “free” gradient estimates appears in the zero-order case, but not in first-order (see our comment on lines 226-228). This suggests the nature of acceleration may be different.

The authors' answers all make sense to me.

One thing is that I think the title could be revised: I'd say a more accurate title would be something like

The limited impact of non-i.i.d. noise on zeroth-order optimization

(or perhaps The limited impact of non-i.i.d. noise on gradient-free optimization or somehting). This would be more accurate than the current title, and more accurately capture the message of the paper than

Gradient-Free Approaches is a Key to an Efficient Interaction with Markovian Stochasticity.

I think this is especially true as (i) there is no "key" to gradient-free approaches, (ii) there is nothing Markovian in the paper or noise, and (iii) the current title is just bad grammar.