Prime factorization using quantum variational imaginary time evolution

The road to computing on quantum devices has been accelerated by the promises that come from using Shor’s algorithm to reduce the complexity of prime factorization. However, this promise hast not yet been realized due to noisy qubits and lack of robust error correction schemes. Here we explore a promising, alternative method for prime factorization that uses well-established techniques from variational imaginary time evolution. We create a Hamiltonian whose ground state encodes the solution to the problem and use variational techniques to evolve a state iteratively towards these prime factors. We show that the number of circuits evaluated in each iteration scales as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(n^{5}d)$$\end{document}O(n5d), where n is the bit-length of the number to be factorized and d is the depth of the circuit. We use a single layer of entangling gates to factorize 36 numbers represented using 7, 8, and 9-qubit Hamiltonians. We also verify the method’s performance by implementing it on the IBMQ Lima hardware to factorize 55, 65, 77 and 91 which are greater than the largest number (21) to have been factorized on IBMQ hardware.

Quantum computation is likely to revolutionize how computation is performed in the field of science, engineering and finance. Computations are performed on quantum states that make use of superposition and entanglement to allow for speedups. Future potential applications include cryptography 1 , search problems 2 , simulation of quantum systems 3 , quantum annealing 4 , machine learning 5 , computation biology 6 , quantum materials 7 , and problems in optimization 8 . Refer 9 for an extensive introduction into the field of quantum computation. Major efforts toward scaling the current algorithm focuses on developing error correction and mitigation schemes, as well as designing operations that make use of fewer ancilla qubits and gate operations 10 . In line with the current major developments, we investigate a more practical near-term scheme for prime factorization that is likely to achieve good results on noisy qubits.
Prime factorization involves expressing a composite number as the product of its prime factors. For generic large numbers that lack any structure, the quadratic sieve is the most commonly used technique. This can be computationally expensive and is exploited in RSA cryptography to guarantee information security over networks. From Shor's 11 work on period finding, it was shown that one could exploit the quantum Fourier transform to compute factors in steps that scaled polynomial in the number of bits. Subsequently, Vandersypen et al. 12 realized this experimentally by factorizing N = 15 using spin-1/2 nuclei as qubits and then Lucero et al. 13 by using superconducting qubits. A simplified version of Shor's algorithm was worked out by Geller et al. 14 for products of Fermat primes (3, 5, 17, 257 and 65537) . Later, Jian et al. 15 proposed an alternative method to compute factors by solving an optimization problem using quantum annealing demonstrated on the D-Wave quantum annealer. The largest experimental realization of general method factorization schemes includes Shor's algorithm to factorize 21 13 and optimization using the D-Wave quantum annealer to factorize 223357 15 . In addition, large numbers with specific properties have also been factorized by exploiting structure contained in the number. In Ref. 16 , the authors use a multiplication table to factorize 56153 using only 4 qubits. Despite being able to factorize large numbers with relatively few qubits by exploiting the structure in the number, these methods do not reveal the power of quantum computers against generic numbers. Studies have not been made with respect to convergence on the solution with increasing the number of qubits.
In this paper we explore how one could use imaginary time evolution to factorize numbers with relatively higher probability. Imaginary time evolution has long been used as a theoretical tool in physics to compute ground-state wave-functions 17 . Shingu et al. 18 show how using imaginary time evolution one could train a  20 showed how one could make use of variational circuits to create states that represent the dynamics of the imaginary time evolution. They use it to compute the ground state wavefunction of hydrogen and lithium hydride, while Yeter-Aydeniz et al. 21 demonstrate that the QITE algorithm serves as a quantum computing benchmark for computational chemistry methods. Herein, we develop an optimization function using the method introduced by Burges 22 and then use it as a Hamiltonian to perform imaginary time evolution on a uniform superposition of all possible considered solutions to the factorization problem. We employ variational circuits to prepare a quantum state that encodes the solution and classically train all parameters. Simulations using Python Numpy packages and IBM-QASM are used to verify the performance. The robustness of the techniques against noise is verified on IBMQ hardware for up to 5 qubits allowing factorization of numbers up to 91 with a probability over 73% . To the best of our knowledge this is the largest number that has been factorized on a quantum circuit using a general purpose algorithm that does not exploit any structural properties of the number. Each iteration involves evaluating a number of circuits that scales as O(n 5 d) where n is the bit-string length of the number to be factorized and d is the circuit depth. A few reasons make factorization an ideal candidate to be solved using Imaginary time evolution with variational circuits. The number of terms in the Hamiltonian expansion scales polynomial with respect to the number of qubits. The amplitude of coefficients in the state vector is real, which simplifies the parameterization of the landscape to be explored and updates to be made. Each iteration only needs to amplify the amplitude of the solution rather than exactly simulate the imaginary time evolution. These factors make Variational QITE a good algorithm to be studied in the context of optimization, and make prime factorization a good test bed to evaluate its working. We are hoping that future rigorous work along these lines, with the availability of more less noisy qubits, might provide for ways in which this method complements with VQE for general optimization problems approached using variational methods.

Background
The Schrödinger equation for a closed system describes dynamics according to some Hamiltonian that governs it. Allowing for time to be complex, we are able to create thermal states of specific temperature starting from a maximally mixed configuration, i.e, ρ T=1/τ = e −Hτ /Tr[e −Hτ ] . Preparing a system at low temperature or alternatively letting the system evolve to large imaginary time values we can more often sample the ground state configuration for a given Hamiltonian. We shall make use of this property to encode the required solution of the problem we intend to solve in the later sections.
The quantum state we intend to prepare using QITE is given by, where N(τ ) = 1/ �ψ(0)|e −2Hτ |ψ(0)� is a normalization factor. Equation (1) satisfies the Wick rotated Schrodinger equation Let |φ(τ )� be a state that is prepared by applying a series of unitary gates as follows We choose the initial parameters to create the state on which the evolution is to be performed. We demand that Eq. (3) is an approximate solution to Eq. (2) by demanding the norm of the variation to vanish, i.e,

Solving the above equation using McLachlan's variational principle (see Supplementary section 1 for proof) we obtain
where The derivatives on the state vectors with respect to the parameters of the circuit can be simplified into a sum of other basic circuit functions. For instance, the derivative of the R x , R y and R z rotation gates with respect to their parameter amounts to a mere shift of the parameter by an angle of π radians. We can further express H as a sum of a string of tensor products of Pauli operators allowing both A and C to be computed as the sum of terms that take the form aR(e iφ �0|U|0�) , where a is some real number and φ is a real phase. Terms of this form can be efficiently computed on quantum hardware using the Hadamard test. www.nature.com/scientificreports/ Evaluating the gradients using Eq. (5), we can then proceed to use gradient descent to update the parameters of the circuit as follows where � θ(τ ) = A −1 (τ )C(τ ) . Using sufficient layers to the variational circuit and small δτ time iterations for updating our parameters, we can prepare states that closely emulate the imaginary time evolution on the state |ψ(0)� . We would like to note that the dynamics of ITE is non unitary and thus we are only able to prepare the output state starting from a fixed initial state using variational methods.

Methods
The problem of prime factorization involves expressing a number as a product of its primes. Here we shall focus specifically on biprimes, product of 2 prime numbers. Shor's algorithm aims to solve this problem by reducing it to a period finding problem of the function a x mod N, where a is any random number and N is the number to be factorized. This algorithm uses a quantum subroutine that exploits the power of quantum Fourier transform as a part of a modified phase estimation algorithm. The largest number factorized using Shor's algorithm is 21.
Being highly sensitive to discrete errors, Shor's algorithm fails to scale without viable error correction schemes 23 .
Here we show how to factorize a given number by first casting it as an optimization problem. Let the number N be given by a product of 2 prime numbers p and q. We model the problem of finding p and q given N as an optimization problem cast as the following cost function Note that in Eq. (8), H(p, q) has global minima when p,q factorize N with a minimum value of 0. The cost function can thus be treated as a Hamiltonian whose ground state encodes the solution to our problem. We express both p and q as a binary string summation which will allow us to transfer to a qubit representation space. Assuming p has an m bit representation and q has an l bit representation, we get Note that Eq. (8) has a trivial solution N = N × 1 . To exclude this solution, we can restrict the search space of our solution by using a lower bound. The smallest prime factor we are looking for is greater than 2 and less than N/2. Thus using N − 1 qubits to represent p and q we can discard the trivial solution from showing up in our simulation. Since both p and q are prime numbers, we are allowed to set the bit 0 to be 1 in either expression making it indivisible by 2. To move to a scaled spin representation we transform all the binary variables b i that represents either the p i 's or q i 's to b i = (s i + 1)/2 . In this representation, s i takes the value ±1 such that s i = 1 maps to b i = 1 and s i = −1 maps to b i = 0 . We thus convert a factorization problem into an optimization problem over the spin variables.
Minimizing this optimization function shall result in the solution we seek. We use imaginary time evolution to evolve the output state starting from a uniform superposition over all possible solutions at T = 0 . The circuit parameters are updated according to Eqs. (5) and (7). With every iteration that is performed, the output state evolves through a time step of δτ . After sufficiently large enough time T, we expect the ground state with the least cost to prepare. For ||H|| 1 T >> 1 we get, where p and q represent the bit string solution of the binary representation of the numbers that multiply to give N. The evolution is achieved by using a variational circuit ansatz to prepares states that mimic the evolution of the system from a given starting state (as described in the previous section). The computational cost in performing one iteration of QITE is O(n 5 d) , where n is the bit-length of the number to be factored and d is the depth of the circuit (see Supplementary section 2 for proof).
Here we show an example of how to factorize N = 15 as Mapping the binary variables to spin variables yields, 10) e −HT |+ + · · · +� ≈ p q We thus obtain the factors of 3 and 5 as expected. In our simulations, we shall make use of no more than the number of qubits needed to represent our solution in the optimization. This is to limit the numbers of qubits used, maximize the computational efficiency and reduce the accumulated errors. We generate the output quantum state prepared using QITE with the ansatz shown in Fig. 1. Notice that the circuit consists of R Y rotation gates applied to each qubit followed by a layer of CNOT gates that help with the entanglement of the qubits. Using only R Y gates ensures that we maintain real amplitudes for the state when expressed in the computation basis. This helps with closely following the dynamics of QITE without introducing any nontrivial phase values.
Equation (5) for such a circuit ansatz can be re-expressed in terms of the quantum Fisher information (QFI) and the gradient of Hamiltonian expectation with respect to the current state. In classical probability, the Fisher information characterizes how much a probability distribution varies by changing a parameter that characterizes a distribution. This can be generalized and extended to talk about quantum states where one characterizes how much a state changes with respect to the parameter that governs it. Refer 24 for a brief introduction to quantum Fisher information. Given a quantum state |ψ(θ 1 , θ 2 . . . θ n )� , where θ i represents the governing parameters, the QFI is given by . For the ansatz being used to generate the state in QITE, we note that �ψ i |ψ� = R(�ψ i |ψ�) . This results in the second term vanishing due to the constant normalization of the state i.e, �ψ|ψ� = 1 . Hence we get F ij = 4M ij where M ij = R(�ψ i (θ(τ ))|ψ j (θ(τ ))�) . Note F ij being real symmetric, avoids the unstability that is likely to occur from inverting a skew symmetric matrix with small imaginary coefficients. Furthermore, since the Hamiltonian is an Hermitian operator, we can express the right hand side of Eq. (5) as follows We thus obtain j F ijθj = −2 d dθ i �H� ψ . The QFI of a given circuit and the gradient of an observable with respect to a state parameter can be directly computed in the IBM Qiskit Aqua framework (see Supplementary section 3 for QFI module).

Results
Simulations. Smaller numbers were factorized on Qiskit using the QASM simulator. For 8 qubits and more, the simulations were run using the Numpy package. The amplitude threshold for the correct solution has been set to 0.85, this amounts to a probability of 73% to get the right solution when measured in the computational basis. We have restricted our ansatz to a single layer in addition to the base layer, in light of exploring shallow circuits for the near term quantum computers. The initial circuit parameters have been randomly sampled to help avoid valleys during the training and allow for faster convergence with shallow circuits employed. We have noticed that it also helps in suppressing one of the solutions in the case of symmetrically equivalent solutions (for instance, 391 can be factorized as 17 × 23 and 23 × 17 ), allowing one solution to quickly reach the threshold. Figure 2 plots the amplitude of the solution in the computational basis against the number of iterations for sev- (16) e −HT |+ + +� ≈ |110� shows that the average number of iterations required to factorize increases with the number of qubits though we have not been able to analytically bound it. A few instances involve the amplitude of the solution flattening before continuing to raise any further, making it significantly slow. In Table 1, we note that despite some numbers being close to each other (for instance 1007 and 1081), the number of iterations for convergence to the solution seem to be widely varying. This can be attributed to the differences in bit string representation and thus the point of minima to the created cost function, and also to the fact that the parameters are being randomly initialized. Results on actual hardware for smaller instances have been presented in the next subsection. We have tested this method against a standard Variational Quantum EigenSolver (VQE), and would like to note that no statistically significant difference with regards to the number of iterations for converging to the solution has been observed.
Results on IBM-Hardware. We have made use of the publicly available IBM hardware to simulate factorization for up to 5 qubits. We start the output state with a uniform superposition over all possible solutions. This is straightforwardly instantiated by setting the parameters of R y gates in the base layer to be π 2 and the remaining parameters to be all 0. We have avoided invoking the Hadamard test to compute the M matrix and C vector as this resulted in significant error accumulation. Instead, we have reformulated the computation as indicated in the methods section and have made use of built-in Qiskit modules. The Gradient module with parameter shift method has been used to compute the C vector, while the QFI module with overlapping block diagonal method has been used to compute the M matrix. The built in QFI module makes 1024 shots on every circuit to be evaluated and includes no measurement error mitigation by default. Please refer to the Qiskit source documentation 25 for details on the implementation of these modules. Figure 3 shows the results of factorizing the numbers N = 55, 65, 77 and 91 on IBMQ-lima hardware that supports 5 qubits. Unlike the simulations that hardly show any oscillation in the amplitude, we see oscillations being introduced due to noise when run on the hardware. The convergence of the solution in the presence of hardware noise makes this method a suitable candidate for solving factorization and other similarly framed optimization problems in the near term quantum computers.

Discussion
We have shown how imaginary time evolution can be used to perform optimization and have demonstrated this for the case of prime factorization. We have shown numerical and experimental results from the factorization of several numbers represented by 7, 8, and 9-qubit Hamiltonians. We have shown that imaginary time evolution significantly populates the probability of measuring the correct solution to the factorization problem when run sufficiently long. The observed performance of the method on the IBM-lima hardware shows that the method is robust to noise and works as a proof of concept in the case of a limited number of qubits. Future work towards bounding the number of iterations analytically can help highlight the performance of this method against well known factorization techniques. We hope that the results explored here promote future work in the direction of exploring Quantum Imaginary time evolution in the context of solving optimization problems.