Variational Iterative Rotation Algorithm: Combinatorial Optimization with Classical Kicked Tops

Abstract

We investigate a classical formulation of the Quantum Approximate Optimization Algorithm (QAOA), realized as a Hamiltonian dynamical system of classical kicked tops, which we call the Variational Iterative Rotation Algorithm (VIRAL). The variational parameters are the transverse and longitudinal rotation angles at each of the p layers of the circuit. We find that VIRAL outperforms QAOA on the canonical Sherrington-Kirkpatrick spin-glass benchmark at all circuit depths, with the energy density converging to the ground state value linearly in 1/p. For large circuit depths, the optimized dynamics follows a Floquet protocol in which a pitchfork bifurcation destabilizes the equatorial fixed point and drives the spins toward polar Ising configurations. Our results demonstrate that the effectiveness of QAOA-like protocols derives primarily from their underlying iterative rotation structure, and that a classical implementation of it outperforms its quantum counterpart. We further elucidate its efficiency by reducing the many-body classical evolution to an effective Landau-Lifshitz dynamics for a single spin in a stochastic magnetic field. In this picture, the covariance matrix of the effective field reveals a nearly rank-one structure in which a single mode dominates the stochastic dynamics. In contrast, quantum fluctuations make the noise covariance of the effective quantum model of higher rank, hampering the control of the system. We propose nanometer-scale magnetic tunnel junctions as a natural physical platform for implementing VIRAL, where spin rotations can be realized using magnetic fields and spin torques.

Figures (5)

• Figure 1: The Variational Iterative Rotation Algorithm. Schematic illustration of the classical optimization protocol based on a dynamical system of classical kicked tops. The dynamics is initialized in a uniform configuration in which all unit spins are aligned at an angle $\alpha$ in the $xz$-plane. Each layer $n$ applies to each spin $\mathbf{s}_i$ a site-dependent rotation about the $z$-axis of angle $\gamma_n h_i^z$ generated by the local field $h_i^z = \sum_{j \neq i} J_{ij}\, s_{j}^z$, followed by a global rotation about the $x$-axis of angle $\beta_n$ acting uniformly on all the spins. This sequence is repeated $p$ times. The variational parameters $\alpha$, ${\gamma_n}$, and ${\beta_n}$ are optimized to minimize the average energy function of the final state $\langle H[\mathbf{s}_p]\rangle_J$. The resulting deterministic evolution steers the spins toward Ising configurations with $s_i^z=\pm1$, corresponding to the minimum of the cost function of the original optimization problem.
• Figure 2: Optimization of the SK energy and universality of the optimized control parameters.. (A) Energy density reached after $p$ layers of either classical (blue circles) or quantum (red diamonds) evolution. The classical results were obtained by optimizing the variational parameters on systems of size $N=2^9$, while the final energy density was evaluated on a larger instance with $N=2^{14}$ spins. The classical protocol can be augmented with a final thresholding layer that projects the spins onto their $z$ axis, $\mathrm{sign}(s_i^z)$, yielding improved performance (green triangles). The blue and green points on the $y$ axis correspond to an idealized protocol with $p=3200$ layers, obtained by fitting the optimized $\beta$ angles with a smooth function, resampling 3200 values, and fixing $\gamma=-1.2$; its performance is consistent with the extrapolated $1/p$ scaling. Quantum results are taken from Refs. farhi2022basso22boulebnane2025evidencequantumapproximateoptimization. Dashed lines denote linear fits to the large-$p$ data. Both classical protocols converge to $E/N \approx -0.761$, outperforming the quantum QAOA extrapolation $E/N \approx -0.757$. The solid black horizontal line represents the exact Parisi ground state energy density of $-0.76321...$ . (B) Optimized variational angles $\gamma_n$ (the local field scaling factor) and $\beta_n/\pi$ (the global rotation angle) as a function of the layer index for circuit depths of $p=20$ and $p=40$. Optimization results for system sizes ranging from $N=64$ to $N=512$ exhibit negligible dependence on $N$.
• Figure 3: Phase diagram of the kicked tops and evolution of spin polarizations.(A), Phase diagram in the $(\gamma, \beta)$ plane showing the stability regions of the equatorial fixed point. In the blue region all normal-mode frequencies $\omega_i$ are real. The red boundaries $\gamma = \pm\tan(\beta/2)$ mark the onset of a pitchfork bifurcation into polar fixed points with $s_i^z\neq0$ (red region), where the frequency of the softest mode vanishes ($\omega\to 0$). The orange boundaries $\gamma = \pm\cot(\beta/2)$ mark the onset of period-2 orbits (orange region), where the frequency of the hardest mode reaches $\omega\to\pi$. In the white regions both instabilities can occur. The four boundaries intersect at two multicritical points $(\beta, \gamma) = (\pi/2, \pm1)$ (red dots). The black curve shows the optimal schedule $\{(\gamma_n, \beta_n)\}_{n=1}^{p}$ for $p=40$, which begins near the critical line $\gamma = -\tan(\beta/2)$ and progressively steers the spins toward the Ising fixed point as $\beta\to0$. (B), Evolution of the $z$-components $s^z_i(n)$ of some representative spins across layers for $p = 40$ and $N=2048$. Spins are colored blue (decided) after their last crossing of $s_i^z = 0$ and red (undecided) while they are still actively crossing it. The intermediate layers show intense mixing, where the large transverse kicks drive the undecided spins repeatedly across zero, while the decided spins progressively consolidate near $\pm1$ as $\beta$ is reduced.
• Figure 4: Correlation, response, and small rank structure of the classical dynamics.(A), Correlation $C(n)=\langle\langle s^z(p)\ s^z(n-1)\rangle\rangle$ (blue dots) and response $R(n)=\langle\langle \partial s^z(p)/\partial\theta_n \rangle\rangle$ (orange dots) as a function of the layer index $n$ for a circuit of depth $p=40$. The response is large in the early layers, while the correlation is large in the later layers. (B), Single layer contribution $\varepsilon(n)=\gamma_nC(n)R(n)$ to the total energy Eq. \ref{['eq:energydecomposition']}, showing how much each layer contributes to the ground state energy. (C), Cumulative energy $e(n)=\sum_{k=1}^n \varepsilon(k)$ computed layer by layer via Eq. \ref{['eq:energydecomposition']}. (D), Eigenvalues of the covariance $\Sigma(t,t')$ of the effective field for the classical (blue circles) and quantum (red diamonds) protocols. The classical spectrum is strongly dominated by a single eigenvalue, reflecting the near rank-one structure of the noise covariance, while the quantum case exhibits a higher rank structure induced by quantum fluctuations.
• Figure S1: Distribution of oscillation frequencies in the interacting system. Probability density $P(\omega)$ of the frequencies of normal modes around the fixed point $\mathbf{s}_i^*=(1,0,0)$ for different values of $\gamma$ (with $\beta=\pi/2$). For small $\gamma$ (e.g., $\gamma=0.5$, blue curve), the frequencies are concentrated in a narrow band around $\omega=\pi/2$. As $\gamma$ increases towards the critical value $\gamma_c=1$, the width of the distribution expands and, at the critical point $\gamma=1$ (red curve), the support of $P(\omega)$ spans the entire range $\omega \in [0, \pi]$. The non-zero density at $\omega \to 0$ indicates the emergence of soft modes and at $\omega \to \pi$ the emergence of period-2 orbits.