%%% andrew childs, amchilds@cs.caltech.edu %%% cs20c final project report \documentstyle[11pt,graphics]{article} \setlength{\oddsidemargin}{0in} \setlength{\evensidemargin}{0in} \setlength{\marginparwidth}{0pt} \setlength{\marginparsep}{0pt} \setlength{\textwidth}{6.5in} \setlength{\textheight}{9in} \setlength{\topmargin}{0pt} \setlength{\headheight}{0pt} \setlength{\headsep}{0pt} \begin{document} \newcommand{\mod}{{\rm~mod~}} %%% HEADER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % some of this might be used by LaTeX2HTML \title{Simulating quantum computation} \author{Andrew Childs} \date{11 June 1998} \leftline{\LARGE\bf Simulating quantum computation} \bigskip \leftline{ Andrew Childs ({\tt amchilds@cs.caltech.edu})} \leftline{ CS 20c Final Project \hfill 11 June 1998} \smallskip \hrule %%% ABSTRACT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \bigskip \begin{abstract} By exploiting the features of quantum mechanics, it is possible to represent and process information in such a way that computation may be done more efficiently on a quantutum computer than on any Turing-equivalent classical information processing device. This paper summarizes the principles of quantum computation and outlines Shor's algorithm for efficient factoring. I describe a library of classes for simulating quantum computation in Java and their application to a simulation of Shor's algorithm. Finally, I present results of the simulation for several small integers. \end{abstract} %%% MAIN PAPER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Introduction} The description of algorithmic complexity in classical computer science relies on the Turing model of computation, a specific description of a computing device. Church's thesis asserts that the Turing model of computation is universal for describing ``reasonable'' computation, which can be taken to mean any computation which proceeds by mechanical application of a set of predefined rules. This assertion, which has been found to hold for all feasible classical methods of computation, allows the establishment of complexity classes that describe how the time to execute an algorithm scales with the input size regardless of the machine used. However, in the Turing model, information is taken as an entirely local quantity. On the other hand, quantum mechanical systems allow the establishment of correlations between distinct parts of a system, so that information can be stored in nonlocal correleations among all parts of a system~\cite{pre}. Fundamentally, information is a physical quantity, so we believe it should be governed by the laws of quantum mechanics. Thus a quantum description of information admits the possibility of new definitions of complexity. As Peter Shor first showed in 1994~\cite{shor}, this is in fact the case: there exist problems which we believe can be solved asymptotically faster on a quantum computer than on any classical computer. This paper describes the design and implementation of a system for simulating quantum computation on a classical computer in the Java programming language. I have written a library of routines for representing quantum information and computational operations. This library was then used to implement a simulation of Shor's efficient factoring algorithm. \section{Foundations of quantum computation} In quantum mechanics, the state of a system is represented by a complex-valued wavefunction whose amplitude squared gives the probability that, upon measurement, the system will be observed in a particular state. Then the general state of a two-state quantum system can be written $$ |\psi\rangle = \alpha |0\rangle + \beta |1\rangle $$ with the normalization condition $|\alpha|^2 + |\beta|^2 = 1$, where $\alpha$ and $\beta$ are complex coefficients. This is only one possible representation of the state; the choice of the orthonormal basis $\{|0\rangle, |1\rangle\}$ is arbitrary, but this standard basis will suffice for present purposes (note that any other representation can be achieved by a unitary transformation, i.e.\ by performing some computation). The normalization condition assures that total measurement probability is unity. With $n$ qubits, the system can be in any arbitary superposition of basis states for a Hilbert space of dimension $2^n$. For example, a three-qubit state has the form $$ |\psi\rangle = \alpha_0 |000\rangle + \alpha_1 |001\rangle + \alpha_2 |010\rangle + \alpha_3 |011\rangle + \alpha_4 |100\rangle + \alpha_5 |101\rangle + \alpha_6 |110\rangle + \alpha_7 |111\rangle, $$ where the $\alpha_j$ are complex numbers with $\sum_{j=0}^7 \alpha_j = 1$. Using base ten shorthand, we can also write this state as $$ |\psi\rangle = \sum_{j=0}^7 \alpha_j |j\rangle. $$ By the measurement postulate of quantum mechanics, measurement of a state yields the basis state $|j\rangle$ with probability $|\alpha_j|^2$ and puts the system in the state $|j\rangle$. In other words, subsequent measurement will yield that state $|j\rangle$ with certainty. Furthermore, by the ``no cloning theorem,'' it is not possible to copy a general quantum state~\cite{ste}. Doing so would allow making many copies of the state, which could then be used to determine the magnitudes of the $\alpha_j$ to any desired accuracy. Although it is possible to prepare states which are in a superposition of basis states, it is not possible to observe a general state to determine all of the coefficients $\alpha_j$. In some sense, it is easier to put information into a quantum state than it is to extract it, a limitation that requires a certain cleverness in the design of quantum algorithms. Although it is possible to formulate a quantum analog of the Turing machine~\cite{eke}, a more useful model of computation for describing real algorithms is the quantum gate array, in which a succession of transformations (gates) are applied to the set of qubits. In quantum mechanics, the time evolution of a system is governed by the action of some Hamiltonian which produces a unitary transformation on the system. This restriction is equivalent to the statement that all operations must be reversible. In quantum computation, these unitary transformations, usually applied to some small subset of a system, are refered to as quantum gates. A general $n$-qubit gate has the form $$ \hat G=\sum_{j=1}^{2^n-1}\sum_{k=1}^{2^n-1} a_{jk} | k \rangle\langle j |, $$ where $a_{jk}$ is a complex coefficient that describes how the $j$th basis state of the input affects the $k$th basis state of the output. If we represent states as complex column vectors of the form $$ |\psi\rangle = \left( \matrix{ \alpha_0 \cr \vdots \cr \alpha_n } \right) $$ then a general gate has the form $$ \hat G=\left( \matrix{ a_{00} & \cdots & a_{n0} \cr \vdots & \ddots & \vdots \cr a_{0n} & \cdots & a_{nn} \cr} \right). $$ The physical restriction to unitary transformations means that this matrix must be unitary, i.e.\ that its conjugate transpose must be its inverse. Some specific examples of gates will be discussed in the description of Shor's algorithm. \section{Simulating quantum computation} The formalisms of quantum computation in the standard basis $\{|0\rangle, |1\rangle\}$ were implemented in Java. Although this system was specifically intended for implementing a simulation of Shor's factoring algorithm, the system was made as general as possible to allow potential future use in simulating other algorithms. All code is available on the web at {\tt http://www.cs.caltech.edu/\~\negthinspace\thinspace amchilds/final}. Basis states $|j\rangle$ are represented by the {\tt BasisState} class. States are represented as arrays of boolean values that tell whether each individual qubit is 0 or 1 for that basis state. Most importantly, the {\tt BasisState} class has methods for returning the representation of a given basis state on some subspace specified by an array of qubit indices. These methods are very useful for performing gate applications, measurents, and projections on a subspace of the system. The core information data structure is the {\tt Qubits} class, which represents a fixed-size set of qubits in terms of the complex coefficients for all possible basis states. For a space of $n$ qubits, $2^n$ complex numbers are stored. In an algorithm, qubits can be created and manipulated, but the state amplitudes cannot be direcly accessed. In keeping with this physical restriction, the array of state coefficients is a private data member. An amplitude can only be extracted by performing a measurement, which modifies the system by projecting the measured subspace onto a definite basis state. Methods are provided for printing all state coefficents or amplitudes, but only for purposes of debugging or algorithm analysis. This information is not physically accessible in a single computation, although it can be thought of as the average result of many identitcally prepared computers. Because only the {\tt Qubits} class has access to the coefficients, it must perform all gate applications. Applying a gate of the same size as the number of qubits is trivial matrix multiplication, but determing how to apply gates to a subspace was one of the more challeing aspects of the project. So that complexity is readily quantifiable, algorithms must employ fixed-size gates; for example, Shor demonstrated the construction of a quantum discrete Fourier transform from one- and two-qubit gates. Therefore, essentially all computation operates on a small subspace of the entire system. Initially, I wrote a method in the {\tt Gate} class which converted a given gate into a gate on a larger space. However, I soon learned that this is not a feasible implementation for large systems. The smallest instance of Shor's algorithm uses twelve qubits represented by $2^{12}=4096$ complex numbers. Thus a transformation of the entire space requires $(2^{12})^2 \approx 1.68 \times 10^7$ complex numbers, creating an unreasonably large structure. Instead, the {\tt Qubits} class contains a method which performs the transformation on a subspace by looking up the appropriate entries in the transformation matrix during the multiplication. The time to perform the operation remains unavoidably long, but the space requirements are drastically reduced. Another difficulty was determing exactly how to map the gate matrix elements onto the larger space. The complete space can be thought of as a tensor product (the Hilbert space version of a Cartesian product) of $n$ two-dimensional spaces; even general works from the literature tend to simply state that gates can be applied to a subspace without exactly stating the effect on the qubits not involved in the transformation~\cite{ste}. Through educated guessing and trial and error, I found that the gate transformation corresponds to a block diagonal expansion when the gate operates on the first qubit. In terms of the basis states, if we let $S[j]$ denote the number formed by selecting the bits of the binary expansion of $j$ that correspond to the qubits that the gate operates on and let $\bar S[j]$ denote the number formed by removing those bits, then $$ \sum_{j=1}^{2^n-1} \sum_{k=1}^{2^n-1} a_{jk} | k \rangle\langle j | \to \sum_{j=1}^{2^N-1} \sum_{k=1}^{2^N-1} \bigg\{ \matrix { a_{jk} | S[k] \rangle\langle S[j] | & \bar S[j] = \bar S[k]\cr 0 & \bar S[j] \ne \bar S[k] }, $$ where $N$ is the number of qubits the gate operates on. The {\tt Qubits} class contains a counter for the number of gate applications. As long as fixed-size gates are used, this provides a measure of the running time of an algorithm that can be used to evaluate computational complexity. In order to perform a special kind of operation which simplifies the simulation of many algorithms, the {\tt Qubits} class also contains a method for reversibly computing a function that maps integers to integers. Assuming that the function can be computed classically in polynomial time, computing it directly via this method will not change whether the simulated algorithm can be performed in polynomial time, and it may vastly simplify the description of the algorithm. Jozsa describes both Shor's algorithm and Simon's algorithm in terms of this kind of general function~\cite{joz}. These functions are represented by the {\tt IntFunction} class, an abstract base class that can be subclassed to implement specific functions. General transformations are represented by the {\tt Gate} class. Because all operations that require access to quantum information must be performed by the {\tt Qubits} class, the {\tt Gate} class is relatively simple: it is primarily a data structure that contains the matrix elements of a transformation. There are also methods to check whether an arbitrary gate is unitary, but this functionality is not often needed. Along with the classes that simulate computation, there are two library classes that provide tools to be used in algorithms. The {\tt GateLib} class is a library of common quantum gates. It has static methods which return {\tt Gate} objects for several common one-, two-, and three-qubit gates. The {\tt AlgoLib} class is a library of algorithms which might be useful as subroutines in simulations of larger algorithms. It contains only the quantum discrete Fourier transform as constructed by Shor~\cite{shor}. Both libraries are easily extensible as new gates or subroutines become useful. \section{Shor's factoring algorithm} Shor's discovery of an algorithm for polynomial-time factoring was the first evidence that quantum computers can be asymptotically faster than classical ones. Factoring is certinaly not the only application of quantum computers; in particular, Grover's algorithm for searching an unstructured database provides a quadratic speedup over classical searching, although it does not change the complexity classification~\cite{ste}. However, Shor's algorithm is both widely studied and relatively simple to understand, so it was chosen for implementation. Shor actually solved the problem of factoring by reduction from another problem, that of finding the order $r$ of $x$ mod $n$, the least $r$ such that $x^r \equiv 1 \mod n$~\cite{shor}. Using some basic number theory, one can show that given a random $x \mod n$ coprime to $n$, the greatest common divisor of $n$ and $x^{r/2}-1$ will be a nontrivial factor of $n$ with a probability that approaces unity as $n$ increases. This condition holds as long as $n$ is odd and not a power of a prime number; in these cases, factors can be found efficiently by classical computation. All of the postprocessing can be done in polynomial time on a classical computer, so an efficient quantum algorithm for finding $r$ is effectively an efficient algorithm for factorizing $n$. In Shor's orignial version of the algorithm, no measurements are performed until the computation is complete, whereupon the whole system is measured and part of the result is used to find the order. In a slight modification used by by Jozsa~\cite{eke,joz}, part of the system is measured during the computation, projecting the system down to the subspace of the other part. This modification makes the algorithm more intuitive and faster to simulate without changing its structure, so it will be used in what follows. The algorithm for finding $r$ uses two quantum registers, the first of size $\lceil \log n^2 \rceil$ and the second of size $\lceil \log n \rceil$. For convenience, we define $q \equiv \lceil n^2 \rceil$. To start the computation, the first register is prepared in an even superposition of all states and the function $x^j \mod n$ is computed in the second register, preparing the state $$ {1 \over \sqrt{q}} \sum_{j=0}^{q-1} | j \rangle | x^j \mod n \rangle. $$ Shor describes how this function could be computed on a quantum computer, noting that keeping the input in a separate register makes the operation reversible and therefore allowed~\cite{shor}. Furthermore, although it is the most computationally intensive step of the algorithm, modular exponentiation can be done in polynomial time on a reversible classical computer, so it is clear that the quantum version can be done in polynomial time (specifically, using the algorithm given by Shor~\cite{shor}, in $O((\log n)^2 (\log \log n) (\log \log \log n))$). After preparing this state, the second register is measured. As described by Ekert and Jozsa~\cite{eke}, this picks out all the values for which $x^j \mod n$ is some arbitrary constant. Since the order of $x \mod n$ is $r$, these values have the form $x^{j+kr} \mod n$ for constant $j$ and all appropriate $k$. For simplicity, consider the case where $r$ divides $q$ exactly. In this case, the state after measurement has the form $$ {\sqrt{r \over q}} \sum_{k=0}^{q/r-1} |j+kr\rangle; $$ in other words, it is periodic with period $r$, but with a randomly chosen offset $j$. Next, the crucial step in Shor's algorithm is to perform the quantum discrete Fourier transform (QFT) $$ |j\rangle \to {1 \over \sqrt{q}} \sum_{l=0}^{q-1} \exp\left(2 \pi i {j l \over q} \right) | l \rangle $$ to pick out the peridoicity $r$ in the state of the second register. Following the treatment of Ekert and Jozsa~\cite{eke}, this operation takes the state of the machine to $$ {\sqrt{r \over q}} \sum_{k=0}^{q/r-1} {1 \over \sqrt{q}} \sum_{l=0}^{q-1} \exp\left(2 \pi i {(j+kr) l \over q}\right) | l \rangle = {\sqrt{r} \over q} \sum_{l=0}^{q-1} \exp\left(2 \pi i {j l \over q}\right) \sum_{k=0}^{q/r-q}\exp\left(2 \pi i {k r l \over q}\right) $$ The inner sum is zero unless $l$ is a multiple of $q/r$, say $l=j q/r$, in which case the sum is equal to $q/r$. Thus the state is $$ {1 \over \sqrt{r}} \sum_{j=0}^{r-1} \exp(2 \pi i j l/r) \left| j q/r \right\rangle $$ As Shor shows~\cite{shor}, the QFT can be constructed using two gates: $$ R_j = {1 \over \sqrt{2}} \left( \matrix{ 1 & 1 \cr 1 & -1} \right), $$ the normalized Hadamard transformation on the $j$th qubit, and $$ S_{j,k} = \left( \matrix{ 1 & 0 & 0 & 0 \cr 0 & 1 & 0 & 0 \cr 0 & 0 & 1 & 0 \cr 0 & 0 & 0 & \exp(i \pi / 2^{k-j})} \right), $$ on the $j$th and $k$th qubits, $j Finding order of 2 mod 15 Initializing quantum computation Allocating 8 qubits for first register, 4 qubits for second register Performing modular exponentiation on second register Measure second register: |8> QFT on first register Measure first register: |192> Finished in 29 quantum gate applications Rational approximation to 192 / 2^8 is 3 / 4 -> Order is 4 -> Possible divisor: 3 => Found a divisor: 3 Factorization complete: 15 = 3 * 5 \end{verbatim} \end{quote} Plotting the state of the first register as the computation proceeds clearly shows how the algorithm works. It is important to note that the information contained in these plots does not correspond to physically accessible information for any single instance of the algorithm, but it can be taken to be an average over many instances for which all measurements yield the same value. The random number coprime to $n=15$ is chosen as $x=2$, and $\sum |2^j \mod 15\rangle$ is computed in the second register. In other words, every basis state for which the second register is the result of modular exponentiation on the first register is given an equal amplitude, and all other basis states have zero amplitude. Measuring the second register chooses a particular value for $j$ (in this case, 8), leading to an offset in an otherwise periodic distribution, as shown in Figure 1. \begin{figure}[bht] \smallskip \includegraphics{measured.eps} \caption{\small Contents of the first register following measurement of the second register in the factorization of 15.} \end{figure} Next, the QFT is performed on the first register. As Figure 2 shows, the distribution of peaks has become completely periodic, and the four peaks correspond to the spacing between peaks in Figure 1. Measurement finds the peak at 192, which gives the proper order of 4. A factor of 15 is given by the GCD of 15 and $2^{4/2}-1=3$, which is 3, a factor of 15. Had the measurement picked out either the peak at 0 or the peak at 128, the order would be incorrectly determined, and the computation would have to be repeated. \begin{figure}[thb] \includegraphics{qft.eps} \caption{\small Contents of the first register following application of the QFT in the factorization of 15.} \end{figure} Unfortunately, due to the exponential increase in the representation of a general quantum system, the program could not be run on very many inputs. The numbers 15, 21, 33, 45, and 55 were successfully factored. The factors of 55 were found with two iterations in about 38.6 hours; larger inputs were not feasible given software, hardware, and time constraints. A logarithmic plot of the time per iteration to factor these five numbers (Figure 3) shows that the simulation time is exponential in input size $\log n$. A good general rule is that the calculation takes roughly five times longer for each additional qubit. The running time for 33 and 45 is essentially the same because both require 11 qubits for the first register. \begin{figure}[bht] \scalebox{.45}{\includegraphics{time.eps}} \hskip 24pt \scalebox{.45}{\includegraphics{qops.eps}} \caption{\small (at left) Time required to factor 15, 21, 33, 45, and 55 on a 233 MHz AMD K6 using the IBM compiler.} \caption{\small (at right) Quantum operations required to perform the QFT when factoring numbers from 15 to 128, inclusive.} \end{figure} Because the algorithm is probabilistic, the number of iterations varies from instance to instance. A more complete description of the running time would include the total time for an average instance, but this requires running each instance many times to acquire good statistics. The few instances which were repeated seem to indicate a fairly weak dependence of the number of instances on the input size. Shor shows that the dependence is logarithmic in input size (though it can be made constant using suitable postprocessing)~\cite{shor}. Therefore, this effect will not change whether the simulation time or the number of quantum operations is exponential or polynomial. Without taking into account the mean number of iterations to factor a give number, the number of quantum in the QFT is completely determined by the number of bits needed to perform the computation. Figure 4 plots this relationship for integers between 15 and 128. The need to use entire qubits introduces discretization, but the relationship is clearly linear. Full simulation would include the time to calculate modular exponentiation, which asymptotically dominates the time to perform the QFT but which is also polynomial. \section{Conclusion} By developing a library of Java classes for performing simulation of general quantum computation, I have learned a lot about how computation can be performed using a quantum mechanical representation of information. I have also learned in detail how this sort of computation differs from classical computation. Simulating Shor's algorithm has given me a firm understanding of how quantum computers can provide a theoretical speedup over classical computers. %%% REFERENCES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \begin{thebibliography}{[10]} \bibitem{pre}Preskill, J. (1997), Chapter 1 of lecture notes for Ph 229, {\tt http://www.theory.caltech.edu/ \allowbreak people/preskill/ph229/notes/chap6.ps}. \bibitem{shor}Shor, P. W. (1994), {\it Algorithms for quantum computation: Discrete logarithms and factoring}, in {\it Proceedings of the 35th Annual Symposium on Foundations of Computer Science}, IEEE Computer Society Press. \bibitem{ste}Steane, A. (1997), {\it Quantum computing}, {\tt http://xxx.lanl.gov/ps/quant-ph/9708022}. \bibitem{eke}Ekert, A. and R. Jozsa (1996), {\it Quantum computation and Shor's factoring algorithm}, Rev. Mod. Phys. {\bf 68} (3) 733--53. \bibitem{joz}Jozsa, R. (1996), {\it Quantum algorithms and the Fourier transform}, {\tt http://xxx.lanl.gov/ \allowbreak ps/quant-ph/9707034}. \bibitem{har}Hardy, G. H., and E. M. Wright (1979), {\it An Introduction to the Theory of Numbers}, Clarendon Press, Oxford. \bibitem{ibm}IBM Corporation, High Performance Compiler for Java, {\tt http://www.alphaworks.ibm.com/\allowbreak alphapreview\_tools/}. \end{thebibliography} \end{document}