Quantum Natural Gradient with Geodesic Corrections for Small Shallow Quantum Circuits
Mourad Halla





CQTA, Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen, Germany


Quantum Natural Gradient with Geodesic Corrections for Small Shallow Quantum Circuits
    Mourad Halla
    September 9, 2024
=====================================================================================



§ ABSTRACT 
The Quantum Natural Gradient (QNG) method enhances optimization in variational quantum algorithms (VQAs) by incorporating geometric insights from the quantum state space through the Fubini-Study metric. In this work, we extend QNG by introducing higher-order integrators and geodesic corrections using the Riemannian Euler update rule and geodesic equations, deriving an updated rule for the Quantum Natural Gradient with Geodesic Correction (QNGGC). QNGGC is specifically designed for small, shallow quantum circuits. We also develop an efficient method for computing the Christoffel symbols necessary for these corrections, leveraging the parameter-shift rule to enable direct measurement from quantum circuits. Through theoretical analysis and practical examples, we demonstrate that QNGGC significantly improves convergence rates over standard QNG, highlighting the benefits of integrating geodesic corrections into quantum optimization processes. Our approach paves the way for more efficient quantum algorithms, leveraging the advantages of geometric methods.




§ INTRODUCTION


Quantum computing represents a significant leap in computational science, enabling the resolution of problems that are fundamentally intractable for classical algorithms by leveraging quantum mechanical principles. Among the leading strategies in near-term quantum computing are Variational Quantum Algorithms (VQAs) <cit.>, which are tailored to harness the power of current noisy intermediate-scale quantum (NISQ) devices. A prominent example of VQAs is the Variational Quantum Eigensolver (VQE) <cit.>, a hybrid quantum-classical algorithm designed to find the ground state energy of quantum systems a task of paramount importance in areas such as quantum chemistry, materials science and condensed matter physics. VQE operates by iteratively optimizing a parameterized quantum state, known as an ansatz, using a combination of quantum and classical computations. Quantum processors are employed to prepare and measure the quantum state, while classical optimization algorithms adjust the parameters of the ansatz to minimize the cost function. This synergy between quantum evaluations and classical optimization enables VQE to efficiently explore the solution space, making it a powerful tool for solving complex problems such as those found in High-Energy Physics applications <cit.>.

Optimization techniques are crucial to the performance of VQAs, as they directly influence convergence rates and the quality of the solutions obtained. While traditional methods like vanilla gradient descent (GD) are commonly used due to their straightforward implementation, the unique challenges of quantum optimization landscapes characterized by non-convexity, noise, and barren plateaus demand more sophisticated techniques to enhance convergence and performance.

The Quantum Natural Gradient (QNG) algorithm <cit.>, a generalization of the natural gradient <cit.>, is an advanced optimization technique that enhances the optimization process by incorporating the geometry of the parameter space into the update rules. Unlike standard gradient descent, which assumes a flat parameter space, QNG leverages a Riemannian metric defined by the Fubini-Study metric or more generally the quantum Fisher information <cit.>, capturing the infinitesimal changes in distances between quantum states and aligning parameter updates with the natural curvature of the quantum state manifold. A key factor in QNG’s superior performance, as shown in <cit.>, is its ability to identify regions of high negative curvature early in the optimization process, significantly accelerating convergence. These regions play a crucial role in guiding the optimization along paths that lead to faster descent.

Following the foundational work by <cit.>, several extensions of the QNG have been developed to broaden its applicability and robustness. For example, the QNG was extended to handle noisy and nonunitary circuits, making it more suitable for realistic quantum devices that are subject to imperfections and decoherence <cit.>. Other works propose using simultaneous perturbation stochastic approximation techniques to approximate the QFIM in QNG <cit.>. Recently, <cit.> introduced two methods that significantly reduce the resources needed for state preparations required for QNG: the Random Natural Gradient and the Stochastic-Coordinate Quantum Natural Gradient.

The design of the ansatz, which serves as the parameterized quantum circuits in VQAs, plays a pivotal role in the algorithm’s success. The choice of ansatz directly influences the expressiveness and trainability of the quantum circuit, as well as its efficiency in hardware simulations, thus significantly impacting the overall performance of the VQA <cit.>. A well-designed ansatz must strike a balance between complexity and the capacity to accurately represent the target state, making the selection of the ansatz a crucial factor in the effective implementation of VQAs.

In QNG, Riemannian geometry <cit.> provides the framework to understand the geometric structure of the parameter space. The Fubini-Study metric, after Tikhonov-Regularization due to its general ill-definition in VQA, is a Riemannian metric that defines the intrinsic geometry of the quantum state space, guiding the optimization process by aligning parameter updates with the manifold’s curvature. Christoffel symbols describe how directions change as vectors are transported along this curved space, while geodesics represent the shortest and most efficient paths. Accurately computing these geodesics is essential for incorporating higher-order curvature corrections in QNG, enhancing optimization precision and overall performance.

The concept of using geodesic corrections to enhance optimization on manifolds was first introduced in classical optimization <cit.> and further developed in subsequent work <cit.>, primarily applied to nonlinear least squares problems under specific curvature assumptions. These methods highlighted the critical role of manifold geometry in boosting optimization performance. More recently, geodesic corrections were incorporated into natural gradient optimization within classical machine learning, demonstrating improved convergence by preserving higher-order invariance properties <cit.>. Building on these foundational ideas, our work integrates geodesic corrections into QNG, specifically tailored for variational quantum algorithms. This approach leverages the unique geometric characteristics of the quantum state space, enabling optimization that respects and exploits the manifold's inherent curvature for enhanced performance.

This manuscript presents a comprehensive approach to incorporating geodesic corrections into the QNG for applications in variational quantum algorithms. In Section <ref>, we provide an overview of differential geometry and geodesic equations, establishing the foundational mathematical context. Section <ref> delves into optimizing idealized variational quantum circuits using Quantum Natural Gradient Descent, emphasizing its effectiveness in quantum optimization. Section <ref> introduces higher-order integrators and derives the update rule for the Quantum Natural Gradient with Geodesic Correction (QNGGC). As the update rule relies on the Christoffel symbols of the second kind, Section <ref> is dedicated to efficiently computing these symbols using the parameter-shift rule, enabling direct measurements from quantum circuits. In Section <ref>, we apply the derived update rule to various examples: Examples 1 and 2 involve analytical calculations for numerical simulations, while Example 3 utilizes quantum software, specifically Qiskit, for practical simulations. Finally, in Section <ref>, we provide an outlook on potential extensions and future research directions, highlighting the integration of geodesic corrections into quantum optimization frameworks.




§ DIFFERENTIAL GEOMETRY AND GEODESIC EQUATIONS



To establish our notation, we will briefly review essential concepts in differential geometry relevant to our work. For readers interested in further details, we recommend <cit.> and <cit.> for introductions tailored to physicists, and <cit.> for more comprehensive technical discussions.

A manifold is a topological space that locally resembles Euclidean space and supports a consistent coordinate system. More formally, an n-dimensional manifold ℳ is a set equipped with a collection of coordinate charts {(U_i, φ_i)}, where each U_i ⊂ℳ is an open subset, and φ_i: U_i →ℝ^n is a homeomorphism, meaning that U_i is locally similar to ℝ^n. For ℳ to be a differentiable manifold, the transition maps between overlapping charts, φ_j ∘φ_i^-1: φ_i(U_i ∩ U_j) →φ_j(U_i ∩ U_j), must be smooth (infinitely differentiable). This smooth structure allows us to perform calculus on the manifold.

A Riemannian manifold is a differentiable manifold ℳ equipped with a metric tensor g_ij, which is a symmetric, positive-definite tensor field assigning an inner product to each tangent space T_pℳ at a point p ∈ℳ. In local coordinates {x^i}, the metric tensor defines the line element:


    ds^2 = g_ij  dx^i   dx^j,


where g_ij = g_ji. The metric allows for measuring lengths and angles on the manifold. The length of a smooth curve γ: [a, b] →ℳ is given by:


    L(γ) = ∫_a^b √(g_ij dγ^i/dτdγ^j/dτ)  dτ.


A geodesic on a manifold is a curve whose velocity vector remains parallel to itself along the curve, representing the straightest possible path given the manifold’s geometry. This property is formally expressed by stating that the geodesic has zero covariant acceleration, which accounts for the curvature of the manifold rather than the usual notion of acceleration in Euclidean space. Mathematically, a curve γ(τ) is a geodesic if it satisfies the geodesic equation:


    d^2 x^i/dτ^2 + Γ^i_jkdx^j/dτdx^k/dτ = 0,


where τ is an affine parameter along the curve, and Γ^i_jk are the Christoffel symbols of the second kind, defined by:


    Γ^i_jk = 1/2 g^il( ∂_j g_lk + ∂_k g_lj - ∂_l g_jk).


The Christoffel symbols define the covariant derivative, which maps tensor fields to other tensor fields, adapting to the curvature of the manifold.

The exponential map is a tool in Riemannian geometry that relates the tangent space at a point to the manifold. Given a tangent vector v ∈ T_pℳ, the exponential map exp_p maps v to a point on the manifold reached by traveling along the geodesic starting at p with initial velocity v. Formally,


    Exp_p(v) = γ(1),

where γ is the geodesic satisfying γ(0) = p and γ̇(0) = v. By rescaling the parameter v by a factor of h, the following relation holds:

    Exp_p(h v) = γ(h).


To approximate Exp_p(hv), we expand γ(h) around h = 0 using a Taylor series:


    γ(h) ≈γ(0) + hγ'(0) + h^2/2γ”(0) + 𝒪(h^3).



This expansion provides insight into the behavior of the geodesic near the starting point p. The first-order term, γ(0) + hγ'(0), corresponds to the initial position and velocity, giving a linear approximation of the geodesic. The second-order term, h^2/2γ”(0), accounts for the curvature effects, refining the approximation by including the acceleration of the geodesic. Since γ”(0) satisfies the geodesic equation, it inherently reflects the manifold's curvature. Thus, for small h, the exponential map can be locally approximated as:


    Exp_p(hv) ≈γ(0) + h γ'(0) + h^2/2γ”(0),



capturing both the initial direction and the curvature of the geodesic path.

In the context of quantum information, particularly for the Quantum Natural Gradient, the manifold of interest is the parameter space of quantum states, where the metric tensor relevant for us is defined by the Fubini-Study metric. This metric provides a Riemannian structure that captures the infinitesimal distance between quantum states, facilitating more efficient optimization in variational quantum algorithms. In the next sections, we review how to apply the Fubini-Study metric to the Quantum Natural Gradient and how to distinguish geodesic corrections to QNG from the perspective of the exponential map.






§ OPTIMIZING IDEALIZED VARIATIONAL QUANTUM CIRCUITS WITH QUANTUM NATURAL GRADIENT


This section provides an overview of VQAs and their optimization, focusing on idealized variational quantum circuits. For a more comprehensive understanding, refer to <cit.> for reviews of VQAs and <cit.> for details on QNG.

Variational quantum circuits are constructed using a family of parameterized unitary transformations. For an n-qubit system, the state space is represented by a Hilbert space of dimension N = 2^n, which can be decomposed as a tensor product of two-dimensional spaces: ℂ^N = (ℂ^2)^⊗ n. The parameterized circuits are typically composed of sequences of unitary transformations:


    U(θ) = U_L(θ_L) U_L-1(θ_L-1) ⋯ U_1(θ_1),


where θ_l represents the parameters for the l-th layer. Each unitary transformation U_l(θ_l) can be decomposed as:


    U_l(θ_l) = ∏_l e^-i θ_l K_l W_l,


where K_l are Hermitian operators, and W_l are fixed entangling unitary operators acting across the qubits.

The goal of VQE is to minimize a cost function, typically defined as the expectation value of an observable Ô with respect to the quantum state |ψ(θ)⟩ = U(θ) |ψ_0⟩:


    ℒ(θ)    = ⟨ψ(θ)|Ô|ψ(θ)⟩,


where |ψ_0⟩ is the initial state. The objective is to find the optimal parameters θ^* that minimize ℒ(θ).

The optimization is typically performed using gradient descent, updating the parameters iteratively:


    θ_t+1   = θ_t - η ∂_j ℒ(θ_t),


where ∂_j:= ∂/∂θ_j. Equation (<ref>) can be interpreted as an approximation of the solution to an ordinary differential equation (ODE) using the Euler method:


    θ̇   = - λ ∂_j ℒ(θ),

where η = h λ is the learning rate, with λ being a time scale constant that affects the speed but not the trajectory of the system, and h is the step size. 

However, this ODE is not invariant under reparameterizations of the parameters θ. For instance, if we rescale the parameters θ→ 2θ, the gradient ∂_j ℒ would scale as 1/2∂_j ℒ, leading to inconsistencies in the optimization process.

The core of this issue lies in the differential geometric nature of the gradient. The parameter update θ̇ transforms as a vector in the tangent space T_θℳ, while the gradient ∂_j ℒ is a covector (or 1-form) in the cotangent space T^*_θℳ. Since the ODE in Eq. (<ref>) attempts to relate objects in different spaces with distinct transformation rules, it is not an invariant relation.


QNG alleviates this issue by approximately solving an invariant ODE. The key idea is to raise the index of the gradient using a metric tensor g_ij. By raising the index of ∂_j ℒ on the right-hand side of the gradient descent ODE, the new ODE becomes:


    θ̇   = -λ g^ij∂_j ℒ(θ),


which is now a vector in T_θℳ, thereby resolving the type mismatch in Eq. (<ref>).

This new ODE is invariant under reparameterizations, ensuring that the forward Euler approximation:


    θ_t+1   = θ_t - η g^ij∂_j ℒ(θ_k),


remains consistent across different parameter spaces. The metric tensor g_ij is recognized as the Fubini-Study metric, derived from the real part of the Quantum Geometric Tensor:


    g_i j   = Re( ⟨∂_i ψ | ∂_j ψ⟩) - ⟨∂_i ψ | ψ⟩⟨ψ | ∂_j ψ⟩.


It is important to note that the metric tensor in VQAs does not always define a legitimate Riemannian metric, as it is often degenerate, meaning it may not be invertible. To address this issue, regularization techniques such as Tikhonov regularization are applied by adding a small constant multiplied by the identity matrix (λ I) to the metric, ensuring it is well-defined and invertible. Additionally, the computation of this metric tensor becomes increasingly intensive as the number of parameters in the circuit grows. Therefore, approximations, such as block-diagonal or diagonal forms, are crucial for practical applications <cit.>.

Thus, QNG utilizes a regularized approximated, often diagonal or block-diagonal, inverse metric tensor to perform parameter updates that respect the intrinsic geometry of the quantum state space, leading to improved convergence speed and accuracy in optimization within VQAs. In the next section, we explore how this approach can be further refined through higher-order integrators and the inclusion of geodesic corrections.




§ HIGHER-ORDER INTEGRATORS AND GEODESIC CORRECTION



The forward Euler method, commonly used in the QNG, provides only a first-order approximation to the exact solution of the natural gradient ordinary differential equation (ODE). For higher accuracy, higher-order integrators should be employed.

To further refine this approach, we can use the Riemannian Euler method, which leverages the exponential map to update the parameters, ensuring that the updates align with the geometry of the manifold. Using (<ref>), the Riemannian Euler update rule is given by:


    γ̇^i_t(h)    =θ_t+1 =  Exp_θ_t(-h λ  g^ij(θ_t) ∂_j ℒ(θ_t)),


where


    γ^i_t(0)    = θ_t, 
    γ̇^i_t(0)    = -λ g^i j∂_j ℒ(θ_t).



Here, the exponential map Exp translates the current parameter θ_t along the geodesic defined by the natural gradient. This approach is effective because it preserves invariance properties under reparameterization, owing to the characteristics of the exponential map. However, directly computing the exponential map is challenging, as it requires solving the geodesic equation exactly. To approximate this computation, we explore using the first and second-order derivatives.

The first derivatives approximate the geodesic as:


    γ^i_t(h)    ≈θ_t + h γ̇^i_t(0) 
    ⇒θ_t+1 = θ_t - η g^ij∂_j  ℒ.


where the learning rate η = h λ. This approximation corresponds to the naive Quantum Natural Gradient update rule, utilizing only the first-order information. For a more precise approximation, we can incorporate second-order information from the geodesic equation:


    γ^i_t(h)    ≈θ_t + h γ̇^i_t(0) + 1/2 h^2 γ̈^i_t(0),


where, according to the geodesic equation (<ref>):


    γ̈^i_t(0)    = - 1/2Γ^i_lmγ̇^l_t(0) γ̇^m_t(0).


The resulting update rule with the geodesic correction is:


    θ_t+1   = θ_t + h γ̇^i_t(0) - 1/2 h^2 Γ^i_l  mγ̇^l_t(0) γ̇^m_t(0),


We now combine all the relevant equations and heuristically allow the correction term in equation (<ref>) to depend on a tunable parameter b rather than being fixed to η^2. In other words, η^2 is too small to effectively capture the geodesic correction effect. This approach is justified because we use an approximate Fubini-Study metric, such as a diagonal metric, and fully capturing the curvature may require going beyond a second-order approximation in (<ref>). Such higher-order corrections would be computationally intensive. To avoid this complexity, we adjust the correction term flexibly with the parameter b, allowing us to capture the geodesic correction effects without the need for exact higher-order terms, which are computationally demanding.

This leads to the following update rule, incorporating the geodesic correction into the QNG:


    θ_t+1   = θ_t - η g^ij(θ_t) ∂_j ℒ(θ_t) - b/2Γ^i_lm(θ_t) 
    ( g^lj(θ_t) ∂_j ℒ(θ_t) )
    ( g^mj(θ_t) ∂_j ℒ(θ_t) ).


The update rule in equation (<ref>) shows that the correction term depends on the Christoffel symbols of the second kind. In the next section, we will discuss how to compute these symbols efficiently.







§ COMPUTING THE CHRISTOFFEL SYMBOLS WITH THE PARAMETER-SHIFT RULE


The Christoffel symbols, Γ^i_jk, are essential components in the update rule (<ref>), as they capture the curvature of the parameter space defined by the Fubini-Study metric. Traditionally, these symbols are computed by differentiating the Fubini-Study metric (<ref>), which involves first calculating the metric tensor and then deriving the Christoffel symbols through classical differentiation. This approach relies heavily on classical post-processing and does not leverage the quantum device directly.

To enable direct estimation of the Christoffel symbols on a quantum device, we reformulate their computation using the parameter-shift rule specifically adapted to the quantum circuit of interest. This method allows the Christoffel symbols to be directly accessed through quantum measurements, bypassing the need for classical differentiation and making the process more efficient and directly linked to the quantum state preparation. For a detailed overview of the parameter-shift rule and its application in deriving the Fubini-Study metric, please refer to the appendix and the related work in <cit.>.

The Fubini-Study metric for a pure variational quantum state |ψ(θ) ⟩ can be represented as a second-order tensor:


    g_j_1 j_2 (θ) = -1/2∂^2/∂θ_j_1∂θ_j_2. | ⟨ψ (θ') | ψ (θ) ⟩|^2 |_θ' = θ


Using the parameter-shift rule, this metric can be reformulated as:


    g_j_1 j_2 (θ) = -1/8   [ 
    | ⟨ψ (θ) | ψ (θ + (𝐞_j_1 + 𝐞_j_2) π / 2) ⟩|^2 
       - | ⟨ψ (θ) | ψ (θ + (𝐞_j_1 - 𝐞_j_2) π / 2) ⟩|^2 
       - | ⟨ψ (θ) | ψ (θ + (-𝐞_j_1 + 𝐞_j_2) π / 2) ⟩|^2 
       + | ⟨ψ (θ) | ψ (θ - (𝐞_j_1 + 𝐞_j_2) π / 2) ⟩|^2 
    ].


where 𝐞_j_1 and 𝐞_j_2 are unit vectors in the parameter space, and the shifted parameter is given by s = π/2.

If the metric tensor simplifies to diagonal elements, we have:


    g_j j(θ) = 1/4[ 1 - | ⟨ψ(θ) | ψ(θ + π e_j) ⟩|^2 ].


In this scenario, the Christoffel symbols can be directly computed from the derivatives of the metric tensor using the higher-order parameter-shift rule, enabling their extraction directly from the quantum circuit. This approach bypasses traditional classical computations, allowing for efficient and accurate estimation directly on the quantum hardware. For the diagonal case, the results can be encapsulated in the following proposition:


The Christoffel symbols of the second kind for the diagonal metric (<ref>) are given by:

    Γ^i_jk   = 1/2g_ii(θ)[1/8δ_ij( 
    - | ⟨ψ(θ) | ψ(θ + π e_i +π2 e_k) ⟩|^2 
    + | ⟨ψ(θ) | ψ(θ + π e_i - π2 e_k) ⟩|^2
    )
       +1/8δ_ik( 
    - | ⟨ψ(θ) | ψ(θ + π e_i +π2 e_j) ⟩|^2 
    + | ⟨ψ(θ) | ψ(θ + π e_i - π2 e_j) ⟩|^2
    ) 
       -1/8δ_jk( 
    - | ⟨ψ(θ) | ψ(θ + π e_j +π2 e_i) ⟩|^2 
    + | ⟨ψ(θ) | ψ(θ + π e_j - π2 e_i) ⟩|^2
    ) 
    ]




See Appendix <ref>.





§ APPLICATION EXAMPLES


In this section, we present examples to illustrate the practical application of the previously developed theories.



 §.§ Example 1 : Single Qubit

In this example, we examine a single qubit scenario. For the VQE, we take the ansatz state defined as:

    |ϕ(θ)⟩ = cosθ_0|0⟩ + e^2iθ_1sinθ_0|1⟩ =|ϕ(θ)⟩ = [        cosθ_0; e^2iθ_1sinθ_0 ].


The Hamiltonian for which we want to find the ground state is H = σ_x. Thus, the cost function for the system is:

    f(θ) = ⟨ϕ(θ)|H|ϕ(θ)⟩ = sin(2θ_0)cos(2θ_1),

and its gradient vector is:

    ∂ f(θ)/∂θ = [  2cos(2θ_0)cos(2θ_1); -2sin(2θ_0)sin(2θ_1) ].


The Fubini-Study metric for this single qubit example, using equation (<ref>), is given by:

    F = [           1           0;           0 sin^2(2θ_0) ].


The Christoffel symbol of this metric, calculated using equation (<ref>), the no zeros terms are:

    Γ^0_  1 1   = - 2 sin(2 θ_0) cos(2 θ_0),  
    Γ^1_  0 1   =Γ^1_  1 0 =   2 cos(2 θ_0)sin(2 θ_0)







Figure <ref> demonstrates the impact of geodesic corrections on optimization performance. The inclusion of geodesic corrections in QNGGC significantly enhances convergence speed (in terms of steps) compared to GD and QNG, as depicted in subplots (a) and (b), where QNGGC rapidly approaches the target energy E^* = -1. Subplot (c) shows the norm of the geodesic correction term, b/2Γ^i_lm(θ_t) 
( g^lj(θ_t) ∂_j ℒ(θ_t) )
( g^mj(θ_t) ∂_j ℒ(θ_t) ), which plays a crucial role in guiding the optimization process, particularly in the initial iterations. The subsequent reduction of this norm towards zero indicates that the convergence of the gradient is being effectively managed, leading to more stable updates. Subplot (d) illustrates the optimization paths in the parameter space, where QNGGC provides a more direct and efficient trajectory compared to GD and QNG. This behavior underscores the advantage of leveraging curvature corrections, resulting in more precise and adaptive parameter updates.






 §.§ Example 2: Two-Qubit Simulation of the Hydrogen Molecule (H_2)


In this example, we focus on finding the ground state of the hydrogen molecule (H_2) using a two-qubit variational quantum eigensolver (VQE) approach. The Hamiltonian for the system is given by <cit.>:


    H = α (σ_z ⊗ I + I ⊗σ_z) + βσ_x ⊗σ_x,


where α = 0.4 and β = 0.2.  In this Hamiltonian, σ_z and σ_x are Pauli matrices that represent spin operators acting on individual qubits. This Hamiltonian has four eigenvalues:


    h_1 = √(4α^2 + β^2),    h_2 = β,    h_3 = -β,    h_4 = -√(4α^2 + β^2),



with the ground state corresponding to the minimum eigenvalue h_4. The corresponding eigenvector is the ground state, |ψ_min⟩, which is expressed as:


    |ψ_min⟩∝ -β|00⟩ + (2α + √(4α^2 + β^2)) |11⟩.



The ansatz for the VQE is designed to approximate this ground state. It is given as follows (refer to Figure <ref> for a visual representation):


    |ψ⟩ = (CRY(θ)_q_0, q_1) (CRX(θ)_q_0, q_1) ( R_y(2θ_0) ⊗ R_y(2θ_1) ) |0⟩⊗|0⟩


where R_y(θ) denotes the single-qubit rotation operator defined as R_y(θ) = e^-iθσ_y/2, and the entangling gates act between two qubits, q_0 and q_1, as CRX(θ)_q_0, q_1 = I ⊗ |0⟩⟨ 0| + X ⊗ |1⟩⟨ 1| and CRY(θ)_q_0, q_1 = I ⊗ |0⟩⟨ 0| + R_y(θ) ⊗ |1⟩⟨ 1|. Please refer to Appendix <ref> for details on the cost function and its gradient.



In this example, the Fubini-Study metric F is calculated as:


    F = [                            1                            0            cos(θ_1) sin(θ_1);                            0                            1           -cos(θ_0) sin(θ_0);            cos(θ_1) sin(θ_1)           -cos(θ_0) sin(θ_0) 1/2(1 - cos(2θ_0) cos(2θ_1)) ]




To avoid singularities in F when updating the parameters in the VQE, we apply Tikhonov regularization by adding a small constant λ as a multiple of the identity matrix I. The inverse Fubini-Study metric F^-1 and the Christoffel symbols for this example are provided in Appendix <ref>.






Figure <ref> demonstrates that QNGGC significantly outperforms GD and QNG, rapidly reducing the energy cost and achieving faster convergence to the target energy. The geodesic correction helps QNGGC navigate the parameter space more effectively, avoiding inefficient trajectories and reaching optimal fidelity faster than the other methods.



Figure <ref> uses the same conditions as Figure <ref>, but with a different value of b = 0.4 and initial parameters set to a more challenging configuration, [π/2, π/2, 0]. This demonstrates how the QNGGC optimizer efficiently identifies a path to escape from plateaus, leading to faster convergence to the ground state compared to GD and QNG.

Figure <ref> presents averaged results over 50 runs with random initializations, highlighting QNGGC’s robust performance across varied starting conditions.







 §.§ Example 3: Transverse Field Ising Model


To enhance our previous results, we extend our approach to simulations involving larger qubit systems, including 4 and 7 qubits, using the quantum software Qiskit <cit.>. We employ the parameter-shift rule to compute both the gradient and Christoffel symbols. The ansatz used in VQE is the EfficientSU2 ansatz <cit.>, with one repetition as the circuit depth and linear entanglement for simplicity, as illustrated in Figure <ref>. Our objective is to determine the ground state energy of the Transverse Field Ising Model <cit.> under open boundary conditions:


    H = -∑_i Z_i Z_i+1 - h ∑_i X_i,




where Z_i and X_i are Pauli operators acting on the i-th qubit, with the transverse field strength fixed at h = 10 for our simulations. This model captures the interplay between nearest-neighbor spin interactions and a transverse magnetic field, driving the system between distinct quantum phases. For smaller values of the transverse field h, the spins tend to align along the z-direction, corresponding to a ferromagnetic ordered phase. As h increases, the system becomes increasingly disordered and transitions into a paramagnetic phase. The critical point at h = 1 marks the quantum phase transition between these phases, where spontaneous symmetry breaking and critical quantum fluctuations occur.


To benchmark the effectiveness of our optimizers, we present in Fig. <ref>, similar to previous examples, log_10(|⟨ H ⟩ - E_0|) to evaluate how quickly and accurately each optimization method converges to the ground state energy. This provides insights into their efficiency and robustness. Parameter tuning was conducted using a grid search to identify the optimal settings for each optimizer.





§ CONCLUSIONS AND OUTLOOK


In this manuscript, we extended the QNG method by incorporating higher-order integration techniques inspired by the Riemannian Euler update rule and a geometric perspective, leading to the development of an update rule based on geodesic equations. We introduced a tunable parameter b to heuristically adjust the effect of the geodesic correction, allowing for flexible integration of higher-order information. Through various examples, we demonstrated that QNGGC offers a more effective optimization strategy compared to traditional QNG by more closely following the direct optimization path. To compute the Christoffel symbols directly from quantum circuits, we proposed a method using the parameter-shift rule, enabling direct evaluation within quantum circuits.

The QNGGC update rule reduces the number of steps needed to reach the target energy but remains time-intensive compared to native QNG due to the computational cost of Christoffel symbol calculations. Currently, QNGGC is mainly applicable to small, shallow quantum circuits.

Future work could explore extending our results to noisy and nonunitary circuits, following approaches such as those presented in <cit.>. Developing more efficient methods for computing Christoffel symbols, such as stochastic parameter-shift rules and simultaneous perturbation stochastic approximation techniques for the Quantum Fisher Information <cit.>, would also be valuable. Additionally, randomness-based methods for Christoffel symbol computation <cit.> represent a promising avenue for further research. These extensions would enhance the broader applicability and efficiency of geodesic corrections in quantum optimization frameworks. Further research will also focus on moving beyond second-order approximations. Our approach could also be applied to time-dependent optimizers, particularly those inspired by the metric tensor and differential geometry.






§ ACKNOWLEDGEMENTS


The author thanks Karl Jansen, Stefan Kühn, Yahui Chai, Yibin Guo, Tim Schwägerl, and Cenk Tüysüz from CQTA, DESY, as well as Tobias Hartung (Northeastern University, London) and Naoki Yamamoto (Keio University, Japan) for their helpful discussions and insightful feedback.












§ PARAMETER-SHIFT RULE AND METRIC TENSOR CALCULATION


This section provides a brief review of the parameter-shift rule and its application in computing both the gradient and the metric tensor for variational quantum algorithms. The parameter-shift rule allows for the direct estimation of gradients and higher-order derivatives on quantum hardware. For further details, we refer the reader to <cit.>.



 §.§ Parameter-Shift Rule for Gradient Estimation


The parameter-shift rule provides a method to calculate the derivative of an expectation value with respect to a parameter θ_i in a quantum circuit. Consider an observable Ô and a parameterized quantum state |ψ(θ)⟩; the expectation value of the observable is given by:


    ⟨Ô⟩(θ) = ⟨ψ(θ) | Ô | ψ(θ) ⟩.


The derivative of this expectation value with respect to the parameter θ_i can be computed using the parameter-shift rule as:


    ∂_i⟨Ô⟩(θ)    = 1/2[⟨Ô⟩(θ + s 𝐞_i) - ⟨Ô⟩(θ - s 𝐞_i)],


where s is the shift in the parameter, typically s = π/2, and 𝐞_i is the unit vector in the direction of θ_i.





 §.§ Fubini-Study Metric


To derive the Fubini-Study metric tensor g_ij(θ), we begin by examining the infinitesimal distance between nearby quantum states in the parameter space. The metric tensor captures the local geometric structure of this space, reflecting how the quantum state |ψ(θ)⟩ evolves under small displacements <cit.> <cit.> <cit.>.

We start by considering the normalization condition of the quantum state:


    ⟨ψ(θ) | ψ(θ) ⟩ = 1.


Differentiating this condition twice with respect to the parameters θ^i and θ^j, we obtain:


    ⟨ψ(θ) | ∂_i ∂_j ψ(θ) ⟩ + ⟨∂_i ψ(θ) | ∂_j ψ(θ) ⟩ + ⟨∂_j ψ(θ) | ∂_i ψ(θ) ⟩ = 0.


Next, consider the quantum state |ψ(θ + δθ)⟩ near |ψ(θ)⟩, where δθ is the displacement vector. Using a Taylor expansion, we express the shifted state as:


    |ψ(θ + δθ)⟩ = |ψ(θ)⟩ + ∂_i |ψ(θ)⟩δθ^i + 1/2∂_i ∂_j |ψ(θ)⟩δθ^i δθ^j + 𝒪(δθ^3).


Taking the inner product of |ψ(θ)⟩ with this expanded state, we get:


    ⟨ψ(θ) | ψ(θ + δθ) ⟩ = 1 + ⟨ψ(θ) | ∂_i ψ(θ) ⟩δθ^i + 1/2⟨ψ(θ) | ∂_i ∂_j ψ(θ) ⟩δθ^i δθ^j + 𝒪(δθ^3).


The fidelity between the states |ψ(θ)⟩ and |ψ(θ + δθ)⟩ is given by:


    |⟨ψ(θ) | ψ(θ + δθ) ⟩|^2    = 1 + ⟨ψ(θ) | ∂_i ψ(θ) ⟩δθ^i + ⟨∂_i ψ(θ) | ψ(θ) ⟩δθ^i 
          + 1/2⟨ψ(θ) | ∂_i ∂_j ψ(θ) ⟩δθ^i δθ^j + 𝒪(δθ^3).


Using the normalization and differentiation results, we find the infinitesimal squared distance between the quantum states as:


    d^2(P_ψ, P_ψ + δθ) = Re[⟨∂_i ψ | ∂_j ψ⟩ - ⟨∂_i ψ | ψ⟩⟨ψ | ∂_j ψ⟩] δθ^i δθ^j.


This form directly defines the Fubini-Study metric tensor:


    g_ij(θ) = Re[⟨∂_i ψ(θ) | ∂_j ψ(θ) ⟩ - ⟨∂_i ψ(θ) | ψ(θ) ⟩⟨ψ(θ) | ∂_j ψ(θ) ⟩].


Thus, the Fubini-Study metric tensor g_ij(θ) quantifies the infinitesimal squared distance between nearby quantum states in the parameter space, providing a geometric interpretation of the state evolution. To express this metric tensor in terms of measurable quantities, we employ the parameter-shift rule, which enables direct computation on quantum hardware. The resulting expressions for the metric tensor components are given in (<ref>) for the full form and (<ref>) for the diagonal approximation.





 §.§ Proof of the Proposition


Using the higher-order derivatives strategy with the parameter-shift rule, the differentiation of the metric (<ref>) with respect to θ_j is given by:


    ∂_k g_j,j(θ) = 1/8[
    - | ⟨ψ(θ) | ψ(θ + π𝐞_j + π/2𝐞_k) ⟩|^2 
    + | ⟨ψ(θ) | ψ(θ + π𝐞_j - π/2𝐞_k) ⟩|^2
    ].


Using the definition of the Christoffel symbols (<ref>), for a diagonal metric tensor, the Christoffel symbols are determined by the following specific cases:

1. For i = j = k:
   
    Γ^i_ii = 1/2 g_ii∂ g_ii/∂θ^i,


2. For i = j ≠ k:
   
    Γ^i_ik = 1/2 g_ii∂ g_ii/∂θ^k,


3. For i = k ≠ j:
   
    Γ^i_ij = 1/2 g_ii∂ g_ii/∂θ^j,


4. For j = k ≠ i:
   
    Γ^i_jj = -1/2 g_ii∂ g_jj/∂θ^i,

   
5. For i ≠ j ≠ k:
   
    Γ^i_jk = 0.

   
The case i = j = k is ignored as shifting the same parameter twice disrupts distinct mixed derivative calculations.

Using the Kronecker delta notation, the expression can be summarized as:


    Γ^i_jk = 1/2 g_ii( δ_ij∂ g_ii/∂θ^k + δ_ik∂ g_ii/∂θ^j - δ_jk∂ g_jj/∂θ^i),


where δ_ij is the Kronecker delta, defined as:


    δ_ij = 
    1    if  i = j, 
    
    0    if  i ≠ j.


By substituting equation (<ref>) into equation (<ref>), we obtain the final expression for the Christoffel symbols for the diagonal metric (<ref>). 

□




§ DETAIL CALCULATION FOR EXAMPLE 2




In this appendix, we provide the detailed expressions for the inverse Fubini-Study metric and the Christoffel symbols for Example 2, which simulates the Hydrogen molecule (H_2).



 §.§ Inverse Fubini-Study Metric

The inverse Fubini-Study metric F^-1, after applying Tikhonov regularization, is given by:





 §.§ Christoffel Symbols

The Christoffel symbols Γ^i_jk for this example are given by:





 §.§ Ansatz, Cost Function, and Gradient


The ansatz given in equation (<ref>) can be expanded as follows:


    |ψ⟩   = (cos(θ_0) cos(θ_1))|00⟩
          + (cos(θ_0) cos(θ_2) sin(θ_1) - cos(θ_1) sin(θ_0) sin(θ_2))|01⟩
          + (sin(θ_0) sin(θ_1))|10⟩
          + (cos(θ_1) cos(θ_2) sin(θ_0) + cos(θ_0) sin(θ_1) sin(θ_2))|11⟩.
 

Using this ansatz and the hydrogen, H_2, Hamiltonian, we derive the following cost function:


    ℒ(θ_0, θ_1, θ_2, θ_3)    = 2 cos^2(θ_0) (-αsin^2(θ_1) sin^2(θ_2) + αcos^2(θ_1) + βsin(θ_1) cos(θ_1) sin(θ_2)) 
          - 2 sin^2(θ_0) cos(θ_1) (αcos(θ_1) cos^2(θ_2) + βsin(θ_1) sin(θ_2)) 
          + 2 sin(θ_0) cos(θ_0) cos(θ_2) (β - 2 αsin(θ_1) cos(θ_1) sin(θ_2))


The gradient of the cost function with respect to the parametersθ_0, θ_1, and θ_2 is expressed as:


    ∂ℒ/∂θ_0   = -sin(2θ_0) (α + 2αcos(2θ_1) + αcos(2θ_2) + 2βsin(2θ_1) sin(θ_2)) 
          + cos(2θ_0) (2βcos(θ_2) - αsin(2θ_1) sin(2θ_2)), 
    ∂ℒ/∂θ_1   = cos(2θ_0) (-2αsin(2θ_1) + 2βcos(2θ_1) sin(θ_2)) 
          - α(2 sin(2θ_1) sin^2(θ_2) + cos(2θ_1) sin(2θ_0) sin(2θ_2)), 
    ∂ℒ/∂θ_2   = βcos(2θ_0) cos(θ_2) sin(2θ_1) 
          - sin(2θ_0) (αcos(2θ_2) sin(2θ_1) + βsin(θ_2)) + α(-cos(2θ_0) + cos(2θ_1)) sin(2θ_2).






99
Cerezo2021
M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, et al., “Variational Quantum Algorithms,” Nature Reviews Physics, 3(9), 625-644 (2021). <https://doi.org/10.1038/s42254-021-00348-9>.

McClean2016
J. R. McClean, J. Romero, R. Babbush and A. Aspuru-Guzik, “The theory of variational hybrid quantum-classical algorithms,” New Journal of Physics, 18(2), 023023 (2016). <https://doi.org/10.1088/1367-2630/18/2/023023>.

Bharti2022
K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, et al., “Noisy intermediate-scale quantum algorithms,” Reviews of Modern Physics, 94(1), 015004(69) (2022). <https://doi.org/10.1103/RevModPhys.94.015004>.

Peruzzo2014
A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O'Brien, “A variational eigenvalue solver on a photonic quantum processor,” Nature Communications 5, 4213 (2014). <https://www.nature.com/articles/ncomms5213>doi:10.1038/ncomms5213.

karl
A. Di Meglio, et al., “Quantum Computing for High-Energy Physics: State of the Art and Challenges,” PRX Quantum 5, 037001 (2024). <https://doi.org/10.1103/PRXQuantum.5.037001>doi:10.1103/PRXQuantum.5.037001.

Stokes2020
J. Stokes, J. Izaac, N. Killoran, and G. Carleo, “Quantum Natural Gradient,” Quantum 4, 269 (2020). <https://doi.org/10.22331/q-2020-05-25-269>.

Amari1998
S.-I. Amari, “Natural Gradient Works Efficiently in Learning,” Neural Computation 10 (2), 251-276 (1998). <https://doi.org/10.1162/089976698300017746>

Meyer
J. Jakob Meyer,, “Fisher Information in Noisy Intermediate-Scale Quantum Applications
,” Quantum 5, 539 (2021). <https://doi.org/10.22331/q-2021-09-09-539>.


Katabarwa2022
A. Katabarwa, S. Sim, D. E. Koh, and P.-L. Dallaire-Demers, “Connecting geometry and performance of two-qubit parameterized quantum circuits,” Quantum 6, 782 (2022). <https://doi.org/10.22331/q-2022-08-23-782>.

Koczor
B. Koczor and  C. Benjami, “Quantum natural gradient generalized to noisy and nonunitary circuits,” Phys. Rev. A 106, 062416 (2022). <https://doi.org/10.1103/PhysRevA.106.062416>.

Gacon
J. Gacon, C. Zoufal, G. Carleo and S. Woerner, “Simultaneous Perturbation Stochastic Approximation of the Quantum Fisher Information
,” Quantum 5, 567 (2021). <https://doi.org/10.22331/q-2021-10-20-567>.

Kolotouros
I. Kolotouros and P. Wallden, “Random Natural Gradient
,” arXiv:2311.04135 (2023). <https://doi.org/10.48550/arXiv.2311.04135>.


Kandala2017
A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, “Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets,” Nature 549, 242-246 (2017). <https://doi.org/10.1038/nature23879>.


Sim2019Advanced
S. Sim, P. D. Johnson, and A. Aspuru-Guzik, “Expressibility and Entangling Capability of Parameterized Quantum Circuits for Hybrid Quantum-Classical Algorithms,” Advanced Quantum Technologies 2, 1900070 (2019). <https://doi.org/10.1002/qute.201900070>.

Tobias
L. Funcke, T. Hartung, K. Jansen, S. Kühn, and P. Stornati, “Dimensional Expressivity Analysis of Parametric Quantum Circuits,” Quantum 5, 422 (2021). <https://doi.org/10.22331/q-2021-03-29-422>.


Lee2013
J. M. Lee, Introduction to Smooth Manifolds, 2nd ed., Springer, 2012. <https://link.springer.com/book/10.1007/978-1-4419-9982-5>.

Wald1984
R. M. Wald, General Relativity, University of Chicago Press, 1984. <https://doi.org/10.7208/chicago/9780226870373.001.0001>

Frankel2011
T. Frankel, The Geometry of Physics: An Introduction, 3rd ed., Cambridge University Press, 2011. <https://doi.org/10.1017/CBO9781139061377>


Transtrum2011
M. K. Transtrum, B. B. Machta, and J. P. Sethna, “Geometry of nonlinear least squares with applications to sloppy models and optimization, ” Physical Review Letters E 83, 036701 (2011). <http://dx.doi.org/10.1103/PhysRevE.83.036701>

Transtrum2012
M. K. Transtrum and J. P. Sethna, “Geodesic acceleration and the small-curvature approximation for nonlinear least squares,” arXiv:1207.4999 (2012). <https://doi.org/10.48550/arXiv.1207.4999>

Song2018
Y. Song, J. Song, and S. Ermon, “Accelerating Natural Gradient with Higher-Order Invariance,”  Proceedings of the 35th International Conference on Machine Learning (ICML), Stockholm, Sweden, PMLR 80 (2018). <https://doi.org/10.48550/arXiv.1803.01273>


Yamamoto
N. Yamamoto, “On the natural gradient for variational quantum eigensolver,” arXiv:1909.05074 (2019). <https://doi.org/10.48550/arXiv.1909.05074>

BravoPrieto2020
C. Bravo-Prieto, J. Lumbreras-Zarapico, L. Tagliacozzo, and J. I. Latorre, “Scaling of variational quantum circuit depth for condensed matter systems,” Quantum 4, 272 (2020). <https://doi.org/10.22331/q-2020-05-28-272>

Mitarai
K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii, “Quantum circuit learning,” Phys. Rev. A 98, 032309 (2018). <https://doi.org/10.1103/PhysRevA.98.032309>

Schuld
M. Schuld, Ville Bergholm, C. Gogolin, J. Izaac, and N. Killoran, “Quantum circuit learning,” Phys. Rev. A 99, 032331 (2018). <https://doi.org/10.1103/PhysRevA.99.032331>

Mari
A. Mari, T. R. Bromley, and N. Killoran, “Estimating the gradient and higher-order derivatives on quantum hardware,” Phys. Rev. A 103, 012405 (2021). <https://doi.org/10.1103/PhysRevA.103.012405>

Wierichs
D. Wierichs, J. Izaac, C. Wang, and C. Yen-Yu Lin, “General parameter-shift rules for quantum gradients,” Quantum 6, 677 (2022). <https://doi.org/10.22331/q-2022-03-30-677>.


sun
R. Wiersema, D. Lewis, D. Wierichs, J. Carrasquilla, and N. Killoran, “Here comes the SU(N): multivariate quantum gates and gradients,” Quantum 8, 1275 (2024). <https://doi.org/10.22331/q-2024-03-07-1275>.

Riemannianflow
R. Wiersema and N. Killoran, “Optimizing quantum circuits with Riemannian gradient flow,” Phys. Rev. A 107, 062421 (2023). <https://doi.org/10.1103/PhysRevA.107.062421>.

Qiskit
H. Abraham, A. Akhalwaya, G. Aleksandrowicz, et al., “Qiskit: An Open-source Framework for Quantum Computing,” 2024. <https://qiskit.org>.

Mathematica
Wolfram Research, Inc., “Mathematica, Version 14.1,” Champaign, IL, 2024. <https://www.wolfram.com/mathematica/>.

Provost1980
J. P. Provost and G. Vallee, “Riemannian Structure on Manifolds of Quantum States,” Commun. Math. Phys. 76, 289-301 (1980). <https://doi.org/10.1007/BF02193559>.

DowlingNielsen2006
M. R. Dowling and M. A. Nielsen, “The geometry of quantum computation,” arXiv preprint arXiv:quant-ph/0701004, 2006. <https://doi.org/10.48550/arXiv.quant-ph/0701004>.

Haug2021
T. Haug, K. Bharti, and M. S. Kim, “Capacity and Quantum Geometry of Parametrized Quantum Circuits,” PRX Quantum 2, 040309 (2021). <https://doi.org/10.1103/PRXQuantum.2.040309>.

Heydari2015
H. Heydari, “Geometric formulation of quantum mechanics,” arXiv preprint arXiv:1503.00238, 2015. <https://doi.org/10.48550/arXiv.1503.00238>.






