Quantum Phase Estimation Algorithm for Finding Polynomial Roots

Theerapat Tansuwannont; Surachate Limkumnerd; Sujin Suwanna; Pruet Kalasuwan

doi:10.1515/phys-2019-0087

Article Open Access

Quantum Phase Estimation Algorithm for Finding Polynomial Roots

Theerapat Tansuwannont , Surachate Limkumnerd , Sujin Suwanna and Pruet Kalasuwan

Published/Copyright: December 31, 2019

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Open Physics Volume 17 Issue 1

Abstract

Quantum algorithm is an algorithm for solving mathematical problems using quantum systems encoded as information, which is found to outperform classical algorithms in some specific cases. The objective of this study is to develop a quantum algorithm for finding the roots of n^th degree polynomials where n is any positive integer. In classical algorithm, the resources required for solving this problem increase drastically when n increases and it would be impossible to practically solve the problem when n is large. It was found that any polynomial can be rearranged into a corresponding companion matrix, whose eigenvalues are roots of the polynomial. This leads to a possibility to perform a quantum algorithm where the number of computational resources increase as a polynomial of n. In this study, we construct a quantum circuit representing the companion matrix and use eigenvalue estimation technique to find roots of polynomial.

Keywords: Quantum algorithm; Quantum phase estimation; Eigenvalue problems; Polynomial roots

PACS: 02.90.+p; 03.67.a

1 Introduction

Roots finding is a centuries-old problem that has continued to attract considerable research interests and efforts due to its relevance in many fields of mathematics and physics involving geometry, number theory, probability and combinatorics. It is well known that for a polynomial of degree 4 or less, there exists a formula or procedure to solve for its roots exactly [1]. However, such a task is impossible for a polynomial of degree 5 or greater [2]. Many root-finding algorithms have been devised for obtain approximated roots of a polynomial of arbitrary degree [3, 4, 5].

The plausibility of a quantum computer—a new type of computation, which embraces quantum mechanics into its information, algorithms, and output measurements—has entailed quantum algorithms. A simple enhancement by rather non-intuitive mathematics of quantum mechanics like superposition, the uncertainty principle, and entanglement, brings quantum algorithms forth to a new level of computation unreachable before by conventional computers with classical algorithms. For example, the Shor’s algorithm, an algorithm to factorize a large integer into a product of primes, is proven in principles to overcome the classical-algorithm limit in terms of speed [6]. A breakthrough in quantum simulation is expected to bring eminent impact into science and technology [7, 8]. Recent progress on the actual quantum devices, such as a successful small-molecule simulation [9, 10, 11, 12] or probing the statistics of quantum systems [13, 14, 15], have yielded high promises and attracted immense interests. The advancement of computation and simulation in the aforementioned examples owes largely to a common underlying method called phase estimation algorithm (PEA) [16]. Still, PEA can be improved, especially at the fundamental algorithm of finding roots of a polynomial.

However, PEA is only applicable to unitary operators which are not always the case for some quantum algorithms; for instance, the phase measurement under the circumstance where decoherence is present in the process [17]. In a measurement process of the quantum algorithm, the non-unitary matrices also play key roles as projective operators. In order to modify existing PEA to be

suitable for eigenvalue problems comprehensively, a programmable circuit, and measurement of the control and ancillary qubits are recently exploited to tailor-made any arbitrary matrix [18]. The great advantage of this proposed scheme is that any matrix can be constructed and the control gate of the respective matrix can be realized, paving ways to build a quantum computer which can calculate eigenvalue of any matrix. However, some drawbacks exist as the algorithm itself may not be efficient for complicated matrices, which quantum complexity arises following the increasing number of non-zero matrix elements [19]. Further investigation on the algorithm in terms of appropriate complexity is still needed.

Our main aim in this paper is to propose a modified quantum phase estimation algorithm for finding polynomial roots, where we present a benchmark implementation of quantum non-unitary eigenvalue calculation scheme for polynomials. This specific task represents the least complex eigenvalue, which the algorithm can be fruitful without too much concern over the complexity.

This article is organized as follows. The remaining subsections of this section will cover key concepts and ideas about the phase estimate algorithm, and the iterated phase estimate algorithm (IPEA) for unitary operators, as well as the quantum algorithm to find complex eigenvalues of a general matrix. In Section II, we present our modified PEA and IPEA, together with the companion matrix approach, and more importantly, the circuit design to estimate roots of a polynomial of degree n. There we focus our presentation of the circuit operation and outputs, leaving the discussion and complexity analysis in Section III. Finally, the conclusions are summarized in Section IV.

1.1 PEA and IPEA for Eigenvalue Problems of Unitary Operators

In the original version of PEA, a phase φ arising after a unitary evolution U with eigenvalue exp(2πφi) is operated on its basis. Because a quantum evolution can be interpreted by a phase factor U = exp(−iHt/~), where H is a Hamiltonian of a finite system, the phase as a result of the phase estimation algorithm is indeed the eigenvalue of the Hamiltonian. PEA has also been introduced as a potential quantum tool to effectively solve various eigenvalue problems involving unitary operators [20, 21, 22, 23]. The unitary operators play a central role in all of the quantum algorithms, as they are required for universal quantum computations [24].

In order to estimate the value of a phase parameter ω_j up to the b bit-precision using PEA, b ancillary qubits in control register are required. In practice, however, the number of qubits which can be implemented is very limited. Iterative Phase Estimation Algorithm (IPEA) is an algorithm improved from the original PEA with an aim to estimate ω_j up to b^th digit while using only one ancillary qubit together with b iterations as a result of scalable inverse quantum Fourier transform in a semi-classical manner [25, 26]. In order to explain the algorithm as illustrated in Figure 1, we first assume that the phase parameter ω_j has a binary expansion no more than b digits (written as ω_j = 0.x₁x₂x₃ . . . xb000 . . .). Initially, all of the ancillary qubits are prepared in state |0〉 and the target register is prepared in the eigenstate | ψ j 〉 of unitary operator U. A Hadamard gate is applied to the control register in order to prepare state (|0〉 + |1〉) / 2 . In the first iteration (k = 1), a c − U 2 b − 1 and

Figure 1

A circuit for the k^th iteration of the IPEA where θ_k is the feedback from prior iterations.

(1) Z θ k = 1 0 0 e − i θ k ,

where θ₁ = 0 are applied. After that, the second Hadamard gate is applied on the control qubit and its state is measured in the computational basis {|0〉, |1〉}. This results in state

(2) 1 2 [ ( 1 + e 2 π i ( 0. x b ) ) | 0 〉 + ( 1 − e 2 π i ( 0. x b ) ) | 1 〉 ] ,

whose measurement gives either 0 or 1, and is determined by the majority probability between |0〉 and |1〉. This measurement result consequently dictates the value of x_b. The next iteration is performed with the c − U 2 b − k and Z (θ_k), where θ_k = 2π (0.0x_b_−k+2x_b_−k+3 . . . x_b) is calculated by the feed-forwarded measurement result of the prior iterations up to x_b_−k+1. The algorithm is finished when the digit x₁ is obtained.

The original IPEA has been used to determine only the phase parameter ω_j of the eigenvalue λ j = e 2 π i ω j of a unitary matrix U. In general, however, an eigenvalue of a non-unitary operator can be written as λ j = | λ j | e 2 π i ω j . The standard IPEA therefore cannot be applied without the knowledge of modulus |λ_j|.

1.2 Quantum Algorithm for Finding Complex Eigenvalues of General Matrices

Recently, Daskin et al. have introduced their technique to find the complex eigenvalues of general matrices [18]. In order to employ the IPEA on the non-unitary operators, first of all, the non-unitary operator O has to be controlled by a phase qubit and the c − U 2 b − k in Figure 1 is replaced by c- O 2 b − k .

In the scheme proposed by Daskin et al., the decomposition of the control gate of the non-unitary operator O of size N = 2^m uses the programmable circuit design, which requires m +1 ancillary qubits and m main qubits [18]. The operator O generally has a eigenvalue of the form | λ j | e 2 π i ω j . In case that ω_j = 0.x₁000 . . ., an operation of c-O followed by the Hadamard gate gives an output state

(3) 1 + | λ j | e 2 π i ( 0. x 1 ) | 0 〉 p + 1 − | λ j | e 2 π i ( 0. x 1 ) | 1 〉 p | ψ j 〉 m ,

where p and m denote phase qubit and main qubits, respectively. As can be seen, the c − O 2 b − k can be realized by the decomposition proposed by Daskin et al. It is also possible to estimate its phase ω_j via the IPEA from the probability shown in (3) in the same fashion as the phase estimation results determined by (2). However, the true novelty of the scheme is in the estimation of |λ_j|—taking the calculation to a complete eigenvalue estimation for any non-unitary matrix. Following the result shown in (3), the value of |λ_j| is related to the probability P₀ or P₁ of finding the phase qubit in states |0〉 or |1〉, respectively. Let P = max{P₀, P₁}, so that we can estimate |λ_j| as

(4) | λ j | = 2 N 2 P − 1 ,

where N is the dimension of matrix. In practice, |λ_j| is determined by the statistics of the measurement. We can also improve the accuracy of the estimation by using the statistics from other iterations. For the k^th iteration after which c − O 2 b − k is operated followed by Z(θ_k) and the Hadamard gate, the relationship between P^(k) and |λ_j| becomes

(5) | λ j | 2 b − k = 2 N 2 P ( k ) − 1.

Since we can estimate both |λ_j| and ω_j, the complex eigenvalue λ_j can be determined.

2 Quantum Algorithm for Finding Polynomial Roots

2.1 Companion Matrix Approach

As the aim of this study is to find roots of a generic polynomial of degree n, we can formulate this problem as the eigenvalue problem of a non-unitary operator. First of all, consider

(6) p ( x ) = x n + a n − 1 x n − 1 + ⋯ + a 1 x + a 0 ;

it can be factorized into the form

(7) p ( x ) = ( x − z 1 ) ( x − z 2 ) ⋯ ( x − z n ) ,

where z₁, z₂, . . . , z_n ∈ C are the roots of p(x). From general linear algebra [27], the roots of polynomial p(x) are eigenvalues of its companion matrix defined as

(8) C p = 0 1 0 0 1 0 ⋱ ⋱ 0 1 − a 0 − a 1 ⋯ − a n − 2 − a n − 1

with respect to the basis {1, x, x², . . . , xⁿ⁻¹}.

Daskin’s algorithm requires that the absolute value of every coefficient a_i must be less than or equal to 1, since rotation gates are used to simulate these coefficients. Therefore, we introduce a scaling method to meet this requirement. Let a_max denote the greatest absolute value of a₀, a₁, . . . , a_n. We choose a basis of circuit in the x-mode or (1/x)-mode depending on whether |a_n| or |a₀| is greater to maximize the success probability of the circuit scheme.

In case |a_n| > |a₀|, the x-mode will be chosen, so the polynomial p(x) can be equivalently expressed in the form

(9) a n a max x n + a n − 1 a max x n − 1 + ⋯ + a 1 a max x + a 0 a max = 0.

Let μ = a_max/a_n be a scaling factor. Then the corresponding eigenvalue equation is written as

(10) 0 a n a max 0 0 a n a max 0 ⋱ ⋱ 0 a n a max − a 0 a max − a 1 a max ⋯ − a n − 2 a max − a n − 1 a max 1 x ⋮ x n − 2 x n − 1 = 0 1 μ 0 0 1 μ 0 ⋱ ⋱ 0 1 μ − a 0 ′ − a 1 ′ ⋯ − a n − 2 ′ − a n − 1 ′ 1 x ⋮ x n − 2 x n − 1 = x μ 1 x ⋮ x n − 2 x n − 1 .

where a i ′ = a i / a max . The eigenvalue of modified companion matrix is x/μ.

On the other hand, if |a₀| > |a_n|, the (1/x)-mode will be used. Dividing the polynomial p(x) by a_maxxⁿ leads to

(11) a 0 a max 1 x n + a 1 a max 1 x n − 1 + ⋯ + a n − 1 a max 1 x + a n a max = 0.

In this case, a scaling factor is μ = a_max/a₀, and the corresponding eigenvalue equation is in the form

(12) 0 a 0 a max 0 0 a 0 a max 0 ⋱ ⋱ 0 a 0 a max − a n a max − a n − 1 a max ⋯ − a 2 a max − a 1 a max 1 1 / x ⋮ 1 / x n − 2 1 / x n − 1 = 0 1 μ 0 1 0 1 μ 0 ⋱ ⋱ 0 1 μ − a 0 ′ − a 1 ′ ⋯ − a n − 2 ′ − a n − 1 ′ 1 1 / x ⋮ 1 / x n − 2 1 / x n − 1 = 1 μ x 1 1 / x ⋮ 1 / x n − 2 1 / x n − 1 ,

where a i ′ = a n − i / a max in this case. Note that the eigenvalue of the modified companion matrix is 1 /μx.

However, the traditional companion matrix as described in (8) has 1’s in the upper diagonal entries but all of such entries of the modified companion matrices as shown in (10) and (12) have absolute values less than 1. To rectify this, we introduce a scaling gate S_m_,μ which will be explained in details later; see Equation (21).

2.2 Quantum Circuit Design

Our design of the respective algorithm relies on Polynomial Representative Circuit (PRC), a circuit to represent this modified companion matrix as illustrated in Figure 2. PRC requires m main qubits and 2 ancillary qubits where 2^m = n is a degree of the polynomial. (Although it is inconvenient, the circuit is also applicable for n ≠ 2^m simply by shifting the degree of polynomial up to the nearest power of 2.) First, let the main qubit be prepared in the initial state :

Figure 2

Polynomial Representative Circuit (PRC) which is used to represent an operation of a companion matrix.

(13) | α 〉 = α 0 α 1 ⋮ α n − 2 α n − 1 ,

and we define |β〉 as a result of C_p operating on |α〉; i.e.

(14) C p | α 〉 = | β 〉 .

Multiplying |α〉 by the modified companion matrix from (10) or (12) gives |β〉 in the form:

(15) | β 〉 = α 1 α 2 ⋮ α n − 1 − a ′ → ⋅ α → ,

where a ′ → ⋅ α → = a 0 ′ α 0 + a 1 ′ α 1 + ⋯ + a n − 1 ′ α n − 1 . Similar to the circuit introduced by Daskin et al., our circuit consists of Input Modification Block, Formation Block, and Combination Block. The main ingredient of the Input Modification Block is a cyclic-swap gate C_s applied on the main qubits. The matrix representation of the gate is

(16) C s = 0 1 0 ⋯ 0 0 0 1 ⋯ 0 ⋱ ⋱ 0 0 ⋯ 0 1 1 0 ⋯ 0 0 ,

which can be implemented by the Toffoli gates as shown in Figure 3.

Figure 3

The cyclic-swap gate can be implemented by a sequence of Toffoli gates.

The aim of the operator C_s is to generate the matrix element of the companion matrix from row 1 to row n − 1. Since the operation is underpinned by the presences of the sequences of Toffoli gates, the algorithm will be plagued by the huge complexity. The complexity of the algorithm is quite large, and yet still smaller than that of Daskin’s scheme, as their matrix elements are generated by the formation block which incurs more complexity.

The formation block in our version plays a role of the controlled gate of an operator F_μ, which represents the components in the last row of the modified companion matrix C_p as in (10) or (12). The rotation gate R_i is represented by a matrix as follows:

(17) R i = a i ′ 1 − a i ′ 2 − 1 − a i ′ 2 a i ′ ,

where i = 0, 1, 2, . . . , n − 1. Accordingly, the array of R_i forms the block matrix F _μ as follows:

(18) F μ = R 1 ⋱ R n − 1 R 0 .

The operation of F _μ can be simulated by a sequence of controlled-rotation gates as in Figure 4. The operation of F _μ will be performed on main qubits and the second ancillary qubit in case that the state of first ancillary qubit is |1〉.

Figure 4

The formation block can be simulated by a sequence of controlled-rotation gates.

Next step, in the combination block, we define the operator C as follows:

(19) C = X H ⊗ m ⊗ I = ∙ ∙ ∙ ∙ ⋯ ∙ ∙ ∙ ∙ ∙ ∙ ⋯ ∙ ∙ ⋮ ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ∙ ∙ ∙ ∙ ⋯ ∙ ∙ 1 0 1 0 ⋯ 1 0 ∙ ∙ ∙ ∙ ⋯ ∙ ∙ .

where "•" represents the elements we can neglected because they will be filter out by post-selection at the final stage of the algorithm.

There are three sub-tasks to undertake in the combination block. First, a controlled-C gate with the operator C is operated on the main qubits conditioning to the state of the first ancillary qubit as |1〉 to create the a ′ → ⋅ α → component. As a result, the operation of C_s, F _μ, and C on main

qubits and second ancillary qubit conditioning to the first ancillary state |1〉 give the following output:

(20) 1 2 m ∙ ∙ ⋮ ∙ a → ⋅ α → ∙

with the probability amplitude 1 / 2 m . This amplitude is required to be balanced with the case where the state of the first ancillary qubit is |0〉. At this stage of the algorithm, the scaling gate S_m_,μ is introduced to balance this probability and generate 1/μ in the upper diagonal entries of the modified companion matrix. It is defined as

(21) S m , μ = 1 2 m μ 1 2 m μ 2 − 1 − 2 m μ 2 − 1 1 .

This gate is, in fact, a rotation gate R_y(θ), where θ = 2 cos − 1 ⁡ ( 1 / 2 m μ ) . Scaling gate S_m_,μ is to be operated on the second ancillary qubit in case that first ancillary qubit is in state |0〉. After the controlled-S_m_,μ gate, the state is transformed into the following,

(22) | Ψ 1 〉 = 1 2 { | 0 〉 ⊗ [ I ⊗ m ⊗ S m , μ C s ⊗ I | α 〉 ⊗ | 0 〉 ] − | 1 〉 ⊗ C F μ C s ⊗ I | α 〉 ⊗ | 0 〉 } = 1 2 m + 1 μ [ | 0 〉 ⊗ [ α 1 ∙ α 2 ∙ ⋯ α n − 1 ∙ α 0 ∙ ] + | 1 〉 ⊗ [ ∙ ∙ ⋯ ∙ − a ′ → ⋅ α → ∙ ] ] .

However, referring to (22), the final state is not exactly |β〉. The last task is just to swap between the coefficient α₀ and − a ′ → ⋅ α → in (22) using the Toffoli gate in Figure 2, which results in

(23) | Ψ 2 〉 = 1 2 m + 1 μ | 0 〉 ⊗ | β 〉 ⊗ | 0 〉 + | ψ ⊥ 〉

where | ψ ⊥ 〉 refers to the case that the ancillary qubits do not all give the result ‘0’. Finally, the post-selection only the results of both ancillary qubits gives ‘0’, the output of the algorithm becomes |β〉 with the success probability of 1 / ( 2 m + 1 μ 2 ) .

2.3 Polynomial Root-finding by Eigenvalue Estimation Technique

In order to find roots of the polynomial, we will use the circuit shown in Figure 5. A controlled operation of the PRC by the phase qubit is denoted by c-C_p in the figure. To describe the operation, we will firstly assume that the main qubits are initially prepared in an eigenstate | ψ j 〉 of the companion matrix. An expected state from the operation will be in the form

Figure 5

A circuit scheme for finding polynomial roots.

(24) | β 〉 = | λ j | e 2 π i ω j | ψ j 〉 .

Here we will assume that ω_j has a binary expansion in the form ω_j = 0.x₁x₂x₃ . . . x_b, where b is bit-precision. As in IPEA, the c-C_p will be operated 2^b^−k times in the k^th iteration as illustrated in Figure 6. The result of the first iteration is given by

Figure 6

In order to estimate polynomial roots up to b bit-precision, c-C_p must be operated 2^b^−k times in the k^th iteration.

(25) 1 2 1 2 m μ 2 b − 1 | 0 〉 | 0 〉 | ψ j 〉 | 0 〉 + | λ j | 2 b − 1 e 2 π i ( 0. x b ) | 1 〉 | 0 〉 | ψ j 〉 | 0 〉 .

Similarly, for the k^th iteration, we have

(26) 1 2 1 2 m μ 2 b − k ( | 0 〉 | 0 〉 | ψ j 〉 | 0 〉 + | λ j | 2 b − k e 2 π i ( 0. x b − k + 1 x b − k + 2 … x b ) | 1 〉 | 0 〉 | ψ j 〉 | 0 〉 ) .

After the operation of Z(θ_k) with θ₁ = 0 and θ_k = 2π (0.0x_b_−k+2x_b_−k+3 . . . x_b) followed by the Hadamard gate, the phase qubit will be in the following state,

(27) 1 2 2 1 2 m μ 2 b − k { [ 1 + | λ j | 2 b − k e 2 π i ( 0. x b − k + 1 ) ] | 0 〉 + [ 1 − | λ j | 2 b − k e 2 π i ( 0. x b − k + 1 ) ] | 1 〉 } .

The value of x_b_−k+1 can be either 0 or 1. Therefore, the probabilities of finding the phase qubit in states |0〉 or |1〉 given that the both ancillary qubits give result 0 depend on the value of x_b_−k+1; namely,

(28) P 0 = 1 + 2 cos ⁡ ( 2 π 0. x b − k + 1 ) | λ j | 2 b − k + | λ j | 2 b − k + 1 8 2 m μ 2 2 b − k ,

(29) P 1 = 1 − 2 cos ⁡ ( 2 π 0. x b − k + 1 ) | λ j | 2 b − k + | λ j | 2 b − k + 1 8 2 m μ 2 2 b − k .

Since x_b_−k+1 can be either 0 or 1, the value of cos(2π0.x_b_−k+1) is either +1 or −1. In practice, x_b can be obtained by comparing P₀ and P₁ [25], i.e., x_b = 0 if and only if P₀ > P₁; and x_b = 1 if and only if P₀ < P₁. In addition, the value of |λ_j| from the k^th iteration can be calculated from the equation

(30) | λ j | 2 b − k + 1 = 4 2 m μ 2 2 b − k | P 0 + P 1 | − 1.

It should be emphasized that in later iterations, the parameter θ_k used in Z(θ_k) is constructed from x_b_−k+1 from the prior iterations as in IPEA. Finally, an estimate of λ_j can now be obtained and the corresponding root of the polynomial can be calculated depending on which mode (x-mode or 1/x-mode) is being used.

However, in general, the eigenstates of the companion matrix are unknown. The following approach is to estimate the greatest eigenvalue |λ_max|. Suppose that an initial state is prepared in a mixed state, the density operator can be expressed as

(31) ρ = ∑ j A j | ψ j 〉 ⟨ ψ j | ,

where A_j is a probability of preparing the initial state in the eigenstate | ψ j 〉 . The operation of c-C_p when the phase qubit is in state |1〉 transforms the density operator as

(32) ρ ↦ C p ρ C p † .

For the k^th iteration after which c − C p 2 b − k is operated and followed by Z(θ_k) and the Hadamard gate, the probabilities of finding phase qubit in states |0〉 or |1〉 given that both ancillary qubits give result 0 are

(33) P 0 = ∑ j A j 1 ± 2 | λ j | 2 b − k + | λ j | 2 b − k + 1 / κ ,

(34) P 1 = ∑ j A j 1 ∓ 2 | λ j | 2 b − k + | λ j | 2 b − k + 1 / κ .

where κ = 8 2 m μ 2 2 b − k . In Equations (33) and (34), the terms with the largest eigenvalue |λ_max| will dominate if the number of iterations is large enough. Hence, the probabilities P₀ and P₁ will be reduced to the following forms:

(35) P 0 ≈ 1 + A max | λ max | 2 b − k + 1 ± 2 | λ max | 2 b − k / κ ,

(36) P 1 ≈ 1 + A max | λ max | 2 b − k + 1 ∓ 2 | λ max | 2 b − k / κ .

The value of x _b_−k+1 can be found by comparing P₀ and P₁. Even without knowing the probability A_max, an estimate of |λ_max| can be obtained from

(37) | λ max | 2 b − k = 2 P 0 + P 1 − 2 / κ | P 0 − P 1 | .

After the largest root is found, it can be factorized from the polynomial, and the same technique and procedure can be repeated to calculate the other roots.

3 Discussion

In order to justify the efficiency of the quantum algorithm for finding roots of a polynomial, we have to compare it with the classical algorithm for solving the same problem. One of the most efficient classical algorithms for finding roots of a polynomial was created by Pan in 2002 [28]. We shall compare these two versions of the algorithm based on (i) resources required for calculation and (ii) algorithmic complexity.

We start by comparing the number of bits and qubits required by the computations. In order to find roots of an n^th degree polynomial, the quantum algorithm requires O(log n) qubits. In contrast, the Pan’s classical root-finding algorithm requires O(n) or O(n log n) bits. This makes it obvious that the quantum algorithm requires significantly fewer bits than its classical counterpart. In this way, especially for larger n, the quantum version of algorithm is capable for finding roots of a much higher order degree than the classical one.

Next, we will compare their algorithmic complexities. In Pan’s algorithm, the number of operations required to find the roots is

(38) O ( ( n log 2 ⁡ n ) ( log 2 ⁡ n + log ⁡ b ) )

where b is the bit precision of the solutions. In contrast, since any m-qubit unitary gates can be simulated using only single-qubit gates and CNOT gates [29], it is appropriate to compute the complexity of a quantum circuit in terms of the number of single-qubit gates and CNOT gates required to construct the circuit. From polynomial root-finding circuit, many unitary operations are controlled by several qubits. We will use the following corollary to calculate complexity of these gates (See Corollary 7.12 in Barenco et al. [24]).

Corollary: For any unitary U, the corresponding c-U gate controlled by (m − 2)-qubit can be simulated by O(m) basic operations in m-qubit network, where the initial value of one qubit is fixed and incurs no net change.

In order to apply this corollary with the root-finding circuit, we need one more ancillary qubit and set its initial value to state |0〉. Note that this ancillary qubit can also be reused to simulate several controlled-unitary operations. In order to compute the overall complexity of the whole circuit, a complexity of each part may be found separately. The cyclic-swap gate can be simulated by m CNOT gates. Number of control qubits of each gate varies from 0 to m−1. For other CNOT gates with multiple control qubits, its operation can be simulated by O(i+2) operations where i > 1 is the number of control qubits. Therefore the overall complexity of the cyclic-swap gate is 2 + ∑ i = 2 m − 1 i ∼ O ( m 2 ) . The formation gate consists of 2^m controlled rotation gates. Each gate, which is further controlled by m qubits, can be simulated by O(m + 2) basic operations. Hence the total complexity of the formation gate is O(2^mm). The combination gate consists of m Hadamard gates and m NOT gates which brings its overall complexity to O(m). Finally, the complexity of the (m+1)-qubit CNOT gate before the application of the last Hadamard gate can be readily computed by the corollary above, which amounts to O(m).

The total complexity of the complete circuit is the sum of complexity of each part as described above. However, in the complexity calculation, only the greatest term is kept. In the case of large n (where n is the degree of a polynomial), we can clearly see that the dominant term comes from the formation gate. Therefore, the complexity of the circuit in terms of the degree of a polynomial is

(39) O ( 2 m m ) ∼ O ( n log ⁡ n ) ,

for one iteration, or

(40) O ( k n log ⁡ n ) ,

for k iterations. Compared with the classical case for b-bit precision, our approach needs k = 2^b and the total complexity of the quantum version is

(41) O ( 2 b n log ⁡ n ) .

Comparing between (38) and (41) (see Figure 7), at certain values of b with increasing n, the quantum algorithm performs faster than the classical one at high n.However, there are some disadvantage influences of b, when the computation is undertaken with higher bit-precision requirements on low degree polynomials.

Figure 7

We compare the number of operations needed for our quantum algorithm (magenta for b = 3, red for b = 8) and the best known classical algorithm (blue for b = 3, green for b = 8), with the increasing n. Our algorithm clearly outperforms the classical one for the calculation with high n for both values of b (upper figure), especially at low bit-precision. For low n, the bit-precision b increases the number of operations exponentially, instead of linearly as in the Pan’s classical algorithm.

In addition, the rising complexity of the algorithm for high degree polynomials can possibly be inflicted by its inscalable success probability. Following (22), the success probability is diminished to 1/2^mμ and will require more iterations to maintain the precision. As a consequence, the increase iterations will effectively raise the complexity.

However, here we can validate that the quantum version of the algorithm is less complex than the classical version in the case where the polynomial has a high degree.On the contrary, the quantum algorithm may be an overkill when finding the roots, with high-bit precision, of a low degree polynomial.

4 Conclusion

Polynomials are the key foundation of various calculations in mathematics and physics. The computational advantages provided by this quantum-based algorithm could be very useful for a broad field of studies [30, 31, 32, 33, 34, 35, 36, 37]. However, to implement this algorithm in further application efficiently, some requirements of the particular task must be considered. In summary, we have provided a quantum algorithm for finding roots of the n^th degree polynomial partially based on Daskin et al.’s circuit for finding complex eigenvalues of a general matrix [18]. To make a comparison with classical version of the algorithm, resources and algorithmic complexities are considered. The quantum version requires significantly fewer number of (quantum) bits than its classical counterpart for a high degree polynomial. In terms of algorithm complexities, the quantum algorithm also trumps the classical algorithms for a high degree polynomial, requiring low bit-precision solutions. The growth in complexity stems from a larger number of iterations needed to achieve the desired precision. Although our result clearly shows that finding the roots of polynomials using quantum information scheme is possible, the most important challenge, however, remains in strengthening the algorithm to overcome the classical algorithm both in the utilized resources and the chosen precisions. Another challenge lies of course in the practical issue of a working quantum computer. It is well known that the current quantum computer technology still falls short of the theoretical requirement of the algorithm, especially in terms of the number of entangled qubits and multiple-qubit quantum operations.

Acknowledgement

We would like to thank Assist. Prof. Dr. Kwan Arayathanitkul for helpful discussions.T.T. acknowledges the support of The Queen Sirikit Scholarship under The Royal Patronage of Her Majesty Queen Sirikit of Thailand. This work is a collaboration of Collaborative Research Unit on Quantum Information, Mahidol University and Optical and Quantum Communication (OQC) Laboratory, National Electronics and Computer Technology Center (NECTEC). Grant No. 035/2557 from the Development and Promotion of Science and Technology Talents project (DPST) scholarship, research fund for DPST graduate with first placement is acknowledged.

References

[1] Weisstein E.W., CRC encyclopedia of mathematics, 2nd ed., CRC Press Taylor & Francis, Boca Raton, Florida, 2009Search in Google Scholar

[2] Jacobson N.,Basic algebra 1, 2nd ed.,Dover, 2009Search in Google Scholar

[3] McNamee J.M., Numerical Methods for Roots of Polynomials, Part 1, Elsevier, Amsterdam, The Netherlands, 2007Search in Google Scholar

[4] Mekwi W.R., Iterative methods for roots of polynomials, Master’s thesis, University of Oxford, Oxford, UK, 2001Search in Google Scholar

[5] Pan V.Y., Solving a Polynomial Equation: Some History and Recent Progress, SIAM Review,1997, 39, 187-22010.1137/S0036144595288554Search in Google Scholar

[6] Shor P.W., Algorithms for quantum computation: discrete logarithms and factoring, In: Goldwasser S.(Ed.), Proceedings 35th Annual Symposium on Foundations of Computer Science (20-22 Nov 1994, Santa Fe, NM, USA), IEEE Comput. Soc. Press, 1994, 124–13410.1109/SFCS.1994.365700Search in Google Scholar

[7] Lloyd S., Universal Quantum Simulators, Science, 1996, 273, 1073–107810.1126/science.273.5278.1073Search in Google Scholar PubMed

[8] Feynman R.P., Simulating Physics with Computers, Int. J. Thy. Phys., 1982, 21, 467–48810.1007/BF02650179Search in Google Scholar

[9] Aspuru-Guzik A., Dutoi A.D., Love P.J., Head-Gordon M., Simulated Quantum Computation of Molecular Energies, Science, 2005, 309, 1704–170710.1126/science.1113479Search in Google Scholar PubMed

[10] Ma X.S., Dakic B., Naylor W., Zeilinger A., Walther P., Quantum simulation of the wavefunction to probe frustrated Heisenberg spin systems, Nature Phys., 2011, 7, 399–40510.1038/nphys1919Search in Google Scholar

[11] Brown K.R., Clark R.J., Chuang I.L.,Limitations of Quantum Simulation Examined by Simulating a Pairing Hamiltonian Using Nuclear Magnetic Resonance, Phys. Rev. Lett., 2006, 97, 05050410.1103/PhysRevLett.97.050504Search in Google Scholar PubMed

[12] Negrevergne C., Somma R., Ortiz G., Knill E., Laflamme R., Liquid-state NMR simulations of quantum many-body problems, Phys. Rev. A, 2005, 71, 03234410.1103/PhysRevA.71.032344Search in Google Scholar

[13] Matthews J.C.F., Poulios K., Meinecke J.D.A., Politi A., Peruzzo A., Ismail N., Wörhoff K., Thompson M.G., O’Brien J.L.,Observing fermionic statistics with photons in arbitrary processes, Sci. Rep., 2013, 3, 153910.1038/srep01539Search in Google Scholar PubMed PubMed Central

[14] Ovidiu I.P., Correlation functions and momentum distribution of one-dimensional hard-core anyons in optical lattices, Journal of Statistical Mechanics: Theory and Experiment, 2015, 1, P0100410.1088/1742-5468/2015/01/P01004Search in Google Scholar

[15] Crespi A., Suppression laws for multiparticle interference in Sylvester interferometers, Phys. Rev. A, 2015, 91, 01381110.1103/PhysRevA.91.013811Search in Google Scholar

[16] A.Y. Kitaev, Quantum measurements and the Abelian Stabilizer Problem, Coll. Comput. Complex., 1996, 3, 3Search in Google Scholar

[17] Kacprowicz M., Demkowics-Dobrzanski R., Wasilewski w., Banszek K., Walmsley I.A., Experimental quantum-enhanced estimation of a lossy phase shift, Nature Photonics, 2010, 4, 357–36010.1038/nphoton.2010.39Search in Google Scholar

[18] Daskin A., Grama A., Kais S., A universal quantum circuit scheme for finding complex eigenvalues, Quantum Information Processing, 2014, 13, 333–353.10.1007/s11128-013-0654-1Search in Google Scholar

[19] Daskin A., Grama A., Kollias G., Kais S., Universal programmable quantum circuit schemes to emulate an operator, J. Chem. Phys., 2012, 137, 23411210.1063/1.4772185Search in Google Scholar PubMed

[20] Abrams D.S., Lloyd S., Quantum Algorithm Providing Exponential Speed Increase for Finding Eigenvalues and Eigenvectors, Phys. Rev. Lett., 1999, 83, 516210.1103/PhysRevLett.83.5162Search in Google Scholar

[21] Zhou X.Q., Kalasuwan P.,Ralph T.C., O’Brien J.L., Calculating unknown eigenvalues with a quantum algorithm, Nature Photonics, 2013, 7, 223–22810.1038/nphoton.2012.360Search in Google Scholar

[22] Wang H., Wu L.A., Liu Y.x., Nori F., Measurement-based quantum phase estimation algorithm for finding eigenvalues of non-unitary matrices, Phys. Rev. A, 2010, 82, 06230310.1103/PhysRevA.82.062303Search in Google Scholar

[23] Barz S., Kassal I., Ringbauer M., Lipp Y.O., Dakic B., Aspuru-Guzik A., Walther P., Univariate Polynomials: Nearly Optimal Algorithms for Numerical Factorization and Root-finding, Scientific Reports, 2014, 4, 611EP10.1038/srep06115Search in Google Scholar PubMed PubMed Central

[24] Barenco A., Bennett C.H., Cleve R., Vincenzo D.P.D., Margolus N., Shor P., Sleator T., Smolin J.A., Weinfurter H., Elementary gates for quantum computation, Phys. Rev. A, 1995, 52,345710.1103/PhysRevA.52.3457Search in Google Scholar PubMed

[25] Griflths R.B., Niu C.S., Semiclassical Fourier Transform for Quantum Computation, Phys. Rev. Lett., 1996, 76, 322810.1103/PhysRevLett.76.3228Search in Google Scholar PubMed

[26] Chiaverini J., Britton J., Leibfried D., Knill E., Barrett M.D., Blakestad R.B., Itano W.M., Jost J.D., Langer C., Ozeri R., Schaetz T., Wineland D.J., Implementation of the Semiclassical Quantum Fourier Transform in a Scalable System, Science, 2005, 308, 997–100010.1126/science.1110335Search in Google Scholar PubMed

[27] Buchberger B., An algorithmic method in polynomial ideal theory, D. Reidel Publishing Company, Dordrecht, Boston Langcaster, 198510.1007/978-94-009-5225-6_6Search in Google Scholar

[28] Pan V.Y., Univariate Polynomials: Nearly Optimal Algorithms for Numerical Factorization and Root-finding, Journal of Symbolic Computation, 2002, 33, 701–73310.1006/jsco.2002.0531Search in Google Scholar

[29] Barenco A., A universal two-bit gate for quantum computation, Proc. R. Soc. Lond. A, 1995, 449, 67910.1098/rspa.1995.0066Search in Google Scholar

[30] Hoffmann C.M., Geometric and Solid Modeling, San Mateo, California, 1989Search in Google Scholar

[31] Go N., Scherga H.A., Ring Closure and Local Conformational Deformations of Chain Molecules, Macromolecules, 1970, 3, 178– 18710.1021/ma60014a012Search in Google Scholar

[32] Faugeras O., Three-Dimensional Computer Vision : A Geometric Viewpoint, MIT Press, Cambridge, Massachusate, 1993Search in Google Scholar

[33] Bishop G., Fuchs H, Research directions in virtual environments, Computer Graphics, 1992, 26, 15310.1145/142413.142416Search in Google Scholar

[34] Baleanu D., Asad J.H., Jajarmi A., New aspects of the motion of a particle in a circular cavity, P ROMANIAN ACAD A, 2018, 19, 361–36710.1186/s13662-018-1684-xSearch in Google Scholar

[35] Baleanu D., Asad J.H., Jajarmi A., Classical and fractional aspects of two coupled pendulums, ROM REP PHYS, 2019, 71, 103Search in Google Scholar

[36] Massey J., Step-by-step Decoding of the Bose-Chadhuri-Hocquenghem Codes, IEEE Trans. Inf. Theory, 1965, 11, 580–58510.1109/TIT.1965.1053833Search in Google Scholar

[37] Hong J., Vetterli M., Simple ALgorithms for BCH Decoding, IEEE Trans. Commun., 1995, 43, 2324–233310.1109/26.403765Search in Google Scholar

Received: 2019-08-14

Accepted: 2019-10-25

Published Online: 2019-12-31

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/phys-2019-0087

Keywords for this article

Quantum algorithm; Quantum phase estimation; Eigenvalue problems; Polynomial roots

Creative Commons

BY 4.0