Variational quantum algorithms rely on the optimization of parameterized quantum circuits in noisy settings. The commonly used back-propagation procedure in classical machine learning is not directly applicable in this setting due to the collapse of quantum states after measurements. Thus, gradient estimations constitute a significant overhead in a gradient-based optimization of such quantum circuits. This paper introduces a random coordinate descent algorithm as a practical and easy-to-implement alternative to the full gradient descent algorithm. This algorithm only requires one partial derivative at each iteration.
Motivated by the behavior of measurement noise in the practical optimization of parameterized quantum circuits, this paper presents an optimization problem setting that is amenable to analysis. Under this setting,
the random coordinate descent algorithm exhibits the same level of stochastic stability as the full gradient approach, making it as resilient to noise. The complexity of the random coordinate descent method is generally no worse than that of the gradient descent and can be much better for various quantum optimization problems with anisotropic Lipschitz constants. Theoretical analysis and extensive numerical experiments validate our findings.

§ INTRODUCTION

Variational quantum algorithms have emerged as a promising application for near-term quantum devices, addressing various computational challenges with enhanced efficiency [46, 15]. These algorithms encompass several notable approaches, such as the variational quantum eigensolver [54], the quantum approximate optimization algorithm [22, 15, 40], and quantum machine learning [38, 8, 63, 14]. They are designed to operate in a hybrid quantum-classical fashion [43, 21].
In these algorithms, the quantum component involves the implementation of parameterized quantum gate operations. By performing measurements, a cost function (and optionally, its gradient) is obtained as the output. The classical computational procedure then utilizes an iterative method to produce updates for the parameters, which are subsequently leveraged to refine and reprogram the quantum circuits. This iterative process continues until convergence is achieved, forming a feedback loop that continues to improve the algorithm's performance.

In variational quantum algorithms, the optimizable parameters are defined within parameterized quantum circuits (PQCs) [7, 52, 64]. A PQC is a sequence of unitary operators represented by parameterized quantum gates that can be readily implemented on a quantum computer. Assuming we are working in an $n$-qubit Hilbert space, a parameterized quantum circuit can be expressed as follows:

\begin{equation}\label{eq: pqc}
U(\bm \theta) = \prod_{j=1}^J U_j(\bm \theta_j) W_j.
\end{equation}

Here, ${\bm \theta} \!= \!\{\bm \theta_j\}_{j=1}^J$ are the parameters that we need to optimize, $U_j(\bm \theta_j)\in \mathbb{C}^{2^n \times 2^n}$ are the parameterized unitary operators, and $W_j \in \mathbb{C}^{2^n \times 2^n}$ are fixed unitary operators. For instance, a simple example of a PQC consisting only of one-qubit Pauli rotation operators takes the form

\[ U_j(\bm \theta_j) = \bigotimes_{m=1}^{M} e^{- i\theta_{j,k_{j,m}} \sigma_{j,k_{j,m}}},  \]

where $\sigma_{j,k_{j,m}}\in\mathbb{C}^{2\times 2}$ is a single-qubit Pauli matrix that acts on $k_{j,m}$-th qubit, $\theta_{j,k_{j,m}}$ represents one of the parameters in $\bm \theta$, and $W_j$'s can be used to represent quantum gates that do not require parameterization, such as the controlled-NOT (CNOT) gate.

Let $d$ be the dimension of the parameters, and we write $\bm\theta=(\theta_1,\theta_2,\cdots,\theta_d)$.
We then optimize the parameter $\bm \theta$ by minimizing a properly chosen cost function $f(\bm \theta)$.
As an example, the variation quantum eigensolvers (VQE) finds the smallest eigenvalue (ground-state energy) and its corresponding eigenvector (ground state) of a given Hamiltonian matrix $H$  by minimizing the energy of the state:

\begin{equation}\label{eq:vqe}
\theta^*=\mathrm{argmin}_{\bm\theta}  f(\bm \theta)=\mathrm{argmin}_{\bm\theta}\bra{U(\bm \theta) \psi_0} H \ket{U(\bm \theta) \psi_0}\,.
\end{equation}

Here, $\ket{\psi_0}\in\mathbb{C}^{2^n}$ is a predetermined initial state that can be easily prepared on a quantum computer. For each given $\bm{\theta}$, $U(\bm{\theta})$ is implemented on a quantum computer to evolve $\ket{\psi_0}$, and the corresponding energy $f(\bm{\theta})$ and its gradient $\nabla_{\bm \theta} f(\bm{\theta})$ can be  consequently obtained with measurements. By solving the optimization problem (<ref>), the minimum value gives an approximation to the smallest eigenvalue of $H$, while $U(\bm{\theta}^*)\ket{\psi_0}$ approximates the corresponding eigenvector.

§.§ Problem setup

Although the problem of optimizing parameters in VQAs resembles classical optimization problems in machine learning, there exist key differences, particularly in how the cost function is evaluated and the level of accuracy that can be obtained for function and gradient evaluations. Firstly, quantum circuits used for estimating partial derivatives in various directions are typically different. This is predominantly because there is no straightforward method (in parallel to backpropagation) to estimate the entire gradient at once, given the inherent nature of quantum states. The predominant method for computing partial derivatives in a PQC is called the parameter-shift rule [18, 71, 5], which can only be applied to evaluate one component of the partial derivatives at a time.
As a result, the estimation of the gradient, $\nabla f(\bm\theta)$, typically incurs a cost that is $d$ times greater than the cost associated with merely estimating a single partial derivative, $\partial_i f(\bm\theta)$.

Secondly, the evaluation of any given quantity, a function value or a partial derivative, requires measurement from quantum computers and is subject to measurement noise. We note that this noise is associated with a finite sampling space. For example, a measurement of the Hamiltonian in (<ref>), which is defined in a finite-dimensional Hilbert space, yields one of its eigenvalues corresponding to the ansatz. Thus, with an increased number of samples or measurements, the central limit theorem suggests that the distribution of the sample average of the function value or the partial derivative can be approximated by a Gaussian distribution, and as a result, the accuracy of function and gradient evaluations can be relatively low. Therefore, the optimization algorithm must be designed to be resilient to measure noise.

In an idealized scenario, we may assume that both the function value and the partial derivatives incorporated into the optimization routine are subject to some Gaussian noise. But the magnitude of corresponding noises can differ up to a constant, especially in situations where the parameter shift rule is applicable (see [68]). With this consideration, the problem of optimizing PQCs can be stated as follows:

[Optimizing parameterized quantum circuits] Finding an efficient algorithm to solve the optimization problem,

\begin{equation}
{\bm \theta}^*=\mathrm{argmin}_{\bm \theta\in\mathbb{R}^d}f(\bm \theta),
\end{equation}

under the following assumptions:

* The cost of evaluating a partial derivative scales linearly with that of a function value.

* Every evaluation of the function and partial derivative is susceptible to Gaussian noise:

\begin{equation}\label{eqn:g_i}
f(\bm \theta) \to f(\bm \theta)+N(0,\sigma^2_1(\bm\theta))\,, \quad
\partial_i f(\bm \theta) \to \partial_i f(\bm \theta)+N(0,\sigma^2_2(\bm\theta))\,.
\end{equation}

Here, $\sigma_1(\bm\theta)$ and $\sigma_2(\bm\theta)$ depend on the real implementation and are not necessarily the same (see [68] for example). For simplicity, in our later analysis, we assume that $\sigma_2(\bm\theta)$ has a uniform upper bound $\sigma_\infty$ (see <ref>).

§.§ Optimization methods

One widely used approach for optimizing VQA is through the application of gradient descent (GD) [67, 79]. The classical gradient descent method involves iteratively updating the parameters $\bm \theta$ by utilizing the gradient of the cost function.

\begin{equation}\label{eq:def_GD}
\bs{\theta}_{n+1} = \bs{\theta}_n - a_n\nabla f(\bs{\theta}_n),
\end{equation}

where $a_n$ denotes the learning rate. In light of the measurement process in quantum computing, we consider the noisy gradient descent: Rather than implementing <ref> with exact $\nabla f(\bs{\theta}_n)$, we apply an unbiased estimator $g(\bs{\theta})$ [$g(\bs{\theta})$ satisfies $\mathbb E[g(\bs{\theta})] = \nabla f(\bs{\theta})$] (for example, (<ref>)). Consequently,  the parameter update  involves the following iteration,

\begin{equation}\label{eq:def_noisy_GD}
\bs{\theta}_{n+1} = \bs{\theta}_n - a_n g(\bs{\theta}_n).
\end{equation}

Since $g(\bs{\theta}_n)$ is an unbiased estimation, Eq. (<ref>) is equivalent to Eq. (<ref>) in the expectation sense. Specifically, by taking the conditional expectation on both sides, we have

\begin{equation}\label{eqn:expect_gd}
\mathbb{E}(\bs{\theta}_{n+1}|\bs{\theta_n})=\bs{\theta}_n - a_n\nabla f(\bs{\theta}_n)
\end{equation}

where $\mathbb{E}(\cdot|\bs{\theta_n})$ denotes the conditional expectation given $\bs{\theta}_n$.

While noisy gradient descent avoids the need for precise gradient information, it still requires the approximated full gradient information at each iteration. As argued before, in the context of VQA, it is often necessary to compute $d$ partial derivatives separately for each direction, which makes the cost of each updating step at least $d$. In this paper, we introduce an alternative optimization method called random coordinate descent (RCD) [73, 48, 60] for addressing <ref>, with the goal of eliminating the cost dependency on $d$ in each step. RCD can be viewed as a variant of gradient descent (GD) where the full gradient in GD is approximated by a randomly selected component of $\nabla f(\bs{\theta}_n)$ in each iteration. Specifically, one RCD iteration can be expressed as:

\begin{equation}\label{eq: def_RCD}
\bs{\theta}_{n+1} = \bs{\theta}_n - a_n\bs{e}_{i_n}\partial_{i_n}f(\bs{\theta}_n)\,.
\end{equation}
Here $\bs{e}_{i_n}$ is the $i_n$-th unit direction,  $f_{i_n}'(\bs{\theta}_n)$ is the corresponding partial derivative of the cost function, and  $i_n$ is  a random index uniformly drawn from $\{1,2,\cdots, d\}$. Similar to Eq. (<ref>),
we can write the noisy RCD as:

\begin{equation}\label{eq:def_noisy_RGD}
\bs{\theta}_{n+1} = \bs{\theta}_n - a_n\bs{e}_{i_n}g_{i_n}(\bs{\theta}_n)\,.
\end{equation}

It is important to emphasize that in each iteration of RCD (<ref>), only one partial derivative information is needed. Consequently, within the scope of VQA (as stated in the first assumption of <ref>), the cost per step of RCD is $d$ times smaller than that of GD.

§.§ Contribution

This paper primarily focuses on the investigation of RCD in the context of noisy gradient evaluation. Our analysis is conducted in a specific comparison with GD, and we illustrate that, under specific conditions, RCD can serve as a favorable alternative for optimizing parameterized quantum circuits. The main contributions of this study can be summarized as follows:

* We show that RCD is theoretically no worse than GD when measuring the complexity by the number of partial derivative calculations (Theorems <ref> and <ref>), assuming the presence of noise and the local PL condition.  A summary of the complexities of the two methods is presented in Table <ref> for comparison.  It is important to highlight that the inequality $L_{\mathrm{avg}}\leq L\leq dL_{\mathrm{avg}}$ always holds. Consequently, when the optimization problem is highly anisotropic, i.e., $L\gg L_{\mathrm{avg}}$, RCD is more cost-effective than GD. In the most extreme case when $L$ is nearly equal to $dL_{\mathrm{avg}}$, RCD can reduce the complexity by a factor of $d$ compared to GD.

* We demonstrate that (noisy) GD and RCD converge with high probability under the local PL condition (Assumption <ref>) and are stable under noisy gradient information. Specifically, if the initial parameter $\bs{\theta}_0$ resides within the basin $\mathcal{N}(\mathcal{X})$ surrounding the global minimum, both noisy methods ensure that the subsequent parameters $\bs{\theta}_n$ will remain consistently within this basin until they converge with the same high probability (Lemmas <ref> and <ref>). To the best of the our knowledge, such stochastic stability has not been established for the optimization methods in variational quantum algorithms.

* We provide extensive empirical evidence demonstrating that RCD consistently delivers superior performance compared to GD (Sections <ref> and <ref>). Our numerical findings support the theoretical observation that RCD can take a larger learning rate than GD, leading to faster convergence.

Algorithm     Iteration cost     Iterations to reach $\epsilon$ tolerance     Total cost
GD      $\Omega(d)$     $\tilde{\mathcal{O}}\left(\max\left\{\frac{L\sigma^2_\infty d}{\mu^2\epsilon},\frac{L}{\mu}\log\left(\frac{1}{\epsilon}\right)\right\}\right)$      $\Omega(\frac{L\sigma^2_\infty d^2}{\epsilon})$

RCD      $\Omega(1)$     $\tilde{\mathcal{O}}\left(\max\left\{\frac{L_{\mathrm{avg}}\sigma^2_\infty d^2}{\mu^2\epsilon},\frac{dL_{\max}}{\mu}\log\left(\frac{1}{\epsilon}\right)\right\}\right)$      $\Omega(\frac{L_{\mathrm{avg}}\sigma^2_\infty d^2}{\epsilon})$

Comparison of the gradient descent and the randomized coordinate descent methods with an unbiased noisy gradient estimation. $d$ is the dimension of the parameter, and the smoothness constants $L$ and $L_{\mathrm{avg}}$ are defined in (<ref>) and (<ref>), respectively. $\sigma^2_\infty$ is a bound for the measurement noise defined in (<ref>). In the table, we limit our attention to the situation where the learning rate is fixed.

§.§ Related works

Gradient descent with noise
The noisy gradient descent (<ref>) is a popular optimization method in the classical machine learning community. Notable examples are the stochastic gradient descent (SGD) [9] or the perturbed gradient descent (PGD)  [30]. The convergence properties of the noisy gradient descent method in (<ref>) have been extensively studied  [47, 50, 58, 51, 68, 41]. For classical machine learning, these previous works except [68] established that when the cost function is $L$ smooth, $\mu$ strong convex (or Polyak-ojasiewicz condition (PL) [56]) and satisfies an additional condition, $f(\bs{\theta}_n)$ converges linearly to an approximation of $f_{\min}$. In the recent work [68], a similar theoretical result was shown for the noisy GD method applied to quantum optimization problems.

Randomized coordinate descent
The randomized coordinate descent method (RCD) has proven its efficiency over GD in many large-scale optimization problems. The convergence properties of RCD have been extensively explored in the fields of machine learning and optimization [73, 48, 60, 49, 39, 17]. For example, it was shown in [48] that when $f$ is strongly convex, the convergence complexity of RCD can be consistently lower than or equal to that of GD. Here, complexity refers to the total number of partial derivative calculations required for convergence. Later, for strongly convex functions, RCD accelerations were achieved with adaptive momentum-based strategies in various regimes [39, 49]. For the non-convex optimization, recent work [17] shows the global convergence behavior of RCD with a focus on saddle point avoidance.
convergence rates of RCD have been scarcely studied for nonconvex optimization problems. More importantly, most related works focused on the case where partial derivatives are computed exactly, while in this work, we deal with the case where partial derivatives are estimated, which is subject to noise, and we will refer to it as noisy RCD (<ref>).

Locally-defined convex conditions for convergence analysis
One limitation of the conventional convergence analysis is its reliance on assumptions of global convex [9] or global PL [68] conditions for the cost function $f(\bm \theta)$. However, we show that such global assumptions are not satisfied in quantum problem applications with PQCs, as elaborated in <ref>.
Thus, one must weaken such a global assumption to a local one in the analysis. Convergence analysis under local assumptions requires more sophisticated techniques (see [23, 34, 53, 45] and therein), but it provides important insights that help to interpret empirical results. In our work, we make a local non-convex condition based on the local PL condition [41]. Under this condition and suitable assumptions for the cost function, we

By employing a stochastic stability argument, we demonstrate that the noisy GD and RCD methods maintain a comparable convergence rate under our local PL condition with high probability (refer to <ref> and <ref>). To the best of the authors' knowledge, this paper is the first to provide a rigorous result for the complexity of noisy GD and RCD under a local PL condition designed for variational quantum algorithms built from PQCs.

Other quantum optimization methods
Another promising direction of research in variational quantum algorithms is zero-order optimization, more commonly known as gradient-free methods. Notably, policy gradient-based techniques have shown their effectiveness in noise robust optimization in the NISQ [75]. Sung et al.[66] construct models based on the previous method and further improve the sample efficiency of the methods. Furthermore, these zero-order optimization methods leverage the strengths of reinforcement learning [76, 24, 11, 12], Monte Carlo tree search [77, 44, 61], and natural evolutionary strategies [2, 80, 27],  Bayesian [70, 69], as well as Gaussian processes [81].

In addition to these zero-order methods, several other optimization methods have been proposed recently [31, 59, 65, 26, 25]. One interesting example is the quantum natural gradient [65] (QNG), an approximate second-order method, that incorporates the quantum geometric tensor, which is similar to the natural gradient in classical machine learning. While an outcome of measurement is used as an estimate of the gradient in the QNG or the noisy gradient (<ref>) from (<ref>),
the Jordan algorithm [31] encodes the partial derivatives as binary numbers in the computational basis.
This algorithm was later  improved by Gilyen et al.  [26] using high-order finite difference approximations,  and applications to VQAs for a certain class of smooth functions were considered. However, the methods [31, 26] require a significant number of ancilla qubits and complex control logics, due to the binary encoding of partial derivatives. Alternatively, [1] proposed a quantum backpropagation algorithm, which uses $\log d$ copies of the quantum state to compute $d$ derivatives. The overhead for computing $d$ derivatives is $\text{polylog}(d)$ times that of function evaluation (therefore mimicking backpropagation). One of the main drawbacks of their algorithm is that there is an exponential classical cost associated with the process.
For a more restrictive class of cost functions (polynomial functions),
[59] proposed a framework to implement the gradient descent and Newton's methods. This method also requires the coherent implementation of the cost function on a quantum computer using e.g., sparse input oracle, and thus can be challenging to implement in near-term devices.

§.§ A numerical
illustration: Variational quantum eigenvalue solver

As a brief illustration of the performance of noisy GD versus RCD methods, we consider the transverse-field Ising model,

\begin{equation}\label{eq:isingmodel}
H = J\sum_{j=1}^{N-1}Z_jZ_{j+1}+\Delta\sum_{j=1}^NX_j,
\end{equation}
with the coefficient $J=1$ and $\Delta = 1.5$.  Here, $N$ denotes the number of qubits, and $X_j,Z_j$ are Pauli operators acting on the $j$-th qubit. In <ref>, we set $N=10$. To implement the quantum circuits, we use Qiskit Aer-simulator [57]
with the command “result.get_counts" that outputs measurement outcomes as classical bitstrings. We utilize the resulting classical bitstrings to compute partial derivatives by applying the parameter shift rule [68]. Thus, the result in <ref> takes into account the measurement noise.

In each experiment, 10 independent simulations are used  with a fixed initialization. The parameterized quantum circuit used for estimating the ground state energy of the Hamiltonian (<ref>) is given in  Figure <ref> (Appendix <ref>).

We compare the optimization performance of the two methods in terms of the number of partial derivative evaluations. The optimization results in <ref> suggest that RCD requires nearly 4 times fewer partial derivative evaluations than GD to converge to an energy ratio of 0.96 and a fidelity of 0.9, both of which are higher than the energy ratio and the fidelity obtained from GD. This observation can be explained by the analysis in <ref>, i.e.,  RCD can be more efficient than GD when the ratio of Lipschitz constants ($L/L_{\mathrm{avg}}$ or $L/L_{\max}$) is significantly larger than $1$. Specifically, the ratio of the total computational cost of GD to RCD can be linked to the Lipschitz ratios, as summarized in <ref>.  For instance, in the lower panels of <ref>, we observe that the ratio $L/L_{\mathrm{avg}}$ and $L/L_{\max}$ remains above 30 and 8 throughout the iterations. The faster convergence of RCD can be attributed to these large Lipschitz ratios.

The comparison of the performance of GD (red) and RCD (blue) for optimizing the Hamiltonian  (<ref>). The unit of the $x$-axis labels the number of partial derivative evaluations as an indication of the computational complexity.
The top  panels show the approximation of the ground state, including the energy ratio (left) and  fidelity (right). In the bottom panels, we show the ratios of Lipschitz constants obtained from the two methods are compared: $\frac{L}{L_{avg}}$ (left) and $\frac{L}{L_{\max}}$ (right).

§ PRELIMINARIES AND MAIN RESULTS

Before we establish results pertinent to the performance of RCD, we first establish consistent notations and assumptions, which are presented in Section <ref>. Following that, we outline our key theoretical findings in Section <ref>.

§.§ Notations and assumptions

Given a vector $\bm v\in\mathbb{R}^d$, we use standard norms for $\bm v$, including the $2$-norm $\|\bm v\|_2:=\sqrt{\sum_iv_i^2}$ and the $\infty$-norm $\|\bm v\|_{\infty}:=\max_i\left|v_i\right|$. In order to ensure the convergence of gradient-based methods, we list several technical assumptions.

We assume the cost function $f$ satisfies the $L$-smoothness. Specifically, it satisfies the following assumption:
The cost function $f$ is $L$-smooth, in that,

\begin{equation}\label{L2}
\|\nabla f(\bs{\theta}) - \nabla f(\bs{\theta}')\|_2\leq L\|\bs{\theta}-\bs{\theta}'\|_2,\quad \text{for all }\bs{\theta},\bs{\theta}'\in\mathbb{R}^d.
\end{equation}

Since the gradient is Lipschitz continuous, the partial derivatives are Lipschitz continuous as well. We define the componentwise Lipschitz constants,
We say that a function $f$ is $L_i$-smooth with respect to the $i$-th component if

\begin{equation}\label{Li}
|\partial_i f(\bs{\theta}+\bs{e}_ih)-\partial_i f(\bs{\theta})|\leq L_i|h| \quad \forall h\in \mathbb{R},
\end{equation}
where $\partial_i f(\bs{\theta})$ denotes the partial derivative in the $i$-th direction.

From these componentwise Lipschitz constants, we denote the maximum and average of those constants as

\begin{equation}\label{Lmax}
L_{\max}:=\max_{1\leq i \leq d} L_i,\quad L_{\mathrm{avg}}=\frac{1}{d}\sum^d_{i=1}L_i\,.
\end{equation}

As shown in [73], in general we have,

\begin{equation} \label{eq:LiL}
L_i \le L_{avg}\le L_{\max} \le L \le d L_{\max}\,.
\end{equation}

Another interpretation is through the hessian:  When $f$ is twice continuously differentiable, the condition (<ref>) is equivalent to $\nabla^2 f(x)\preceq L I_d$, and similiarly, the condition (<ref>) is equivalent to $\sup_{\bs{\theta}}\left|\partial_{i}^2 f(\bs{\theta})\right|\leq L_i$. We note that both the upper and lower bounds of $L$ in terms of $L_{\max}$ in (<ref>) are tight. If $\nabla^2 f$ is a diagonal matrix, then $L_{\max}=L$, both being the largest diagonal element of $\nabla^2 f$.
(This is the case in which all coordinates are independent of each other, for example, $f = \sum_i\lambda_ix_i^2$.)
On the other hand, if $\nabla^2f = \mathsf{e}\cdot\mathsf{e}^\top$
where $\mathsf{e}\in\mathbb{R}^d$ satisfies $\mathsf{e}_i=1$ for all $i$, then $L = dL_{\max}$.
This is a situation where $f$ is highly anisotropic, e.g.,  $f = (\sum_i x_i)^2/2$, where $L=d$ and $L_{\max}=1$.
In addition, when $L_{avg}=L$, we see that $L_{avg}=L_{\max}=L_i$ for all $i$.

Next, it is important to note that the estimation of the gradients in quantum computing can be susceptible to noise, which stems from the inherent nature of quantum measurements. Consequently, in our analysis and comparative studies of different optimization methods, we will take into account the presence of noise. To facilitate such analysis, we make the following assumption:

[Bounds of the noise with respect to the $2$-norm]
Given any $\bs{\theta}\in\mathbb{R}^d$, we assume that we can find an unbiased random estimate $\bs{g}(\bs{\theta})$ for the gradient $\nabla f(\bs{\theta})$, meaning that

\[
\mathbb{E}\left[\bs{g}(\bs{\theta})\right]=\nabla f(\bs{\theta})\,.
\]

Furthermore, we assume that there exists a constant $\sigma^2_\infty>0$ such that

\begin{eqnarray}  \sigma^2_\infty>\sup_{\bs{\theta}\in\mathbb{R}^d}\max_{1\leq i \leq d}\mathbb{E}\left[|\partial_i f(\bs{\theta}) - g_i(\bs{\theta})|^2\right]\,.\label{eqn:sigma_infinite}
\end{eqnarray}

Here, we also assume $g(\bs{\theta})$ is independent for different $\bs{\theta}$.

Additionally, we assume the existence of a basin encompassing the global minimum, within which $f$ satisfies the Polyak-Łojasiewicz condition (PL) [56], equivalently, the local PŁ condition [41].

[Local PL condition] Define $\mathcal{X}$ as the set of global minima and $f_{\min}$ as the global minimum value evaluated over $\mathcal{X}$. Then there exists a $\delta_f,\mu>0$ such that for any $\bs\theta\in\mathcal{N}(\mathcal{X}):=f^{-1}([f_{\min},\delta_f))$,

\begin{equation*}
\|\nabla f(\bs{\theta})\|^2\geq 2\mu \left(f(\bs{\theta})-f_{\min}\right).
\end{equation*}

It is worthwhile to highlight that the PL condition is defined not on the entire space $\mathbb{R}^d$ but $\mathcal{N}(\mathcal{X})$, which is reasonable in the context of the variational quantum algorithm. We support this argument with the following remark.

Let $f(\bs\theta)$ be a cost function defined by some parameterized quantum circuit (<ref>). Note that $f$ is  periodic and smooth, due to its specialized form. By the extreme value theorem, we see that there exist global maximum and minimum of $f$, denoted by $\bs\theta_{\max}$ and $\bs\theta_{\min}$. In general, $f$ is not constant, which means that $f_{\max}>f_{\min}$. Had $f$ satisfied the global PL condition, it would have followed that at the global maximum $\bs\theta_{\max}$,

\begin{equation}
0=\|\nabla f(\bs\theta_{\max})\|^2\geq 2\mu \left(f_{\max}-f_{\min}\right)\geq 0,
\end{equation}
which gives a contradiction to the general case that $f_{\max}>f_{\min}$. As another case study, if $f$ is assumed to be convex, namely,

\begin{equation}
f(\bs \theta')\geq f(\bs \theta) + (\nabla f(\bs\theta),\bs\theta'-\bs\theta)\text{ for all }\bs\theta,\bs\theta'\in\mathbb{R}^d,
\end{equation}
then setting $\bs\theta = \bs\theta_{\max}$ and $\bs\theta'=\bs\theta_{\min}$ results in a contradiction. Therefore, the cost function $f$ that is constructed from an ansatz similar to (<ref>), will not satisfy global PL or convex conditions in general.

§.§ Main result: complexity comparison of GD and RCD

In this study, our main focus is to compare the complexity of noisy gradient descent (GD) and randomized coordinate descent (RCD) under the assumptions of a local Polyak-Łojasiewicz (PL) condition <ref>. For the sake of simplicity, in the remaining part of this paper, we will refer to “noisy gradient descent" and “noisy randomized coordinate descent" as “GD" and “RCD", respectively, without explicitly mentioning the term “noisy".

The main theoretical results are summarized in the following two theorems:

Assume $f$ is a $L$-smooth function that satisfies assumption <ref> and $g$ satisfies assumption <ref>. Given $\epsilon>0$ small enough, if $f(\bs{\theta}_1)\leq \delta_f$ and $a_n=\Theta(\min\left\{\mu \epsilon/(L\sigma^2_\infty d),1/L\right\})$ in GD (<ref>), then with probability $1-f(\bs{\theta}_1)/\delta_f-o(1)$, there  exists  at least one

\begin{equation}
n<N=\tilde{\Theta}(\max\left\{L\sigma^2_\infty d/(\mu^2\epsilon),L/\mu\right\})
\end{equation}
such that $f(\bs{\theta}_n)\leq f_{\min}+\epsilon$.

Assume $f$ is a $L$-smooth function that satisfies assumption <ref> and $g$ satisfies assumption <ref>. Given $\epsilon>0$ small enough, if $f(\bs{\theta}_1)\leq \delta_f$ and $a_n=\Theta(\max\left\{\mu \epsilon/(L_{\mathrm{avg}}\sigma^2_\infty d),1/L_{\max}\right\})$ in RCD (<ref>), then with probability $1-f(\bs{\theta}_1)/\delta_f-o(1)$, there  exists  at least one

\begin{equation}
n<N=\tilde{\Theta}(\max\left\{L_{\mathrm{avg}}\sigma^2_\infty d^2/(\mu^2\epsilon),L_{\max}d/\mu\right\})
\end{equation}

such that $f(\bs{\theta}_n)\leq f_{\min}+\epsilon$.

Based on the theorem mentioned above, to achieve $f(\boldsymbol{\theta}_n)-f_{\min} \leq \epsilon$, we can select the learning rate $a_n = \frac{\mu \epsilon}{L\sigma^2_\infty d}$ for GD and $a_n = \frac{\mu \epsilon}{L_{\mathrm{avg}}\sigma^2_\infty d}$ for RCD. Recalling equation (<ref>), we observe that $L_{\mathrm{avg}} \leq L$, which means that we could use a larger learning rate for RCD. This choice aligns with the learning rates utilized in the numerical experiments presented in Section <ref> as well as those in Section <ref>.

We compare the complexity of the noisy GD and RCD methods with the estimates of the number of iterations. First, according to the above result, we conclude that the number of iterations required for GD is $N=\tilde{\Theta}\left(\frac{L\sigma^2_\infty d}{\mu^2\epsilon}\right)$[This complexity aligns with the classical theoretical results for gradient descent (GD), which typically assume the presence of strong convexity or a local PL condition for the function $f$.], while for RCD,  we have $N=\tilde{\mathcal{O}}\left(\frac{L_{\mathrm{avg}}\sigma^2_\infty d^2}{\mu^2\epsilon}\right)$. Notably, in RCD, there is an additional factor of $d$, which can be understood in the expectation sense: During each iteration of the noisy RCD, the randomness arises from two sources: the random direction $i_n$ and the noisy partial derivative $g_{i_n}(\bm \theta_n)$. By taking the conditional expectation with respect to $\bs{\theta}_n$, we obtain:

\begin{equation}\label{eqn:expect_rcd}
\mathbb{E}(\bs{\theta}_{n+1}|\bs{\theta_n})=\bs{\theta}_n - \frac{a_n}{d}\nabla f(\bs{\theta}_n)\,.
\end{equation}

Compared with (<ref>), there is an extra $1/d$ factor in the expectation of RCD. Consequently, in each iteration, the rate of decay of the cost function is smaller in RCD compared to GD.
Consequently, we anticipate that RCD would necessitate more iteration steps to achieve convergence. On the other hand, it is also important to note that in certain scenarios where $L_{\mathrm{avg}}d$ is comparable to $L$, the number of iterations required for RCD is comparable  to that of GD.

Meanwhile, it is important to point out that a more practical criterion for comparing the two methods is the cumulative cost of each method, which is represented by the number of partial derivative calculations from the quantum circuits. This is because quantum algorithms for estimating the gradient have a cost proportional to $d$.  Since each iteration of GD needs to calculate the full gradient ($d$ partial derivatives), the total number of partial derivative estimations in GD is

\[
N_{\text{partial,GD}}=\tilde{\Theta}\left(\frac{L\sigma^2_\infty d^2}{\mu^2\epsilon}\right)\,.
\]

In contrast, the number of partial derivative estimations in RCD is:

\[
N_{\text{partial,RCD}}=\tilde{\mathcal{O}}\left(\frac{L_{\mathrm{avg}}\sigma^2_\infty d^2}{\mu^2\epsilon}\right)\,.
\]

From equation (<ref>), we can deduce that:

\[
\tilde{\Omega}(N_{\text{partial,RCD}})=N_{\text{partial,GD}}=d\tilde{\mathcal{O}}(N_{\text{partial,RCD}})\,.
\]

This suggests that the computational cost of RCD is $L/L_{\mathrm{avg}}$ times cheaper than that of GD. In an extreme case where $f$ is highly skewed, i.e., $L/L_{\mathrm{avg}} \approx d$,  RCD can reduce the computational cost by a factor of the dimension $d$, which will be a significant reduction for large quantum circuits.

In addition to the complexity result, it is worth noting that the two methods exhibit similar success probability, which is approximately $1-f(\bs{\theta}_1)/\delta_f$, as indicated by the two aforementioned theorems. This observation is quite surprising, as each iteration of RCD appears noisier due to the random selection of the updating direction $i_n$.  Intuitively, this suggests that we might need to choose a smaller learning rate $a_n$ to ensure stability in RCD, which would consequently increase its complexity. However, our theory unveils that choosing a similar learning rate $a_n$ is adequate to stabilize RCD. To elucidate this point, it's important to recognize that, on average, RCD behaves equivalently to GD. By conducting more iterations, RCD can approximate its average behavior (expectation), effectively mitigating the extra randomness introduced by $i_n$. This compensation mechanism ensures that the success probabilities remain consistent between the two methods.

§ PROOF OF MAIN RESULTS

In this section, we provide the proofs for Theorems <ref> and <ref>. We will start by showing the stochastic stability of the two methods in Section <ref>. This will guarantee that the parameter $\bm\theta$ is likely to stay close to the global minimum until attaining a small loss. Following that, in Section <ref>, we utilize the local Polyak-Łojasiewicz (PL) condition around the global minimum to establish the convergence of $f(\bs{\theta}_n)$. In all of the following theoretical results and the corresponding proofs in the appendices, we assume $f_{\min}=0$ without loss of generality by modifying the original function as

\begin{eqnarray}\label{eq:shiftfmin}
f(\bs\theta)\leftarrow f(\bs\theta) - f_{\min}.
\end{eqnarray}
Thus, all results in this section can be reformulated for the original cost function by the substitution (<ref>), which will yield  <ref>.

§.§ Stochastic stability

In the context of optimization, stability and convergence are not separate properties. In a deterministic algorithm, convergence immediately guarantees stability. However, this connection does not hold for stochastic processes in general. For instance, when optimization methods such as noisy GD, SGD, or noisy RCD are applied, discrete-time stochastic processes are generated. In such cases, a convergence theory must be developed for a collection of random paths, which can exhibit different convergence behaviors among themselves.

In our specific case, we anticipate that when $\bs{\theta}_n$ remains within the basin $\mathcal{N}(\mathcal{X})$ and the learning rate is correctly chosen, both the GD and the RCD methods, when the gradient is exactly calculated, converge to a global minimum due to the local PL condition stated in assumption <ref>. However, in the presence of noise in the gradient and the use of a constant learning rate, it is generally impossible to ensure that $\bs{\theta}_n\in \mathcal{N}(\mathcal{X})$ almost surely, unless a different strategy is adopted such as the decreasing learning rates [53, 23, 34]. On the other hand, the purpose of the optimization algorithm is to minimize the loss function, which means that it suffices to ensure stability until small loss is achieved. To quantify such a likelihood, in this section, we demonstrate that when $\bs{\theta}_0\in \mathcal{N}(\mathcal{X})$, there exists a finite probability that $\bs{\theta}_n$ obtained from GD and RCD remain within $\mathcal{N}(\mathcal{X})$ until achieving a small loss. This provides a high probability of convergence for the two methods.

We summarize the result for noisy GD in the following lemma.

Assume that $f$ is a $L$ smooth function that satisfies the assumption <ref> and $g$ satisfies the assumption <ref>. If $f(\bs{\theta}_1)\leq \delta_f$ and the learning rate is chosen as follows,

\[
a_n=a<\min\left\{\frac{1}{L}\;,\frac{2\mu \delta_f}{L\sigma^2_\infty d}\right\}\,,
\]

then,  with high probability, iterations of noisy GD (<ref>) remain in $f^{-1}([0,\delta_f))$ until a small loss is achieved. Specifically,

\begin{equation}\label{eqn:p_gd}
\mathbb{P}\left\{\exists N>0\ \text{such that}\ f(\bs\theta_N)\notin\mathcal{N}\ \text{and}\  f(\bs\theta_n)>\frac{La\sigma^2_\infty d}{\mu},\ \forall n<N\right\}\leq \frac{f(\bs{\theta}_{1})}{\delta_f}.
\end{equation}

In light of equation (<ref>), if we select the learning rate $a_n$ to be sufficiently small, then with a probability of $1-\frac{f(\bs{\theta}_{1})}{\delta_f}$, the parameters are guaranteed to achieve a small loss before escaping the basin.

Despite infrequent updates of the gradient components, RCD still demonstrates a similar level of stochastic stability. This key observation is summarized in the following lemma:

Assume that $f$ is a $L$-smooth function that satisfies assumption <ref> and $g$ satisfies assumption <ref>. Given any $f(\bs{\theta}_1)<\delta_f$, if one chooses the learning rate

\[
a_n=a<\min\left\{\frac{1}{L_{\max}},\; \frac{d}{\mu}, \;\frac{2\mu \delta_f}{L_{\mathrm{avg}}\sigma^2_\infty d}\right\}\,,
\]

then,  with high probability, iterations from the noisy RCD (<ref>) stay at $f^{-1}([0,\delta_f))$ until achieving a small loss. Specifically,

\[
\mathbb{P}\left\{\exists N>0\ \text{such that}\ f(\bs\theta_N)\notin\mathcal{N}\ \text{and}\  f(\bs\theta_n)>\frac{L_{\mathrm{avg}}a\sigma^2_\infty d}{\mu},\ \forall n<N\right\}\leq \frac{f(\bs{\theta}_{1})}{\delta_f}.
\]

The proofs of Lemma <ref> and <ref> are provided in Appendices <ref> and <ref>, respectively. The core concept of these proofs is based on the construction of a specialized  supermartingale and the utilization of Markov's inequality. For example, to prove  Lemma <ref>, we define a stochastic process

\[
\begin{aligned}
&f(\bs\theta_n)\mathbb{I}_n,\quad n<\tau\\
&f(\bs\theta_\tau)\mathbb{I}_\tau,\quad n\geq\tau
\end{aligned}\right.\,.
\]

where the indicator random variable is given by,

\[
\begin{aligned}
&1,\quad \text{if}\quad \{\bs{\theta}_k\}^{n-1}_{k=1}\subset f^{-1}([0,\delta_f))\\
&0,\quad \text{otherwise}.
\end{aligned}\right.\,,
\]

and the stopping time
\[
\tau=\inf \left\{k:f(\bs\theta_k)\leq \frac{La\sigma^2_\infty d}{\mu}\right\}\,.
\]

We observe that $V_n$ is a meticulously crafted supermartingale, allowing us to distinguish between stable and unstable events. In particular, we demonstrate that if $\bs\theta_n$ exits the basin before it reaches $f(\bs\theta_n)=\frac{La\sigma^2_\infty d}{\mu}$ (an unstable event), then $\sup_n V_n\geq \delta_f$.
Therefore, we can employ $V_n$ as a categorizer and the probability of failure of GD can be characterized by the value of $V_n$. More specifically,

\[
\mathbb{P}\left\{\exists N>0\ \text{such that}\ f(\bs\theta_N)\notin\mathcal{N}\ \text{and}\  f(\bs\theta_n)>\frac{La\sigma^2_\infty d}{\mu},\ \forall n<N\right\}\leq \mathbb{P}\left\{\sup_{n}V_n\geq \delta_f\right\}\,.
\]

Except for its use as a categorizer, we have designed $V_n$ in such a way that it is a supermartingale, meaning $\mathbb{E}(V_{n+1}|\bs{\theta}_{k\leq n})\leq V_{n}$. Therefore, we can use Markov's inequality for supermartingales to bound the supremum of $V_n$ and achieve the desired result.

§.§ Convergence analysis

In this section, we present the convergence properties of noisy GD and RCD methods.
It is important to note that Theorems <ref> and <ref> directly follow from Theorems <ref> and <ref>, respectively.

Our first theorem shows the convergence performance of the noisy GD method,
Assume $f$ is a $L$-smooth function that satisfies Assumption <ref> and $g$ satisfies
<ref>. Given any precision $0<\epsilon<\delta_f$, the initial guess $f(\bs{\theta}_1)<\delta_f$, and the probability of failure $\eta\in\left(\frac{f(\bs\theta_1)}{\delta^{-1}_f},1\right)$, we choose the learning rate in (<ref>)

\[
a_n=a=\mathcal{O}\left(\min\left\{\frac{1}{L},\frac{\mu \epsilon}{L\sigma^2_\infty d}\right\}\right),
\]

and the total number of iterations

\[
N=\Omega\left(\frac{1}{\mu a\eta}\log\left(\frac{f(\bs{\theta_1})}{\left(\eta-\frac{f(\bs\theta_1)}{\delta_f}\right)\epsilon}\right)\right)\,.
\]

Then, with probability $1-\eta$, we can find at least one $\bs\theta_m$ with $1\leq m\leq N$ such that $f(\bs\theta_m)\leq \epsilon$. In particular,

\[
\mathbb{P}\left\{\exists m\leq N,\ f(\bs{\theta}_m)\leq\epsilon \right\}\geq 1-\eta,
\]

Next, we state the convergence property of the noisy RCD method in the following theorem,

Assume $f$ is a $L$-smooth function that satisfies Assumption <ref> and $g$ satisfies Assumption <ref>. Given any precision $0<\epsilon<\delta_f$, the initial guess $f(\bs{\theta}_1)<\delta_f$, and the probability of failure $\eta\in\left(\frac{f(\bs\theta_1)}{\delta^{-1}_f},1\right)$, we choose the learning rate in (<ref>)

\[
a_n=a=\mathcal{O}\left(\min\left\{\frac{1}{L_{\max}},\frac{d}{\mu},\frac{\mu \epsilon}{L_{\mathrm{avg}}\sigma^2_\infty d}\right\}\right),
\]

and the total number of iterations

\[
N=\Omega\left(\frac{d}{\mu a\eta}\log\left(\frac{f(\bs{\theta_1})}{\left(\eta-\frac{f(\bs\theta_1)}{\delta_f}\right)\epsilon}\right)\right)\,.
\]

Then, with probability $1-\eta$, we can find at least one $\bs\theta_m$ with $1\leq m\leq N$ such that $f(\bs\theta_m)\leq \epsilon$. In particular,

\[
\mathbb{P}\left\{\exists m\leq N,\ f(\bs{\theta}_m)\leq\epsilon \right\}\geq 1-\eta,
\]

The proofs of these theorems can be found in the appendix <ref>.

We emphasize that <ref> and <ref> are general convergence results that require only mild conditions. Specifically, <ref> can be used to demonstrate the stability and convergence of the traditional SGD algorithm when the right assumptions are in place. A convergence result analogous to the one previously discussed has been investigated in <cit.>, where the authors impose a more stringent requirement on the cost function [6]. In our work, we demonstrate the convergence of noisy GD using more sophisticated techniques in probability theory and adopt a weak version of probabilistic convergence [33]. In addition, our approach can be directly extended to show the convergence of noisy RCD as in <ref>, which to the best of our knowledge, has not been established before. These two theorems suggest that with a high probability, the loss function can achieve small loss during the training process. In other words, it is likely that the parameter $\bs\theta$ remains in the basin $\mathcal{N}$ until the precision $\epsilon$ is attained at some point. After that, the optimization algorithm could diverge unless a certain strategy is applied, for example, a schedule of decreasing learning rates or an early stopping criterion.

Our theoretical result clarifies a relation between the learning rate and the desired precision in optimization. For example, the precision $\epsilon$ is manifested in the upper bounds of the learning rates in <ref> and <ref>. Thus, to reach precision $\epsilon$, it is suggested to use an $\mathcal{O}(\epsilon)$ learning rate. Otherwise, due to the stability issue, the trajectory is no longer guaranteed to converge to the precision with positive probability.

We present the roadmap for proving Theorem <ref> as follows: Define the stopping time

\[
\tau=\inf \left\{k:f(\bs\theta_k)\leq \epsilon\right\}\,.
\]

To prove Theorem <ref>, it suffices to demonstrate that the probability of failure $\mathbb{P}(\tau>N)$ is small. Since the learning rate $a_n$ is selected to be sufficiently small and, according to the lemma <ref>, it is likely that $\bs\theta_n$ will remain within the basin until the $\epsilon$ loss is achieved[Rigorously, we must also take into account the possibility that the optimization algorithm does not reach $\epsilon$ loss in a finite number of iterations.]. Thus, informally, it suffices for us to assume $\bs\theta_n\in\mathcal{N}$. The next step is to find an upper bound for the probability of failure $p_{\mathrm{fail}}=\mathbb{P}(\tau>N)$. Using the local PL condition, we can show that when $\epsilon<f(\bs\theta_n)<\delta_f$,

\[
\mathbb{E}(f(\bs{\theta}_{n+1})|\bs{\theta}_n)\leq \left(1-\frac{\mu a}{2}\right)f(\bs{\theta}_n)\,,
\]

meaning that the conditional expectation of $f(\bs{\theta}_{n+1})$ decays to zero with rate $\left(1-\frac{\mu a}{2}\right)$. Inspired by this observation, we can construct a supermartingale to show that, if $\tau>N$, with high probability, we have $\inf_{1\leq n\leq N}f(\bs\theta_n)\leq \epsilon$. We note that this event is complementary to the failure event $\{\tau>N\}$. Consequently, we obtain an upper bound for $p_{\mathrm{fail}}$.

§ NUMERICAL RESULTS

In <ref>, depicted in <ref>, we have demonstrated that the noisy RCD leads to faster convergence than the noisy GD for VQE problems. In this section, we extend our investigation to gauge the efficiency of noisy RCD applied to various other variational quantum algorithms, especially those involving non-convex optimization problems. The implementation of these algorithms is executed on classical computers. To emulate quantum measurement noise, the partial derivatives undergo perturbation through additive Gaussian noise as outlined in <ref> [
The derivative with noise is computed by adding Gaussian noise to the original derivative: $\partial_i f(x) \leftarrow \partial_i f(x) + \epsilon$, where $\epsilon$ follows a Gaussian distribution, denoted as $\mathcal N(0, \sigma)$. In this notation, $\sigma$ signifies the standard deviation, defining the intensity of the Gaussian noise.]. Subsequently, we substantiate this approximation through a numerical experiment on a quantum simulator. This experiment further also proposes suitable values for the strength of the Gaussian noise that we will introduce in the upcoming numerical tests to appropriately mimic the measurement noise.

In the experiment presented in <ref>, we utilize Qiskit-0.44.1 [57].

The algorithms for subsequent examples are implemented using Numpy [29] and Jax [10]. We conducted each experiment ten times, employing different random initializations for each run. All tests are executed on an Intel(R) Xeon(R) CPU @ 2.20GHz, complemented by a T4 GPU.

§.§ Analyzing the noise distribution

Building on the numerical experiment detailed in <ref> and executed in Qiskit, we analyze the statistics of the partial derivatives derived from the quantum circuit.
<ref> showcases the histograms representing 10,000 estimates of partial derivatives with respect to the initial 12 directions, while the histograms for the remaining directions are presented in <ref>. Each estimate of the partial derivatives is averaged over 1000 shots.
From all histograms, we can clearly see that the distribution is closely approximated by a Gaussian distribution.
In addition, the magnitude of the standard deviation of partial derivative estimates across all directions is comparable. These observations support assumptions of the noise model in <ref>. For simplicity, we will employ the Gaussian noise model in our subsequent investigations to compare the performance of the noisy GD and RCD methods.

Histograms of the estimated partial derivatives:  Each panel displays the histogram of 10000 partial derivative estimates in one of the first 12 directions, which are obtained by applying the parameter-shift rule for the ansatz in <ref>. The sampling of the partial derivatives is carried out at a suboptimal point chosen from one simulation used in <ref>, where the fidelity is about 0.889.

In the next two sections, we conduct a comprehensive comparison between noisy RCD and GD across a broad spectrum of variational quantum algorithms and applications.

§.§ VQE with a varied circuit structure

In <ref>, we utilize the VQE for the TFIM (<ref>) employing both the noisy GD and the noisy RCD. In this section, we tackle the same optimization task but with a modified setup.
Specifically,  <ref> depicts the PQC [32] utilized in the experiments showcased in <ref>, distinct from that presented in <ref>.

In the experiments illustrated in <ref>, each optimization outcome derives from 10 identical simulations with the same initial condition. We set the learning rates for the RCD and GD at $0.3$ and $0.05$, respectively. Each experiment utilizes 10,000 shots, with $18$ trainable parameters. Results shown in <ref> demonstrate that, compared to GD, RCD requires nearly three times fewer partial derivative evaluations to converge.

@C=1.0em @R=0.2em @!R

q_0 :       q_0 :       R_Y ([0])     R_Z ([3])     1     2     R_Y ([6])     R_Z ([9])         1     2     R_Y ([12])     R_Z ([15])

q_1 :       q_1 :       R_Y ([1])     R_Z ([4])             1     R_Y ([7])     R_Z ([10])             1     R_Y ([13])     R_Z ([16])

q_2 :       q_2 :       R_Y ([2])     R_Z ([5])                 R_Y ([8])     R_Z ([11])                 R_Y ([14])     R_Z ([17])

A variational circuit ansatz is employed for the Transverse-Field Ising Model expressed in Equation (<ref>), utilizing $3$ qubits. This circuit is a parameterized construct comprised of alternating rotation and entanglement layers. Each rotation layer involves the application of single qubit gates, specifically Rotation-y and Rotation-z gates, to all qubits. In contrast, the entanglement layer employs two-qubit gates, namely the controlled-X gate, to facilitate entanglement among the qubits. The ansatz is designated with $18$ parameters.

Performance comparison between GD (red) and RCD (blue) in terms of energy ratio and Lipschitz constant ratios for optimizing the Hamiltonian  (<ref>). The energy ratio $ E / E_{GS}$ is presented in the left panel, while the Lipschitz constant ratios, denoted as $\frac{L}{L_{\text{avg}}}$ and $\frac{L}{L_{\max}}$, are shown in the middle and right panels respectively. The shaded areas in each panel represent variations observed across multiple trials.

§.§ Quantum Approximate Optimization Algorithm (QAOA) for quantum Hamiltonians

The Quantum Approximate Optimization Algorithm (QAOA)
[22], originally devised for solving combinatorial problems, is a leading example for demonstrating quantum advantage
on near-term quantum computers. As introduced in [22], the QAOA utilizes a parametrized quantum circuit (PQC), which naturally enables optimization through the variational quantum algorithm.

In a generalized QAOA model, we begin with an initial quantum state $\ket{\psi_i}$, which can be easily prepared in experiments, and let it evolve by a parameterized unitary transformation,

\begin{equation}
\ket{\psi(\bs\alpha,\bs\beta)} = U(\{\alpha_j, \beta_j\}_{j=1}^p) \ket{\psi_i} = e^{-i H_2 \beta_p} e^{-i H_1 \alpha_p} \cdots e^{-i H_2 \beta_1} e^{-i H_1 \alpha_1} \ket{\psi_i}
\label{eqn:qaoawf},
\end{equation}
where the vector $\bs\alpha$ (or $\bs\beta$) enumerates the parameters $\alpha_j$ (or $\beta_j$), and thus the total number of parameters is $2p$ and the unitary transformation alternates between two kinds of parameterized unitary transformations.
With this ansatz, the optimization is performed with the parameters $\{\alpha_j,\beta_j\}$  associated with the application-dependent Hamiltonian matrices $H_1$ and $H_2$, respectively.

In the subsequent sections, we will consider optimization problems based on the QAOA (<ref>). We will conduct a comparative analysis of the noisy GD and RCD for various QAOA models that will span a range of systems, including the Ising model (refer to Section <ref>), the Heisenberg model (refer to Section <ref>), and Variational Quantum Factoring (refer to Section <ref>).

§.§.§ QAOA – Ising Model

In this section, we parameterize the transverse-field Ising model by a  Hamiltonian

\begin{equation}\label{eq: ising qaoa}
H[h]= \sum_{j=1}^{N-1}Z_{j+1}Z_j + \sum_{j=1}^{N}(Z_j +  h X_j),
\end{equation}

where $N$ denotes the total number of qubits. The global control field $h\in\{\pm 4\}$ takes two discrete values, corresponding to the two alternating QAOA generators $H_1=H[-4]$ and $H_2=H[+4]$ [13, 75]. The initial state \( \ket{\psi_i} \) corresponds to the ground state of \( H[-2] \), while the desired target state \( \ket{\psi_\ast} \) is selected as the ground state of \( H[+2] \). The variational problem aims to optimize the fidelity [Fidelity serves as a metric for optimization. However, one caveat of utilizing fidelity is its reliance on the ground state. In this context, we assume the presence of an oracle capable of producing the fidelity value. Subsequently, we also employ energy as an observable metric for optimization purposes.],

\begin{eqnarray}
\max_{\{\alpha_i, \beta_i\}_{i=1}^p} \mathcal{F}(\{\alpha_i, \beta_i\}_{i=1}^p) = \max_{\{\alpha_i, \beta_i\}_{i=1}^p}\abs{ \mel{\psi_{*}}{ U(\{\alpha_i, \beta_i\}_{i=1}^p)}{ \psi_i} }^2,
\label{eqn:qaoa_HM}
\end{eqnarray}

\begin{equation}
U(\{\alpha_i, \beta_i\}_{i=1}^p) \ket{\psi_i}= e^{-i H_2 \beta_p} e^{-i H_1 \alpha_p} \cdots e^{-i H_2 \beta_1} e^{-i H_1 \alpha_1}\ket{\psi_i}.
\label{eqn:qaoa_HM_psi}
\end{equation}
We note that the fidelity optimization (<ref>) is equivalent to the optimization of the form (<ref>) by letting the Hamiltonian be $\ketbra{\psi_*}$.

In the numerical test, we choose a system from \eqref{eq: ising qaoa} with three qubits ($N=3$), and  then apply both GD and RCD methods in the optimization.
Figure~\ref{fig:qaoa} shows the optimization results obtained from the noisy GD and RCD with the respective learning rates of 0.0045 and 0.015 by using an ansatz defined with 20 parameters. By adjusting the learning rate and tracking the stability, We observe that RCD permits a larger learning rate in comparison to GD, while maintaining the stability. Similar to the results presented in \cref{fig:vqe}, we compare the performance of the two methods in terms of the number of partial derivative evaluations. From \cref{fig:qaoa}, We observe that noisy RCD converges much faster than noisy GD. While RCD achieves a fidelity near 1 with 500 partial derivative evaluations, GD only attains a fidelity below 0.25 with an equivalent number of evaluations.
This computational effectiveness of RCD can be attributed to the large ratios of Lipschitz constants shown in \cref{fig:qaoa}, which are obtained along the optimization trajectories.

\begin{figure}[htbp]
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/tfim2_training_new.png}
\end{subfigure}
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/tfim2_training_aux_new.png}
\end{subfigure}

\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/tfim2_lips_new.png}
\end{subfigure}
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/tfim2_lips_max_new.png}
\end{subfigure}
\caption{ \label{fig:qaoa} Performance comparison between the noisy GD and RCD for the Ising model \eqref{eq: ising qaoa}. The corresponding Lipschitz constant ratios, denoted as $\frac{L}{L_{\text{avg}}}$ and $\frac{L}{L_{\max}}$, are presented in the bottom figures. The shaded areas within the figures represent variations that have been observed across ten random realizations. The optimization is performed for parameters with dimension equal to 20.
\end{figure}

\subsubsection{QAOA -- Heisenberg Model} \label{sec:heisenberg}

Our second test problem with QAOA is the (anisotropic) spin-$1$ Heisenberg model, $H =H_1+H_2$, with the alternating Hamiltonians given by, \[H_1\!\!=\!\! J\sum_{j=1}^N (X_{j+1}X_j \!+\! Y_{j+1}Y_j), \quad H_2= \Delta \sum_{j=1}^N Z_{j+1}Z_j, \nonumber\] with anisotropic parameter $Δ/J=0.5$ (topological/Haldane~[16, 55, 37, 78]). For the Heisenberg model, we consider a system consisting of eight qubits ($N=8$) and choose the fidelity as a measure for optimization, similar to the setup for the results in \cref{fig:qaoa}. We set the antiferromagnetic initial state to \(\ket{\psi_i} =  \ket{10101010} \). The target state  is the ground state of the Hamiltonian \(H \!=\! H_1 + H_2\). We employ the QAOA ansatz represented by Eqn.~\eqref{eqn:qaoa_HM_psi} and carry out the fidelity optimization detailed in Eqn.~\eqref{eqn:qaoa_HM}.

Figure~\ref{fig:heisenberg} showcases the performance outcomes  from noisy GD and RCD simulations with learning rates set to 0.01 and 0.1, respectively. This QAOA model involves  28 parameters.
The fidelity result shows that RCD converges to the target state much faster than GD. This phenomenon can be elucidated by noting that the ratios of Lipschitz constants derived from both noisy methods, $L/L_avg$ and $L/L_max$, average around 10 and 6 along the trajectories, respectively. Especially, the magnitude of the ratio $L/L_max$ is similar to that of the ratio of the numbers of partial derivative evaluations to reach a high fidelity $>0.8$ from both noisy methods, as shown in \cref{fig:heisenberg}. Based on the observed numerical results, a high ratio of $L/L_max$ is responsible for the  efficiency of RCD in this optimization problem.

\begin{figure}[htbp]
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/heisenberg2_training.png}
\end{subfigure}
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/heisenberg2_training_aux.png}
\end{subfigure}

\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/heisenberg_lips.png}
\end{subfigure}
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/heisenberg_lips_max.png}
\end{subfigure}
\caption{ \label{fig:heisenberg} Performance comparison between noisy GD and RCD for the Heisenberg model. The corresponding Lipschitz constant ratios, denoted as $\frac{L}{L_{\text{avg}}}$ and $\frac{L}{L_{\max}}$, are presented in the middle and right. The shaded areas within the figure represent variations that have been observed across ten random realizations. The optimization is performed in dimensions of 28. }
\end{figure}

\subsection{QAOA for classical combinatorial optimization problems}

Quadratic Unconstrained Binary Optimization (QUBO) problems have significant applications in fields such as finance, logistics, and machine learning, etc. Recognized as one prominent optimization model in quantum computing, QUBO consolidates a wide range of combinatorial optimization problems~[20, 35, 4, 28] and translates them into identifying the ground state of classical Ising models~[42].

The goal of QUBO is to identify a sequence of binary variables ($0$ or $1$) that minimize a quadratic function. Specifically, a cost function $f_Q$ is constructed over the set of binary vectors, $𝔹^n$:
\begin{equation}\label{eq:qubof}
f_{Q}(x)=x^{\top} Q x=\sum_{i,j=1}^n Q_{i j} x_{i} x_{j}.
\end{equation}
In this context, $𝔹={0,1}$ signifies the set of binary values (or bits), and $𝔹^n$ represents the collection of binary vectors with length $n>0$. A symmetric, real-valued matrix $Q ∈ℝ^n ×n$ is introduced, with each element $Q_i j$ determining the weight for the corresponding pair of indices $i, j ∈1, …, n$. For example, if $i=j$, the term $Q_ii x_i^2$ contributes $Q_ii$ to the function value when $x_i=1$. On the other hand, if $i≠j$, the term $Q_ij x_i x_j$ contributes $Q_ij$ to the function value when both $x_i=1$ and $x_j=1$.

Overall, QUBO seeks to minimize the function $f_Q$ over the set of binary vectors by determining an optimal minimizer $x^*$,
\begin{equation}
x^{*}=\underset{x \in \mathbb{B}^{n}}{\arg \min } f_{Q}(x).
\end{equation}
Incorporating the variational quantum algorithm into QUBO, we reformulate the cost function using the substitution:
\begin{equation}
x_i = \frac{1-Z_i}{2}\text{ or }\frac{1+Z_i}{2},
\end{equation}where the variable $x_i$ is supplanted by the Pauli Z matrix operating on the $i$-th qubit. This replacement facilitates the formulation of a model Hamiltonian. Its ground state can be approximated by minimizing the expected energy via the variational quantum algorithm, as elaborated in \cref{sec:maxcut}.

In the following sections, we evaluate the performance of the noisy GD and RCD across various QUBO applications, focusing on the ground state energy estimation. These applications encompass Max-Cut in \cref{sec:maxcut}, the Traveling Salesman Problem in \cref{sec:tsp}, and Variational Quantum Factoring in \cref{sec:factor}.

\subsubsection{Max-Cut} \label{sec:maxcut}

For the Max-Cut problem, the graph employed in our numerical experiments is presented as follows:

\begin{center}
\begin{tikzpicture}
[scale=.8,auto=left,every node/.style={circle,fill=blue!20}]

\node (n0) at (2,4) {0};
\node (n1) at (0,2) {1};
\node (n2) at (2,0) {2};
\node (n3) at (4,2) {3};

\foreach \from/\to in {n0/n1,n0/n2,n0/n3,n1/n2,n2/n3}
\draw (\from) -- (\to);
\end{tikzpicture}
\end{center}
The global cost function is designed to maximize $C = ∑_(i,j) ∈E x_i (1 - x_j)$, where $E$ represents the edges in the graph. For the given graph, the QUBO problem can be formulated as:
\[
\min_{x_i\in\{0,1\}}-3x_{0}^{2} + 2x_{0}x_{1} + 2x_{0}x_{2} + 2x_{0}x_{3} - 2x_{1}^{2} + 2x_{1}x_{2} - 3x_{2}^{2} + 2x_{2}x_{3} - 2x_{3}^{2}\,.
\]
In order to construct the corresponding Hamiltonian, we associate the binary variables \( x_i \) with the Pauli \( Z \) matrices, denoted as \( Z_i \), which act on individual qubits. Taking into account the relationship between the binary variables \( x_i \) and the Pauli matrices \( Z_i \), defined by the equation \( x_i = \frac{{1 - Z_i}}{2} \), the cost function is articulated by the Hamiltonian:
\begin{equation}
H = \frac{1}{2} I - 3 Z_0 + \frac{1}{2} Z_0 Z_1 + \frac{1}{2} Z_0 Z_2 + \frac{1}{2} Z_0 Z_3 + \frac{1}{2} Z_1 Z_2 + \frac{1}{2} Z_2 Z_3.
\end{equation}

Using this Hamiltonian, we construct a parameterized quantum circuit with four qubits ($N=4$) and 20 parameters. The circuit consists of alternating single-gate rotations, denoted as $U_single (θ)=∏_i=1^n RY(θ_i)$~\footnote{Each layer of rotation gates includes a rotation-Y gate applied to every qubit.} and entangler gate $U_entangler $~\footnote{The entanglement layer incorporates two-qubit gates for qubit entanglement without tunable parameters. In this experiment, the entangler gate employs controlled-Z gates. For a comprehensive explanation, refer to the circuit architecture in Appendix~\ref{sec:circuit_qubo}}. The configuration of the parametrized quantum circuit is illustrated in Figure~\ref{fig:circuit_qubo} in Appendix~\ref{sec:circuit_qubo}. This structure resembles the variational quantum circuit of the QAOA, with the ansatz given by $|ψ(θ)⟩=[U_single (θ) U_entangler ]^m|+⟩$. For the optimization process, we assign a learning rate of 0.1 for GD and 3.0 for RCD and select energy as the optimization metric.

As illustrated in \cref{fig:maxcut}, the RCD also outperforms GD in this case, as it converges to an energy ratio of 1 with roughly 200 partial derivative evaluations. In contrast, the GD achieves only an average of 0.75 with 1000 derivative evaluations. The superior performance of RCD in \cref{fig:maxcut} can again be attributed to the significant values of $L/L_avg$ and $L/L_max$,  both exceeding an order of magnitude of 3. As observed from the optimization result, a high ratio of $L/L_avg$ is indicative of the rapid convergence of RCD in this application.

\begin{figure}[htbp]
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/maxcut_training.png}
\end{subfigure}
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/maxcut_lips.png}
\end{subfigure}
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/maxcut_lips_max.png}
\end{subfigure}
\caption{ \label{fig:maxcut} Performance comparison between noisy GD and RCD for the Max-cut problem. The corresponding Lipschitz constant ratios, denoted as $\frac{L}{L_{\text{avg}}}$ and $\frac{L}{L_{\max}}$, are presented in the middle and right panels. The shaded areas within the figure represent variations that have been observed across ten random realizations. The optimization process has been performed in 20 dimensions.}

\end{figure}

\subsubsection{Traveling Salesman Problem (TSP)} \label{sec:tsp}

We have designed a numerical test for the Traveling Salesman Problem (TSP) using three cities as an example. The intercity costs for these cities are 48, 63, and 91 respectively.  The cost of TSP is defined as
\[
C(\mathbf{x})=\sum_{i, j} w_{i j} \sum_{p} x_{i, p} x_{j, p+1}+A \sum_{p}\left(1-\sum_{i} x_{i, p}\right)^{2}+A \sum_{i}\left(1-\sum_{p} x_{i, p}\right)^{2}\,,
\]
where $i$ labels the node, $p$ indicates its order, and $x_i, p$ is in the set ${0, 1}$ and the penalty parameter $A$ is set sufficiently large to effectively enforce constraints. More details regarding the expansion of $C(x)$ can be found in \cref{sec:cost}.

Utilizing the defined cost function, we establish a model Hamiltonian in the same manner as presented in \cref{sec:maxcut}. We aim to prepare its ground state to address the QUBO problem. A detailed representation of the Hamiltonian is available in \cref{sec:cost}. We construct a parameterized quantum circuit comprising alternating single-gate rotations, represented by $U_single (θ)=∏_i=1^n R_Y(θ_i)$ and entangler gate $U_entangler$. This circuit resembles the one depicted in Figure~\ref{fig:circuit_qubo} in Appendix~\ref{sec:circuit_qubo}, albeit with a greater number of qubits. The total number of trainable parameters is 90, which requires nine qubits ($N=9$) and ten  alternating layers. We employ energy as the measure for the optimization cost function.

In the left panel in \cref{fig:tsp}, the optimization results obtained from the noisy RCD and GD are plotted. Notably, GD exhibits slower convergence compared to RCD in achieving an energy ratio of 1. The employment of 90 parameters in the optimization, a number markedly greater than those in prior applications, might account for this disparity. This increased parameter count likely requires additional iterations and partial derivative evaluations when applying GD. Similar to previous results, the two types of  Lipschitz constant ratios are obtained and shown along with the iterations in \cref{fig:tsp}. Again, we can see that the values of the ratios are considerably large, especially during the initial stage of the optimization, underlining the efficiency of RCD in the optimization process.

\begin{figure}[htbp]
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/tsp_training_shift.png}
\end{subfigure}
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/tsp_lips.png}
\end{subfigure}
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/tsp_lips_max.png}
\end{subfigure}
\caption{ \label{fig:tsp} Performance comparison between noisy GD and RCD for the TSP problem. The corresponding Lipschitz constant ratios, denoted as $\frac{L}{L_{\text{avg}}}$ and $\frac{L}{L_{\max}}$, are presented in the middle and right panels. The shaded areas within the figure represent variations that have been observed across ten random realizations. The optimization process has been performed in 90 dimensions. In the first panel, $\tilde E / \tilde E_{GS}$ is defined as $(E - c) / (E_{GS} - c)$, where $c/E_{GS} = 3000$. For clarity in the presentation, we adjust the energy by a constant. }
\end{figure}

\subsubsection{Variational Quantum Factoring}\label{sec:factor}

Our next QUBO problem is designed as a variational quantum factoring task. For this task, we formulated the optimization problem within the framework of quantum adiabatic computation~[19, 3].
For example, to factorize $143$ into the product of two prime numbers, let $143 = pq$, where
\begin{align*}
p &= 8 + 4p_2 + 2p_1 + 1, \\
q &= 8 + 4q_2 + 2q_1 + 1.
\end{align*}
Upon direct computation, the relations are simplified to
\begin{align}
p_1 + q_1 - 1 &= 0, \\
p_2 + q_2 - 1 &= 0, \\
p_2q_1 + p_1q_2 - 1&=0.
\end{align}
To solve this system of equations, we introduce a cost function
\begin{equation}
c(p_1, q_1, p_2, q_2) = (p_1+q_1-1)^2+(p_2+q_2-1)^2+(p_2 q_1+p_1 q_2-1)^2.
\label{eqn:ss_transform}
\end{equation}
By borrowing techniques (see Appendix.~\ref{sec:tech_qfac} for more details) from~[74, 62], the cost function can be reduced to
\begin{equation}
c(p_1, q_1, p_2, q_2) = 5-3 p_1-p_2-q_1+2 p_1 q_1-3 p_2 q_1+2 p_1 p_2 q_1-3 q_2+p_1 q_2+2 p_2 q_2+2 p_2 q_1 q_2.
\end{equation}
Following the methods detailed in the QUBO, we treat $(p_1, q_1, p_2, q_2)$ as boolean functions and substitute each boolean with $1/2 (1 - Z_i)$ as we did in previous sections. Then, the problem can be reformulated into the Ising Hamiltonian,
\begin{equation}
\begin{aligned}
H =& -3 \mathbb{I} + \frac{1}{2} Z_0 + \frac{1}{4} Z_1 + \frac{3}{4} Z_0Z_2 + \frac{1}{4} Z_2 - \frac{1}{4} Z_1Z_2 + \frac{1}{4} Z_0Z_1 \\
&- \frac{1}{4} Z_0Z_1Z_2 + \frac{1}{2} Z_3 + \frac{1}{4} Z_0Z_3 + \frac{3}{4} Z_1Z_3 + \frac{1}{4} Z_2Z_3 - \frac{1}{4} Z_1Z_2Z_3.
\end{aligned}
\label{eqn:factor}
\end{equation}
The ground states of this Hamiltonian are $|0110⟩$ and $|1001⟩$, which respectively correspond to the solutions for the factorization of the number $143$. We summarize it as follows,
\begin{eqnarray}
& (p_{1}, p_{2}, q_{1}, q_{2})=(0,1,1,0) \longleftrightarrow (p,q) = (13,11)\\
& (p_{1}, p_{2}, q_{1}, q_{2})=(1,0,0,1)  \longleftrightarrow (p,q) = (11,13)\\
& p = 8 + 4p_2 + 2p_1 + 1 \text{ and }q = 8 + 4q_2 + 2q_1 + 1\quad \text{Boolean functions}.
\end{eqnarray}

In our numerical experiment, we select the mixer Hamiltonian $H_2 = ∑X_i$ and set up a 20-layer QAOA, which corresponds to 40 parameters~\footnote{The QAOA ansatz builds the variational circuit by alternating between the parametrized unitary evolution associated with the problem Hamiltonian $H$ and the mixer Hamiltonian $H_2$.}. We set the learning rates to $0.0001$ for GD and $0.005$ for RCD and choose the energy as a measure for optimization. Even with a small step size, the variance of GD is notably large. Employing a larger step size for GD further exacerbates the results.

In \cref{fig:factor}, the optimization results of the Hamiltonian \eqref{eqn:factor} are depicted, showing that the number of partial derivative evaluations for the RCD to reach an energy ratio of 1 is about 400 whereas the GD seems to require more than 1000 to the same tolerance. As discussed previously, this observation aligns with prior discussions, particularly given the pronounced magnitude of the Lipschitz constant ratios evident in  \cref{fig:factor}.

\begin{figure}[htbp]
\centering
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/factor_training.png}
\end{subfigure}
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/factor_lips.png}
\end{subfigure}
\begin{subfigure}[b]{0.32\textwidth}
\centering
\includegraphics[width=\textwidth]{figs/factor_lips_max.png}
\end{subfigure}
\caption{ \label{fig:factor}  Performance comparison between noisy GD and RCD for the quantum factoring problem. The corresponding Lipschitz constant ratios, denoted as $\frac{L}{L_{\text{avg}}}$ and $\frac{L}{L_{\max}}$, are presented in the middle and right panels. The shaded areas within the figure represent variations that have been observed across ten random realizations. The optimization process has been performed in 40 dimensions.  }

\end{figure}

\section{Conclusion}

We considered the use of a noisy random coordinate descent method to analyze its potential advantage over the noisy gradient descent, which evaluates all partial derivatives at each step, in the context of variational quantum optimization.
Most previous works on randomized coordinate descent algorithms studied the case of convex cost functions, which do not fit into most variational quantum applications that involve non-convex cost functions. In this work, we generalized the conventional convergence analysis of randomized coordinate descent to local convergence analysis under a local-PL condition that can capture a large class of non-convex optimization. In particular, we proved that noisy randomized coordinate descent can converge faster than noisy gradient descent in terms of the total cost, measured in terms of the total number of partial derivative estimations. In addition, we conducted extensive numerical experiments implementing both methods for many interesting quantum optimization problems. We observed that noisy randomized coordinate descent typically demands less measurement cost than noisy gradient descent, thereby demonstrating its efficiency in many non-convex quantum applications.
From an optimization standpoint, variational quantum optimization as outlined in \cref{problem_1} raises many interesting questions. For example, can second order, or zeroth order optimization methods (i.e., methods using only function evaluation) be more efficient compared to the current gradient-based algorithms? In a technical viewpoint, another question is whether the stability result \cref{lemma: stability of GD} can be generalized so that the event covers the case that the iteration diverges at some timepoint, but it remains in the entire basin until then, $f^-1[0,δ_f)$, not necessarily in the region above the set of global minima, $f^-1(Laσ_∞^2d/μ,δ_f)$. If this is possible to show, then it will provide stronger result such as the stability of the noisy GD and RCD within the entire basin as the stability of Markov Chain in [36].

\section*{Acknowledgements.   } TK and XL's research is supported by the National Science Foundation Grants DMS-2111221 and CCF-2312456. TK is supported by a KIAS Individual Grant CG096001 at Korea Institute for Advanced Study. This material is based upon work supported by the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Quantum Systems Accelerator (ZD).  Additional funding is provided by the SciAI Center funded by the Office of Naval Research (ONR), under Grant Number N00014-23-1-2729 (LL). LL is a Simons investigator in Mathematics.

\appendix

\section{Stochastic stability of noisy GD}\label{sec: appendix B}

In this section, we prove Lemma \ref{lemma: stability of GD}.

\begin{proof}[Proof of Lemma \ref{lemma: stability of GD}]
Define the probability filtration: $\mathcal{F}_n=\sigma\left(\bs{\theta}_k\middle|k\leq n\right)$ and the stopping time\footnote{It is straightforward to see
\[
\{\tau\leq n\}\in \mathcal{F}_n,\quad \{\tau>n\}\in \mathcal{F}_n\,.
\]}
\[
\tau=\inf \left\{k:f(\bs\theta_k)\leq \frac{La\sigma^2_\infty d}{\mu}\right\}\,,
\]
which is the smallest timepoint where the noisy GD achieves $f(\bs\theta_k)\leq\frac{La\sigma^2_\infty d}{\mu}$.

Define the indicator function $\mathbb{I}_n$:
\begin{equation}\label{eq: characteristic ft2}
\mathbb{I}_n=\left\{
\begin{aligned}
&1,\quad \text{if}\quad \{\bs{\theta}_k\}^{n-1}_{k=1}\subset f^{-1}\left(\left[0,\delta_f\right)\right)\\
&0,\quad \text{otherwise}
\end{aligned}\right.\,,
\end{equation}
and the stochastic process
\[
\begin{aligned}
&f(\bs\theta_n)\mathbb{I}_n,\quad n<\tau\\
&f(\bs\theta_\tau)\mathbb{I}_\tau,\quad n\geq\tau
\end{aligned}\right.\,.
\]
According to the definition of $V_n$, there are complementary and exclusive events (cases):
\begin{itemize}
\item Case 1: If there exists $0<n<\infty$ such that: 1. $\bs\theta_n\notin\mathcal{N}$; 2. For any $m<n$, $\bs\theta_m\in\mathcal{N}$ and $f(\bs\theta_m)>\frac{La\sigma^2_\infty d}{\mu}$. Then
\[
V_n\geq\delta_f\Rightarrow \sup_n V_n\geq\delta_f\,.
\]
\item Case 2: For any $n<\tau$, $f(\bs\theta_n)\in\mathcal{N}$.

\end{itemize}
We observe that Case 2 is the stable situation, indicating that $f(\bs\theta_n)$ remains in the basin of the global minimum until it achieves a small loss\footnote{We emphasize that Case 2 also includes the situation where $f(\bs\theta_n)$ remains in the basin and never achieves the small loss.}.
To prove \eqref{eqn:p_gd}, it suffices to show that
\begin{equation}\label{eqn:P_omega_1_bound}
\mathbb{P}(\Omega_1)\leq \frac{f(\bs{\theta}_{1})}{\delta_f},
\end{equation}where $\Omega_1$ denotes the event associated with Case 1.

Now, we show that $V_n$ is a supermartingale to bound $\sup_n V_n$. Taking the conditional expectation, we obtain
\[
\mathbb{E}(V_{n+1}|\mathcal{F}_{n})=\mathbb{E}(V_{n+1}|\mathcal{F}_{n},\mathbb{I}_n=1,\tau\leq n)\mathbb{P}(\tau\leq n)+\mathbb{E}(V_{n+1}|\mathcal{F}_{n},\mathbb{I}_n=1,\tau>n)\mathbb{P}(\tau>n)\,,
\]
where we use $\mathbb{I}_{n+1}\leq \mathbb{I}_n$. There are two terms in the above equation:
\begin{itemize}
\item For the first term, when $\tau\leq n$, we obtain $V_{n+1}=V_\tau=V_n$. This implies
\begin{equation}\label{eqn:first_term_estimation}
\mathbb{E}(V_{n+1}|\mathcal{F}_{n},\mathbb{I}_n=1,\tau\leq n)=V_n\,.
\end{equation}
\item For the second term, when $\tau>n$, we have $f(\bs\theta_{n})> \frac{La\sigma^2_\infty d}{\mu}$. Then, taking the conditional expectation  yields
\begin{equation}\label{eqn:second_term_estimation}
\begin{split}
\leq &f(\bs{\theta}_n)-a\|\nabla f(\bs{\theta}_n)\|^2+\frac{La^2}{2}(\|\nabla f(\bs{\theta}_n)\|^2+\sigma^2_\infty d)\\
\leq &(1-\mu a)f(\bs{\theta}_n)+\frac{La^2\sigma^2_\infty d}{2}\\
< &(1-\mu a)f(\bs{\theta}_n) +\frac{\mu a}{2}f(\bs{\theta}_n)\\
\leq &\left(1-\frac{\mu a}{2}\right)f(\bs{\theta}_n)\mathbb{I}_n=\left(1-\frac{\mu a}{2}\right)V_n,
\end{split}
\end{equation}
where we use assumption \ref{as: local PL} and $a<\frac{1}{L}$ in the second inequality, $\tau>n$ in the third inequality.
\end{itemize}
Combining \eqref{eqn:first_term_estimation} and \eqref{eqn:second_term_estimation}, we obtain
\begin{equation}\label{eqn:iteration}
\mathbb{E}(V_{n+1}|\mathcal{F}_{n})=V_n\mathbb{P}(\tau\leq n)+\left(1-\frac{\mu a}{2}\right)V_n\mathbb{P}(\tau>n)\leq V_n\,.
\end{equation}
Thus, $V_{n}$ is a supermartingale.

Now, we consider the Case 1 event:
\[
\Omega_1=\left\{\exists n>1, \bs\theta_n\notin\mathcal{N}\ \text{and}\ f(\bs\theta_m)>\epsilon\text{ with }\bs\theta_m\in\mathcal{N},\ \forall 1\leq m<n\right\}\subset\left\{\sup_n V_n\geq \delta_f\right\}\,.
\]
Because $V_{n}$ is a supermartingale, we obtain Case 1 happens with small probability:
\[
\mathbb{P}(\Omega_1)\leq \frac{V_1}{\delta_f}=\frac{f(\bs\theta_1)}{\delta_f}\,.
\]
This concludes the proof.
\end{proof}

\section{Stochastic stability of noisy RCD}\label{sec: appendix C}

In this section, we prove \cref{lemma: stability of RCD} with a slight modification of the proof in \cref{sec: appendix B}. From a theoretical viewpoint, the difference between the noisy GD and RCD methods is made by the construction of gradient estimate (e.g. see \eqref{eq:def_GD} and \eqref{eq: def_RCD}). Compared to GD,  the additional randomness of RCD comes with the random selection of a component as in \eqref{eq: def_RCD}. This difference affects the recursive inequality \eqref{eqn:second_term_estimation} mainly in the previous proof, where we considered the properties of the gradient estimator. From this observation, it suffices to derive a recursive inequality similar to \eqref{eqn:second_term_estimation} to prove \cref{lemma: stability of RCD}.

Note that the sampling of a component within RCD is performed before estimating a partial derivative. Thus, the first step is to take expectation on the partial derivative estimate,
\begin{equation}\label{eq: expectation on xi_RCD}
\mathbb{E}_{\xi_{i_n}}[f(\bs{\theta}_{n+1})]\leq f(\bs{\theta}_n)-a\mathbb{E}_{\xi_{i_n}}[\partial_{i_n}f(\bs{\theta}_n)g_{i_n}(\bs{\theta}_n)]+\frac{L_{i_n}a^2}{2}\mathbb{E}_{\xi_{i_n}}[|g_{i_n}(\bs{\theta}_n)|^2]\,,
\end{equation}where $i_n$ is uniformly sampled index and $g_i_n$ is the corresponding unbiased estimate for the partial derivative. Let $ℱ_n$, $τ$, $𝕀_n$, and $V_n$ be as defined in the previous section. By considering the inequality \eqref{eq: expectation on xi_RCD} and the conditional expectation in \eqref{eqn:second_term_estimation}, we achieve the following result by taking expectations with respect to the random index $i_n$,
\begin{equation}\label{eqn:second_term_estimation_RCD}
\begin{split}
\leq &\left(f(\bs{\theta}_n)-\frac{a}{d}\|\nabla f(\bs{\theta}_n)\|^2+\frac{L_{\max}a^2}{2d}\|\nabla f(\bs{\theta}_n)\|^2+\frac{L_{\mathrm{avg}}\sigma^2_\infty da^2}{2d}\right)\mathbb{I}_{n+1}\\
\leq &\left(\left(1-\frac{\mu a}{d}\right)f(\bs{\theta}_n)+\frac{L_{\mathrm{avg}}a^2\sigma^2_\infty }{2}\right)\mathbb{I}_{n+1}\\
< &\left(\left(1-\frac{\mu a}{d}\right)f(\bs{\theta}_n) +\frac{\mu a}{2d}f(\bs{\theta}_n)\right)\mathbb{I}_{n+1}\\
= &\left(1-\frac{\mu a}{2d}\right)f(\bs{\theta}_n)\mathbb{I}_n=\left(1-\frac{\mu a}{2d}\right)V_n,
\end{split}
\end{equation}provided that $f(θ_n)>L_avgaσ_∞^2d/μ$ and $a_n=a<min{1/L_max,  d/μ,  2μδ_f/L_avgσ^2_∞d}$.

Similar to \eqref{eqn:iteration}, in the case of RCD, \eqref{eqn:second_term_estimation_RCD} implies
\begin{equation}\label{eqn:iteration_RCD}
\mathbb{E}(V_{n+1}|\mathcal{F}_{n})=V_n\mathbb{P}(\tau\leq n)+\left(1-\frac{\mu a}{2d}\right)V_n\mathbb{P}(\tau>n)\leq V_n\,,
\end{equation}
which implies $V_n$ forms a supermartingale. The remaining proof of \cref{lemma: stability of RCD} follows the same steps as the proof of \cref{lemma: stability of GD}, so we will not include them here.

\section{The proofs of \cref{thm:GD_noisy} and \cref{thm:RCD_noisy}}\label{sec: appendix D}

We first show the convergence rate of the noisy GD method, followed by a similar analysis for the noisy RCD method. The following proofs are similar to those in \cref{sec: appendix B} and \cref{sec: appendix C} with minor differences.

\begin{proof}[Proof of Theorem \ref{thm:GD_noisy}]
Define the probability filtration: $\mathcal{F}_n=\sigma\left(\bs{\theta}_k\middle|k\leq n\right)$ and the stopping time
\[
\tau=\inf \left\{k:f(\bs\theta_k)\leq \epsilon\right\}\,,
\]
which is the smallest timepoint where the noisy GD achieves $f(\bs\theta_k)\leq\epsilon$. Then, our ultimate goal is to show that the event that $\inf_{1\leq n\leq N}f(\bs\theta_n)\leq \epsilon$ occurs with high probability, say, at least $1-\eta$.
Then, our goal is to show that for any $\eta\in\left(\frac{f(\bs\theta_1)}{\delta_f},1\right)$, there exists a sufficiently large $N$ such that
\begin{equation}\label{pfail}
p_{\mathrm{fail}}:=\mathbb{P}(\tau>N)\leq \eta.
\end{equation}
Define the indicator function $\mathbb{I}_n$:
\[
\mathbb{I}_n=\left\{
\begin{aligned}
&1,\quad \text{if}\quad \{\bs{\theta}_k\}^{n-1}_{k=1}\subset f^{-1}\left(\left[0,\delta_f\right)\right)\\
&0,\quad \text{otherwise}
\end{aligned}\right.\,,
\]
and the stochastic process
\[
\begin{aligned}
&f(\bs\theta_n)\mathbb{I}_n,\quad n<\tau\\
&f(\bs\theta_\tau)\mathbb{I}_\tau,\quad n\geq\tau
\end{aligned}\right.\,.
\]
Define the unstable event:
\[
\Omega=\left\{\exists n>1, \bs\theta_n\notin\mathcal{N}\ \text{and}\ f(\bs\theta_m)>\epsilon,\ \forall 1\leq m<n\right\}\subset\left\{\sup_n V_n\geq \delta_f\right\}\,.
\]
According to Lemma \ref{lemma: stability of GD} and the proof in \cref{sec: appendix B}, for learning rate $a$ with $\frac{La\sigma^2_\infty d}{\mu}<\epsilon$, we obtain $\Omega$ happens with small probability:
\begin{equation}\label{eqn:failure_1}
\mathbb{P}(\Omega)\leq \frac{V_1}{\delta_f}=\frac{f(\bs\theta_1)}{\delta_f}\,.
\end{equation}
Recalling \eqref{pfail}, we note that, for any $n\leq N$,
\[
\mathbb{P}(\tau>n)\geq p_{\mathrm{fail}}\,.
\]
Plugging this into \eqref{eqn:iteration}, we obtain that
\begin{equation}
\begin{split}
\mathbb{E}(V_{n+1}|\mathcal{F}_n)&= \left(1-\mathbb{P}(\tau>n)+\left(1-\frac{\mu a}{2}\right)\mathbb{P}(\tau>n)\right)V_n\\
& = \left(1-\frac{\mu a \mathbb{P}(\tau>n)}{2}\right)V_n\\
&\leq \left(1-\frac{\mu a p_{\mathrm{fail}}}{2}\right)V_n.
\end{split}
\end{equation}By taking the total expectation on both sides and using a telescoping trick, we achieve that
\begin{equation}
\mathbb{E}(V_{n+1})\leq \left(1-\frac{\mu a p_{\mathrm{fail}}}{2}\right)^n V_1=\left(1-\frac{\mu a p_{\mathrm{fail}}}{2}\right)^n f(\bs\theta_1).
\end{equation}
This means that if the probability of failure,
$p_{\mathrm{fail}}$, is large, the expectation of $V_{n+1}$ decreases quickly. By Markov's inequality, we have
\[
\mathbb{P}\left(V_N> \epsilon\right)\leq \frac{\left(1-\frac{\mu a p_{\mathrm{fail}}}{2}\right)^{N-1} f(\bs\theta_1)}{\epsilon}\,,
\]
\[
\mathbb{P}\left(V_N\leq \epsilon\right)\geq 1-\frac{\left(1-\frac{\mu a p_{\mathrm{fail}}}{2}\right)^{N-1} f(\bs\theta_1)}{\epsilon}\,.
\]
Now, if we consider the event $\{V_N\leq \epsilon\}$, it is the union of the following two events (not necessarily exclusive and complementary), which are slightly different from the ones in \cref{sec: appendix B}:
\begin{itemize}
\item $\Omega_1$: There exists $n\leq N$ such that $f(\bs\theta_n)\leq \epsilon$ and $\bs\theta_n\in\mathcal{N}$. This means \[\inf_{1\leq n\leq N}f(\bs\theta_n)\leq \epsilon\,.\]

We want to show that $\Omega_1$ happens with high probability.

\item $\Omega_2$: There exists $n< N$ such that $f(\bs \theta_n)>\delta_f$ and $f(\bs\theta_m)>\epsilon$ for any $m<n$.

We note that, when $\Omega_2$ happens, we have $V_{n+1}=0$ with $f(\bs\theta_n)>\delta_f$, which implies $\Omega_2\subset \Omega$. According to \eqref{eqn:failure_1}, we obtain
\[
\mathbb{P}\left(\Omega_2\right)\leq \mathbb{P}\left(\Omega\right)\leq \frac{f(\bs\theta_1)}{\delta_f}\,.
\]
\end{itemize}
Now, we give a lower bound for the event $\Omega_1$:
\begin{equation}\label{lowerbound}
\mathbb{P}\left(\inf_{1\leq n\leq N}f(\bs\theta_n)\leq \epsilon\right)=\mathbb{P}\left(\Omega_1\right)\geq \mathbb{P}\left(V_N\leq \epsilon\right)-\mathbb{P}(\Omega_2)\geq 1-\frac{\left(1-\frac{\mu a p_{\mathrm{fail}}}{2}\right)^N f(\bs\theta_1)}{\epsilon}-\frac{f(\bs\theta_1)}{\delta_f}.
\end{equation}Notice
\[
\mathbb{P}\left(\inf_{1\leq n\leq N}f(\bs\theta_n)\leq \epsilon\right)\leq \mathbb{P}(\tau\leq N)=1-p_{\mathrm{fail}}\,.
\]
Combining the above two inequalities, we have
\begin{equation}\label{eqn:must_hold}
p_{\mathrm{fail}}\leq \frac{\left(1-\frac{\mu a p_{\mathrm{fail}}}{2}\right)^N f(\bs\theta_1)}{\epsilon}+\frac{f(\bs\theta_1)}{\delta_f}\,.
\end{equation}
Next, we show \eqref{pfail} using the proof by contradiction. Assume that the conclusion of the theorem is not true, meaning that for some $\eta\in\left(\frac{f(\bs\theta_1)}{\delta_f},1\right)$ and every $N$,
\[
p_{\mathrm{fail}}> \eta\,.
\]
When $p_{\mathrm{fail}}>\eta$ and $N=\frac{2}{\mu a\eta}\log\left(\frac{f(\bs{\theta_1})}{\left(\eta-\frac{f(\bs\theta_1)}{\delta_f}\right)\epsilon}\right)$, then
\[
\begin{aligned}
&\frac{\left(1-\frac{\mu a p_{\mathrm{fail}}}{2}\right)^N f(\bs\theta_1)}{\epsilon}+\frac{f(\bs\theta_1)}{\delta_f}\\
<&\frac{\left(1-\frac{\mu a \eta}{2}\right)^N f(\bs\theta_1)}{\epsilon}+\frac{f(\bs\theta_1)}{\delta_f}\\
\leq&\frac{\exp\left(-\frac{\mu a \eta N}{2}\right) f(\bs\theta_1)}{\epsilon}+\frac{f(\bs\theta_1)}{\delta_f}=\eta<p_{\mathrm{fail}}
\end{aligned}
\]
where we use $p_{\mathrm{fail}}>\eta$ in the first inequality and $(1-x)^N\leq \exp(-xN)$ in the second inequality. This contradicts to \eqref{eqn:must_hold}. Thus, \eqref{pfail} must be true and we conclude the proof.

\end{proof}

\begin{proof}[Proof of Theorem \ref{thm:RCD_noisy}]
Denote the probability of failure
\[
\]

Similar to the calculation in the previous proof, from \eqref{eqn:iteration_RCD}, we have
\begin{equation}\label{lowerbound_RCD}
\mathbb{P}\left(\inf_{1\leq n\leq N}f(\bs\theta_n)\leq \epsilon\right)\geq 1-\frac{\left(1-\frac{\mu a p_{\mathrm{fail}}}{2d}\right)^N f(\bs\theta_1)}{\epsilon}-\frac{f(\bs\theta_1)}{\delta_f}.
\end{equation}
With the same logic below \eqref{lowerbound}, we conclude the proof of \cref{thm:RCD_noisy}.
\end{proof}

\newpage

\section{Parameterized Circuit for the VQE}\label{sec:circuit}

\begin{figure}[htbp]
\centering
\Qcircuit @C=0.9em @R=.7em {
& \gate{R_Y(\frac{3\pi}{2})} & \multigate{1}{R_{ZZ}(\theta_1)} & \qw
& \gate{R_X(\theta_2)} & \qw &
\cdots &   & \multigate{1}{R_{ZZ}(\theta_{39})} & \qw
& \gate{R_X(\theta_{40})} & \qw \\
& \gate{R_Y(\frac{3\pi}{2})} & \ghost{R_{ZZ}(\theta_1)} &  \multigate{1}{R_{ZZ}(\theta_1)} & \gate{R_X(\theta_2)} & \qw & \cdots &  & \ghost{R_{ZZ}(\theta_{39})} &  \multigate{1}{R_{ZZ}(\theta_{39})} & \gate{R_X(\theta_{40})} & \qw \\
& \gate{R_Y(\frac{3\pi}{2})} & \multigate{1}{R_{ZZ}(\theta_1)} & \ghost{R_{ZZ}(\theta_1)} & \gate{R_X(\theta_2)} & \qw & \cdots &  & \multigate{1}{R_{ZZ}(\theta_{39})} & \ghost{R_{ZZ}(\theta_{39})} & \gate{R_X(\theta_{40})} & \qw\\
& \gate{R_Y(\frac{3\pi}{2})} & \ghost{R_{ZZ}(\theta_1)} & \multigate{1}{R_{ZZ}(\theta_1)} & \gate{R_X(\theta_2)} & \qw & \cdots &  & \ghost{R_{ZZ}(\theta_{39})} & \multigate{1}{R_{ZZ}(\theta_{39})} & \gate{R_X(\theta_{40})} & \qw\\
& \gate{R_Y(\frac{3\pi}{2})} & \multigate{1}{R_{ZZ}(\theta_1)} & \ghost{R_{ZZ}(\theta_1)} & \gate{R_X(\theta_2)} & \qw & \cdots &  & \multigate{1}{R_{ZZ}(\theta_{39})} & \ghost{R_{ZZ}(\theta_{39})} & \gate{R_X(\theta_{40})} & \qw\\
& \gate{R_Y(\frac{3\pi}{2})} & \ghost{R_{ZZ}(\theta_1)} & \multigate{1}{R_{ZZ}(\theta_1)} & \gate{R_X(\theta_2)} & \qw & \cdots &  & \ghost{R_{ZZ}(\theta_{39})} & \multigate{1}{R_{ZZ}(\theta_{39})} & \gate{R_X(\theta_{40})} & \qw\\
& \gate{R_Y(\frac{3\pi}{2})} & \multigate{1}{R_{ZZ}(\theta_1)} & \ghost{R_{ZZ}(\theta_1)} & \gate{R_X(\theta_2)} & \qw & \cdots & & \multigate{1}{R_{ZZ}(\theta_{39})} & \ghost{R_{ZZ}(\theta_{39})} & \gate{R_X(\theta_{40})} & \qw\\
& \gate{R_Y(\frac{3\pi}{2})} & \ghost{R_{ZZ}(\theta_1)} & \multigate{1}{R_{ZZ}(\theta_1)} & \gate{R_X(\theta_2)} & \qw & \cdots &  & \ghost{R_{ZZ}(\theta_{39})} & \multigate{1}{R_{ZZ}(\theta_{39})} & \gate{R_X(\theta_{40})}
& \qw \\
& \gate{R_Y(\frac{3\pi}{2})} & \multigate{1}{R_{ZZ}(\theta_1)} & \ghost{R_{ZZ}(\theta_1)} & \gate{R_X(\theta_2)} & \qw & \cdots & & \multigate{1}{R_{ZZ}(\theta_{39})} & \ghost{R_{ZZ}(\theta_{39})} & \gate{R_X(\theta_{40})} & \qw \\
& \gate{R_Y(\frac{3\pi}{2})} & \ghost{R_{ZZ}(\theta_1)} & \qw & \gate{R_X(\theta_2)} & \qw & \cdots &  & \ghost{R_{ZZ}(\theta_{39})} & \qw & \gate{R_X(\theta_{40})}
& \qw }
\caption{A QAOA-like ansatz motivated by [72] is used for the TFIM model \eqref{eq:isingmodel} with $10$ qubits.
For the result in \cref{fig:vqe}, 40 parameters are assigned with 20 layers of alternating rotation ZZ gates and rotation X gates. The circuit depth is 61.}
\label{fig:circuit}
\end{figure}

\newpage

\section{Parameterized Circuit for the VQE in QUBO experiments}\label{sec:circuit_qubo}

The quantum circuit described below is utilized in the QUBO experiments.

\begin{figure}[htbp]
\centering
\Qcircuit @C=1em @R=1em {
& \gate{R_Y(\theta_1)} & \ctrl{1} & \qw & \qw & \qw  & \gate{R_Y(\theta_5)} & \ctrl{1} & \qw & \qw & \qw  & \cdots & & \gate{R_Y(\theta_{n-3})} & \ctrl{1} & \qw & \qw & \qw  \\
& \gate{R_Y(\theta_2)} & \gate{Z} & \ctrl{1} & \qw & \qw  & \gate{R_Y(\theta_6)} & \gate{Z} & \ctrl{1} & \qw & \qw  & \cdots & & \gate{R_Y(\theta_{n-2})} & \gate{Z} & \ctrl{1} & \qw & \qw  \\
& \gate{R_Y(\theta_3)} & \qw & \gate{Z} & \ctrl{1} & \qw  & \gate{R_Y(\theta_7)} & \qw & \gate{Z} & \ctrl{1} & \qw  & \cdots & & \gate{R_Y(\theta_{n-1})} & \qw & \gate{Z} & \ctrl{1} & \qw  \\
& \gate{R_Y(\theta_4)} & \qw & \qw & \gate{Z} & \qw  & \gate{R_Y(\theta_8)} & \qw & \qw & \gate{Z} & \qw  & \cdots & & \gate{R_Y(\theta_n)} & \qw & \qw & \gate{Z} & \qw  \\

\caption{A parametrized quantum circuit is employed in the QUBO experiments. This circuit features alternating layers of single rotation gates and entangling controlled-z gates. The adjustable parameters are exclusively found in the single rotation gates, and these parameters vary across different layers and qubits.}
\label{fig:circuit_qubo}
\end{figure}

\section{Additional histograms of partial derivative estimates}\label{sec:noise}

\cref{fig:vqe_noise} plots the histograms with respect to the first 12 parameters among 40. The rest of 28 histograms are shown in the following figures. It is observed in all figures that the variances of partial derivative estimates in all directions are a similar magnitude of value.

\begin{figure}[htbp]
\centering
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d13.jpg}
\end{subfigure}
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d14.jpg}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d15.jpg}
\end{subfigure}

\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d16.jpg}
\end{subfigure}
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d17.jpg}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d18.jpg}
\end{subfigure}

\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d19.jpg}
\end{subfigure}
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d20.jpg}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d21.jpg}
\end{subfigure}

\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d22.jpg}
\end{subfigure}
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d23.jpg}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d24.jpg}
\end{subfigure}

\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d25.jpg}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d26.jpg}
\end{subfigure}

\caption{The histogram of partial derivative estimates with respect to the 13-th to the 26-th parameters are plotted with the same setup as in \cref{fig:vqe_noise}.}
\end{figure}

\begin{figure}[htbp]
\centering
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d27.jpg}
\end{subfigure}
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d28.jpg}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d29.jpg}
\end{subfigure}

\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d30.jpg}
\end{subfigure}
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d31.jpg}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d32.jpg}
\end{subfigure}

\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d33.jpg}
\end{subfigure}
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d34.jpg}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d35.jpg}
\end{subfigure}

\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d36.jpg}
\end{subfigure}
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d37.jpg}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d38.jpg}
\end{subfigure}

\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d39.jpg}
\end{subfigure}
\centering
\begin{subfigure}[b]{0.3\textwidth}
\centering
\includegraphics[width=\textwidth]{ExperimentData/noise/d40.jpg}
\end{subfigure}

\caption{The histogram of partial derivative estimates with respect to the 27-th to the 40-th parameters are plotted with the same setup as in \cref{fig:vqe_noise}.}
\end{figure}

\section{Cost function for the TSP}\label{sec:cost}

First, the cost function is defined as

\[
C(\mathbf{x})=\sum_{i, j} w_{i j} \sum_{p} x_{i, p} x_{j, p+1}+A \sum_{p}\left(1-\sum_{i} x_{i, p}\right)^{2}+A \sum_{i}\left(1-\sum_{p} x_{i, p}\right)^{2}\,,
\]
where $A=10000$ , \( w_{12} = w_{21} = 48 \), \( w_{13} = w_{31} = 91 \), and \( w_{23} = w_{32} = 63 \), and \(w_{ii} = 0, i = 1, 2, 3\).

We can introduce a new Boolean variable, denoted by \( \tilde{x}_{3i + j - 4} = x_{i,j} \), where \( i, j = 1, 2, 3 \). For simplicity, in the following formula, we will use \( x_0, \ldots, x_8 \) to represent \( \tilde{x}_0, \ldots, \tilde{x}_8 \). With this notation, the expanded form of the cost function can be expressed as:

\begin{equation*}
\begin{split}
C({\bm x}) = -200000 x_0 & - 200000 x_1 - 200000 x_2 - 200000 x_3 - 200000 x_4 - 200000 x_5 \\
& - 200000 x_6 - 200000 x_7 - 200000 x_8 \\
& + \left[ 200000 x_0x_1 + 200000 x_0x_2 + 200000 x_0x_3 + 48 x_0x_4 + 48 x_0x_5 \right.\\
& + 200000 x_0x_6 + 91 x_0x_7 + 91 x_0x_8 + 200000 x_1x_2 + 48 x_1x_3 \\
& + 200000 x_1x_4 + 48 x_1x_5 + 91 x_1x_6 + 200000 x_1x_7 + 91 x_1x_8 \\
& + 48 x_2x_3 + 48 x_2x_4 + 200000 x_2x_5 + 91 x_2x_6 + 91 x_2x_7 \\
& + 200000 x_2x_8 + 200000 x_3x_4 + 200000 x_3x_5 + 200000 x_3x_6 \\
& + 63 x_3x_7 + 63 x_3x_8 + 200000 x_4x_5 + 63 x_4x_6 + 200000 x_4x_7 \\
& + 63 x_4x_8 + 63 x_5x_6 + 63 x_5x_7 + 200000 x_5x_8 + 200000 x_6x_7 \\
& \left.+ 200000 x_6x_8 + 200000 x_7x_8 \right] + 600000
\end{split}
\end{equation*}

In order to build the corresponding Hamiltonian, we align the binary variables \( x_i \) with the Pauli \( Z \) matrices, which operate on individual qubits, and are represented by \( Z_i \). Taking into account the relationship between the binary variables \( x_i \) and the Pauli \( Z \) matrices, defined by the equation \( x_i = \frac{{1 - Z_i}}{2} \), we can express the Hamiltonian for QUBO as follows,
\begin{equation*}
\begin{split}
H_{\text{TSP}} = 600303.0 &-100069.5 Z_0 -100055.5 Z_4 + 12.0 Z_4 Z_0 -100069.5 Z_1\\
&-100055.5 Z_5 + 12.0 Z_5 Z_1 -100069.5 Z_2 -100055.5 Z_3 + 12.0 Z_3 Z_2\\
&-100077.0 Z_7 + 22.75 Z_7 Z_0 -100077.0 Z_8 + 22.75 Z_8 Z_1\\
&-100077.0 Z_6 + 22.75 Z_6 Z_2 + 12.0 Z_3 Z_1 + 12.0 Z_4 Z_2\\
&+ 12.0 Z_5 Z_0 + 15.75 Z_7 Z_3 + 15.75 Z_8 Z_4 + 15.75 Z_6 Z_5\\
&+ 22.75 Z_6 Z_1 + 22.75 Z_7 Z_2 + 22.75 Z_8 Z_0\\
&+ 15.75 Z_6 Z_4 + 15.75 Z_7 Z_5 + 15.75 Z_8 Z_3\\
&+ 50000.0 Z_3 Z_0 + 50000.0 Z_6 Z_0 + 50000.0 Z_6 Z_3\\
&+ 50000.0 Z_4 Z_1 + 50000.0 Z_7 Z_1 + 50000.0 Z_7 Z_4\\
&+ 50000.0 Z_5 Z_2 + 50000.0 Z_8 Z_2 + 50000.0 Z_8 Z_5\\
&+ 50000.0 Z_1 Z_0 + 50000.0 Z_2 Z_0 + 50000.0 Z_2 Z_1\\
&+ 50000.0 Z_4 Z_3 + 50000.0 Z_5 Z_3 + 50000.0 Z_5 Z_4\\
&+ 50000.0 Z_7 Z_6 + 50000.0 Z_8 Z_6 + 50000.0 Z_8 Z_7
\end{split}
\end{equation*}

\section{Technique used in quantum factoring}\label{sec:tech_qfac}

The introduced technique proposes an alternative formulation for equations of the type \(A B + S = 0\). Here, \(A\) and \(B\) represent Boolean variables, while \(S\) denotes integers with \(S \in \mathbb{Z}\). The optimization algorithm targets the minimization of the quadratic version of this equation.

Given the problem Hamiltonian, defined as \( H = (A B + S)^2 \), it can be restructured as:
\begin{equation}
H = 2 \left[ \frac{1}{2} \left( A + B - \frac{1}{2} \right) + S \right]^2 - \frac{1}{8}.
\end{equation}

While the two Hamiltonians are not generally equivalent, they do share the same minimizer due to their underlying Boolean function properties. For instance:

\begin{itemize}
\item When \(AB=1\): The minimizer for the first Hamiltonian dictates \(S=-1\). In the reformulated version, we get
\begin{align*}
H &= 2 \left[ \frac{1}{2} \left( 1 + 1 - \frac{1}{2} \right) - 1 \right]^2 - \frac{1}{8} \\
&= 0.
\end{align*}

\item When \(AB=0\): At least one of \(A\) or \(B\) is zero. Assuming \(A=0\) (without loss of generality) and due to the minimizer, we get \(S=0\). This also minimizes the reformulated Hamiltonian since, regardless of whether \(B\) is 0 or 1, the result remains 0.
\end{itemize}

Thus, the reformulated version can be employed interchangeably in certain scenarios. However, this updated representation leads to a significant reduction in the many-body interactions observed experimentally. Specifically, the quartic terms in the Ising Hamiltonian are eliminated, simplifying experimental realizations.
As a result, the third Hamiltonian term $(p_2q_1 + p_1q_2 - 1 ) ^ 2$ in Eqn.(\ref{eqn:ss_transform}) is reformulated as:
\[ H^{\prime} = 2 \left[ \frac{1}{2} \left( p_{1} + q_{2} - \frac{1}{2} \right) + p_{2} q_{1} - 1 \right]^2 - \frac{1}{8}. \]

\bibliographystyle{plain}
\begin{thebibliography}{10}

[1]
Amira Abbas, Robbie King, Hsin-Yuan Huang, William~J Huggins, Ramis Movassagh,
Dar Gilboa, and Jarrod~R McClean.
\newblock On quantum backpropagation, information reuse, and cheating
measurement collapse.
\newblock {\em arXiv preprint arXiv:2305.13362}, 2023.

[2]
Abhinav Anand, Matthias Degroote, and Alán Aspuru-Guzik.
\newblock Natural evolutionary strategies for variational quantum computation.
\newblock {\em Machine Learning: Science and Technology}, 2(4):045012, jul

[3]
Eric Anschuetz, Jonathan Olson, Al{\'a}n Aspuru-Guzik, and Yudong Cao.
\newblock Variational quantum factoring.
\newblock In {\em Quantum Technology and Optimization Problems: First
International Workshop, QTOP 2019, Munich, Germany, March 18, 2019,
Proceedings 1}, pages 74--85. Springer, 2019.

[4]
Martin Anthony, Endre Boros, Yves Crama, and Aritanan Gruber.
\newblock Quadratic reformulations of nonlinear binary optimization problems.
\newblock {\em Mathematical Programming}, 162:115--144, 2017.

[5]
Leonardo Banchi and Gavin~E Crooks.
\newblock Measuring analytic gradients of general quantum evolution with the
stochastic parameter shift rule.
\newblock {\em Quantum}, 5:386, 2021.

[6]
Raef Bassily, Mikhail Belkin, and Siyuan Ma.
\newblock On exponential convergence of sgd in non-convex over-parametrized
\newblock {\em arXiv preprint arXiv:1811.02564}, 2018.

[7]
Marcello Benedetti, Erika Lloyd, Stefan Sack, and Mattia Fiorentini.
\newblock Parameterized quantum circuits as machine learning models.
\newblock {\em Quantum Science and Technology}, 4(4):043001, 2019.

[8]
Jacob Biamonte, Peter Wittek, Nicola Pancotti, Patrick Rebentrost, Nathan
Wiebe, and Seth Lloyd.
\newblock Quantum machine learning.
\newblock {\em Nature}, 549(7671):195--202, 2017.

[9]
L{\'e}on Bottou, Frank~E Curtis, and Jorge Nocedal.
\newblock Optimization methods for large-scale machine learning.
\newblock {\em Siam Review}, 60(2):223--311, 2018.

[10]
James Bradbury, Roy Frostig, Peter Hawkins, Matthew~James Johnson, Chris Leary,
Dougal Maclaurin, George Necula, Adam Paszke, Jake Vander{P}las, Skye
Wanderman-{M}ilne, and Qiao Zhang.
\newblock {JAX}: composable transformations of {P}ython+{N}um{P}y programs,

[11]
Marin Bukov.
\newblock Reinforcement learning for autonomous preparation of
floquet-engineered states: Inverting the quantum kapitza oscillator.
\newblock {\em Phys. Rev. B}, 98:224305, Dec 2018.

[12]
Marin Bukov, Alexandre G.~R. Day, Dries Sels, Phillip Weinberg, Anatoli
Polkovnikov, and Pankaj Mehta.
\newblock Reinforcement learning in different phases of quantum control.
\newblock {\em Phys. Rev. X}, 8:031086, Sep 2018.

[13]
Marin Bukov, Alexandre~GR Day, Dries Sels, Phillip Weinberg, Anatoli
Polkovnikov, and Pankaj Mehta.
\newblock Reinforcement learning in different phases of quantum control.
\newblock {\em Physical Review X}, 8(3):031086, 2018.

[14]
M~Cerezo, Guillaume Verdon, Hsin-Yuan Huang, Lukasz Cincio, and Patrick~J
\newblock Challenges and opportunities in quantum machine learning.
\newblock {\em Nature Computational Science}, 2(9):567--576, 2022.

[15]
Marco Cerezo, Andrew Arrasmith, Ryan Babbush, Simon~C Benjamin, Suguru Endo,
Keisuke Fujii, Jarrod~R McClean, Kosuke Mitarai, Xiao Yuan, Lukasz Cincio,
\newblock Variational quantum algorithms.
\newblock {\em Nature Reviews Physics}, 3(9):625--644, 2021.

[16]
Wei Chen, Kazuo Hida, and B.~C. Sanctuary.
\newblock Ground-state phase diagram of $s=1$ $\mathrm{XXZ}$ chains with
uniaxial single-ion-type anisotropy.
\newblock {\em Phys. Rev. B}, 67:104401, Mar 2003.

[17]
Ziang Chen, Yingzhou Li, and Jianfeng Lu.
\newblock On the global convergence of randomized coordinate gradient descent
for nonconvex optimization.
\newblock {\em SIAM Journal on Optimization}, 33(2):713--738, 2023.

[18]
Gavin~E Crooks.
\newblock Gradients of parameterized quantum gates using the parameter-shift
rule and gate decomposition.
\newblock {\em arXiv preprint arXiv:1905.13311}, 2019.

[19]
Nikesh~S Dattani and Nathaniel Bryans.
\newblock Quantum factorization of 56153 with only 4 qubits.
\newblock {\em arXiv preprint arXiv:1411.6758}, 2014.

[20]
Pablo D\'{\i}ez-Valle, Diego Porras, and Juan~Jos\'e Garc\'{\i}a-Ripoll.
\newblock Quantum variational optimization: The role of entanglement and
problem hardness.
\newblock {\em Phys. Rev. A}, 104:062426, Dec 2021.

[21]
Suguru Endo, Zhenyu Cai, Simon~C Benjamin, and Xiao Yuan.
\newblock Hybrid quantum-classical algorithms and quantum error mitigation.
\newblock {\em Journal of the Physical Society of Japan}, 90(3):032001, 2021.

[22]
Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.
\newblock A quantum approximate optimization algorithm.
\newblock {\em arXiv preprint arXiv:1411.4028v1}, Nov 2014.

[23]
Benjamin Fehrman, Benjamin Gess, and Arnulf Jentzen.
\newblock Convergence rates for the stochastic gradient descent method for
non-convex objective functions.
\newblock {\em The Journal of Machine Learning Research}, 21(1):5354--5401,

[24]
Thomas F{\"o}sel, Murphy~Yuezhen Niu, Florian Marquardt, and Li~Li.
\newblock Quantum circuit optimization with deep reinforcement learning.
\newblock {\em arXiv preprint arXiv:2103.07585}, 2021.

[25]
Pan Gao, Keren Li, Shijie Wei, Jiancun Gao, and Guilu Long.
\newblock Quantum gradient algorithm for general polynomials.
\newblock {\em Physical Review A}, 103(4):042403, 2021.

[26]
Andr{\'a}s Gily{\'e}n, Srinivasan Arunachalam, and Nathan Wiebe.
\newblock Optimizing quantum optimization algorithms via faster quantum
gradient computation.
\newblock In {\em Proceedings of the Thirtieth Annual ACM-SIAM Symposium on
Discrete Algorithms}, pages 1425--1444. SIAM, 2019.

[27]
Alessandro Giovagnoli, Yunpu Ma, and Volker Tresp.
\newblock Qneat: Natural evolution of variational quantum circuit architecture.
\newblock {\em arXiv preprint arXiv:2304.06981v1}, Apr 2023.

[28]
Fred Glover, Gary Kochenberger, and Yu~Du.
\newblock A tutorial on formulating and using qubo models.
\newblock {\em arXiv preprint arXiv:1811.11538v6}, Nov 2018.

[29]
Charles~R. Harris, K.~Jarrod Millman, St{\'{e}}fan~J. van~der Walt, Ralf
Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor,
Sebastian Berg, Nathaniel~J. Smith, Robert Kern, Matti Picus, Stephan Hoyer,
Marten~H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime~Fern{\'{a}}ndez
del R{\'{i}}o, Mark Wiebe, Pearu Peterson, Pierre G{\'{e}}rard-Marchant,
Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph
Gohlke, and Travis~E. Oliphant.
\newblock Array programming with {NumPy}.
\newblock {\em Nature}, 585(7825):357--362, September 2020.

[30]
Chi Jin, Praneeth Netrapalli, Rong Ge, Sham~M Kakade, and Michael~I Jordan.
\newblock On nonconvex optimization for machine learning: Gradients,
stochasticity, and saddle points.
\newblock {\em Journal of the ACM (JACM)}, 68(2):1--29, 2021.

[31]
Stephen~P Jordan.
\newblock Fast quantum algorithm for numerical gradient estimation.
\newblock {\em Physical Review Letters}, 95(5):050501, 2005.

[32]
Abhinav Kandala, Antonio Mezzacapo, Kristan Temme, Maika Takita, Markus Brink,
Jerry~M. Chow, and Jay~M. Gambetta.
\newblock Hardware-efficient variational quantum eigensolver for small
molecules and quantum magnets.
\newblock {\em Nature}, 549(7671):242--246, 2017.

[33]
Bobby Kleinberg, Yuanzhi Li, and Yang Yuan.
\newblock An alternative view: When does sgd escape local minima?
\newblock In {\em International conference on machine learning}, pages
2698--2707. PMLR, 2018.

[34]
Taehee Ko and Xiantao Li.
\newblock A local convergence theory for the stochastic gradient descent method
in non-convex optimization with non-isolated local minima.
\newblock {\em Journal of Machine Learning}, 2(2):138--160, 2023.

[35]
Gary Kochenberger, Jin-Kao Hao, Fred Glover, Mark Lewis, Zhipeng L{\"u}, Haibo
Wang, and Yang Wang.
\newblock The unconstrained binary quadratic programming problem: a survey.
\newblock {\em Journal of Combinatorial Optimization}, 28(1):58--81, Jul 2014.

[36]
Harold~J Kushner and G~George Yin.
\newblock Applications to learning, state dependent noise, and queueing.
\newblock {\em Stochastic Approximation Algorithms and Applications}, pages
25--46, 1997.

[37]
A~Langari, F~Pollmann, and M~Siahatgar.
\newblock Ground-state fidelity of the spin-1 heisenberg chain with single ion
anisotropy: quantum renormalization group and exact diagonalization
\newblock {\em Journal of Physics: Condensed Matter}, 25(40):406002, sep 2013.

[38]
Martin Larocca, Nathan Ju, Diego Garc{\'\i}a-Mart{\'\i}n, Patrick~J Coles, and
Marco Cerezo.
\newblock Theory of overparametrization in quantum neural networks.
\newblock {\em arXiv preprint arXiv:2109.11676}, 2021.

[39]
Yin~Tat Lee and Aaron Sidford.
\newblock Efficient accelerated coordinate descent methods and faster
algorithms for solving linear systems.
\newblock In {\em 2013 ieee 54th annual symposium on foundations of computer
science}, pages 147--156. IEEE, 2013.

[40]
Jiaqi Leng, Yuxiang Peng, Yi-Ling Qiao, Ming Lin, and Xiaodi Wu.
\newblock Differentiable analog quantum computing for optimization and control.
\newblock {\em arXiv preprint arXiv:2210.15812}, 2022.

[41]
Chaoyue Liu, Libin Zhu, and Mikhail Belkin.
\newblock Loss landscapes and optimization in over-parameterized non-linear
systems and neural networks.
\newblock {\em Applied and Computational Harmonic Analysis}, 59:85--116, 2022.

[42]
Andrew Lucas.
\newblock Ising formulations of many np problems.
\newblock {\em Frontiers in physics}, 2:5, 2014.

[43]
Jarrod~R McClean, Jonathan Romero, Ryan Babbush, and Al{\'a}n Aspuru-Guzik.
\newblock The theory of variational hybrid quantum-classical algorithms.
\newblock {\em New Journal of Physics}, 18(2):023023, 2016.

[44]
Fan-Xu Meng, Ze-Tong Li, Xu-Tao Yu, and Zai-Chen Zhang.
\newblock Quantum circuit architecture optimization for variational quantum
eigensolver via monto carlo tree search.
\newblock {\em IEEE Transactions on Quantum Engineering}, 2:1--10, 2021.

[45]
Panayotis Mertikopoulos, Nadav Hallak, Ali Kavis, and Volkan Cevher.
\newblock On the almost sure convergence of stochastic gradient descent in
non-convex problems.
\newblock {\em Advances in Neural Information Processing Systems},
33:1117--1128, 2020.

[46]
Nikolaj Moll, Panagiotis Barkoutsos, Lev~S Bishop, Jerry~M Chow, Andrew Cross,
Daniel~J Egger, Stefan Filipp, Andreas Fuhrer, Jay~M Gambetta, Marc Ganzhorn,
\newblock Quantum optimization using variational algorithms on near-term
quantum devices.
\newblock {\em Quantum Science and Technology}, 3(3):030503, 2018.

[47]
Eric Moulines and Francis~R Bach.
\newblock Non-asymptotic analysis of stochastic approximation algorithms for
machine learning.
\newblock In {\em Advances in Neural Information Processing Systems}, pages
451--459, 2011.

[48]
\newblock Efficiency of coordinate descent methods on huge-scale optimization
\newblock {\em SIAM Journal on Optimization}, 22(2):341--362, 2012.

[49]
Yurii Nesterov and Sebastian~U Stich.
\newblock Efficiency of the accelerated coordinate descent method on structured
optimization problems.
\newblock {\em SIAM Journal on Optimization}, 27(1):110--123, 2017.

[50]
Lam~M Nguyen, Phuong~Ha Nguyen, Marten van Dijk, Peter Richt{\'a}rik, Katya
Scheinberg, and Martin Tak{\'a}{\v{c}}.
\newblock Sgd and hogwild! convergence without the bounded gradients
\newblock In {\em Proceedings of the 35th International Conference on Machine
Learning}, 2018.

[51]
Phuong~Ha Nguyen, Lam~M Nguyen, and Marten van Dijk.
\newblock Tight dimension independent lower bound on optimal expected
convergence rate for diminishing step sizes in sgd.
\newblock In {\em 33rd Conference on Neural Information Processing Systems
(NeurIPS 2019)}, 2019.

[52]
Mateusz Ostaszewski, Edward Grant, and Marcello Benedetti.
\newblock Structure optimization for parameterized quantum circuits.
\newblock {\em Quantum}, 5:391, 2021.

[53]
Vivak Patel, Shushu Zhang, and Bowen Tian.
\newblock Global convergence and stability of stochastic gradient descent.
\newblock {\em Advances in Neural Information Processing Systems},
35:36014--36025, 2022.

[54]
Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou,
Peter~J. Love, Al{\'a}n Aspuru-Guzik, and Jeremy~L. O'Brien.
\newblock A variational eigenvalue solver on a photonic quantum processor.
\newblock {\em Nature Communications}, 5:4213, 2014.

[55]
Frank Pollmann, Ari~M. Turner, Erez Berg, and Masaki Oshikawa.
\newblock Entanglement spectrum of a topological phase in one dimension.
\newblock {\em Phys. Rev. B}, 81:064439, Feb 2010.

[56]
Boris Polyak.
\newblock Gradient methods for the minimisation of functionals.
\newblock {\em Ussr Computational Mathematics and Mathematical Physics},
3:864--878, 12 1963.

[57]
{Qiskit contributors}.
\newblock Qiskit: An open-source framework for quantum computing, 2023.

[58]
Alexander Rakhlin, Ohad Shamir, and Karthik Sridharan.
\newblock Making gradient descent optimal for strongly convex stochastic
\newblock In {\em Proceedings of the 29th International Conference on Machine
Learning}, 2012.

[59]
Patrick Rebentrost, Maria Schuld, Leonard Wossnig, Francesco Petruccione, and
Seth Lloyd.
\newblock Quantum gradient descent and newton’s method for constrained
polynomial optimization.
\newblock {\em New Journal of Physics}, 21(7):073023, 2019.

[60]
Peter Richt\'{a}rik and Martin Tak\'{a}\u{a}\'{z}.
\newblock Iteration complexity of randomized block-coordinate descent methods
for minimizing a composite function.
\newblock {\em Math. Program.}, 144(1–2):1–38, 2014.

[61]
Bodo Rosenhahn and Tobias~J Osborne.
\newblock Monte carlo graph search for quantum circuit optimization.
\newblock {\em arXiv preprint arXiv:2307.07353}, 2023.

[62]
Gernot Schaller and Ralf Sch{\"u}tzhold.
\newblock The role of symmetries in adiabatic quantum algorithms.
\newblock {\em arXiv preprint arXiv:0708.1882}, 2007.

[63]
Maria Schuld, Ilya Sinayskiy, and Francesco Petruccione.
\newblock An introduction to quantum machine learning.
\newblock {\em Contemporary Physics}, 56(2):172--185, 2015.

[64]
Sukin Sim, Peter~D Johnson, and Al{\'a}n Aspuru-Guzik.
\newblock Expressibility and entangling capability of parameterized quantum
circuits for hybrid quantum-classical algorithms.
\newblock {\em Advanced Quantum Technologies}, 2(12):1900070, 2019.

[65]
James Stokes, Josh Izaac, Nathan Killoran, and Giuseppe Carleo.
\newblock Quantum natural gradient.
\newblock {\em Quantum}, 4:269, 2020.

[66]
Kevin~J Sung, Jiahao Yao, Matthew~P Harrigan, Nicholas~C Rubin, Zhang Jiang,
Lin Lin, Ryan Babbush, and Jarrod~R McClean.
\newblock Using models to improve optimizers for variational quantum
\newblock {\em Quantum Science and Technology}, 5(4):044008, sep 2020.

[67]
Ryan Sweke, Frederik Wilde, Johannes Meyer, Maria Schuld, Paul~K. Faehrmann,
Barth{\'{e}}l{\'{e}}my Meynard-Piganeau, and Jens Eisert.
\newblock Stochastic gradient descent for hybrid quantum-classical
\newblock {\em {Quantum}}, 4:314, 2020.

[68]
Ryan Sweke, Frederik Wilde, Johannes Meyer, Maria Schuld, Paul~K F{\"a}hrmann,
Barth{\'e}l{\'e}my Meynard-Piganeau, and Jens Eisert.
\newblock Stochastic gradient descent for hybrid quantum-classical
\newblock {\em Quantum}, 4:314, 2020.

[69]
Shiro Tamiya and Hayata Yamasaki.
\newblock Stochastic gradient line bayesian optimization for efficient
noise-robust optimization of parameterized quantum circuits.
\newblock {\em npj Quantum Information}, 8(1):90, 2022.

[70]
Simone Tibaldi, Davide Vodola, Edoardo Tignone, and Elisa Ercolessi.
\newblock Bayesian optimization for qaoa.
\newblock {\em arXiv preprint arXiv:2209.03824}, 2022.

[71]
David Wierichs, Josh Izaac, Cody Wang, and Cedric Yen-Yu Lin.
\newblock General parameter-shift rules for quantum gradients.
\newblock {\em Quantum}, 6:677, 2022.

[72]
Roeland Wiersema, Cunlu Zhou, Yvette de~Sereville, Juan~Felipe Carrasquilla,
Yong~Baek Kim, and Henry Yuen.
\newblock Exploring entanglement and optimization within the hamiltonian
variational ansatz.
\newblock {\em PRX Quantum}, 1(2):020319, 2020.

[73]
S.~J. Wright.
\newblock Coordinate descent algorithms.
\newblock {\em Mathematical Programming, Series B}, 151(1):3--34, 2015.

[74]
Nanyang Xu, Jing Zhu, Dawei Lu, Xianyi Zhou, Xinhua Peng, and Jiangfeng Du.
\newblock Quantum factorization of 143 on a dipolar-coupling nuclear magnetic
resonance system.
\newblock {\em Physical review letters}, 108(13):130501, 2012.

[75]
Jiahao Yao, Marin Bukov, and Lin Lin.
\newblock Policy gradient based quantum approximate optimization algorithm.
\newblock In {\em Mathematical and Scientific Machine Learning}, pages
605--634. PMLR, 2020.

[76]
Jiahao Yao, Paul Kottering, Hans Gundlach, Lin Lin, and Marin Bukov.
\newblock Noise-robust end-to-end quantum control using deep autoregressive
policy networks.
\newblock In Joan Bruna, Jan Hesthaven, and Lenka Zdeborova, editors, {\em
Proceedings of the 2nd Mathematical and Scientific Machine Learning
Conference}, volume 145 of {\em Proceedings of Machine Learning Research},
pages 1044--1081. PMLR, 16--19 Aug 2022.

[77]
Jiahao Yao, Haoya Li, Marin Bukov, Lin Lin, and Lexing Ying.
\newblock Monte carlo tree search based hybrid optimization of variational
quantum circuits.
\newblock In Bin Dong, Qianxiao Li, Lei Wang, and Zhi-Qin~John Xu, editors,
{\em Proceedings of Mathematical and Scientific Machine Learning}, volume 190
of {\em Proceedings of Machine Learning Research}, pages 49--64. PMLR, 15--17
Aug 2022.

[78]
Jiahao Yao, Lin Lin, and Marin Bukov.
\newblock Reinforcement learning for many-body ground-state preparation
inspired by counterdiabatic driving.
\newblock {\em Phys. Rev. X}, 11:031070, Sep 2021.

[79]
Xiao Yuan, Suguru Endo, Qi~Zhao, Ying Li, and Simon~C. Benjamin.
\newblock Theory of variational quantum simulation.
\newblock {\em {Quantum}}, 3:191, 2019.

[80]
Tianchen Zhao, Giuseppe Carleo, James Stokes, and Shravan Veerapaneni.
\newblock Natural evolution strategies and variational monte carlo.
\newblock {\em Machine Learning: Science and Technology}, 2(2):02LT01, dec

[81]
Daiwei Zhu, Norbert~M Linke, Marcello Benedetti, Kevin~A Landsman, Nhung~H
Nguyen, C~Huerta Alderete, Alejandro Perdomo-Ortiz, Nathan Korda, A~Garfoot,
Charles Brecque, et~al.
\newblock Training of quantum circuits on a hybrid quantum computer.
\newblock {\em Science advances}, 5(10):eaaw9918, 2019.

\end{thebibliography}

\end{document}