Fermionic Quantum Approximate Optimization Algorithm

Quantum computers are expected to accelerate solving combinatorial optimization problems, including algorithms such as Grover adaptive search and quantum approximate optimization algorithm (QAOA). However, many combinatorial optimization problems involve constraints which, when imposed as soft constraints in the cost function, can negatively impact the performance of the optimization algorithm. In this paper, we propose fermionic quantum approximate optimization algorithm (FQAOA) for solving combinatorial optimization problems with constraints. Specifically FQAOA tackle the constrains issue by using fermion particle number preservation to intrinsically impose them throughout QAOA. We provide a systematic guideline for designing the driver Hamiltonian for a given problem Hamiltonian with constraints. The initial state can be chosen to be a superposition of states satisfying the constraint and the ground state of the driver Hamiltonian. This property is important since FQAOA reduced to quantum adiabatic computation in the large limit of circuit depth p and improved performance, even for shallow circuits with optimizing the parameters starting from the fixed-angle determined by Trotterized quantum adiabatic evolution. We perform an extensive numerical simulation and demonstrates that proposed FQAOA provides substantial performance advantage against existing approaches in portfolio optimization problems. Furthermore, the Hamiltonian design guideline is useful not only for QAOA, but also Grover adaptive search and quantum phase estimation to solve combinatorial optimization problems with constraints. Since software tools for fermionic systems have been developed in quantum computational chemistry both for noisy intermediate-scale quantum computers and fault-tolerant quantum computers, FQAOA allows us to apply these tools for constrained combinatorial optimization problems.

Quantum computers are expected to accelerate solving combinatorial optimization problems, including algorithms such as Grover adaptive search and quantum approximate optimization algorithm (QAOA).However, many combinatorial optimization problems involve constraints which, when imposed as soft constraints in the cost function, can negatively impact the performance of the optimization algorithm.In this paper, we propose fermionic quantum approximate optimization algorithm (FQAOA) for solving combinatorial optimization problems with constraints.Specifically FQAOA tackle the constrains issue by using fermion particle number preservation to intrinsically impose them throughout QAOA.We provide a systematic guideline for designing the driver Hamiltonian for a given problem Hamiltonian with constraints.The initial state can be chosen to be a superposition of states satisfying the constraint and the ground state of the driver Hamiltonian.This property is important since FQAOA reduced to quantum adiabatic computation in the large limit of circuit depth p and improved performance, even for shallow circuits with optimizing the parameters starting from the fixed-angle determined by Trotterized quantum adiabatic evolution.We perform an extensive numerical simulation and demonstrates that proposed FQAOA provides substantial performance advantage against existing approaches in portfolio optimization problems.Furthermore, the Hamiltonian design guideline is useful not only for QAOA, but also Grover adaptive search and quantum phase estimation to solve combinatorial optimization problems with constraints.Since software tools for fermionic systems have been developed in quantum computational chemistry both for noisy intermediate-scale quantum computers and fault-tolerant quantum computers, FQAOA allows us to apply these tools for constrained combinatorial optimization problems.

I. INTRODUCTION
Quantum optimization algorithms have attracted attention because of the potential for quantum computation to establish advantages in practical problems.The targets are combinatorial optimization problems in industry, such as financial optimization [1,2], logistics and distribution optimization [3,4], and energy optimization [5].In these practical problems, it is necessary to find a combination that gives the minimization cost while satisfying the constraints.
A standard quantum approach to solving these optimization problems are quantum annealing based on the adiabatic theorem [6][7][8].This approach slowly transforms a system Hamiltonian from a driver Hamiltonian Ĥd to a cost function-based problem Hamiltonian Ĥp , which leads to optimal solutions of the cost function from the ground state (g.s.) of Ĥd .However, the problems that can be solved with quantum annealing are limited to quadratic unconstrained binary optimization (QUBO) problems, whereas in practical industrial problems some constraints are imposed on general variables that are not restricted to binary values.Thus, in quantum annealing, a penalty function of quadratic form must be incorporated into the Ĥp as a soft constraint.Since this constraint is treated as an approximation, it may occasionally encounter states that are out of the solution space, making the computation inaccurate [2,9,10].Furthermore, for problems with large realistic sizes and complex interactions, the first excitation gap in the adiabatic time evolution becomes smaller and the adiabatic dynamics is time consuming [11].
Next, we will introduce the quantum approximate optimization algorithm (QAOA) [12], a variational algorithm for solving combinatorial optimization problems utilizing the controllability of a universal quantum computer.In principle, there are no hardware connectivity limitations, and quantum gates can handle higher-order interactions as well as second-order [13], which allows efficient encoding of general variables [14].This method covers the quantum adiabatic algorithm (QAA) by taking a large circuit depth p [7,12].Owing to this property, the variational parameters can be presumed to some extent, which may be an advantage even in shallow circuits [15][16][17].However, the same problems as in quantum annealing exist here for constrained optimization problems.
To solve these problems, Hadfield et al. proposed a hard constraint approach: a quantum alternating operator ansatz [18,19], which enforces the generated quantum states into the feasible subspace.In particular, XY -QAOA using a mixer with an XY Hamiltonian [20] has been applied to graph-coloring problems [10], portfolio optimization problems [2], extractive summarization problems [9], etc.In these problems, constraints are im-posed on the bit strings that keep the hamming weight constant, and the initial states are based on Dicke states [21] that satisfy the constraints.
According to these previous studies, the algorithms improve the computational accuracy by restricting the combinatorial search space, however, they have the following common problems.Due to generality of the mixer, there is no systematic description of an ansatz incorporating the hard constraints, and the relations between the initial states and the driver Hamiltonian are not explored enough.Therefore, the algorithm does not guarantee that it will result in QAA in the limit of large p.In addition, although the number of variational parameters remains the same as in conventional methods, no special effort has been made to handle the parameter optimization.Therefore, a systematic and efficient method of solving constrained optimization problems with guaranteed convergence is required.
In this paper, we developed a systematic framework called fermionic QAOA (FQAOA) for constrained optimization problems using fermionic formulation.The FQAOA translates the problem Hamiltonian and the driver Hamiltonian of a conventional spin system into the representation of fermion systems, and the equality constraint is naturally incorporated as a constant number of particles condition, and hence no penalty function is needed.We also show that the g.s. of the driver Hamiltonian can be prepared as an initial state, thus FQAOA will be QAA in the large p limit, allowing us to execute the time-discretized QAA.The parameters used here can then be set as initial parameters for variational quantum algorithms and efficiently computed using gradient descent methods, including parameter shifting methods [22,23].This FQAOA framework is reduced to QAA in the limit of large p.
As a representative example, we will take a portfolio optimization problem.This is a second-order, equalityconstrained, three-variable problem.Numerical simulations demonstrate that the proposed FQAOA provides a significant performance improvement.In particular, the computational accuracy at p = 1 of FQAOA with fixedangles determined by Trotterized QAA outperforms the results of the XY -QAOA using parameter optimization at p = 4 in the previous study [2] by about half of gate operations.In addition, the FQAOA with parameter optimization at p = 4 improves the probability of achieving low-energy states by a factor of 40 compared to the results of XY -QAOA at p = 4 in Ref. [2].Our algorithm enables us to simulate a wide variety of constraint optimization problems with high accuracy by using the platforms of quantum chemical computation both for noisy intermediate-scale quantum computers (NISQ) and faulttolerant quantum computers (FTQC).
This paper is organized as follows.In Sec.II, we introduce the constrained polynomial optimization problems.We formulate the FQAOA for these problems in Sec.III and introduce the quantization of the problems by the fermionic representation.In Sec.[1,24] and the required number of bits D to represent the largest integer I of z l , where ⌈x⌉ is the ceiling function of x.

Encoding
One-hot IV, we present the general framework for portfolio optimization problems and evaluate the numerical results of the FQAOA, where we compare our results with the computational results of other QAOAs, including the previous study.A summary is provided in Sec.V.

II. CONSTRAINED OPTIMIZATION PROBLEMS
The cost function E(z) for polynominal optimization problems with N integer variables z l ∈ {0, 1, 2, • • • , I} takes the following form: where α l1,l2,...,l k represent k-body interactions and V j is the j-th subset of vertex l. z l can be embedded in binary x l,d ∈ {0, 1} as [1,24] where encoding function f d is shown in Table I.In this paper, we do not employ one-hot encoding since it introduces an extra constraint originating from the encoding shown in Table I.For the remaining three encodings using f d , the cost function is transformed into the following Our goal is to find a bit string x * = arg min x E(x) under the constraints in Eq. (4).

III. FRAMEWORK OF FERMIONIC QAOA
We construct a framework of FQAOA shown in Fig. 1, which is a hybrid quantum-classical algorithm for solving constraint optimization problems expressed in fermion form within the unary encoding.The cost function and constraints in Eqs.(3) and ( 4) are mapped to the problem Hamiltonian Ĥp with a fixed number of localized fermions.A parametrized quantum circuit with (γ, β) is used to compute the expectation value E p (γ, β) of Ĥp , which is minimized by the outer parameter update loop.An FQAOA ansatz |φ p (γ, β) consists of an initial state preparation unitary Ûinit and a phase rotation unitary Ûp (γ) and a mixing unitary Ûm (β), where Ĥd is the driver Hamiltonian describing the nonlocal fermions.By carefully designing Ĥd and Ûinit , the constraints are satisfied at any steps without additional penalty function.As will be explained below, the framework is constructed so it comes down to QAA in the large p limit.The components of the FQAOA ansatz are described in detail in the following subsections.

A. Mapping to fermionic formulation
In this paper, the binary x l,d are mapped to the number operators of fermions on the (l, d) site as where ĉ † l,d (ĉ l,d ) is the creation (annihilation) operator on (l, d) site, which satisfies the anti-commutation relations: The computational basis corresponding to the bit string x ∈ {0, 1} N D can be written where the subscript (l, d) is collectively denoted as i and |vac is a vacuum satisfying ĉi |vac = 0. Note that the operators have the anti-commutation relations in Eq. ( 8), so the sequence of operators in Eq. ( 9) follows a fixed order.In this paper, we apply the snake type Jordan-Wigner (JW) ordering shown in Fig. 2, in which the relation between i and (l, d) can be explicitly denoted by

B. Problem Hamiltonian with linear constraints
The cost function with constraints in Eq. ( 3) with (4) are mapped to the following eigenvalue problems: where Ĥp and Ĉj is problem Hamiltonian and j-th constraint operator, respectively, which can be explicitly expressed as follows: where Ĥp is a model with k = 1, 2 . . ., K-body interactions of localized fermions.Thus, the constrained optimization problems have been replaced by the eigenvalue problems of obtaining the g.s.|φ x * and its energy E min = E(x * ) = min x E(x) under the constraints of Eq. (14).

C. Driver Hamiltonian and initial states
We address the design of the driver Hamiltonian Ĥd and initial states |φ 0 = Ûinit |vac .First, we assume that [ Ĥd , Ĥp ] = 0 to introduce hybridization between different basis states |φ x , Thus Ĥd represent non-local fermions represented by hopping terms, while Ĥp denote a localized ones.Conditions that Ĥd and |φ 0 must satisfy are as follows The following conditions are imposed as necessary conditions to satisfy the constraints at any approximation level p all times during the time evolution process.
• Condition II A series of hopping terms allows transitions between any two localized states |φ x ′ and |φ x that satisfy the constraints Framework of FQAOA, where (γ (0) , β (0) ) are the parameters in the time-discretized QAA and (γ * , β * ) are the parameters resulting from the parameter update, and Ûp(γ) = exp(−iγ Ĥp), Ûm(β) = exp(−iβ Ĥd ), and Ûinit is the unitary operator for preparing the ground state of the driver Hamiltonian Ĥd .
The orderings of operators in the Jordan-Wigner (JW) encoding [26] used in this paper, so-called snake type JW ordering, where sequential numbers i = 1, 2, ..., N D shown in blue are assigned in the order according to Eqs. (11) and (12).

• Condition III
The initial state |φ 0 = Ûinit is the g.s. of Ĥd and at the same time satisfies the constraints where E 0 is the g.s.energy of Ĥd .
By using the designed Ĥd and its g.s. as the initial state |φ 0 , the constrained optimization problems can be solved by performing adiabatic time evolution from the non-local initial state to the desired localized state.It is expected that FQAOA in this setting will allow more efficient solving of constrained optimization problems with shallow circuits.

D. Fermionic QAOA ansatz
The FQAOA ansatz satisfying the hard constraint is written in the following The Ûm (γ) is the mixer that changes the fermion configurations, the generator of which is the driver Hamiltonian Ĥd satisfying condition I and II in Sec.III C. The initial state |φ 0 = Ûinit |vac satisfies condition III in Sec.III C. For efficient quantum computation, it is important that Ûinit can be implemented in polynomial time.We will see this is actually the case for a linear order in terms of the number of qubits in the next section with a concrete example of the portfolio optimization problem.
The FQAOA contains the QAA [7] in the limit of large p. Taking the time-dependent Hamiltonian as Ĥ(t) = (1 − t/T ) Ĥd + (t/T ) Ĥp , discretizing time t by t j = (2j − 1)T /2p (j = 1, 2, . . ., p) with the execution time T = p∆t, the parameters assigned are as follows: the derivation of which by QAA is shown in Appendix A. As there, the key condition for performing the above QAA is that the initial state must be the g.s. of Ĥd under the constraint, which is satisfied by condition III.Previous studies on constrained optimization problems [9,10] often use the Dicke state as the initial state, and some of them do not satisfy this condition.
E. Parameter optimization for fermionic QAOA The objective of FQAOA is to obtain a bitstring that gives the lowest energy.To this end, the energy expectation value E p (γ, β) needs to be computed, which is obtained as the output of the FQAOA ansatz in Eq. ( 21) as The minimum energy at the approximation level p is determined by the following parameter optimization [12]: where the parameter set (γ * , β * ) determines the approximate variational wave function of the g.s. as |ψ p (γ * , β * ) .
The wave function |ψ p (γ, β) determines a probability P x (γ, β) that a bit string x is observed as In this paper, to avoid the difficulty of parameter optimization, the parameters in Eq. ( 22) are used as initial values for the parameter optimization.The above mentioned calculation procedure is summarized in Fig. 1.

IV. PORTFOLIO OPTIMIZATION PROBLEM
As an example of the application of FQAOA, we take a portfolio optimization problem and evaluate its performance by comparing with a previous study [2].The cost function for the optimization problem under the constraint of the total number of stock holdings being M is defined by [1,2,25] where the subscripts l and l ′ are the indices of the stocks and σ l,l ′ and µ l denote the asset covariance and average return, respectively.An integer w l is the discrete portfolio asset position to be held, representing long (w l ≥ 1), short (w l ≤ −1), or no-hold (w l = 0) position.The parameter λ for 0 ≤ λ ≤ 1 adjust the asset manager's risk tolerance.The first and second terms represent risk and return, respectively.Thus, if λ is large, risk is reduced.There are two cases in the portfolio optimization problems: one in which short position is not considered, where w l ∈ {0, 1, • • • , I}, and the other in which it is, where without short positions, In this paper, the latter is investigated with unary encoding f d = 1 and I = D, where D is defined to be even.The cost function and constraint in Eqs. ( 26) and ( 27) can be rewritten explicitly in the binary form as follows where To convert these problems to the fermionic formulation, we can simply let x l,d −→ nl,d for both cases in Eq. (29).

A. Problem Hamiltonian with particle number conservation
We present the problem Hamiltonian Ĥp and the eigenvalue problems for solving the portfolio optimization problems using quantum computation.The eigenvalue problems with Ĥp are as follows with where E(x) is energy eigenvalue defined in Eq. ( 30), M ′ = N D/2 − M , and |φ x is the same as Eq. ( 9).The correspondence between the stock positions, variables, and quantum states for a specific stock is shown in Table II, where the mapping to the spin-1/2 quantum states are also shown.Details of mapping to spin systems are described in Sec.IV D 1.
Hopping terms of Ĥt in Eq. (37), and that of Ûδ (β) for δ =I, II, III, IV, and BC in Eq. ( 56) and (b) FSWAP gate network of Fδ for δ =I and II in Eq. ( 62) on the D-leg ladder lattice.

B. Driver Hamiltonian and initial states on D-leg ladder lattice
Following the three guidelines in Sec.III C, we propose a driver Hamiltonian Ĥd and an initial state |φ 0 for unary encoding.First, from condition I, Ĥd is assumed to be a tight-binding model consisting of hopping terms.Next, from the condition II, the hopping terms are ĉ † l,d ĉl±1,d and ĉ † l,d ĉl,d±1 .The Ĥd designed in this way is the tight-binding model on a D-leg ladder lattice.Finally, from the codition III, the |φ 0 is set to the g.s. of Ĥd with particle number M ′ .
The tight binding Hamiltonian on the D-leg ladder lattice, which is adopted as Ĥd in this paper, is defined as where t and t ⊥ are the longitudinal and transverse hopping integrals in the tight-binding model, respectively.The correspondence between lattice labels and hopping is shown by the colored bonds in Fig. 3 (a), where the periodic boundary condition ĉN+1,d = ĉ1,d is imposed.
The energy eigenvalues ε k,m and the quasiparticles αk,m in Eq. ( 37) have the following form The initial state (the g.s. of Ĥd = Ĥt ) in FQAOA is as follows where the product for j is taken such that

C. Implementation on quantum circuits
Here we describe how to implement the ansatz shown in Eq. ( 21) on quantum circuits as Fig. 1.The ansatz can be divided into two parts: initial state preparation Ûinit and phase rotation (mixing) unitary Ûp (γ) [ Ûm (β)].The JW transformation [26] is needed to implement Fermionic operators.In this paper, we follow the transformation shown by Babbush et al [27].The number of gates required for the calculation is summarized in Table V of Appendix B.

Initial states preperation
The initial states |φ 0 of the FQAOA on the D-leg ladder lattice are explicitly expressed by the Slater determinant in Eq. (41).The states can be prepared by quantum circuits by using Givens rotations.Algorithms to prepare this initial state examined using quantum devices [28,29].In this paper, the algorithm of Jiang et al. [30] is used to prepare the initial state on quantum circuits, where the general Slater determinant for n qubits can be computed efficiently at a circuit depth of at most n for particle number n/2 [30].In the following, we will first present the general concepts, and then demonstrate their implementation on actual quantum circuit through a specific example.
Following Jiang et al. [30], the general construction of Ûinit is explained below.First, to reduce the number of operations in the quantum circuit, a unitary matrix V acts as follows: where The unitary matrix V is determined so that the components of the upper right triangular region of the matrix φ 0 are set to zero and the contribution of det(V ) can be ignored since it only changes the global phase.Next, with unitary matrix U satisfying the following equation: the following transformation is performed on Eq. ( 42) where the influence of Λ = M ′ j=1 λ j can be ignored.The U satisfying Eq. ( 43) can be constructed using a sequence of Givens rotations as [30], which can be written by which is a generalization of the usual Givens rotation to a complex matrix [30].The transformation of the basis ĉ † by G k can be written using the operator Ĝk as follows: Ĝk where (47) Therefore, the transformation of the basis ĉ † by U in Eq. ( 44) can be written as using Û = Ĝ1 Ĝ2 • • • ĜNG .Thus, the initial state can be written in the following form where the influence of Φ = NG k=1 ϕ k /2 can be ignored.Finally the unitary operator used for the initial state preparation can be written as The following is an example of a circuit for initial state preparation with M ′ = 4 particles on a N D = 4 × 2 ladder lattice.First, as shown in Eq. ( 42), to shorten the circuit depth, the components in the upper right triangular region of the matrix φ 0 in Eq. ( 40) are set to zero by the unphysical unitary transformation V as The matrix elements assigned numbers j (j = 1, 2, • • • 16) in the above matrix can be zeroed by acting G † k in Eq. (45) to the right in numeric order.Putting a series of Givens rotations as e iλ1 0 0 0 0 0 0 0 0 e iλ2 0 0 0 0 0 0 0 0 e iλ3 0 0 0 0 0 0 0 0 e iλ4 0 0 0 0 which corresponds to Eq. ( 43), and the initial state can be written by Eq. ( 44) , which can be rewritten in the form of Eq. ( 49).Finally, the quantum circuit for preparing the initial state |φ 0 = Ûinit |vac can be written as where The correspondence between i and (l, d) is denoted in Eqs. ( 11) and (12).Note that this operation is only valid between adjacent sites along the JW ordering in Fig. 2.

Phase rotation and mixing unitary
This section describes how to implement unitary operators of the phase rotation Ûp and the mixing Ûm on the D-leg ladder lattice in Fig. 3 in a quantum circuit.
First, we show the implementation of Ûp (γ).In portfolio optimization problems, the polynominal interactions are two-body terms.In the fermionic representation, the correspondence nl,d ↔ (1 − Ẑl,d )/2 allows us to implement the interaction term as Next, we describe the implementation of Ûm (β) = exp(−iβ Ĥt ).In this paper, Ûm (β) is implemented using the Trotter decomposition to keep the level of approximation the same as in the previous studies being compared.The effectiveness of this method has been demonstrated in quantum devices [29].Another technique to deal with this has been proposed using fermionic fast fourier tansformation [10,31,32], which transform Ĥt into a diagonal representation and performs the computation efficiently.
The followings are specific decomposition of Ûm (β) = exp(−iβ Ĥt ) used in this paper, Ûm (β) = ÛIV (β) ÛIII (β) ÛBC (β) ÛII (β) ÛI (β), (56) where Ûδ (δ = I, II, III, IV, and BC) consists of a set of commutable hopping pairs, which are shown as colorcoded bonds in Fig. 3 (a).ÛI,II causes hopping in the leg direction along the JW ordering shown in Fig. 2. which can be written as where i l,d is defined by Eq. ( 12) and where t is the hopping integral with t = t or t ⊥ .The hopping terms in Ûδ (β) (δ =III and IV) in Fig. 3 (a) that are out of the JW ordering in Fig. 2 can be efficiently computed using the network of fermionic swap (FSWAP) gates introduced in the literature [33].The implementation in this paper followed the method in Ref. [34].ÛIII (β) ÛIV (β) can be written with V (β) to the N -th power as ÛIV ÛIII = V (β) N with where Fδ (δ = I and II) are depicted in Fig. 3 (b) and f i swap is the FSWAP operator [27,35] which exchanges fermionic quantum states between i and i + 1 sites in the leg direction, taking into account the anti-commutation relation.
The hopping terms ÛBC (β) in Fig. 3 (a) can be embedded into ÛIV ÛIII where the quantum states at l = 1 and N are adjacent in the FSWAP network as where p = ⌊(N + 1)/2⌋ and q = p + 1.
As an example, the first part of the Ûm in Eq. ( 56) in the case of M ′ = 4 particles on a N D = 4 × 2 ladder lattice with t = t ⊥ ≡ t is shown below where and The correspondence between i and (l, d) is denoted in Eqs.(11) and (12).Note that Eqs. ( 67) and ( 68) are only valid between adjacent sites along the JW ordering in Fig. 3 (b).

D. Numerical evaluation of FQAOA
In this section, the validity of our proposed computational framework FQAOA is investigated through numerical simulations.First, we examine the proposed driver Hamiltonian in the conventional adiabatic time-evolution framework.Next, by comparing the statistical data obtained by FQAOA with other methods, we evaluate the performance of FQAOA assuming the use of NISQ devices.

Computational frameworks to be compared
Our proposed FQAOA is compared with two computational frameworks: the conventional X-QAOA [7] with soft constraint and the XY -QAOA with hard constraint [2].These are based on spin representations with spin-1/2 operator ŝα i for directions α = x, y, and z at the i site.In the X-QAOA, the initial state |φ 0 , which does not satisfy the constraints, is transferred to the approximate optimal solution indicated by the problem Hamiltonian Ĥp with a penalty term Ĥpen added.In the XY -QAOA, on the other hand, the Ĥpen is not required, because an eigenstate of the constraint operator Ĉ is chosen as an |φ 0 and a mixer Û XY m = exp(−i Ĥd ) satisfying the commutation relation [ Ĥd , Ĉ] = 0 in Eq. ( 17) is applied.
The problem Hamiltonian Ĥp and the constraint operator Ĉ are obtained by applying the mapping nl,d ↔ 1/2 − ŝz l,d to Eq. ( 34) and (35) as respectively.The driver Hamiltonian Ĥd and initial state |φ 0 explained below are summarized in Table III.
of ĤX is simply expressed as , which is used as the initial state |φ 0 .Since [ Ĥd , Ĉ] = 0 a penalty Hamiltonian Ĥpen has to be added to the Ĥp as where A is a parameter that adjusts the strength of the penalty.Another comparison (referred as XY -QAOA-I and II) is based on the quantum alternating operator ansatz [18], which Hodson et al. applied to the portfolio optimization problem [2].We here introduce the following XY Hamiltonian on a D circle lattice shown in Fig. 4 as where the condition [ ĤXY , Ĉ] = 0 is satisfied.In this case, the hard constraint can be imposed on the ansatz in Eq (21) by preparing an initial state that satisfies the constraints.In Ref. [2], the special case of D = 2 is applied.For implementation of the mixer Û XY m (β) = exp(−iβ ĤXY ), the following Trotter decomposition is applied where U XY δ (δ= I, II, and BC) consists of sets of commutable exchange pairs, which are shown as color-coded bonds in left panel of Fig. 4.
Next, we explain the initial state |φ 0 used in the methods XY -QAOA-I and II.Hodson et al. [2] adopted |φ 0 with long positions for M specific stocks and no-hold positions for the others as where the correspondence between spin states and positions is shown in Table II.We refer to their method using this initial state as XY -QAOA-I, in which the simulation results strongly depend on the initial stock holdings.Therefore, we introduce another method XY -QAOA-II, in which the initial stock position of XY -QAOA-I are superposed in all cases as where P is permutation of N symbols for the index l.
Note that since neither |φ I nor |φ II is a g.s. of the ĤXY , the QAA calculation cannot be applied.

Computational details
In the following, we compare our results with a previous study by Hodson et al [2].To compare with the previous study, the case where short positions can be taken (D = 2) for eight stocks (N = 8) and the total number of stocks held (M = 4) are used in the calculations.In Eq. ( 26), the parameters σ l,l ′ and µ l are the values in Fig. 2 and Table IV of Ref. [2], respectively.We also set λ = 0.9 and A = 0.003 according to the paper.The energy is scaled by W , which is the energy range of Ĥp in Eq. ( 34) under the constraints of Eq. (33).The hopping parameters in Ĥt are set to t = t ⊥ = W/W t , where the W t is the energy range of Ĥt at t = t ⊥ = 1.Numerical simulations have been performed using the fast quantum simulator qulacs [36].For the parameter optimization to obtain E p (γ * , β * ) in Eq. ( 24), the Broyden-Fletcher-Goldfarb-Shanno (BFGS) and conjugate gradient (CG) algorithm has been applied in our FQAOA.On the other hand, in X-QAOA and XY -QAOA-II, the BFGS with 100 basin hopping is applied to avoid trapping in the local minimum.The Nelder-Mead method has been applied in the previous studies (XY -QAOA-I) [2].The methods used for the parameter optimization are summarized in the fourth column of the Table III.

Simulation results using fixed-angle FQAOA
Hereafter we will show the simulation results of our proposed FQAOA.First, we compare the results of the fixed-angle FQAOA with that of the X-QAOA, where the variational parameters are fixed to (γ, β) = (γ (0) , β (0) ) in Eq. ( 22) for various W ∆t. We then show that FQAOA is reduced to QAA in the large p limit.
The simulation results of the residual energy ∆E = E p (γ (0) , β (0) ) − E min for various p by using fixed-angle QAOA are shown in Fig. 5.This corresponds to a timediscretized QAA with execution time T = p∆t and the discretized time width ∆t is adjusted for verification.These T and ∆t are scaled by W , which is the energy range of Ĥ′ p in Eq. (72) for X-QAOA and W = W for FQAOA.
According to the results in Fig. 5, the computational error of fixed-angle FQAOA is smaller than that of the conventional method X-QAOA for all p W ∆t region.This is because the FQAOA ansats strictly satisfies the constraint at all approximation levels p, whereas that of X-QAOA treats the constraints approximately and the penalty imposed on states out of the solution space increases the energy.The data in FQAOA for sufficiently Comparison of the residual energy ∆E = Ep(γ (0) , β (0) ) − Emin for fixed-angle FQAOA and X-QAOA, where (γ (0) , β (0) ) is obtained by the discrete-time QAA in Eq. ( 22) and Emin is the ground state energy of Ĥp under the constraint.As a reference, the data for p = 1, 2, . . .small W ∆t = 0.1 lie on the power law ∆E(T ) ∝ T −1/2 , which confirms the asymptotic approach to the g.s. by QAA [37][38][39].Also in X-QAOA, the asymptotic lines are parallel and their pre-factors are 1,000 orders of magnitude larger.Therefore, we can conclude that the FQAOA method, which intrinsically treats the constraints as particle number conservation laws, is superior to the conventional X-QAOA method.As a point of reference, open circles represent the results of fixed-angle FQAOA for p = 1, 2, . . ., 10 at W ∆t = 10.The results obtained here at p = 1 are superior to the previous study of XY -QAOA-II as will be seen next.

Simulation results using FQAOA with paremeter optimization
Here we present simulation results from the framework of FQAOA.First, to see the effect of parameter optimization, the results of the cumulative probability F (E) obtained by (a) fixed-angle FQAOA and (b) FQAOA with parameter optimization are shown in Fig. 6.The scaled variational parameters used in each calculation are presented in Fig. 7.
The F (E) is expressed by the following equations with probability distribution P (E) = ∂F (E)/∂E, which is derived from where P x (γ, β) is the probability of observing bit string x defined by Eq. (25).By a simple calculation, the error ∆E can be obtained by the following equation Since the area enclosed by the rectangle and the resulting curve in Fig. 6 corresponds to the ∆E, the smaller this area is, the smaller the expected value of the error.First we investigate the results of fixed-angle FQAOA shown in Fig. 6 (a).The F (E) at p = 0 is derived from the initial state |φ 0 , where the probability distribution is close to that of the randomly sampled under the constraint, so the |φ 0 is represented by a superposition of all classical states satisfying the constraint.The F (E) of fixed-angle FQAOA at p = 1 is larger than the value of the previous study XY -QAOA-I at p = 4 shown in Fig. 8 of the Ref. [2], so the computational error ∆E in the former is smaller than that in the latter.As a result, the fixed angle FQAOA at p = 1 outperforms the XY -QAOA-I at p = 4 [2] in the following three respects: the computational accuracy, the number of gate operations, and no need for variational parameter optimization.The actual number of gate operations for fixed angle FQAOA at p = 1 is about half the number of those for XY -QAOA-I at p = 4; specifically, the number of operations for single-and two-qubits is 556 and 656 for the former, respectively, while the latter is 936 and 1092.These values are derived from the Table VI in Appendix B.
We now turn to the results in (b).We see that the area between the rectangle and the curves decreases at each p compared to that in (a), indicating that the ∆E in Eq. (80) decreases for all p.In particular, P (E) = ∂F (E)/∂E increases in the low energy region.This means that the parameter optimization of FQAOA increases the probability of realization of low-energy (low-cost) states.Indeed, it can be seen that the probability of obtaining a solution in the low-energy region 0 ≤ E ≤ W/100, represented by F (E/100), increases from 0.357 to 0.521 by parameter optimization at p = 10.
We would like to make mention of the parameter optimization in Eq. ( 24).In the previous study of XY -QAOA-I [2], the Nelder-Mead method was applied, which is not realistic to implement in real devices because it requires multiple trials for initial parameter settings.In addition, there is no guarantee that the accuracy will improve with an increase in p for a small number of trials [2].On the other hand, in the FQAOA of this study, by setting (γ (0) , β (0) ) as initial parameters, the accuracy of the fixed-angle FQAOA is at least guaranteed.This appropriate initial parameter settings also allow the parameter to be updated efficiently in a single local search using the BFGS and CG methods.Figure 7 shows the results of the updated variational parameters.The optimized parameters scale well for increasing p, suggesting that (γ (0) , β (0) ) is a good starting point for any p.The value of E(γ * , β * ) is checked to be in agreement with the  The FQAOA also shows superiority in variational parameter optimization calculations.In the XY -QAOA methods, the fixed-angle QAOA reduced to QAA is not applicable.Therefore, stochastic methods have to be applied for parameter optimization due to the lack of a policy for setting initial parameters, which increases the computation time by the number of initial parameter sampling.On the other hand, as we have seen above, FQAOA provides good accuracy without parameter optimization, and moreover, the parameter-optimized results as variational algorithm can be obtained without initial parameter sampling.This property is a significant advantage in performing calculations on real devices.

V. SUMMARY AND DISCUSSION
In this paper, we proposed fermionic QAOA (FQAOA) for constrained optimization problems.This algorithm obtains the optimized solution by changing the physical system from the non-local fermions described by the driver Hamiltonian to the localized fermions corresponding to the optimization problem, where the constraint condition is regarded as the number of fermions.In this paper, we express the design guideline of the driver Hamiltonian in a fermionic formulation and propose a strategy to set the ground state of the driver Hamiltonian as the initial state.We also proposed to use the parameters of the fixed-angle QAOA corresponding to the time-discretized QAA as the initial variational parameters.As a demonstration of this algorithm, a portfolio optimization problem was taken up.The fixed-angle QAOA calculation confirms that the residual energy cal-culated by the FQAOA and the X-QAOA ansatz decays at similar powers, however, the FQAOA method, which essentially treats the constraint as a particle number conservation law, is superior to the conventional X-QAOA method.The resulting computational accuracy of the fixed-angle FQAOA at p = 1 exceeded that of the parameter-optimized XY -QAOA at p = 4 in the previous study [2] in about half the number of gate operations.Furthermore, the parameter-optimized FQAOA at p = 4 improved the probability of achieving low-energy states by a factor of 40 compared to that in the XY -QAOA at p = 4 [2].
In the present paper, non-interacting fermions were set as the initial state.On the other hand, the initial state preparation using Givens rotation allows the implementation of any Slater determinant.Thus, for example, the Hartree-Fock solution can be set as the initial state.In the most promising field of quantum chemical calculations, such calculations have been realized on quantum devices [28].Therefore, our algorithm is expected to significantly improve the computational accuracy of constrained combinatorial optimization problems by using the quantum chemical computing tools such as Open-Fermion [40], both for NISQs and FTQCs.
We note here that the unary encoding was applied in the portfolio optimization problem.For other encodings, FQAOA can be performed by preparing the initial state in a way that satisfies Eq. (31).However, since f d = 1, the constraint l x l,d = m d is imposed independently for each d, where d m d = M ′ .Correspondingly, since there are multiple initial states that depend on the assignment of m d , multiple FQAOAs must be attempted to search for the optimal solution.In addition, one-hot encoding and other efficient encodings may impose restrictions on embedding general variables into binary values.In this case, FQAOA can be performed in the same way by imposing constraints in encoding.
Finally, it is necessary to mention the quantum advantage.This is not limited to QAOA, but is an open problem in the field of quantum computing.In the context of QAOA, for example, warm-start techniques have been shown to improve the performance of the algorithm for some optimization problems [41,42], where the quantum algorithm can inherit the performance guarantees of classical algorithms.However, the effectiveness of this strategy is a future challenge as it depends on the available resources.The number of gate operations for FQAOA obtained by counting the number of gates in section IV C is shown in Table V.As a comparison, the counts for the four methods (X-QAOA, XY -QAOA-I and II, and FQAOA) at D = 2 are also shown in Table VI.For Ûp = exp(−iγ Ĥp ), all methods are implemented with the same quantum circuit, and the number of single-and two-qubit gate operations are as follows Also, since N l=1 ŝz l,1 is commutative with both Ĥp and ĤXY , the R z gate in the first line of Eq. (B1) can be ignored.
FIG. 5.Comparison of the residual energy ∆E = Ep(γ (0) , β (0) ) − Emin for fixed-angle FQAOA and X-QAOA, where (γ (0) , β (0) ) is obtained by the discrete-time QAA in Eq. (22) and Emin is the ground state energy of Ĥp under the constraint.As a reference, the data for p = 1, 2, . . ., 10 at W ∆t = 10 are shown by the open circles.Two bold gray lines represent the power law of ∆E ∝ T −1/2 for execution time W T = p W ∆t.
FIG. 5.Comparison of the residual energy ∆E = Ep(γ (0) , β (0) ) − Emin for fixed-angle FQAOA and X-QAOA, where (γ (0) , β (0) ) is obtained by the discrete-time QAA in Eq. (22) and Emin is the ground state energy of Ĥp under the constraint.As a reference, the data for p = 1, 2, . . ., 10 at W ∆t = 10 are shown by the open circles.Two bold gray lines represent the power law of ∆E ∝ T −1/2 for execution time W T = p W ∆t.

FIG. 8 .
FIG. 8. Statistical analysis of computational errors ∆E in Eq. (81) depending on the QAOA level p.In each box plot, the box shows quartiles, the divider is the median and the ends of the whisker indicate a range where the cumulative probability is between 0.05 to 0.95.Circles with broken lines represent the expectation values ∆E/W = [Ep(γ * , β * ) − Emin]/W , where Emin is the g.s.energy of Ĥp under the constraint.The inset in (c) shows expectation values ∆E/W in log scale, where ∆E using Ep(γ (0) , β (0) ) at W ∆t = 10 are shown by the open circles.

2 D 2 −
+ 10N D − 6N 2N 2 D + 2N D − 2NTABLE VI.Comparison of the gate counts of Ûinit and Ûm at D = 2.The numbers of single-and two-qubit gates for each method are shown in the upper and lower rows, respectively.M 2 ) 4N 2 + 2N Appendix B: Comparison of Gate Counts

N 2 ×
(2N + 1), 2N (2N − 1), respectively.The only difference between XY -QAOA I and II is the initial state.In XY -QAOA-I, N − M triplet states are generated for specific unowned shares as shown in Eq. (76).On the other hand, XY -QAOA-II requires the generation of a superposition of unowned shares as shown in Eq. (77).This initial state can be prepared by the following unitary operationÛinit = exp i π 2 ĈN,N−M ,(B1)where ĈN,N−M is the operator that generates a Dicke state of hamming weight N − M on N qubits with d = 2.The number of gate operations for single-and two-qubits in ĈN,N−M , respectively are as follows[43] 4N M − 4M 2 − 2N + 1, 5N M − 5M 2 − 2N.

TABLE I .
Encoding function f d

TABLE II .
Mapping of stock positions to variables and quantum states when short, no-hold, and long positions are considered (D=2).|φ l (|s l,1 , s l,2 ) is the state ket of the occupation (spin) representation for l-th stock, where |vac represents a vacuum.

TABLE III .
Comparison of calculation methods, where the driver Hamiltonian Ĥd , initial states |φ0 , and optimization methods of variational parameters are shown.In FQAOA, Broyden-Fletcher-Goldfarb-Shanno (BFGS) and conjugate gradient (CG) method give the same results. ).
[2]h those of previous study XY -QAOA-I[2], the calculation error is reduced to 1/6 and the probability of realization of low-energy states is increased by a factor of 40.The largest low-energy probability is F (W/100) = 0.521 for the parameter-optimized FQAOA at p = 10.

TABLE V .
Numbers of single-and two-qubit gates of Ûinit, Ûp, and Ûm in FQAOA ansatz.