A tree tensor network approach to simulating Shor's algorithm

TL;DR

This paper introduces a novel tensor network method to simulate and analyze the entanglement structure of Shor's algorithm, revealing how entanglement scales with system size and modular periodicity.

Contribution

It presents a tree tensor network approach that efficiently maps to a matrix product state, enabling detailed study of entanglement in Shor's wavefunction.

Findings

Entanglement grows exponentially with qubits before saturation.

Saturation occurs at a scale proportional to modular periodicity.

The method captures non-local entanglement structure explicitly.

Abstract

Simulating quantum systems constructively furthers our understanding of qualitative and quantitative features which may be analytically intractable. In this letter, we directly simulate and explore the entanglement structure present in a paradigmatic example of quantum information: Shor's wavefunction. The methodology employed is a dynamical tensor network which is initially constructed as a tree tensor network, inspired by the modular exponentiation quantum circuit, and later efficiently mapped to a matrix product state. Utilizing the Schmidt number as a local entanglement metric, our construction explicitly captures the wavefunction's non-local entanglement structure and an entanglement scaling relation is discovered. Specifically, we see that entanglement across a bipartition grows exponentially in the number of qubits before saturating at a critical scale which is proportional to the modular periodicity.

Equations

$$|\Psi ⟩=\frac{1}{\sqrt{r}}\sum _{i=0}^{r-1}\left(\sum _{j=0}^{\lceil (2^{2l}-1)/r\rceil }|jr+i⟩\right)\otimes |x^i\mathrm{mod}N⟩.$$

$$S=\{\begin{array}{ll}{2}^n,& \text{if}n<l_{r(\tilde{r})}\\ r(\tilde{r}),& \text{otherwise}\end{array}$$

Full text

00footnotetext: This manuscript has been authored by UT-Battelle, LLC, under Contract No. DE-AC0500OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for the United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan.

Introduction.— Tensor networks are graphical data structures consisting of nodal tensors with indexed edges which represent physical and virtual degrees of freedom. The graphical connectivity and dimensionality of the virtual degrees of freedom encode entanglement in a localized manner. Arising naturally from the need to efficiently decompose entangled many-body wavefunctions, tensor networks have been tremendously successful in numerically finding ground states of local Hamiltonians in low-dimensional condensed matter systems White_DMRG ; Orus_review and also have applications in simulating some quantum circuits Shi_2006 ; Markov_Shi ; Biamonte_2013 ; Wang_2015 .

Along with the ability to calculate local observables, tensor network representations of many-body wavefunctions provide explicit insight into entanglement structure, which differs vastly based on system dimensionality and criticalityHastings_2011 ; Schollwoeck_2011 ; Pirvu_2012 . Aside from finding ground states, tensor networks have recently successfully simulated the evolution of open quantum systemsWerner_2014 ; Schroder_2016 . Despite these benefits, finding an appropriate and efficient tensor network representation for a physical system is not a simple task. This is especially true for higher dimensional systems, for which no general efficient methods are known Orus_review .

In contrast to condensed matter states residing on a lattice, dimensionality or local geometry are ill-defined quantities for quantum information theoretic states which are generated by logical quantum circuits. Motivated by the algorithmic structure and permutational invariance of qubits involved in the modular exponentiation step, we develop a bipartite tree tensor network (TTN) structure which naturally simulates states generated through Shor’s algorithm. Inspecting the entanglement via the TTN decomposition demonstrates a clear and manifest area law violation which we discussed in detail. Additionally, our TTN state is efficiently converted into a matrix product state (MPS) tensor network so that the output state can be subjected to the quantum Fourier transform component of the algorithm. The tensor networks we develop, although not fully efficient, are useful for representing moderate size systems, with up to 39 qubit wavefunctions constructed on a laptop computer. Our dynamical tree tensor network methodology is generally applicable to many familiar information theoretic bipartite quantum states, e.g. those generated in dual register systems by the quantum phase estimation and hidden subgroup algorithmsNielsen .

Shor’s Wavefunction.— We first outline the logical operations comprising Shor’s algorithm Shor in order to develop an intuition for the form of an appropriate tensor network. To factor a natural number $N=pq$, with $p,q$ large prime numbers, we draw a random integer $x\in\mathbb{Z}_{N}$ and use Shor’s algorithm Shor to find the characteristic modular periodicity $r$ given by $x^{r}\mod N=1$. Assuming that $gcd(x,N)=1$ – in the unlikely case that $x=p$ or $q$, $N$ is trivially factored – one initializes $2l$ ($l$) qubit top (bottom) register, for a total of $3l$ qubits, where $l=\log_{2}(N)$ is the number of bits required to represent $N$. The top register is initialized into the product state $|+\rangle^{\otimes 2l}=1/(\sqrt{2})^{2l}\sum_{i=0}^{2^{2l}-1}|i\rangle$ while the bottom register is initialized as $|{1}\rangle$ which is the integer basis representation of the computational basis state $|{0,0,...,0,1}\rangle$. We represent the bottom register in the integer basis for the remainder of this letter. The composite initial product state is thus $|\Psi_{i}\rangle=|+\rangle^{\otimes 2l}_{top}\otimes|1\rangle_{bot}$.

The modular exponentiation (ME) unit (see Fig. 1 yellow box) entangles each top register qubit with the entire bottom register via controlled modular multiplication operators $U\equiv U(x,N)$. The operator $U$ is a rank $2$ tensor with dimensions $2^{l}\times 2^{l}$ and matrix elements satisfying $U|{b}\rangle=|{xb\mod{N}}\rangle$. Powers of $U^{2^{i}}$ are generated by $i$ iterative matrix multiplications. Upon application of the last controlled operator, the state reads $|{\Psi}\rangle=2^{-l}\sum_{i=0}^{2^{2l}-1}|i\rangle\otimes|x^{i}\mod N\rangle$. Because $x^{r}\mod N=1$ we may group together like bottom register basis vectors and write the state as

[TABLE]

Bipartite tensor network complexity.— Eq. 1 is by definition a bipartite Schmidt decomposition between the two registers and reveals several interesting features. We see $|\Psi\rangle$ is $r$-entangled across the bipartition, that is, the Schmidt coefficient $1/\sqrt{r}$ appears $r$ times. In the worst case $r\sim\mathcal{O}(N)$ so the inter-register entanglement scales exponentially in the number of qubits $l$ Orus_2004 . The equality of all Schmidt coefficients also foreshadows different correlations scaling compared to ground states of local Hamiltonians, which have exponentially (or power-law) decaying correlations. Eq. 1 also demonstrates that it is natural to decompose $|\Psi\rangle$ across the inter-register bipartition, and we shall retain this feature in our tree network.

The quantum Fourier transformation is known to be efficiently simulable Aharonov_2006 ; Yoran_2007 , suggesting that the non-trivial part of the computation occurs during the modular exponentiation step. We thus pose the following question: what are the entanglement properties of the basis state $\sum_{j=0}^{\lceil(2^{2l}-1)/r\rceil}|jr+i\rangle$? The top register qubits are clearly entangled with one another via their interaction with the bottom register. We therefore know from Eq. 1 that $r$ sets an upper bound on the amount of entanglement. Below we elucidate the entanglement properties of the top register basis states by developing a tensor network representation whose geometry is consistent with the inter-register bipartition and, more importantly, by the permutational invariance of the top register qubits.

A first attempt at a tensor network was performed in Ref. 6, which treats the bottom register as a qudit lying at one end of a MPS. The ME algorithm was performed by contracting two-local controlled modular multiplication gates along with a series of swap gates. In doing so, the complexity of storing the state is reduced from $\mathcal{O}(2^{3l})\rightarrow\mathcal{O}(2^{l}r)+const$ with the a constant given by $\sum_{j}2*d^{(j)}_{l}d^{(j)}_{r}$ where $d^{(j)}_{l(r)}$ are the virtual bond dimensions to the right and left of the $j$th top register qubit. While this approach was successful in simulating Shor’s algorithm, artificially large virtual bond dimensions were generated by successive swap operations. This leads to a situation where $d^{(j)}_{l(r)}=r$ for many bonds when, as we shall see, few virtual bonds of that size are necessary.

Tree Generation Algorithm.— We now introduce a natural tensor network which maintains the inter-register bipartition and distributes entanglement in an unbiased manner. This network is dynamically constructed by following the ME sub-circuit, as shown in Fig. 1, with intermediate virtual updates performed, as shown in Fig. 2 and discussed below, in between circuit operations. Our construction algorithm goes as follows. (i) Contract the $i$th controlled $U^{2^{i}}$ operator with the $i$th single qubit $|{+}\rangle$ tensor and the current bottom qudit state as shown in panel $(a)$. (ii) Perform internal operations updating and generating virtual indices. This cascades the $i$th qubit from the tree root (i.e. directly connected to the bottom register) to a new bottom tree branch as illustrated in panels $(b-e)$. Repeat the procedure for all qubits indexed by $i\in\mathbb{Z}_{2l}$.

The internal updates consists of the following steps: (i) After applying a controlled $U^{2^{i}}$ gate, a SVD separates the $i$th qubit from the root qudit (Fig. 2 panel $(b)$). The $i$th qubit is maximally entangled with the bottom register via a $\chi=2$ dimensional auxiliary edge with singular values $(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})$ (ii) Generate the new inter-register entanglement bond by performing an SVD between the bottom register and its local complement formed by the union of the new qubit and the previous tree root as indicated by the dashed purple line in panel $(b)$. Recall that this bond’s dimensionality eventually saturates at $r$. (iii) The tensors encircled by the purple box in panel $(c)$ are contracted in order to ‘lower’ the qubit through the tree before, (iv) another SVD along a bipartition, which is chosen to direct the qubit through a specific path, is performed. Step (iv) is identical to step (ii) but occurs further down the tree. Repeat steps (iii,iv) until each qubit settles into its final location at the bottom of the tree. An example of a final tree configuration is illustrated in Fig. 3.

Note that the choice of a binary tree is arbitrary and that the entanglement features discussed in the next section hold for aribtrary tree data structure. Also note that the number of virtual updates cascading the $i$th qubit is clearly upper bounded by final tree depth. Since the final tree depth is logarithmic in the number of qubits, that is with depth $\lceil\log(2l)\rceil$, the tree construction procedure is efficient. This trade-off can be compared to that in an MPS based simulation for which at least $2l$ swaps are performed. Thus, a logarithmic number of updates to generate an unbiased and natural representation of the state is well justified and we now discuss the emergent entanglement properties.

Entanglement Features.— After cascading all qubits, the tree tensor network exactly encodes the wavefunction appearing in Eq. 1. The entanglement structure for the top register, not apparent in Eq. 1, is now revealed by the holographic dimension Orus_2014 constructed by our virtual updates. In our analysis we use the Schmidt number, given by the bond dimension of the virtual index connecting bipartitions, as the metric for shared entanglement. This is an appropriate metric because, unlike ground states of local Hamiltonians which have exponential or power-law decaying Schmidt coefficients, the Schmidt coefficients are equal in magnitude. Thus the Schmidt number completely describes the entanglement which can therefore be visualized as done in Fig. 3, where the drawn bonds are weighted as $\log_{2}{(d_{i})}+1$, where $d_{i}$ is the local bond dimension which is also labeled.

At the bottom tree level where qubits first connect to their parent branches all qubits are maximally entangled to the rest of the network, with equal Schmidt coefficients $\lambda_{0}=\lambda_{1}=\frac{1}{\sqrt{2}}$. At the next level, all pairs of qubits (e.g. Fig. 3 blue highlighted qubits) are still maximally entangled with their complement, i.e. with degenerate Schmidt coefficients $\lambda_{i}=0.5$ for $i=(0,1,2,3)$. This trend, with clusters of $2^{n}$ qubits maximally entangled to the rest of the state by $2^{n}$ identical Schmidt coefficients (e.g. Fig. 3 cluster of 4 green highlighted qubits, and so on) continues up to a critical size, at which point the entanglement rapidly saturates.

The qubit cluster size at which the entanglement scaling saturates depends on which qubits are selected and is either $l_{r}=\lceil\log_{2}(r)\rceil$ or $l_{\tilde{r}}=\lceil\log_{2}(\tilde{r})\rceil$, where $\tilde{r}=r/2^{m}$ and $m$ is the largest integer such that $2^{m}$ divides $r$. The critical length scale, specific to the choice of $N$ and $x$ determining $r$, can be understood by the following arguments. The $r$ dimensional tree root bond mediates the $r$-fold entanglement across the register bipartition as per Eq. 1. Further, $r$ constrains the intra-register entanglement because all entanglement between qubits, stored in the bulk, was generated by controlled modular multiplication gates acting solely on the qudit. Descending from the tree root, bond dimensions decrease from $r$ to either $\tilde{r}$, or powers of 2 less than $r,\tilde{r}$. An example is provided by Fig. 3 (a) where we have illustrated the final tree generated for $N=3403$. Thus the entanglement scales as:

[TABLE]

where $n$ refers to the number of qubits and the saturation dimension, $r$ vs $\tilde{r}$, depends on whether qubits belong to the first $l_{r}$ qubits (as seen from left to right in Fig. 3) or to the remainder of the register. Eq. 2 defines a class states whose entanglement is reminiscent of quantum error correcting codes with entanglement set by a distance $d$ for $[[n,k,d]]$ codes Gottesman_97 .

We now address the seemingly strange feature of why the first $l_{r}$ qubits are more entangled than the others. Note modular exponentiation with at least $l_{r}$ qubits is needed to generate the $r$-fold basis vectors on either side of the register bipartition. Modular exponentiation with the remaining qubits extends the orthonormal basis vectors $|jr+i\rangle$ as the Hilbert space grows, leaving the bipartite entanglement at order $r$. To understand why the latter qubits are less entangled than the former, consider a single orthonormal Schmidt vector $\sum_{j=0}^{\lceil(2^{2l}-1)/r\rceil}|jr+i\rangle$. The qubits are projected into such a state upon the measurement of the bottom register qudit. Since $r=2^{m}\tilde{r}$ (if r is odd the algorithm restarts with a different $x$) and $i<r$, the last $m$ bits for each $jr+i\rangle$ are the identical. Qubits $1-m$ are therefore disentangled from the remaining state, and the remaining entanglement now follows the scaling law in Eq.2 saturating at $\tilde{r}$. Fig 3 (b) illustrates this point by re-plotting the tree after a qudit measurement and bond updates are performed. Note the changes in the bond dimensions along the left side of the tree.

MPS conversion and interpretation.— We now briefly comment on the conversion of a tree network to a MPS network, which is useful in simulating Shor’s algorithm in its entirety, due to an efficient representation of the QFT component for MPS systems Aharonov_2006 ; Yoran_2007 ; Wang_2015 . A computational basis state measurement $\mathcal{M}_{i}$ (illustrated by the nodes indexed $(i_{1},...,i_{l})$ in Fig. 1) projects the qudit register onto a basis state $|{i}\rangle$ from Eq. 1. We can probabilistically simulate this measurement by contracting the bottom register tensor with basis states $|{i^{\prime}}\rangle$ to find non-zero matrix elements and then choose one randomly, with each measurement outcome being equally probable.

A series of virtual updates now reduce the tree network into the MPS form. The first (leftmost) first qubit in Fig 3 is already in the MPS form because a single tensor connects a physical and virtual degree of freedom. We proceed by selecting the next qubit and contracting its parent bonds until it connects to the virtual degree of freedom to the right of qubit 1. After contracting the parent bonds of qubit 3, a decomposition is performed to its right in order to generate a virtual MPS bond not connected to qubit 4. This procedure is iteratively repeated for all remaining qubits, at which point the resulting network is an MPS. Again, the number of updates is bounded by the number of qubits and the logarithmic tree depth.

Inspecting the MPS virtual Schmidt coefficients provides a complimentary perspective to that provided by Fig. 3. In Fig. 4 we plot the virtual bond dimensions across the MPS network for 24 simulations involving the same $N=1763$ with $x$ randomly chosen. Many distinct $x$ balues share the same order $r_{x}$, so their entanglement spectrums are superimposed. The MPS spectrum in Fig. 4 verifies the entanglement scaling described by Eq. 2, namely, (i) the first $m$ qubits are disentangled, (ii) entanglement grows exponentially up to a critical length scale, and (ii) the entanglement saturates at the scale $\tilde{r}\equiv r/2^{m}$.

Discussion and conclusion.— In this letter we have explicitly constructed a bipartite tree tensor network naturally representing Shor’s wavefunction and used the resulting decompositions to elucidate the multi-partite entanglement properties. Our representation is generated in a circuit-like manner, i.e. by contracting the modular multiplication gates in series, along with a logarithmic number of decompositions between successive gates. The tensor network developed is fully compatible with the MPS formalism, allowing for an efficient translation between the two tensor network representations. Our method allows us to examine the entanglement scaling, as quantified by the Schmidt coefficients, for specific instances of Shor’s algorithm. One distinguishing feature is that the Schmidt coefficients are constant across all decompositions and irrespective of their number and magnitude. The Schmidt rank across system bipartitions involving regions smaller than a critical threshold scale exponentially in the number of qubits enclosed, i.e. regions containing less than $\log_{2}(\tilde{r})$, where $\tilde{r}\propto r$, qubits are maximally entangled to their compliment. In contrast, the entanglement for bipartitions involving larger regions saturates at $\tilde{r}$.

Our explicit entanglement analysis, for wavefunctions leading to quantum speedups, offers new insight into quantum complexity theory. We have also reiterated the difficulty in classically simulating Shor’s algorithm by demonstrating that conventional tensor network approximations cannot be made due to the equality of Schmidt coefficients. This contrasts condensed matter simulations in which truncation is a crucial numerical approximation, for example, underlying the efficiency of the 1D DMRG algorithm. It is clear that similar truncation techniques cannot apply to our system, so it is interesting to consider alternative approximations. As the computational bottleneck was provided by the full representation of the bottom register qudit, we briefly consider approximations and generalizations of our TTN to reduce this cost and allow the simulation of larger systems. A first approximation could be performed by simply eliminating a set ratio of the bottom register basis vectors. This will deform the modular multiplication operators and lead to a state similar to Eq. 1, except coefficients having support on, and being generated via the eliminated basis vectors will vanish. The modular periodicity $r$ will be partially encoded in this state and may still be extracted after the QFT.

A completely different approach would be to explicitly expand the bottom register qudit into its own tree such that the bipartite tree tensor network consists of two opposing trees. The structure of the bottom register tree would be of tremendous interest, and would likely provide insight into the complexity of the ME protocol. For this to work, the application of a modular multiplication gate should be broken down into its primitive components, i.e. modular addition, etc. vanMeter_2005 ; Pavlidis_2013 and both trees updated dynamically. Assuming the entire algorithm is not classically efficiently stimulable, an additional computational bottleneck may then arise from the presumably large (possibly exponential) number of virtual updates propagating the entanglement onto the inter-register bond.

Acknowledgements.— E. D. would like to thank R. Bennink for careful reading of the manuscript. Research sponsored by the Intelligence Community Postdoctoral Research Fellowship and the Laboratory Directed Research and Development Program of Oak Ridge National Laboratory, managed by UT-Battelle, LLC, for the U.S. Department of Energy. This manuscript has been authored by UT-Battelle, LLC, under Contract No. DE-AC0500OR22725 with the U.S. Department of Energy.