Proof. Since \(F({\boldsymbol{\theta}}) = \mathbb{E}_{\xi}[f({\boldsymbol{\theta}}; \xi)]\), the Lipschitz continuity of \(F({\boldsymbol{\theta}})\) follows: \[|F({\boldsymbol{\theta}}_1) - F({\boldsymbol{\theta}}_2)| = |\mathbb{E}_{\xi}[f({\boldsymbol{\theta}}_1; \xi) - f({\boldsymbol{\theta}}_2; \xi)]| \overset{(a)}{\leq} \mathbb{E}_{\xi}[|f({\boldsymbol{\theta}}_1; \xi) - f({\boldsymbol{\theta}}_2; \xi)|] \overset{(b)}{\leq} L_0 \norm{{\boldsymbol{\theta}}_1 - {\boldsymbol{\theta}}_2} \;,\] where \((a)\) follows from Jensen’s inequality, and \((b)\) comes from Lem. 4.

The Lipschitz continuity for the gradient and Hessian of \(F({\boldsymbol{\theta}})\) can be proved similarly. ◻

Proof. The unbiasedness of the gradient estimator \(\hat{\nabla} F({\boldsymbol{\theta}})\) follows from Stein’s identity: \[\label{eq:bias-true-grad-deriv} \begin{align} \mathbb{E}[\hat{\nabla} F({\boldsymbol{\theta}})] =& \mathbb{E}\left[\frac{1}{K - 1} \sum_{k=1}^{K} \frac{f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}_k; \xi) - \frac{1}{K} \sum_{k'=1}^{K} f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}_{k'}; \xi)}{\mu} {\mathbf{u}}_k\right] \\ =& \mathbb{E}\left[\frac{1}{K - 1} \sum_{k=1}^{K} \frac{\frac{K - 1}{K} f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}_k; \xi) - \frac{1}{K} \sum_{k'=1, k' \neq k}^{K} f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}_{k'}; \xi)}{\mu} {\mathbf{u}}_k \right] \\ =& \mathbb{E}\left[\frac{1}{K} \sum_{k=1}^{K} \frac{f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}_k; \xi)}{\mu} {\mathbf{u}}_k \right] - \mathbb{E}\left[\frac{1}{K(K - 1)} \sum_{k=1}^{K} \sum_{k'=1, k' \neq k}^{K} \frac{f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}_{k'}; \xi)}{\mu} {\mathbf{u}}_k \right] \\ \overset{(a)}{=}& \mathbb{E}_{\mathbf{u}}\left[\frac{F({\boldsymbol{\theta}}+ \mu {\mathbf{u}})}{\mu} {\mathbf{u}}\right] \overset{(b)}{=} \nabla F_{\mu}({\boldsymbol{\theta}}) \;, \end{align}\tag{15}\] where \((a)\) comes from the dependence of different random directions \(\{{\mathbf{u}}_k\}_{k=1}^K\), and \((b)\) follows from Stein’s identity.

For the bias between the smoothed gradient and the true gradient, we have: \[\norm{\nabla F_{\mu}({\boldsymbol{\theta}}) - \nabla F({\boldsymbol{\theta}})}^2 \overset{(a)}{\le} \mathbb{E}_{{\mathbf{u}}} \left[\norm{\nabla F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) - \nabla F({\boldsymbol{\theta}})}^2\right] \overset{(b)}{\le} L_1^2 \mu^2 \mathbb{E}_{{\mathbf{u}}} \left[\norm{{\mathbf{u}}}^2\right] \overset{(c)}{=} L_1^2 \mu^2 d \;,\] where \((a)\) follows from Jensen’s inequality, \((b)\) comes from Lem. 4, and \((c)\) results from the fact that \(\mathbb{E}_{{\mathbf{u}}}[\norm{{\mathbf{u}}}^2] = d\) for \({\mathbf{u}}\sim {\mathcal{N}}(0, {\mathbf{I}}_d)\).

Moreover, since the gradient estimator \(\hat{\nabla} F({\boldsymbol{\theta}})\) is constructed from \(K\) dependent random directions, we have the variance bound: \[\require{physics} \begin{align} &\mathbb{E}\left[\norm{\hat{\nabla} F({\boldsymbol{\theta}}) - \nabla F_\mu({\boldsymbol{\theta}})}^2\right] \overset{(a)}{=} \mathbb{E}\qty[\norm{\hat{\nabla} F({\boldsymbol{\theta}})}^2] - \norm{\nabla F_\mu({\boldsymbol{\theta}})}^2 \overset{(b)}{\le} \mathbb{E}\qty[\norm{\frac{1}{K - 1} \sum_{k=1}^{K} \frac{f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}_k; \xi) - F_\mu({\boldsymbol{\theta}})}{\mu} {\mathbf{u}}_k}^2] \\ \overset{(c)}{\le}& \frac{K}{(K - 1)^2 \mu^2} \mathbb{E}\qty[\left|f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi) - F_\mu({\boldsymbol{\theta}})\right|^2 \norm{{\mathbf{u}}}^2] \overset{(d)}{\le} \frac{K d \qty(\sigma_\xi^2 + 4 L_0^2 \mu^2 (d + 1))}{(K - 1)^2 \mu^2} \;, \end{align}\] where \((a)\) follows from 15 , \((b)\) comes from getting rid of \(\norm{F_\mu({\boldsymbol{\theta}})}^2\), \((c)\) follows from the dependence of random directions \(\{{\mathbf{u}}_k\}_{k=1}^K\), and \((d)\) results from the following inequality: \[\require{physics} \begin{align} \mathbb{E}\qty[\left|f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi) - F_\mu({\boldsymbol{\theta}})\right|^2 \norm{{\mathbf{u}}}^2] \overset{(a)}{=}& \mathbb{E}\qty[\left|f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi) - F({\boldsymbol{\theta}}+ \mu {\mathbf{u}})\right|^2 \norm{{\mathbf{u}}}^2] + \mathbb{E}\qty[\left|F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) - F_\mu({\boldsymbol{\theta}})\right|^2 \norm{{\mathbf{u}}}^2] \\ \overset{(b)}{\le}& \sigma_\xi^2 \mathbb{E}[\norm{{\mathbf{u}}}^2] + \mathbb{E}\qty[\left|F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) - F_\mu({\boldsymbol{\theta}})\right|^2 \norm{{\mathbf{u}}}^2] \\ \overset{(c)}{\le}& \sigma_\xi^2 d + \mathbb{E}\qty[2 \qty(\left|F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) - F({\boldsymbol{\theta}})\right|^2 + 2 \left|F({\boldsymbol{\theta}}) - \mathbb{E}_{{\mathbf{u}}'} \qty[F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}')]\right|^2) \norm{{\mathbf{u}}}^2] \\ \overset{(d)}{\le}& \sigma_\xi^2 d + \mathbb{E}\qty[2 \qty(\left|F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) - F({\boldsymbol{\theta}})\right|^2 + 2 \mathbb{E}_{{\mathbf{u}}'} \qty[\left|F({\boldsymbol{\theta}}) - F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}')\right|^2]) \norm{{\mathbf{u}}}^2] \\ \overset{(e)}{\le}& \sigma_\xi^2 d + 2 L_0^2 \mu^2 \mathbb{E}\qty[\qty(\norm{{\mathbf{u}}}^2 + \mathbb{E}_{{\mathbf{u}}'} \qty[\norm{{\mathbf{u}}'}^2]) \norm{{\mathbf{u}}}^2] \\ \overset{(f)}{=}& \sigma_\xi^2 d + 2 L_0^2 \mu^2 \mathbb{E}\qty[\qty(\norm{{\mathbf{u}}}^2 + d) \norm{{\mathbf{u}}}^2] \overset{(g)}{=} d \qty(\sigma_\xi^2 + 4 L_0^2 \mu^2 (d + 1)) \;, \end{align}\] where \((a)\) follows from 1 , \((b)\) comes from Assump. 2, \((c)\) results from the fact that \((a + b)^2 \le 2(a^2 + b^2)\), \((d)\) follows from Jensen’s inequality, \((e)\) comes from Lem. 4, and \((c)\), \((f)\), \((g)\) results from the fact that \(\mathbb{E}[\norm{{\mathbf{u}}}^2] = d\) and \(\mathbb{E}[\norm{{\mathbf{u}}}^4] = d(d + 2)\) for \({\mathbf{u}}\sim {\mathcal{N}}(0, {\mathbf{I}}_d)\).

Overall, combining the bias between the smoothed gradient and the true gradient with the variance of the gradient estimator, we have: \[\require{physics} \begin{align} \mathbb{E}\left[\norm{\hat{\nabla} F({\boldsymbol{\theta}}) - \nabla F({\boldsymbol{\theta}})}^2\right] \overset{(a)}{\le}& \mathbb{E}\left[\norm{\hat{\nabla} F({\boldsymbol{\theta}}) - F_\mu({\boldsymbol{\theta}})}^2\right] + \norm{F_\mu({\boldsymbol{\theta}}) - \nabla F({\boldsymbol{\theta}})}^2 \le \frac{K d \qty(\sigma_\xi^2 + 4 L_0^2 \mu^2 (d + 1))}{(K - 1)^2 \mu^2} + L_1^2 \mu^2 d \;, \end{align}\] where \((a)\) follows from 15 . ◻

Proof. \[\norm{\nabla^2 F_{\mu}({\boldsymbol{\theta}}) - \nabla^2 F({\boldsymbol{\theta}})}_F^2 \overset{(a)}{\le} \mathbb{E}_{{\mathbf{u}}} \left[\norm{\nabla^2 F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) - \nabla^2 F({\boldsymbol{\theta}})}_F^2\right] \overset{(b)}{\le} L_2^2 \mu^2 \mathbb{E}_{{\mathbf{u}}} \left[\norm{{\mathbf{u}}}^2\right] \overset{(c)}{=} L_2^2 \mu^2 d \;,\] where \((a)\) follows from Jensen’s inequality, \((b)\) comes from Lem. 4, and \((c)\) results from the fact that \(\mathbb{E}_{{\mathbf{u}}}[\norm{{\mathbf{u}}}^2] = d\) for \({\mathbf{u}}\sim {\mathcal{N}}(0, {\mathbf{I}}_d)\). ◻

Proof. First, \[\require{physics} \label{eq:second-moment-nu-deriv} \begin{align} \mathbb{E}_{{\mathbf{u}}}[\nu^2 \norm{{\mathbf{u}}}^2] =& \frac{1}{\mu^4} \mathbb{E}_{{\mathbf{u}}} \left[\left|f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi) - F_\mu({\boldsymbol{\theta}})\right|^2 \norm{{\mathbf{u}}}^2\right] \\ \overset{(a)}{=}& \frac{1}{\mu^4} \mathbb{E}_{{\mathbf{u}}} \left[\qty(\left|f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi) - F({\boldsymbol{\theta}}+ \mu {\mathbf{u}})\right|^2 + \left|F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) - F_\mu({\boldsymbol{\theta}})\right|^2) \norm{{\mathbf{u}}}^2\right] \\ \overset{(b)}{\le}& \frac{1}{\mu^4} \qty(\sigma_\xi^2 \mathbb{E}_{{\mathbf{u}}}[\norm{{\mathbf{u}}}^2] + \mathbb{E}_{{\mathbf{u}}} \left[\left|F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) - F_\mu({\boldsymbol{\theta}})\right|^2 \norm{{\mathbf{u}}}^2\right]) \\ \overset{(c)}{=}& \frac{1}{\mu^4} \qty(\sigma_\xi^2 d + \mathbb{E}_{{\mathbf{u}}} \left[\left|F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) - F_\mu({\boldsymbol{\theta}})\right|^2 \norm{{\mathbf{u}}}^2\right]) \;, \end{align}\tag{16}\] where \((a)\) follows from 1 , i.e. \(\mathbb{E}_{\xi}[f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi) - F({\boldsymbol{\theta}}+ \mu {\mathbf{u}})] = 0\), \((b)\) comes from Assump. 2, and \((c)\) results from the fact that \(\mathbb{E}_{{\mathbf{u}}}[\norm{{\mathbf{u}}}^2] = d\) for \({\mathbf{u}}\sim {\mathcal{N}}(0, {\mathbf{I}}_d)\).

Furthermore, the second term in 16 can be bounded as: \[\require{physics} \label{eq:second-moment-nu-deriv-2} \begin{align} \mathbb{E}_{{\mathbf{u}}} \left[\left|F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) - F_\mu({\boldsymbol{\theta}})\right|^2 \norm{{\mathbf{u}}}^2\right] \overset{(a)}{\le}& 2 \mathbb{E}_{{\mathbf{u}}} \left[\qty(\left|F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) - F({\boldsymbol{\theta}})\right|^2 + \left|F({\boldsymbol{\theta}}) - \mathbb{E}_{{\mathbf{u}}'} \qty[F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}')]\right|^2) \norm{{\mathbf{u}}}^2\right] \\ \overset{(b)}{\le}& 2 \mathbb{E}_{{\mathbf{u}}} \left[\qty(\left|F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) - F({\boldsymbol{\theta}})\right|^2 + \mathbb{E}_{{\mathbf{u}}'} \qty[\left|F({\boldsymbol{\theta}}) - F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}')\right|^2]) \norm{{\mathbf{u}}}^2\right] \\ \overset{(c)}{\le}& 2 L_0^2 \mu^2 \mathbb{E}_{{\mathbf{u}}} \left[\qty(\norm{{\mathbf{u}}}^2 + \mathbb{E}_{{\mathbf{u}}'} \qty[\norm{{\mathbf{u}}'}^2]) \norm{{\mathbf{u}}}^2\right] \\ \overset{(d)}{=}& 2 L_0^2 \mu^2 \mathbb{E}_{{\mathbf{u}}} \left[\qty(\norm{{\mathbf{u}}}^2 + d) \norm{{\mathbf{u}}}^2\right] \overset{(e)}{=} 4 L_0^2 \mu^2 d (d + 1) \;, \end{align}\tag{17}\] where \((a)\) follows from the fact that \((a + b)^2 \le 2(a^2 + b^2)\), \((b)\) comes from Jensen’s inequality, \((c)\) results from Lem. 4, and \((d)\), \((e)\) results from the fact that \(\mathbb{E}[\norm{{\mathbf{u}}}^2] = d\) and \(\mathbb{E}[\norm{{\mathbf{u}}}^4] = d(d + 2)\) for \({\mathbf{u}}\sim {\mathcal{N}}(0, {\mathbf{I}}_d)\).

Overall, substituting 17 into 16 , we complete the proof: \[\require{physics} \mathbb{E}_{{\mathbf{u}}}[\nu^2 \norm{{\mathbf{u}}}^2] \le \frac{1}{\mu^4} \qty(\sigma_\xi^2 d + 4 L_0^2 \mu^2 d (d + 1)) = \frac{d}{\mu^4} \qty(\sigma_\xi^2 + 4 L_0^2 \mu^2 (d + 1)) \;.\]

Similarly, we have: \[\require{physics} \mathbb{E}_{{\mathbf{u}}}[\nu^2 \norm{{\mathbf{u}}}^4] \le \frac{d (d + 2)}{\mu^4} \qty(\sigma_\xi^2 + 4 L_0^2 \mu^2 (d + 2)) \;.\] \[\require{physics} \mathbb{E}_{{\mathbf{u}}}[\nu^2 \norm{{\mathbf{u}}}^6] \le \frac{d (d + 2)(d + 4)}{\mu^4} \qty(\sigma_\xi^2 + 4 L_0^2 \mu^2 (d + 3)) \;.\] ◻

Proof. The bias of the gradient estimator \(\hat{\nabla} F({\boldsymbol{\theta}})\) follows from Stein’s identity: \[\label{eq:bias-true-grad-deriv-N-greater-than-1} \begin{align} \mathbb{E}[\hat{\nabla} F({\boldsymbol{\theta}}_{t-1})] =& \mathbb{E}\left[\frac{1}{N K - 1} \sum_{n, k = 1}^{N, K} \frac{f({\boldsymbol{\theta}}_{t - n} + \mu {\mathbf{u}}_{t-n, k}; \xi) - \frac{1}{N K} \sum_{n, k = 1}^{N, K} f({\boldsymbol{\theta}}_{t - n} + \mu {\mathbf{u}}_{t-n, k}; \xi)}{\mu} {\mathbf{u}}_{t-n, k}\right] \\ =& \mathbb{E}\left[\frac{1}{N K - 1} \sum_{n, k=1}^{N, K} \frac{\frac{N K - 1}{N K} f({\boldsymbol{\theta}}_{t - n} + \mu {\mathbf{u}}_{t-n, k}; \xi) - \frac{1}{N K} \sum_{\substack{n', k'=1 \\ (n', k') \neq (n, k)}}^{N, K} f({\boldsymbol{\theta}}_{t - n} + \mu {\mathbf{u}}_{t-n, k'}; \xi)}{\mu} {\mathbf{u}}_{t-n, k} \right] \\ =& \mathbb{E}\left[\frac{1}{N K} \sum_{n, k=1}^{N, K} \frac{f({\boldsymbol{\theta}}_{t - n} + \mu {\mathbf{u}}_{t-n, k}; \xi)}{\mu} {\mathbf{u}}_{t-n, k} \right] \\ & \quad - \mathbb{E}\left[\frac{1}{N K(N K - 1)} \sum_{n, k=1}^{N, K} \sum_{\substack{n', k'=1 \\ (n', k') \neq (n, k)}}^{N, K} \frac{f({\boldsymbol{\theta}}_{t - n} + \mu {\mathbf{u}}_{t-n, k'}; \xi)}{\mu} {\mathbf{u}}_{t-n, k} \right] \\ \overset{(a)}{=}& \mathbb{E}_{\mathbf{u}}\left[\frac{1}{N} \sum_{n=1}^N \frac{F({\boldsymbol{\theta}}_{t - n} + \mu {\mathbf{u}}_{t-n})}{\mu} {\mathbf{u}}_{t-n}\right] \overset{(b)}{=} \frac{1}{N} \sum_{n=1}^N \nabla F_{\mu}({\boldsymbol{\theta}}_{t - n}) \;, \end{align}\tag{19}\] where \((a)\) comes from the dependence of different random directions \(\{{\mathbf{u}}_{t-n, k}\}_{k=1}^K\), and \((b)\) follows from Stein’s identity. ◻

Proof. First, at the iteration \(t - n\), we have: \[\require{physics} \label{eq:second-moment-nu-deriv-N-greater-than-1} \begin{align} \mathbb{E}_{{\mathbf{u}}}\qty[\nu^2_{t - n} \norm{{\mathbf{u}}_{t - n}}^2] =& \frac{1}{\mu^4} \mathbb{E}_{{\mathbf{u}}} \left[\left|f({\boldsymbol{\theta}}_{t - n} + \mu {\mathbf{u}}_{t - n}; \xi) - \frac{1}{N} \sum_{n'=1}^N F_\mu({\boldsymbol{\theta}}_{t - n'})\right|^2 \norm{{\mathbf{u}}_{t - n}}^2\right] \\ \overset{(a)}{=}& \frac{1}{\mu^4} \mathbb{E}_{{\mathbf{u}}} \left[\qty(\left|f({\boldsymbol{\theta}}_{t - n} + \mu {\mathbf{u}}_{t - n}; \xi) - F({\boldsymbol{\theta}}_{t - n} + \mu {\mathbf{u}}_{t - n})\right|^2 + \left|F({\boldsymbol{\theta}}_{t - n} + \mu {\mathbf{u}}_{t - n}) - \frac{1}{N} \sum_{n'=1}^N F_\mu({\boldsymbol{\theta}}_{t - n'})\right|^2) \norm{{\mathbf{u}}_{t - n}}^2\right] \\ \overset{(b)}{\le}& \frac{1}{\mu^4} \qty(\sigma_\xi^2 \mathbb{E}_{{\mathbf{u}}}[\norm{{\mathbf{u}}_{t - n}}^2] + \mathbb{E}_{{\mathbf{u}}} \left[\left|F({\boldsymbol{\theta}}_{t - n} + \mu {\mathbf{u}}_{t - n}) - \frac{1}{N} \sum_{n'=1}^N F_\mu({\boldsymbol{\theta}}_{t - n'})\right|^2 \norm{{\mathbf{u}}_{t - n}}^2\right]) \\ \overset{(c)}{=}& \frac{1}{\mu^4} \qty(\sigma_\xi^2 d + \mathbb{E}_{{\mathbf{u}}} \left[\left|F({\boldsymbol{\theta}}_{t - n} + \mu {\mathbf{u}}_{t - n}) - \frac{1}{N} \sum_{n'=1}^N F_\mu({\boldsymbol{\theta}}_{t - n'})\right|^2 \norm{{\mathbf{u}}_{t - n}}^2\right]) \;, \end{align}\tag{20}\] where \((a)\) follows from 1 , i.e. \(\mathbb{E}_{\xi}[f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi) - F({\boldsymbol{\theta}}+ \mu {\mathbf{u}})] = 0\), \((b)\) comes from Assump. 2, and \((c)\) results from the fact that \(\mathbb{E}_{{\mathbf{u}}}[\norm{{\mathbf{u}}}^2] = d\) for \({\mathbf{u}}\sim {\mathcal{N}}(0, {\mathbf{I}}_d)\).

Furthermore, the second term in 20 can be bounded as: \[\require{physics} \label{eq:second-moment-nu-deriv-2-N-greater-than-1} \begin{align} &\mathbb{E}_{{\mathbf{u}}} \left[\left|F({\boldsymbol{\theta}}_{t - n} + \mu {\mathbf{u}}_{t - n}) - \frac{1}{N} \sum_{n'=1}^N F_\mu({\boldsymbol{\theta}}_{t - n'})\right|^2 \norm{{\mathbf{u}}_{t - n}}^2\right] \\ \overset{(a)}{=}& \mathbb{E}_{{\mathbf{u}}} \left[\qty(\left|F({\boldsymbol{\theta}}_{t - n} + \mu {\mathbf{u}}_{t - n}) - F_\mu({\boldsymbol{\theta}}_{t - n})\right|^2 + \left|F_\mu({\boldsymbol{\theta}}_{t - n}) - \frac{1}{N} \sum_{n'=1}^N F_\mu({\boldsymbol{\theta}}_{t - n'})\right|^2) \norm{{\mathbf{u}}_{t - n}}^2\right] \\ \overset{(b)}{\le}& 4 L_0^2 \mu^2 d (d + 1) + \left|F_\mu({\boldsymbol{\theta}}_{t - n}) - \frac{1}{N} \sum_{n'=1}^N F_\mu({\boldsymbol{\theta}}_{t - n'})\right|^2 \cdot \mathbb{E}_{{\mathbf{u}}} \left[\norm{{\mathbf{u}}_{t - n}}^2\right] \\ =& 4 L_0^2 \mu^2 d (d + 1) + \left|\frac{1}{N} \sum_{n'=1}^N \qty(F_\mu({\boldsymbol{\theta}}_{t - n}) - F_\mu({\boldsymbol{\theta}}_{t - n'}))\right|^2 \cdot d \\ \overset{(c)}{\le}& 4 L_0^2 \mu^2 d (d + 1) + \frac{d L_0^2}{N} \sum_{n'=1}^N \norm{{\boldsymbol{\theta}}_{t - n} - {\boldsymbol{\theta}}_{t - n'}}^2 \\ \overset{(d)}{\le}& 4 L_0^2 \mu^2 d (d + 1) + \frac{d L_0^2 \eta^2 \mu^2 K C_\nu}{\lambda^2 N (N K - 1)^2} \cdot \sum_{n'=1}^N |n - n'| \;, \end{align}\tag{21}\] where \((a)\) follows from the fact that \(\mathbb{E}_{{\mathbf{u}}}[F({\boldsymbol{\theta}}+ \mu {\mathbf{u}})] = F_\mu({\boldsymbol{\theta}})\), \((b)\) comes from 17 , \((c)\) results from Lem. 4, and \((d)\) results from: \[\begin{align} \norm{{\boldsymbol{\theta}}_{t - n} - {\boldsymbol{\theta}}_{t - n'}}^2 \overset{(a)}{\le}& \sum_{n'' = \min(n, n')}^{\max(n, n')} \norm{{\boldsymbol{\theta}}_{t - n''} - {\boldsymbol{\theta}}_{t - n'' + 1}}^2 \le \eta^2 \sum_{n'' = \min(n, n')}^{\max(n, n')} \norm{\widetilde{H}^{-1}({\boldsymbol{\theta}}_{t - n''}) \hat{\nabla} F({\boldsymbol{\theta}}_{t - n''})}^2 \\ \overset{(b)}{\le}& \frac{\eta^2 \mu^2 K C_\nu}{\lambda^2 (N K - 1)^2} |n - n'| \;, \end{align}\] where \((a)\) follows from the triangle inequality, and \((b)\) follows from the result in 36 , and we define the constant upper bound \(C_\nu\) such that \(\require{physics} \sum_n \mathbb{E}\qty[\nu_{t-n}^2 \norm{{\mathbf{u}}_{t-n}}^2] \le C\) for all \(n\).

Overall, substituting 21 into 20 , and summing over \(n\), we have: \[\require{physics} \begin{align} \sum_{n=1}^N \mathbb{E}_{{\mathbf{u}}}[\nu_{t-n}^2 \norm{{\mathbf{u}}_{t-n}}^2] \le& \frac{1}{\mu^4} \qty(\sigma_\xi^2 d N + 4 L_0^2 \mu^2 d (d + 1) N + \frac{d L_0^2 \eta^2 \mu^2 K C_\nu}{\lambda^2 N (N K - 1)^2} \cdot \sum_{n, n' = 1}^N |n - n'|) \\ =& \frac{1}{\mu^4} \qty(\sigma_\xi^2 d N + 4 L_0^2 \mu^2 d (d + 1) N + \frac{d L_0^2 \eta^2 \mu^2 K C_\nu (N^2 - 1)}{3 \lambda^2 (N K - 1)^2}) \;. \end{align}\]

Note that the definition of \(C_\nu\) implies that: \[\require{physics} C_\nu = \frac{1}{\mu^4} \qty(\sigma_\xi^2 d N + 4 L_0^2 \mu^2 d (d + 1) N + \frac{d L_0^2 \eta^2 \mu^2 K C_\nu (N^2 - 1)}{3 \lambda^2 (N K - 1)^2}) \;,\] \[\require{physics} C_\nu = \frac{1}{\mu^4} \qty(\sigma_\xi^2 d (d + 2) N + 4 L_0^2 \mu^2 d (d + 2)^2 N + \frac{d (d + 2) L_0^2 \eta^2 \mu^2 K C_\nu (N^2 - 1)}{3 \lambda^2 (N K - 1)^2}) \;,\] which leads to: \[\require{physics} \sum_{n=1}^N \mathbb{E}\qty[\nu_{t-n}^2 \norm{{\mathbf{u}}_{t-n}}^2] \le \frac{\lambda^2 N (N K - 1)^2 d }{\lambda^2 \mu^4 (N K - 1)^2 - L_0^2 \eta^2 \mu^2 K (N^2 - 1) d / 3} \qty(\sigma_\xi^2 + 4 L_0^2 \mu^2 (d + 1)) \;.\]

Similarly, we have: \[\require{physics} \sum_{n=1}^N \mathbb{E}_{{\mathbf{u}}}\qty[\nu_{t-n}^2 \norm{{\mathbf{u}}_{t-n}}^4] \le \frac{\lambda^2 N (N K - 1)^2 d (d + 2)}{\lambda^2 \mu^4 (N K - 1)^2 - L_0^2 \eta^2 \mu^2 K (N^2 - 1) d (d + 2) / 3} \qty(\sigma_\xi^2 + 4 L_0^2 \mu^2 (d + 2)) \;,\] \[\require{physics} \sum_{n=1}^N \mathbb{E}_{{\mathbf{u}}}\qty[\nu_{t-n}^2 \norm{{\mathbf{u}}_{t-n}}^6] \le \frac{\lambda^2 N (N K - 1)^2 d (d + 2) (d + 4)}{\lambda^2 \mu^4 (N K - 1)^2 - L_0^2 \eta^2 \mu^2 K (N^2 - 1) d (d + 2) (d + 4) / 3} \qty(\sigma_\xi^2 + 4 L_0^2 \mu^2 (d + 3)) \;.\] ◻

Proof. We first derive the first derivative of the single-step policy optimization objective 7 by applying the Policy Gradient Thm. [@pg-thm]: \[\nabla F_\mu({\boldsymbol{\theta}}) = \mathbb{E}_{{\mathbf{x}}\sim \pi_{{\boldsymbol{\theta}}}({\mathbf{x}})}\left[ \nabla \ln \pi_{{\boldsymbol{\theta}}}({\mathbf{x}}) \mathbb{E}_{\xi}[f({\mathbf{x}}; \xi)] \right] \;.\]

Taking the second derivative, we have: \[\require{physics} \label{eq:hess-zoo-deriv-1} \begin{align} \nabla^2 F_\mu({\boldsymbol{\theta}}) =& \nabla \int \pi_{{\boldsymbol{\theta}}}({\mathbf{x}}) \nabla \ln \pi_{{\boldsymbol{\theta}}}({\mathbf{x}}) \mathbb{E}_{\xi}[f({\mathbf{x}}; \xi)] d{\mathbf{x}}\overset{(a)}{=} \int \qty(\nabla \pi_{{\boldsymbol{\theta}}}({\mathbf{x}}) \qty(\nabla \ln \pi_{{\boldsymbol{\theta}}}({\mathbf{x}}))^{\top} + \pi_{{\boldsymbol{\theta}}}({\mathbf{x}}) \nabla^2 \ln \pi_{{\boldsymbol{\theta}}}({\mathbf{x}})) \mathbb{E}_{\xi}[f({\mathbf{x}}; \xi)] d{\mathbf{x}}\\ \overset{(b)}{=}& \int \pi_{{\boldsymbol{\theta}}}({\mathbf{x}}) \qty(\qty(\nabla \ln \pi_{{\boldsymbol{\theta}}}({\mathbf{x}})) \qty(\nabla \ln \pi_{{\boldsymbol{\theta}}}({\mathbf{x}}))^{\top} + \nabla^2 \ln \pi_{{\boldsymbol{\theta}}}({\mathbf{x}})) \mathbb{E}_{\xi}[f({\mathbf{x}}; \xi)] d{\mathbf{x}}\\ =& \mathbb{E}_{{\mathbf{x}}\sim \pi_{{\boldsymbol{\theta}}}({\mathbf{x}})}\left[ \left(\left(\nabla \ln \pi_{{\boldsymbol{\theta}}}({\mathbf{x}})\right)\left(\nabla \ln \pi_{{\boldsymbol{\theta}}}({\mathbf{x}})\right)^{\top} + \nabla^2 \ln \pi_{{\boldsymbol{\theta}}}({\mathbf{x}})\right)\mathbb{E}_{\xi}[f({\mathbf{x}}; \xi)] \right] \;, \end{align}\tag{22}\] where \((a)\) comes from the product rule of differentiation, and \((b)\) follows from the fact that \(\nabla \pi_{{\boldsymbol{\theta}}}({\mathbf{x}}) = \pi_{{\boldsymbol{\theta}}}({\mathbf{x}}) \nabla \ln \pi_{{\boldsymbol{\theta}}}({\mathbf{x}})\). ◻

Proof. Taking the reparametrization trick \({\mathbf{x}}= {\boldsymbol{\theta}}+ \mu {\mathbf{u}}\) with \({\mathbf{u}}\sim {\mathcal{N}}(0, {\mathbf{I}}_d)\), we have: \[\label{eq:hess-zoo-deriv-2} \nabla^2 F_\mu({\boldsymbol{\theta}}) \overset{(a)}{=} \mathbb{E}_{\mathbf{u}}\left[ \left({\mathbf{u}}{\mathbf{u}}^{\top} - {\mathbf{I}}_d\right) \frac{\mathbb{E}_{\xi}[f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi)]}{\mu^2} \right] \;,\tag{23}\] where \((a)\) comes from \(\nabla \ln \pi_{{\boldsymbol{\theta}}}({\mathbf{x}}) = \frac{1}{\mu^2}({\mathbf{x}}- {\boldsymbol{\theta}}) = \frac{1}{\mu}{\mathbf{u}}\) and \(\nabla^2 \ln \pi_{{\boldsymbol{\theta}}}({\mathbf{x}}) = -\frac{1}{\mu^2}{\mathbf{I}}_d\) for Gaussian policy.

Since \(\mathbb{E}_{\mathbf{u}}[{\mathbf{u}}{\mathbf{u}}^{\top}] = {\mathbf{I}}_d\), for any constant baseline \(b\) that is independent of \({\mathbf{u}}\), 23 can be further transformed as: \[\label{eq:hess-zoo-deriv-3} \nabla^2 F_\mu({\boldsymbol{\theta}}) = \mathbb{E}_{\mathbf{u}}\left[ \left({\mathbf{u}}{\mathbf{u}}^{\top} - {\mathbf{I}}_d\right) \frac{\mathbb{E}_{\xi}[f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi)]}{\mu^2} \right] + b \cdot \mathbb{E}_{\mathbf{u}}\left[ \left({\mathbf{u}}{\mathbf{u}}^{\top} - {\mathbf{I}}_d\right)\right] = \mathbb{E}_{\mathbf{u}}\left[ \left({\mathbf{u}}{\mathbf{u}}^{\top} - {\mathbf{I}}_d\right) \frac{\mathbb{E}_{\xi}[f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi) - b]}{\mu^2} \right] \;,\tag{24}\] which completes the proof. ◻

Proof. Regardless of the choice of baseline \(b\), since \(\mathbb{E}_{\mathbf{u}}[{\mathbf{u}}{\mathbf{u}}^{\top}] = {\mathbf{I}}_d\), Prop. 1 can always be transformed as: \[\label{eq:hess-zoo-deriv-4} \begin{align} \nabla^2 F_\mu({\boldsymbol{\theta}}) =& \mathbb{E}_{\mathbf{u}}\left[ \frac{f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi) - b}{\mu^2} {\mathbf{u}}{\mathbf{u}}^{\top} \right] - \mathbb{E}_{\mathbf{u}}\left[ \mathbb{E}_{\mathbf{u}}\left[\frac{f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi) - b}{\mu^2}\right] {\mathbf{u}}{\mathbf{u}}^{\top}\right] \\ =& \mathbb{E}_{\mathbf{u}}\left[ \frac{f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi) - b - \mathbb{E}_{\mathbf{u}}\left[f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi) - b\right]}{\mu^2} {\mathbf{u}}{\mathbf{u}}^{\top} \right] \\ =& \mathbb{E}_{\mathbf{u}}\left[ \frac{f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi) - \mathbb{E}_{\mathbf{u}}\left[f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi)\right]}{\mu^2} {\mathbf{u}}{\mathbf{u}}^{\top} \right] \;, \end{align}\tag{25}\] which completes the proof. ◻

Proof. Suppose there exists another baseline \(b' \neq \mathbb{E}_{{\mathbf{u}}\sim {\mathcal{N}}(0, {\mathbf{I}}_d)}\left[\mathbb{E}_{\xi}\left[f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}; \xi)\right]\right]\) that also satisfies the unbiasedness: \[\require{physics} \mathbb{E}_{\mathbf{u}}\left[ \frac{\mathbb{E}_{\xi}\left[f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}; \xi) - b'\right]}{\mu^2} {\mathbf{u}}{\mathbf{u}}^{\top} \right] = \nabla^2 F_\mu({\boldsymbol{\theta}}) = \mathbb{E}_{\mathbf{u}}\left[ \frac{\mathbb{E}_{\xi}\left[f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}; \xi)\right]}{\mu^2} \qty({\mathbf{u}}{\mathbf{u}}^{\top} - {\mathbf{I}}_d)\right] \;.\]

Take the difference on both sides, we have: \[\require{physics} \begin{align} 0 =& \mathbb{E}_{\mathbf{u}}\left[\frac{\mathbb{E}_{\xi}\left[f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}; \xi) - b'\right]}{\mu^2} {\mathbf{u}}{\mathbf{u}}^{\top}\right] - \mathbb{E}_{\mathbf{u}}\left[ \frac{\mathbb{E}_{\xi}\left[f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}; \xi)\right]}{\mu^2} \qty({\mathbf{u}}{\mathbf{u}}^{\top} - {\mathbf{I}}_d)\right] \\ =& \mathbb{E}_{\mathbf{u}}\left[-\frac{b'}{\mu^2} {\mathbf{u}}{\mathbf{u}}^{\top}\right] + \mathbb{E}_{\mathbf{u}}\left[\frac{\mathbb{E}_{\xi}\left[f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}; \xi)\right]}{\mu^2} {\mathbf{I}}_d\right] \\ =& \frac{\mathbb{E}_{\mathbf{u}}\left[\mathbb{E}_{\xi}\left[f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}; \xi)\right]\right] - b'}{\mu^2} {\mathbf{I}}_d \;, \end{align}\] which implies that \(b' = \mathbb{E}_{\mathbf{u}}\left[\mathbb{E}_{\xi}\left[f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}; \xi)\right]\right]\), contradicting the assumption. Therefore, the average baseline is the unique choice that guarantees the unbiasedness of the Hessian estimator in Prop. 2. ◻

Proof. We rewrite the smoothed objective using the importance weight \(w({\mathbf{x}}) = \pi_{{\boldsymbol{\theta}}}({\mathbf{x}})/\rho({\mathbf{x}})\). The Hessian is \(\nabla^2 \int f({\mathbf{x}}) \frac{\pi_{{\boldsymbol{\theta}}}({\mathbf{x}})}{\rho({\mathbf{x}})} \rho({\mathbf{x}}) d{\mathbf{x}}\). Since \(\rho\) is independent of \({\boldsymbol{\theta}}\), the derivatives apply solely to \(\pi_{{\boldsymbol{\theta}}}\). Interchanging differentiation and integration yields \(\mathbb{E}_{{\mathbf{x}}\sim \rho}[f({\mathbf{x}}) \frac{1}{\rho({\mathbf{x}})} \nabla^2 \pi_{{\boldsymbol{\theta}}}({\mathbf{x}})]\). The result follows by expanding \(\nabla^2 \pi_{{\boldsymbol{\theta}}} = \pi_{{\boldsymbol{\theta}}} (\nabla^2 \ln \pi_{{\boldsymbol{\theta}}} + \nabla \ln \pi_{{\boldsymbol{\theta}}} \nabla \ln \pi_{{\boldsymbol{\theta}}}^\top)\).

We start from the Hessian of the smoothed objective in 23 : \[\nabla^2 F_\mu({\boldsymbol{\theta}}) = \mathbb{E}_{\mathbf{u}}\left[ \left({\mathbf{u}}{\mathbf{u}}^{\top} - {\mathbf{I}}_d\right) \frac{\mathbb{E}_{\xi}[f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi)]}{\mu^2} \right] = \int \left({\mathbf{u}}{\mathbf{u}}^{\top} - {\mathbf{I}}_d\right) \frac{\mathbb{E}_{\xi}[f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi)]}{\mu^2} \pi_{{\boldsymbol{\theta}}}({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) d{\mathbf{u}}\;,\] where we rewrite the expectation over \({\mathbf{u}}\) as an integral over the Gaussian policy \(\pi_{{\boldsymbol{\theta}}}({\mathbf{x}})\).

Under the assumptions \(\text{supp}(\pi_{{\boldsymbol{\theta}}}) \subseteq \text{supp}(\rho)\) that ensures the importance sampling ratio \(\mathcal{W}_{\pi/\rho}({\mathbf{u}}) \triangleq \frac{\pi_{{\boldsymbol{\theta}}}({\boldsymbol{\theta}}+ \mu {\mathbf{u}})}{\rho({\mathbf{u}})}\) is well-defined, we can express the Hessian of smoothed objective as: \[\nabla^2 F_\mu({\boldsymbol{\theta}}) = \int \left({\mathbf{u}}{\mathbf{u}}^{\top} - {\mathbf{I}}_d\right) \frac{\mathbb{E}_{\xi}[f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi)]}{\mu^2} \mathcal{W}_{\pi/\rho}({\mathbf{u}}) \rho({\mathbf{u}}) d{\mathbf{u}}= \mathbb{E}_{{\mathbf{u}}\sim \rho}\left[ \left({\mathbf{u}}{\mathbf{u}}^{\top} - {\mathbf{I}}_d\right) \frac{\mathbb{E}_{\xi}[f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi)]}{\mu^2} \mathcal{W}_{\pi/\rho}({\mathbf{u}}) \right] \;.\]

Finally, similar to the derivation in 23 , we have: \[\nabla^2 F_\mu({\boldsymbol{\theta}}) = \mathbb{E}_{{\mathbf{u}}\sim \rho}\left[ \frac{\mathbb{E}_{\xi}[f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi)] - \mathbb{E}_{{\mathbf{u}}\sim \rho}\left[\mathbb{E}_{\xi}[f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi)] \mathcal{W}_{\pi/\rho}({\mathbf{u}})\right]}{\mu^2} {\mathbf{u}}{\mathbf{u}}^{\top} \cdot \mathcal{W}_{\pi/\rho}({\mathbf{u}}) \right] \;.\] ◻

Proof. The gradient of the smoothed objective \(F_\mu({\boldsymbol{\theta}})\) is a direct result of Stein’s identity: \[\nabla F_\mu({\boldsymbol{\theta}}) = \mathbb{E}_{\mathbf{u}}\left[\nabla F({\boldsymbol{\theta}}+ \mu {\mathbf{u}})\right] = \frac{1}{\mu} \mathbb{E}_{\mathbf{u}}\left[\nabla_{\mathbf{u}}F({\boldsymbol{\theta}}+ \mu {\mathbf{u}})\right] \overset{(a)}{=} \frac{1}{\mu} \mathbb{E}_{\mathbf{u}}\left[F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) {\mathbf{u}}\right] \;,\] where \((a)\) follows from Stein’s identity.

Then, the Hessian of the smoothed objective \(F_\mu({\boldsymbol{\theta}})\) can be derived by taking another derivative: \[\require{physics} \nabla^2 F_\mu({\boldsymbol{\theta}}) = \qty(\frac{1}{\mu} \mathbb{E}_{\mathbf{u}}\left[\qty(\nabla F({\boldsymbol{\theta}}+ \mu {\mathbf{u}})) {\mathbf{u}}\right]) = \frac{1}{\mu^2} \mathbb{E}_{\mathbf{u}}\left[\qty(\nabla_{\mathbf{u}}F({\boldsymbol{\theta}}+ \mu {\mathbf{u}})) {\mathbf{u}}^{\top}\right] \;.\] Note that \(\require{physics} \nabla_{\mathbf{u}}\qty(F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) {\mathbf{u}}) = \qty(\nabla_{\mathbf{u}}F({\boldsymbol{\theta}}+ \mu {\mathbf{u}})) {\mathbf{u}}^{\top} + F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) {\mathbf{I}}_d\), we have: \[\require{physics} \nabla^2 F_\mu({\boldsymbol{\theta}}) = \frac{1}{\mu^2} \mathbb{E}_{\mathbf{u}}\left[\nabla_{\mathbf{u}}\qty(F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) {\mathbf{u}}) - F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) {\mathbf{I}}_d\right] \overset{(a)}{=} \frac{1}{\mu^2} \mathbb{E}_{\mathbf{u}}\left[F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) \qty({\mathbf{u}}{\mathbf{u}}^{\top} - {\mathbf{I}}_d)\right] \;,\] where \((a)\) again follows from Stein’s identity. ◻

Proof. We first simplify the Hessian estimator in 8 by: \[\begin{align} \widehat{{\mathbf{H}}}({\boldsymbol{\theta}}) =& \frac{1}{K - 1} \sum_{k=1}^{K} \frac{f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_k; \xi) - \frac{1}{K} \sum_{k'=1}^{K} f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_{k'}; \xi)}{\mu^2} {\mathbf{u}}_k{\mathbf{u}}_k^{\top} \\ =& \frac{1}{K - 1} \sum_{k=1}^{K} \frac{\frac{K - 1}{K} f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_k; \xi) - \frac{1}{K} \sum_{\substack{k'=1 \\ k' \neq k}}^{K} f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_{k'}; \xi)}{\mu^2} {\mathbf{u}}_k{\mathbf{u}}_k^{\top} \;. \end{align}\]

Take the expectation over \(\{{\mathbf{u}}_k\}_{k=1}^{K}\) and \(\xi\) on both sides, we have: \[\require{physics} \begin{align} \mathbb{E}\qty[\widehat{{\mathbf{H}}}({\boldsymbol{\theta}})] =& \frac{1}{K - 1} \sum_{k=1}^{K} \mathbb{E}_{{\mathbf{u}}} \mathbb{E}_{\xi} \left[\frac{\frac{K - 1}{K} f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_k; \xi) - \frac{1}{K} \sum_{\substack{k'=1 \\ k' \neq k}}^{K} f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_{k'}; \xi)}{\mu^2} {\mathbf{u}}_k{\mathbf{u}}_k^{\top} \right] \\ \overset{(a)}{=}& \frac{1}{K - 1} \sum_{k=1}^{K} \mathbb{E}_{{\mathbf{u}}} \left[\frac{\frac{K - 1}{K} F({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_k) - \frac{1}{K} \sum_{\substack{k'=1 \\ k' \neq k}}^{K} F({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_{k'})}{\mu^2} {\mathbf{u}}_k{\mathbf{u}}_k^{\top} \right] \\ =& \frac{1}{K - 1} \sum_{k=1}^{K} \frac{\frac{K - 1}{K} \mathbb{E}_{{\mathbf{u}}} \qty[F({\boldsymbol{\theta}}+ \mu{\mathbf{u}}) {\mathbf{u}}{\mathbf{u}}^{\top}] - \frac{1}{K} \sum_{\substack{k'=1 \\ k' \neq k}}^{K} \mathbb{E}_{{\mathbf{u}}} \qty[F({\boldsymbol{\theta}}+ \mu{\mathbf{u}}) {\mathbf{I}}_d]}{\mu^2} \\ =& \frac{1}{K - 1} \sum_{k=1}^{K} \frac{\frac{K - 1}{K} \mathbb{E}_{{\mathbf{u}}} \qty[F({\boldsymbol{\theta}}+ \mu{\mathbf{u}}) {\mathbf{u}}{\mathbf{u}}^{\top}] - \frac{K - 1}{K} \mathbb{E}_{{\mathbf{u}}} \qty[F({\boldsymbol{\theta}}+ \mu{\mathbf{u}}) {\mathbf{I}}_d]}{\mu^2} \\ =& \mathbb{E}_{{\mathbf{u}}} \left[\frac{F({\boldsymbol{\theta}}+ \mu{\mathbf{u}}) \qty({\mathbf{u}}{\mathbf{u}}^{\top} - {\mathbf{I}}_d)}{\mu^2}\right] \\ \overset{(b)}{=}& \nabla^2 F_\mu({\boldsymbol{\theta}}) \;, \end{align}\] where \((a)\) follows from the independence of \(\frac{1}{K} \sum_{\substack{k'=1 \\ k' \neq k}}^{K} f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_{k'}; \xi)\) and \({\mathbf{u}}_k{\mathbf{u}}_k^{\top}\), and \(\mathbb{E}_{{\mathbf{u}}} [u u^{\top}] = {\mathbf{I}}_d\). Step \((b)\) follows from Lem. 3. ◻

Proof. Since the baseline \(b\) is a free parameter in the Hessian estimator, to ensure unbiasedness, we consider the Monte Carlo implementation of Prop. 1: \[\require{physics} \label{eq:hess-est-gaussian} \widehat{{\mathbf{H}}}({\boldsymbol{\theta}}) = \frac{1}{K} \sum_{k=1}^{K} \frac{f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}_k; \xi) - b}{\mu^2} \qty({\mathbf{u}}_k{\mathbf{u}}_k^{\top} - {\mathbf{I}}_d) \;,\tag{26}\] where the baseline \(b\) is a constant independent of the random directions \(\{{\mathbf{u}}_k\}_{k=1}^{K}\).

We verify that the Hessian estimator in 26 is unbiased: \[\require{physics} \label{eq:hess-est-gaussian-bias} \mathbb{E}\qty[\widehat{{\mathbf{H}}}({\boldsymbol{\theta}})] \overset{(a)}{=} \mathbb{E}_{{\mathbf{u}}} \left[\frac{F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) - b}{\mu^2} \qty({\mathbf{u}}{\mathbf{u}}^{\top} - {\mathbf{I}}_d) \right] \overset{(b)}{=} \mathbb{E}_{{\mathbf{u}}} \left[\frac{F({\boldsymbol{\theta}}+ \mu {\mathbf{u}})}{\mu^2} \qty({\mathbf{u}}{\mathbf{u}}^{\top} - {\mathbf{I}}_d) \right] \overset{(c)}{=} \nabla^2 F_\mu({\boldsymbol{\theta}}) \;,\tag{27}\] where \((a)\) follows from the independence of \(\{{\mathbf{u}}_k\}\) across different queries \(k\), \((b)\) uses the property \(\mathbb{E}_{{\mathbf{u}}}[{\mathbf{u}}{\mathbf{u}}^{\top}] = {\mathbf{I}}_d\), and \((c)\) follows from Stein’s identity.

Next, the variance of the Hessian estimator in 26 is given by: \[\require{physics} \label{eq:var-decomp} \begin{align} &\mathrm{Var}\qty(\widehat{{\mathbf{H}}}({\boldsymbol{\theta}})) = \mathbb{E}\qty[\norm{\widehat{{\mathbf{H}}}({\boldsymbol{\theta}}) - \nabla^2 F_{\mu}({\boldsymbol{\theta}})}_F^2] = \mathbb{E}\left[ \norm{\widehat{{\mathbf{H}}}({\boldsymbol{\theta}}) }_F^2 \right] - 2 \Tr\left( \mathbb{E}\left[\widehat{{\mathbf{H}}}({\boldsymbol{\theta}})\right]^{\top} \nabla^2 F_{\mu}({\boldsymbol{\theta}}) \right) + \norm{\nabla^2 F_{\mu}({\boldsymbol{\theta}})}_F^2 \\ \overset{(a)}{=}& \mathbb{E}\left[ \norm{\frac{1}{K} \sum_{k=1}^{K} \frac{f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_k; \xi) - b}{\mu^2} \qty({\mathbf{u}}_k{\mathbf{u}}_k^{\top} - {\mathbf{I}}_d) }_F^2 \right] - \norm{\nabla^2 F_{\mu}({\boldsymbol{\theta}})}_F^2 \;, \end{align}\tag{28}\] where \((a)\) follows from 27 .

Differentiating the variance with respect to the baseline \(b\) yields: \[\require{physics} \begin{align} \frac{\partial}{\partial b} \mathrm{Var}\qty(\widehat{{\mathbf{H}}}({\boldsymbol{\theta}})) =& 2 \mathbb{E}\left[ \left\langle \frac{1}{K} \sum_{k=1}^{K} \frac{f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_k; \xi) - b}{\mu^2} \qty({\mathbf{u}}_k{\mathbf{u}}_k^{\top} - {\mathbf{I}}_d), \frac{1}{K} \sum_{k=1}^{K} \qty(-\frac{1}{\mu^2}) \qty({\mathbf{u}}_k{\mathbf{u}}_k^{\top} - {\mathbf{I}}_d) \right\rangle_F \right] \\ \overset{(a)}{=}& -\frac{2}{K^2 \mu^4} \sum_{k = 1}^{K} \mathbb{E}\left[ \left\langle \qty(f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_k; \xi) - b) \qty({\mathbf{u}}_k{\mathbf{u}}_k^{\top} - {\mathbf{I}}_d), \qty({\mathbf{u}}_k{\mathbf{u}}_k^{\top} - {\mathbf{I}}_d) \right\rangle_F \right] \\ =& -\frac{2}{K^2 \mu^4} \sum_{k = 1}^{K} \mathbb{E}_{\mathbf{u}}\left[ \left\langle \qty(F({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_k) - b) \qty({\mathbf{u}}_k{\mathbf{u}}_k^{\top} - {\mathbf{I}}_d), \qty({\mathbf{u}}_k{\mathbf{u}}_k^{\top} - {\mathbf{I}}_d) \right\rangle_F \right] \\ =& -\frac{2}{K^2 \mu^4} \sum_{k = 1}^{K} \mathbb{E}_{\mathbf{u}}\left[ \left(F({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_k) - b\right) \norm{\qty({\mathbf{u}}_k{\mathbf{u}}_k^{\top} - {\mathbf{I}}_d)}_F^2 \right] \;, \end{align}\] where \((a)\) follows from the independence of \(\{{\mathbf{u}}_k\}\) across different queries \(k\).

Setting \(\require{physics} \frac{\partial}{\partial b} \mathrm{Var}\qty(\widehat{{\mathbf{H}}}({\boldsymbol{\theta}})) = 0\), we have: \[\require{physics} b^* = \frac{\sum_{k = 1}^{K} \mathbb{E}_{\mathbf{u}}\left[ F({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_k) \norm{\qty({\mathbf{u}}_k{\mathbf{u}}_k^{\top} - {\mathbf{I}}_d)}_F^2 \right]}{\sum_{k = 1}^{K} \mathbb{E}_{\mathbf{u}}\left[ \norm{\qty({\mathbf{u}}_k{\mathbf{u}}_k^{\top} - {\mathbf{I}}_d)}_F^2 \right]} = \frac{\mathbb{E}_{\mathbf{u}}\left[ F({\boldsymbol{\theta}}+ \mu{\mathbf{u}}) \norm{\qty({\mathbf{u}}{\mathbf{u}}^{\top} - {\mathbf{I}}_d)}_F^2 \right]}{\mathbb{E}_{\mathbf{u}}\left[ \norm{\qty({\mathbf{u}}{\mathbf{u}}^{\top} - {\mathbf{I}}_d)}_F^2 \right]} \;.\]

In particular, when \({\mathbf{u}}\sim {\mathcal{N}}(0, {\mathbf{I}}_d)\) as \(d \to \infty\), such that \(\require{physics} \norm{\qty({\mathbf{u}}_k{\mathbf{u}}_k^{\top} - {\mathbf{I}}_d)}_F^2\) is effectively constant, the optimal baseline reduces to the expected function value: \[b^* = \mathbb{E}_{\mathbf{u}}\left[ F({\boldsymbol{\theta}}+ \mu{\mathbf{u}}) \right] = F_\mu({\boldsymbol{\theta}}) \;.\]

Finally, following the discussion in Appx. 11.1, which establishes that Prop. 1 with an averaged baseline is equivalent to Def. 1, we conclude that the Hessian estimator in Def. 1 also attains the minimum variance among all constant baselines. ◻

Proof. We adopt the decomposition of the variance of the Hessian estimator from 28 , and let \(b\) be the averaged baseline in Thm. 5. Then we have: \[\require{physics} \label{eq:var-decomp-1} \begin{align} \mathrm{Var}\qty(\widehat{{\mathbf{H}}}({\boldsymbol{\theta}})) =& \mathbb{E}\left[ \norm{\frac{1}{K - 1} \sum_{k=1}^{K} \frac{f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_k; \xi) - b}{\mu^2} {\mathbf{u}}_k{\mathbf{u}}_k^{\top} }_F^2 \right] - \norm{\nabla^2 F_{\mu}({\boldsymbol{\theta}})}_F^2 \\ \le& \mathbb{E}\left[ \norm{\frac{1}{K - 1} \sum_{k=1}^{K} \frac{f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_k; \xi) - b}{\mu^2} {\mathbf{u}}_k{\mathbf{u}}_k^{\top} }_F^2 \right] \\ \overset{(a)}{=}& \frac{1}{(K - 1)^2 \mu^4} \sum_{k = 1}^{K} \mathbb{E}\left[ \left|f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_k; \xi) - F_\mu({\boldsymbol{\theta}})\right|^2 \norm{{\mathbf{u}}_k{\mathbf{u}}_k^{\top}}_F^2 \right] \\ =& \frac{K}{(K - 1)^2 \mu^4} \mathbb{E}\left[ \left|f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}; \xi) - F_\mu({\boldsymbol{\theta}})\right|^2 \norm{{\mathbf{u}}}^4 \right] \\ \overset{(b)}{\le}& \frac{K d (d + 2) V}{(K - 1)^2 \mu^4} \;, \end{align}\tag{29}\] where \((a)\) follows from the independence of \(\{{\mathbf{u}}_k\}\) across different queries \(k\), \((b)\) comes from Lem. 7. Here, we denote \(V \triangleq \sigma_\xi^2 + 4 L_0^2 \mu^2 (d + 2)\) for simplicity. ◻

Proof. The total error is bounded as: \[\require{physics} \label{eq:bias-variance-decomp} \begin{align} &\mathbb{E}\qty[\norm{\widehat{{\mathbf{H}}}({\boldsymbol{\theta}}) - \nabla^2 F({\boldsymbol{\theta}})}_F^2] \overset{(a)}{=} \mathbb{E}\qty[\norm{\widehat{{\mathbf{H}}}({\boldsymbol{\theta}}) - \nabla^2 F_{\mu}({\boldsymbol{\theta}})}_F^2] + \norm{\nabla F_\mu({\boldsymbol{\theta}}) - \nabla F({\boldsymbol{\theta}})}_F^2 \overset{(b)}{\le} \frac{K d (d + 2) V}{(K - 1)^2 \mu^4} + L_2^2 \mu^2 d \;, \end{align}\tag{30}\] where \((a)\) follows \(\require{physics} \mathbb{E}\qty[\widehat{{\mathbf{H}}}({\boldsymbol{\theta}}) - \nabla^2 F_{\mu}({\boldsymbol{\theta}})] = 0\) from Thm. 4, and \((b)\) comes from Thm. 10 and Lem. 6. Here we denote \(V \triangleq \sigma_\xi^2 + 4 L_0^2 \mu^2 (d + 2)\) in the last step for simplicity. ◻

Proof. The Hessian estimator in Def. 1 can be rewritten as: \[\widehat{{\mathbf{H}}}({\boldsymbol{\theta}}) = \frac{1}{K - 1} \sum_{k=1}^{K} \frac{f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_k; \xi) - \frac{1}{K} \sum_{k'=1}^{K} f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_{k'}; \xi)}{\mu^2} {\mathbf{u}}_k{\mathbf{u}}_k^{\top} \cdot \mathcal{W}({\boldsymbol{\theta}}+ \mu {\mathbf{u}}_k) = U D U^{\top} \;,\] where we denote \(\require{physics} {\mathbf{D}}= \operatorname{diag}(\nu_{1} / \qty(K - 1), \ldots, \nu_{K} / \qty(K - 1)) \in {\mathbb{R}}^{K \times K}\) with \(\nu_{k} = \frac{f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_{k}; \xi) - \frac{1}{K} \sum_{k'=1}^{K} f({\boldsymbol{\theta}}+ \mu{\mathbf{u}}_{k'}; \xi)}{\mu^2} \mathcal{W}({\boldsymbol{\theta}}+ \mu {\mathbf{u}}_k)\), and \({\mathbf{U}}= [{\mathbf{u}}_1, {\mathbf{u}}_2, \ldots, {\mathbf{u}}_{K}] \in {\mathbb{R}}^{d \times K}\).

According to the Woodbury matrix identity [@woodbury-matrix-identity], we have: \[\require{physics} \qty(\widehat{{\mathbf{H}}}({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d)^{-1} = \qty(\lambda {\mathbf{I}}_d + {\mathbf{U}}{\mathbf{D}}{\mathbf{U}}^{\top})^{-1} = \frac{1}{\lambda} {\mathbf{I}}_d - \frac{1}{\lambda^2} {\mathbf{U}}\qty({\mathbf{D}}^{-1} + \frac{1}{\lambda} {\mathbf{U}}^{\top} {\mathbf{U}})^{-1} {\mathbf{U}}^{\top} \;.\]

If we consider the diagonal approximation of \({\mathbf{U}}\) (i.e. assuming \(\{{\mathbf{u}}_k\}_{k=1}^K\) are orthogonal), then \({\mathbf{U}}^{\top} {\mathbf{U}}\) is a diagonal matrix with the \(k\)-th diagonal element being \(\norm{{\mathbf{u}}_k}^2\). The inverse term can be calculated as: \[\require{physics} {\mathbf{D}}^{-1} + \frac{1}{\lambda} {\mathbf{U}}^{\top} {\mathbf{U}}= \operatorname{diag}\qty(\frac{K - 1}{\nu_{k}} + \frac{1}{\lambda} \norm{{\mathbf{u}}_{k}}^2)_{k=1}^{K} \;.\]

Therefore, the approximated inverse Hessian \(\widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}})\) is: \[\require{physics} \widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) = \frac{1}{\lambda} {\mathbf{I}}_d - \frac{1}{\lambda^2} \sum_{k=1}^{K} \frac{{\mathbf{u}}_{k} {\mathbf{u}}_{k}^{\top}}{\frac{K - 1}{\nu_{k}} + \frac{1}{\lambda} \norm{{\mathbf{u}}_{k}}^2} = \frac{1}{\lambda} {\mathbf{I}}_d - \sum_{k=1}^{K} \frac{\nu_{k}}{\lambda \qty(\lambda (K - 1) + \nu_{k} \norm{{\mathbf{u}}_{k}}^2)} {\mathbf{u}}_{k} {\mathbf{u}}_{k}^{\top} \;.\] ◻

Proof. Adopt from inverse Hessian estimator in 13 , we have: \[\require{physics} \widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) \hat{\nabla} F({\boldsymbol{\theta}}) = \left(\frac{1}{\lambda} {\mathbf{I}}_d - \sum_{k = 1}^{K} \frac{\mu \nu_{k}}{\lambda \qty(\lambda (K - 1) + \nu_{k} \norm{{\mathbf{u}}_{k}}^2)} {\mathbf{u}}_{k} {\mathbf{u}}_{k}^{\top}\right) \hat{\nabla} F \;.\]

The first term is simple: \[\frac{1}{\lambda} {\mathbf{I}}_d \hat{\nabla} F({\boldsymbol{\theta}}) = \frac{1}{\lambda (K - 1)} \sum_{k = 1}^{K} \mu \nu_k {\mathbf{u}}_k \;.\]

While for the second term, to reduce the bias introduced by the correlation between the same random directions used in Hessian and gradient estimation, we remove the random direction that is used in the Hessian estimation from the gradient estimation, i.e., we have: \[\require{physics} \begin{align} &\qty(\sum_{k = 1}^{K} \frac{\mu \nu_{k}}{\lambda \qty(\lambda (K - 1) + \nu_{k} \norm{{\mathbf{u}}_{k}}^2)} {\mathbf{u}}_{k} {\mathbf{u}}_{k}^{\top}) \hat{\nabla} F({\boldsymbol{\theta}}) = \qty(\sum_{k = 1}^{K} \frac{\mu \nu_{k}}{\lambda \qty(\lambda (K - 1) + \nu_{k} \norm{{\mathbf{u}}_{k}}^2)} {\mathbf{u}}_{k} {\mathbf{u}}_{k}^{\top}) \qty(\frac{1}{K - 2} \sum_{\substack{k' = 1 \\ k' \neq k}}^{K} \mu \nu_{k'} {\mathbf{u}}_{k'}) \\ =& \frac{1}{K - 2} \sum_{k = 1}^{K} \sum_{\substack{k' = 1 \\ k' \neq k}}^{K} \frac{\mu \nu_{k} \mu \nu_{k'}}{\lambda \qty(\lambda (K - 1) + \nu_{k} \norm{{\mathbf{u}}_{k}}^2)} {\mathbf{u}}_{k} \qty({\mathbf{u}}_{k}^{\top} {\mathbf{u}}_{k'}) = \sum_{k = 1}^{K} \mu \nu_{k} \frac{{\mathbf{u}}_k^\top \qty(\sum_{k' = 1, k' \neq k}^{K} \nu_{k'} {\mathbf{u}}_{k'}) / \qty(K - 2)}{\lambda \qty(\lambda (K - 1) + \nu_{k} \norm{{\mathbf{u}}_{k}}^2)} {\mathbf{u}}_{k} \;. \end{align}\]

Combining the two terms together, we have: \[\require{physics} \widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) \hat{\nabla} F({\boldsymbol{\theta}}) = \sum_{k = 1}^{K} \mu \nu_{k} \left(\frac{1}{\lambda (K - 1)} - \frac{{\mathbf{u}}_k^\top \qty(\sum_{k' = 1, k' \neq k}^{K} \nu_{k'} {\mathbf{u}}_{k'}) / \qty(K - 2)}{\lambda \qty(\lambda (K - 1) + \nu_{k} \norm{{\mathbf{u}}_{k}}^2)}\right) {\mathbf{u}}_{k} \;.\] ◻

Proof. We first use Resolvent identity to decompose the error between two inverse Hessians into: \[\require{physics} \label{eq:inv-hess-bias-decomp} \begin{align} &\mathbb{E}\qty[\norm{\widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) - \qty(\nabla^2 F_{\mu}({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d)^{-1}}_F^2] \\ \overset{(a)}{=}& \mathbb{E}\qty[\norm{ \widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) \qty(\nabla^2 F_{\mu}({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d - \widetilde{{\mathbf{H}}}({\boldsymbol{\theta}})) \qty(\nabla^2 F_{\mu}({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d)^{-1} }_F^2] \\ \overset{(b)}{\le}& \mathbb{E}\qty[\norm{ \widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) }_2^2 \norm{ \nabla^2 F_{\mu}({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d - \widetilde{{\mathbf{H}}}({\boldsymbol{\theta}}) }_F^2] \norm{(\nabla^2 F_{\mu}({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d)^{-1} }_2^2 \;, \end{align}\tag{31}\] where \((a)\) follows from the Resolvent identity \({\mathbf{A}}^{-1} - {\mathbf{B}}^{-1} = {\mathbf{A}}^{-1} ({\mathbf{B}}- {\mathbf{A}}) {\mathbf{B}}^{-1}\) for any invertible matrices \({\mathbf{A}}\) and \({\mathbf{B}}\), and \((b)\) follows from the Cauchy-Schwarz inequality.

We next bound these three terms in 31 respectively. For the first term, we have: \[\require{physics} \norm{ \widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) }_2^2 = \norm{\frac{1}{\lambda} {\mathbf{I}}_d - \sum_{k=1}^{K} \frac{\nu_{k}}{\lambda \qty(\lambda (K - 1) + \nu_{k} \norm{{\mathbf{u}}_{k}}^2)} {\mathbf{u}}_{k} {\mathbf{u}}_{k}^{\top}}_2^2 \;.\]

Note that the summation \(\require{physics} \sum_{k=1}^{K} \frac{\nu_{k}}{\lambda \qty(\lambda (K - 1) + \nu_{k} \norm{{\mathbf{u}}_{k}}^2)} {\mathbf{u}}_{k} {\mathbf{u}}_{k}^{\top}\) is a rank-\(K\) matrix. Under the orthogonality assumption of random directions \(\{u_k\}_{k=1}^K\), the eigenvalues of \(\widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}})\) related to the subspace spanned by \(\{{\mathbf{u}}_k\}_{k=1}^K\) can be computed as: \[\require{physics} \lambda_k = \frac{1}{\lambda} - \frac{\nu_{k} \norm{{\mathbf{u}}_{k}}^2}{\lambda \qty(\lambda (K - 1) + \nu_{k} \norm{{\mathbf{u}}_{k}}^2)} = \frac{K - 1}{\lambda (K - 1) + \nu_{k} \norm{{\mathbf{u}}_{k}}^2} = \frac{1}{\lambda \qty(1 + t_k)} \;,\] where we denote \(\require{physics} t_k \triangleq \nu_k \norm{{\mathbf{u}}_k}^2 / \qty(\lambda (K - 1))\) for simplicity.

Meanwhile, in the orthogonal complement subspace, the eigenvalues are all \(1 / \lambda\). Therefore, we have: \[\norm{ \widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) }_2^2 = \max \left\{\frac{1}{\lambda^2}, \max_{k = 1, \dots, K} \frac{1}{\lambda^2 (1 + t_k)^2} \right\} = \frac{1}{\lambda^2} \max \left\{1, \max_{k = 1, \dots, K} \frac{1}{(1 + t_k)^2} \right\} \;.\]

To avoid singularity, we assume the existence of a positive constant \(\rho \in (0, 1]\) such that \(|t_k + 1| \ge \rho\) for all \(k\). There are two cases to consider:

• For case \(t_k \ge 0\) or \(t_k \le -2\), we have \(1 / (1 + t_k)^2 \le 1\), thus \(\norm{ \widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) }_2^2 = 1 / \lambda^2\).

• For case \(\rho - 1 \le t_k < 0\) or \(-2 < t_k \le -1 - \rho\), we have \(1 / (1 + t_k)^2 \ge 1\), thus: \[\norm{ \widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) }_2^2 = \max_{k = 1, \dots, K} \frac{1}{\lambda^2 (1 + t_k)^2} \le \frac{1}{\lambda^2 \rho^2} \;.\]

Overall, since \(\rho \in (0, 1]\), we have: \[\norm{ \widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) }_2^2 \le \frac{1}{\lambda^2 \rho^2}\]

Then, the second term in 31 can be bounded by Thm. 6: \[\require{physics} \mathbb{E}\qty[\norm{ \nabla^2 F_{\mu}({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d - \widetilde{{\mathbf{H}}}({\boldsymbol{\theta}}) }_F^2] \le \frac{K d (d + 2) V}{(K - 1)^2 \mu^4} + L_2^2 \mu^2 d\]

For the third term in 31 : \[\require{physics} \label{eq:inv-hess-bias-decomp-3} \begin{align} \norm{(\nabla^2 F({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d)^{-1} }_2^2 \le \norm{(\nabla^2 F({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d)^{-1} }_F^2 = \sum_{i=1}^{d} \frac{1}{\qty(\sigma_i + \lambda)^2} \le \frac{d}{\qty(\sigma_{\min} + \lambda)^2} \;, \end{align}\tag{32}\] where \(\sigma_i\) denotes the \(i\)-th eigenvalue of \(\nabla^2 F({\boldsymbol{\theta}})\), and \(\sigma_{\min}\) denotes the smallest eigenvalue respectively.

Overall, we have: \[\require{physics} \mathbb{E}\qty[\norm{\widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) - \qty(\nabla^2 F_{\mu}({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d)^{-1}}_F^2] \le \frac{d}{\lambda^2 \rho^2 \qty(\sigma_{\min} + \lambda)^2} \qty(\frac{K d (d + 2) V}{(K - 1)^2 \mu^4} + L_2^2 \mu^2 d) \;.\] ◻

Proof. We first decompose the error of the inverse Hessian-gradient product approximation into: \[\require{physics} \label{eq:inv-hess-grad-prod-bias-decomp} \begin{align} &\mathbb{E}\qty[\norm{\widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) \hat{\nabla} F({\boldsymbol{\theta}}) - \qty(\nabla^2 F({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d)^{-1} \nabla F({\boldsymbol{\theta}})}^2] \\ \overset{(a)}{\le}& 2 \underbrace{\mathbb{E}\qty[\norm{\widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) \hat{\nabla} F({\boldsymbol{\theta}}) - \qty(\nabla^2 F({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d)^{-1} \hat{\nabla} F({\boldsymbol{\theta}})}^2]}_{\raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {1}}}} \\ & \quad + 2 \underbrace{\mathbb{E}\qty[\norm{\qty(\nabla^2 F({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d)^{-1} \hat{\nabla} F({\boldsymbol{\theta}}) - \qty(\nabla^2 F({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d)^{-1} \nabla F({\boldsymbol{\theta}})}^2]}_{\raisebox{.5pt}{\textcircled{\raisebox{-.9pt} {2}}}} \;, \end{align}\tag{33}\] where \((a)\) follows from the inequality \((a + b)^2 \le 2 (a^2 + b^2)\).

For the first term in 33 , we have: \[\require{physics} \label{eq:inv-hess-grad-prod-bias-decomp-1} \begin{align} &\mathbb{E}\qty[\norm{\widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) \hat{\nabla} F({\boldsymbol{\theta}}) - \qty(\nabla^2 F({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d)^{-1} \hat{\nabla} F({\boldsymbol{\theta}})}^2] \\ \overset{(a)}{=}& \mathbb{E}\qty[\norm{\qty(\nabla^2 F({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d)^{-1} \qty(\widetilde{{\mathbf{H}}}({\boldsymbol{\theta}}) - \nabla^2 F({\boldsymbol{\theta}}) - \lambda {\mathbf{I}}_d) \widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) \hat{\nabla} F({\boldsymbol{\theta}})}^2] \\ \overset{(b)}{\le}& \mathbb{E}\qty[\norm{\widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) \hat{\nabla} F({\boldsymbol{\theta}})}^2 \norm{\widetilde{{\mathbf{H}}}({\boldsymbol{\theta}}) - \nabla^2 F({\boldsymbol{\theta}}) - \lambda {\mathbf{I}}_d}_F^2] \norm{\qty(\nabla^2 F({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d)^{-1}}_2^2 \;, \end{align}\tag{34}\] where \((a)\) follows from the Resolvent identity \({\mathbf{A}}^{-1} - {\mathbf{B}}^{-1} = {\mathbf{B}}^{-1} ({\mathbf{B}}- {\mathbf{A}}) {\mathbf{A}}^{-1}\) for any invertible matrices \({\mathbf{A}}\) and \({\mathbf{B}}\), and \((b)\) follows from \(\norm{{\mathbf{A}}{\mathbf{v}}}^2 \le \norm{{\mathbf{A}}}_2^2 \norm{{\mathbf{v}}}^2\), \(\norm{{\mathbf{A}}{\mathbf{B}}}_2 \le \norm{{\mathbf{A}}}_2 \norm{{\mathbf{B}}}_2\), and \(\norm{{\mathbf{A}}}_2 \le \norm{{\mathbf{A}}}_F\) for any matrices \({\mathbf{A}}, {\mathbf{B}}\) and vector \({\mathbf{v}}\).

We next bound the three terms in 34 respectively. For the first term, we have: \[\require{physics} \label{eq:inv-hess-grad-prod-norm-1} \norm{\widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) \hat{\nabla} F({\boldsymbol{\theta}})}^2 = \norm{\sum_{k = 1}^{K} \mu \nu_{k} \left(\frac{1}{\lambda (K - 1)} - \frac{{\mathbf{u}}_k^\top \qty(\sum_{k' = 1, k' \neq k}^{K} \nu_{k'} {\mathbf{u}}_{k'}) / \qty(K - 2)}{\lambda \qty(\lambda (K - 1) + \nu_{k} \norm{{\mathbf{u}}_{k}}^2)}\right) {\mathbf{u}}_{k}}^2 \;.\tag{35}\]

Under the assumption that \(\norm{{\mathbf{u}}_k} = 1\), and \(\{{\mathbf{u}}_k\}_{k=1}^K\) are orthogonal, the second term in 35 reduces to zero, and we then have: \[\label{eq:inv-hess-grad-prod-norm-2} \norm{\widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) \hat{\nabla} F({\boldsymbol{\theta}})}^2 = \norm{\sum_{k = 1}^{K} \frac{\mu \nu_{k}}{\lambda (K - 1)} {\mathbf{u}}_{k}}^2 \overset{(a)}{=} \frac{\mu^2}{\lambda^2 (K - 1)^2}\sum_{k = 1}^{K} \nu_{k}^2 \norm{{\mathbf{u}}_k}^2 = \frac{K \mu^2}{\lambda^2 (K - 1)^2} \nu^2 \norm{{\mathbf{u}}}^2 \;,\tag{36}\] where \((a)\) follows from the independence of \(\{{\mathbf{u}}_k\}_{k=1}^K\).

Next, the second term in 34 can be bounded by: \[\label{eq:inv-hess-grad-prod-norm-3} \norm{\widetilde{{\mathbf{H}}}({\boldsymbol{\theta}}) - \nabla^2 F({\boldsymbol{\theta}}) - \lambda {\mathbf{I}}_d }_F^2 \overset{(a)}{\le} \norm{\widehat{{\mathbf{H}}}({\boldsymbol{\theta}}) - \nabla^2 F_{\mu}({\boldsymbol{\theta}})}_F^2 + L_2^2 \mu^2 d \overset{(b)}{\le} \frac{K}{(K - 1)^2} \nu^2 \norm{{\mathbf{u}}}^4 + L_2^2 \mu^2 d \;.\tag{37}\] where \((a)\) results from 30 , and \((b)\) follows from 29 .

Combine 36 and 37 , we have: \[\require{physics} \label{eq:inv-hess-grad-prod-norm-4} \begin{align} &\mathbb{E}\qty[\norm{\widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) \hat{\nabla} F({\boldsymbol{\theta}})}^2 \norm{\widetilde{{\mathbf{H}}}({\boldsymbol{\theta}}) - \nabla^2 F({\boldsymbol{\theta}}) - \lambda {\mathbf{I}}_d}_F^2] \le \frac{K \mu^2}{\lambda^2 (K - 1)^2} \mathbb{E}\qty[\nu^2 \norm{{\mathbf{u}}}^2 \left(\frac{K}{(K - 1)^2} \nu^2 \norm{{\mathbf{u}}}^4 + L_2^2 \mu^2 d\right)] \\ =& \frac{K \mu^2}{\lambda^2 (K - 1)^2} \qty(\frac{K}{(K - 1)^2} \mathbb{E}[\nu^4 \norm{{\mathbf{u}}}^6] + L_2^2 \mu^2 d \mathbb{E}[\nu^2 \norm{{\mathbf{u}}}^2]) \;. \end{align}\tag{38}\]

Given Lem. 7 and the following inequality: \[\require{physics} \mathbb{E}[\nu^4 \norm{{\mathbf{u}}}^6] \overset{(a)}{\le} \mathbb{E}[\nu^2 \norm{{\mathbf{u}}}^6] \cdot \max \left|\nu^2\right| \overset{(b)}{\le} \qty(\frac{2 G}{\mu^2})^2 \frac{d (d + 2) (d + 4)}{\mu^4} \left(\sigma_\xi^2 + 4 L_0^2 \mu^2 (d + 3)\right) \;,\] where \((a)\) follows from Holder’s inequality \(\mathbb{E}[|XY|] \le \mathbb{E}[|X|^p]^{1/p} \mathbb{E}[|Y|^q]^{1/q}\) for \(p = 1\) and \(q = \infty\), and \((b)\) comes from Lem. 7 and the bound \(\max |\nu| = \max |f({\boldsymbol{\theta}}+ \mu {\mathbf{u}}; \xi) - F_\mu({\boldsymbol{\theta}})| / \mu^2 = 2 G / \mu^2\).

Therefore, 38 can be further bounded by: \[\require{physics} \mathbb{E}\qty[\norm{\widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) \hat{\nabla} F({\boldsymbol{\theta}})}^2 \norm{\widetilde{{\mathbf{H}}}({\boldsymbol{\theta}}) - \nabla^2 F({\boldsymbol{\theta}}) - \lambda {\mathbf{I}}_d}_F^2] \le \frac{K d}{\lambda^2 (K - 1)^2}\qty(B_1 \sigma_\xi^2 + 4 B_2 L_0^2 \mu^2) \;,\] where we denote \(\require{physics} B \triangleq \frac{4 K G^2}{\qty(K - 1)^2 \mu^6} \qty(d + 2) (d + 4)\), \(B_1 \triangleq \frac{B}{\mu^6} + L_2^2 d\), and \(\require{physics} B_2 \triangleq \frac{B}{\mu^6} \qty(d + 3) + L_2^2 d (d + 1)\) for simplicity.

Adopting the results from 32 , the third term in 34 can be bounded as: \[\require{physics} \norm{(\nabla^2 F({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d)^{-1} }_2^2 \le \norm{(\nabla^2 F({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d)^{-1} }_F^2 = \sum_{i=1}^{d} \frac{1}{\qty(\sigma_i + \lambda)^2} \le \frac{d}{\qty(\sigma_{\min} + \lambda)^2} \;,\] where \(\sigma_i\) denotes the \(i\)-th eigenvalue of \(\nabla^2 F({\boldsymbol{\theta}})\), and \(\sigma_{\min}\) denotes the smallest eigenvalue respectively.

Overall, we have: \[\require{physics} \mathbb{E}\qty[\norm{\widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) \hat{\nabla} F({\boldsymbol{\theta}}) - \qty(\nabla^2 F({\boldsymbol{\theta}}) + \lambda {\mathbf{I}}_d)^{-1} \hat{\nabla} F({\boldsymbol{\theta}})}^2] \le \frac{K d^2 \qty(B_1 \sigma_\xi^2 + 4 B_2 L_0^2 \mu^2)}{\lambda^2 \qty(\sigma_{\min} + \lambda)^2 (K - 1)^2} \;.\] ◻

Proof. Due to the \(L_1\)-Lipschitz smoothness of \(F({\boldsymbol{\theta}})\) given in Lem. 4, the smoothed objective \(F_\mu({\boldsymbol{\theta}})\) is also \(L_1\)-Lipschitz smooth: \[\norm{\nabla F_\mu({\boldsymbol{\theta}}) - \nabla F_\mu({\boldsymbol{\theta}}')} \overset{(a)}{\le} \mathbb{E}_{{\mathbf{u}}} \left[\norm{\nabla F({\boldsymbol{\theta}}+ \mu {\mathbf{u}}) - \nabla F({\boldsymbol{\theta}}' + \mu {\mathbf{u}})}\right] \le L_1 \norm{{\boldsymbol{\theta}}- {\boldsymbol{\theta}}'} \;,\] where \((a)\) follows from Jensen’s inequality.

Therefore, we have: \[\label{eq:F95mu-smoothness} F_\mu({\boldsymbol{\theta}}_{t+1}) - F_\mu({\boldsymbol{\theta}}_t) \le \left\langle\nabla F_\mu({\boldsymbol{\theta}}_t), {\boldsymbol{\theta}}_{t+1} - {\boldsymbol{\theta}}_t \right\rangle + \frac{L_1}{2} \norm{{\boldsymbol{\theta}}_{t+1} - {\boldsymbol{\theta}}_t}^2 \;.\tag{39}\]

For the first term in the RHS of 39 , we have: \[\require{physics} \begin{align} &\mathbb{E}\qty[\left\langle\nabla F_\mu({\boldsymbol{\theta}}_t), {\boldsymbol{\theta}}_{t+1} - {\boldsymbol{\theta}}_t \right\rangle] = - \eta \mathbb{E}\qty[\left\langle\nabla F_\mu({\boldsymbol{\theta}}_t), \widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}_t) \hat{\nabla} F({\boldsymbol{\theta}}_t) \right\rangle] \\ \overset{(a)}{=}& - \eta \mathbb{E}\qty[\left\langle\nabla F_\mu({\boldsymbol{\theta}}_t), \sum_{k = 1}^{K} \mu \nu_{k} \left(\frac{1}{\lambda (K - 1)} - \frac{{\mathbf{u}}_k^\top \qty(\sum_{k' = 1, k' \neq k}^{K} \nu_{k'} {\mathbf{u}}_{k'}) / \qty(K - 2)}{\lambda \qty(\lambda (K - 1) + \nu_{k} \norm{{\mathbf{u}}_{k}}^2)}\right) {\mathbf{u}}_{k} \right\rangle] \\ \overset{(b)}{=}& - \eta \mathbb{E}\qty[\left\langle\nabla F_\mu({\boldsymbol{\theta}}_t), \sum_{k = 1}^{K} \frac{\mu \nu_{k}}{\lambda (K - 1)} {\mathbf{u}}_{k} \right\rangle] = - \frac{\eta}{\lambda} \mathbb{E}\qty[\left\langle\nabla F_\mu({\boldsymbol{\theta}}_t), \hat{\nabla} F({\boldsymbol{\theta}}_t) \right\rangle] \overset{(c)}{=} - \frac{\eta}{\lambda} \norm{\nabla F_\mu({\boldsymbol{\theta}}_t)}^2 \;, \end{align}\] where \((a)\) follows from the definition of \(\widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}_t)\), \((b)\) follows from the orthogonality of \(\{{\mathbf{u}}_k\}_{k=1}^K\), and \((c)\) follows from Lem. 5.

For the second term in the RHS of 39 , we have: \[\require{physics} \frac{L_1}{2} \mathbb{E}\qty[\norm{{\boldsymbol{\theta}}_{t+1} - {\boldsymbol{\theta}}_t}^2] = \frac{L_1 \eta^2}{2} \mathbb{E}\qty[\norm{\widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}_t) \hat{\nabla} F({\boldsymbol{\theta}}_t)}^2] \overset{(a)}{\le} \frac{L_1 \eta^2 K d^2 \qty(B_1 \sigma_\xi^2 + 4 B_2 L_0^2 \mu^2)}{2 \lambda^2 \qty(\sigma_{\min} + \lambda)^2 (K - 1)^2} \;,\] where \((a)\) follows from the result in Thm. 8.

Summing all iterations from \(0\) to \(T - 1\) to both sides of 39 , and denote \(\Delta = F_\mu({\boldsymbol{\theta}}_0) - F_\mu({\boldsymbol{\theta}}^*)\), we have: \[\require{physics} \frac{1}{T} \sum_{t=0}^{T-1} \mathbb{E}\qty[\norm{\nabla F_\mu({\boldsymbol{\theta}}_t)}^2] \le \frac{\lambda \Delta}{\eta T} + \frac{L_1 \eta K d^2 \qty(B_1 \sigma_\xi^2 + 4 B_2 L_0^2 \mu^2)}{2 \lambda \qty(\sigma_{\min} + \lambda)^2 (K - 1)^2} \;.\]

To simplify the expression, we choose \(\require{physics} \eta \sim {\mathcal{O}}\qty(\frac{4 \lambda \qty(\sigma_{\min} + \lambda)^2 (K - 1)^2 \epsilon^2}{L_1 K d^2 \qty(B_1 \sigma_\xi^2 + 4 B_2 L_0^2 \mu^2)}) \sim {\mathcal{O}}\qty(\epsilon^2)\), and \(\require{physics} T \sim {\mathcal{O}}\qty(\frac{2 \lambda \Delta}{\eta \epsilon^2}) \sim {\mathcal{O}}\qty(\epsilon^{-4})\), such that: \[\require{physics} \frac{1}{T} \sum_{t=0}^{T-1} \mathbb{E}\qty[\norm{\nabla F_\mu({\boldsymbol{\theta}}_t)}^2] \le \frac{1}{2} \epsilon^2 + \frac{1}{2} \epsilon^2 = \epsilon^2 \;.\]

Finally, since \(\norm{\nabla F({\boldsymbol{\theta}}_t) - \nabla F_\mu({\boldsymbol{\theta}}_t)} \le \frac{\mu L_1 d}{2}\), we have: \[\require{physics} \begin{align} &\frac{1}{T} \sum_{t=0}^{T-1} \mathbb{E}\qty[\norm{\nabla F({\boldsymbol{\theta}}_t)}] \le \frac{1}{T} \sum_{t=0}^{T-1} \mathbb{E}\qty[\norm{\nabla F_\mu({\boldsymbol{\theta}}_t)}] + \frac{1}{T} \sum_{t=0}^{T-1} \mathbb{E}\qty[\norm{\nabla F({\boldsymbol{\theta}}_t) - \nabla F_\mu({\boldsymbol{\theta}}_t)}] \\ \overset{(a)}{\le}& \sqrt{\frac{1}{T} \sum_{t=0}^{T-1} \mathbb{E}\qty[\norm{\nabla F_\mu({\boldsymbol{\theta}}_t)}^2]} + \frac{\mu L_1 d}{2} \le \epsilon + \frac{\mu L_1 d}{2} \;, \end{align}\] where \((a)\) follows from Jensen’s inequality. ◻

Proof. We adopt the decomposition of the variance of the Hessian estimator with query reuse from 9 , and let \(b\) be the averaged baseline in Thm. 5. Then we have: \[\require{physics} \label{eq:var-decomp-1-query-reuse} \begin{align} \mathrm{Var}\qty(\widehat{{\mathbf{H}}}({\boldsymbol{\theta}}_{t-1})) =& \mathbb{E}\left[ \norm{\frac{1}{N K - 1} \sum_{n,k=1}^{N,K} \frac{f({\boldsymbol{\theta}}_{t-n} + \mu{\mathbf{u}}_{t-n,k}; \xi) - b}{\mu^2} {\mathbf{u}}_{t-n,k}{\mathbf{u}}_{t-n,k}^{\top} }_F^2 \right] - \norm{\nabla^2 F_{\mu}({\boldsymbol{\theta}}_{t-1})}_F^2 \\ \le& \mathbb{E}\left[ \norm{\frac{1}{N K - 1} \sum_{n,k=1}^{N,K} \frac{f({\boldsymbol{\theta}}_{t-n} + \mu{\mathbf{u}}_{t-n,k}; \xi) - b}{\mu^2} {\mathbf{u}}_{t-n,k}{\mathbf{u}}_{t-n,k}^{\top} }_F^2 \right] \\ \overset{(a)}{=}& \frac{1}{(N K - 1)^2 \mu^4} \sum_{n,k = 1}^{N,K} \mathbb{E}\left[ \left|f({\boldsymbol{\theta}}_{t-n} + \mu{\mathbf{u}}_{t-n,k}; \xi) - \frac{1}{N} \sum_{n'=1}^N F_\mu({\boldsymbol{\theta}}_{t-n'})\right|^2 \norm{{\mathbf{u}}_{t-n,k}{\mathbf{u}}_{t-n,k}^{\top}}_F^2 \right] \\ =& \frac{K}{(NK - 1)^2 \mu^4} \sum_{n=1}^N \mathbb{E}\left[ \left|f({\boldsymbol{\theta}}_{t-n} + \mu{\mathbf{u}}_{t-n,k}; \xi) - \frac{1}{N} \sum_{n'=1}^N F_\mu({\boldsymbol{\theta}}_{t-n'})\right|^2 \norm{{\mathbf{u}}_{t-n,k}{\mathbf{u}}_{t-n,k}^{\top}}_F^2 \right] \\ \overset{(b)}{\le}& \frac{N K d (d + 2) \lambda^2 V}{\lambda^2 (N K - 1)^2 \mu^4 - L_0^2 \eta^2 \mu^2 K (N^2 - 1) d (d + 2) / 3} \;, \end{align}\tag{40}\] where \((a)\) follows from the independence of \(\{{\mathbf{u}}_k\}\) across different queries \(k\), \((b)\) comes from Lem. 9. Here, we denote \(V \triangleq \sigma_\xi^2 + 4 L_0^2 \mu^2 (d + 2)\) for simplicity. ◻

Proof. The total error is bounded as: \[\require{physics} \label{eq:bias-variance-decomp-query-reuse} \begin{align} \mathbb{E}\qty[\norm{\widehat{{\mathbf{H}}}({\boldsymbol{\theta}}) - \nabla^2 F({\boldsymbol{\theta}})}_F^2] \overset{(a)}{=}& \mathbb{E}\qty[\norm{\widehat{{\mathbf{H}}}({\boldsymbol{\theta}}) - \nabla^2 F_{\mu}({\boldsymbol{\theta}})}_F^2] + \norm{\nabla F_\mu({\boldsymbol{\theta}}) - \nabla F({\boldsymbol{\theta}})}_F^2 \\ \overset{(b)}{\le}& \frac{N K d (d + 2) \lambda^2 V}{\lambda^2 (N K - 1)^2 \mu^4 - L_0^2 \eta^2 \mu^2 K (N^2 - 1) d (d + 2) / 3} + L_2^2 \mu^2 d \;, \end{align}\tag{41}\] where \((a)\) follows \(\require{physics} \mathbb{E}\qty[\widehat{{\mathbf{H}}}({\boldsymbol{\theta}}) - \nabla^2 F_{\mu}({\boldsymbol{\theta}})] = 0\) from Thm. 4, and \((b)\) comes from Thm. 11 and Lem. 6. Here we denote \(V \triangleq \sigma_\xi^2 + 4 L_0^2 \mu^2 (d + 2)\) in the last step for simplicity. ◻

Proof. From 12 and 13 , the only difference between the exact inverse and the approximation is whether \({\mathbf{M}}^{-1}\) or \(\widetilde{{\mathbf{M}}}^{-1}\) is used in the low-dimensional Woodbury correction: \[\require{physics} \qty(\widehat{{\mathbf{H}}}({\boldsymbol{\theta}})+\lambda {\mathbf{I}}_d)^{-1}-\widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}}) = -\frac{1}{\lambda^2}{\mathbf{U}}\qty({\mathbf{M}}^{-1}-\widetilde{{\mathbf{M}}}^{-1}){\mathbf{U}}^\top \;.\] For Gaussian directions, standard concentration gives \(\norm{{\mathbf{U}}}_2^2=\norm{{\mathbf{G}}}_2=\mathcal{O}_p(d)\), \(\norm{\Delta}_2=\Theta_p(d)\), and \(\norm{{\mathbf{E}}}_F=\mathcal{O}_p(K\sqrt d)\), since the off-diagonal entries \({\mathbf{u}}_i^\top{\mathbf{u}}_j\) have variance \(d\) for \(i\neq j\). Because \(\norm{{\mathbf{D}}^{-1}}_2=\mathcal{O}_p(1)\) and \(\lambda\) is fixed, the diagonal matrix \(\widetilde{{\mathbf{M}}}\) has entries of order \(d\), hence \(\norm{\widetilde{{\mathbf{M}}}^{-1}}_2=\mathcal{O}_p(d^{-1})\).

By the assumed stability of the exact Woodbury inverse, \(\norm{{\mathbf{M}}^{-1}}_2=\mathcal{O}_p(d^{-1})\). Then, by the resolvent identity, \[{\mathbf{M}}^{-1}-\widetilde{{\mathbf{M}}}^{-1} = -\frac{1}{\lambda}{\mathbf{M}}^{-1}{\mathbf{E}}\widetilde{{\mathbf{M}}}^{-1}.\] Therefore, \[\require{physics} \begin{align} &\norm{\qty(\widehat{{\mathbf{H}}}({\boldsymbol{\theta}})+\lambda {\mathbf{I}}_d)^{-1}-\widetilde{{\mathbf{H}}}^{-1}({\boldsymbol{\theta}})}_F \le \frac{1}{\lambda^3}\norm{{\mathbf{U}}}_2^2\norm{{\mathbf{M}}^{-1}}_2\norm{{\mathbf{E}}}_F\norm{\widetilde{{\mathbf{M}}}^{-1}}_2 \\ =&\;\mathcal{O}_p(d)\cdot \mathcal{O}_p(d^{-1})\cdot \mathcal{O}_p(K\sqrt d)\cdot \mathcal{O}_p(d^{-1}) = \mathcal{O}_p\qty(\frac{K}{\sqrt d}) \;, \end{align}\] which completes the proof. ◻