Limitations on the simulation of non-sparse Hamiltonians

Introduction

One of the primary applications of quantum computers is the simulation of quantum systems. Indeed, it was the apparent exponential time complexity of simulating quantum systems on a classical computer that led Feynman to propose the idea of quantum computation [Fey82].

In addition to predicting the behavior of physical systems, Hamiltonian simulation has algorithmic applications. For example, the implementation of a continuous-time quantum walk algorithm is a Hamiltonian simulation problem. Examples of algorithms that can be implemented using Hamiltonian simulation methods include unstructured search [FG96], adiabatic optimization [FGGS00], a quantum walk with exponential speedup over classical computation [CCDFGS03], and the recent NAND tree evaluation algorithm [FGG07].

where $d$ is the maximum number of nonzero entries in any row and $\epsilon$ is the maximum error permitted in the final state (quantified in terms of trace distance).

The dependence of (1) on the simulation time is nearly optimal, since it is not possible to simulate a general sparse Hamiltonian for time $t$ using $o(t)$ queries. Intuitively, there is no generic way to fast-forward through the time evolution of quantum systems. More formally,

For any positive integer $N$ there exists a row-computable sparse Hamiltonian $H$ with $\|{H}\|=1$ such that simulating the evolution of $H$ for time $t=\pi N/2$ within precision $1/4$ requires at least $N/4$ queries to $H$.

More recently, methods have been presented for simulating a Hamiltonian $H$ that is not necessarily sparse. Of course, we do not expect to efficiently simulate a general Hamiltonian, simply because there are too many Hamiltonians to consider (just as we cannot hope to efficiently implement a general unitary operation [Kni95]). However, we can conceivably efficiently simulate non-sparse Hamiltonians with a suitable concise description. In particular, by applying phase estimation to a discrete-time quantum walk derived from $H$, one can simulate $H$ for time $t$ in a number of walk steps that grows only linearly with $t$ [Chi08]. More precisely, we have

Measures of simulation complexity

These properties are reminiscent of the axioms for matrix norms, suggesting that it may be reasonable to quantify the complexity of simulating $H$ in terms of some matrix norm $\nu(H)$. Indeed, results on the simulation of sparse Hamiltonians are typically stated in terms of the spectral norm $\|{H}\|$, and Theorem 2 also involves matrix norms. We now introduce various matrix norms relevant to Hamiltonian simulation.

The spectral norm of a matrix $H$ is defined as

where $\|{v}\|$ is the standard Euclidean vector norm defined as \|{v}\|\mathrel{\mathchoice{\vbox{\hbox{\displaystyle:}}}{\vbox{\hbox{\textstyle:}}}{\vbox{\hbox{\scriptstyle:}}}{\vbox{\hbox{\scriptscriptstyle:}}}{=}}\sqrt{\sum_{i}|v_{i}|^{2}}.

The induced 1-norm of a matrix $H$ is defined as

where $\|{v}\|_{1}$ is the vector 1-norm defined as \|{v}\|_{1}\mathrel{\mathchoice{\vbox{\hbox{\displaystyle:}}}{\vbox{\hbox{\textstyle:}}}{\vbox{\hbox{\scriptstyle:}}}{\vbox{\hbox{\scriptscriptstyle:}}}{=}}{\sum_{i}|v_{i}|}.

The maximum column norm of a matrix $H$ is defined as

The maximum column norm is the maximum Euclidean norm of the columns of $H$. This norm appears in the complexity of an algorithm for simulating Hamiltonians whose graphs are trees [Chi08, Theorem 4] and in the related Proposition 2 in Section 5.

The max norm of a matrix $H$ is defined as

The max norm is just the largest entry of $H$ in absolute value. It is a matrix norm, and is typically much smaller than the other norms mentioned.

The following lemma relates the various norms introduced above.

Furthermore, each of these inequalities is the best possible.

The first inequality follows from the fact that the maximum element in any column cannot be greater than the Euclidean norm of that column. We have

Using the triangle inequality with $\|{H}\|=\max_{\|{v}\|=1}(\sum_{i}|\sum_{j}H_{ij}v_{j}|^{2})^{\frac{1}{2}}$, we get

Now by maximizing over all $v$ with $\|{v}\|=1$ instead of only those with $v_{j}\geq 0$, we get

The last inequality is actually an equality due to the Perron–Frobenius theorem.

For the next inequality, we use the fact that $\|{v}\|_{1}\leq\sqrt{N}\|{v}\|$ for all vectors $v$. This can be proved using the Cauchy–Schwarz inequality, $|\langle u,v\rangle|\leq\|{u}\|\|{v}\|$, by taking $u_{i}=v_{i}/|v_{i}|$. Let $j_{\max}$ be the index $j$ that maximizes $\sum_{i}|H_{ij}|$. Thus $\|{H}\|_{1}=\sum_{i}|H_{ij_{\max}}|=\|{He_{j_{\max}}}\|_{1}$. Using these two inequalities, it follows that

The last inequality is proved using the fact that for any $j$, $H_{ij}\leq\max_{i}H_{ij}$; thus

In Section 5, we discuss some more examples in which Lemma 1 can be strengthened, with emphasis on the implications for simulations.

A no–fast-forwarding theorem for dense Hamiltonians

The no–fast-forwarding theorem (Theorem 1 above) establishes a lower bound for the simulation of sparse Hamiltonians. Although we stated the theorem with $\|{H}\|=1$, any of the norms in Lemma 1 could have been used, since the Hamiltonian used in the proof of the no–fast-forwarding theorem is $2$-sparse, and by (15) the norms differ at most by a factor of $2$. In particular, the theorem could be restated with $\max(H)\leq 1$ or $\|{H}\|_{1}\leq 2$.

In terms of this black-box model, we have the following:

The main idea, as in the proof of Theorem 1 [BACS05], is to construct a Hamiltonian whose simulation for time $t=\pi N/2$ determines the parity of $N$ bits. Since we know that computing the parity of $N$ bits requires at least $N/2$ queries [BBCMW01, FGGS98], this Hamiltonian cannot be simulated with $o(N)$ queries. Moreover, we want this Hamiltonian to be non-sparse.

We start with a simple Hamiltonian $H_{1}$ whose graph is just a line with $N+1$ vertices. Consider the Hamiltonian acting on vectors $|i\rangle$ with $i\in\{0,\ldots,N\}$. The nonzero matrix entries of $H_{1}$ are $\langle i\left|H_{1}\right|i+1\rangle=\langle i+1\left|H_{1}\right|i\rangle=\sqrt{(N-i)(i+1)}/N$ for $i\in\{0,1,\ldots,N-1\}$. This Hamiltonian has $\|{H_{1}}\|=1$, and simulating $H_{1}$ for $t=\pi N/2$ starting with the state $|0\rangle$ gives the state $|N\rangle$ (i.e., $e^{-iH_{1}t}|0\rangle=|N\rangle$).

Now, as in Ref. [BACS05], consider a Hamiltonian $H_{2}$ generated from an $N$-bit string $S_{0}S_{1}\ldots S_{N-1}$. $H$ acts on vertices $|i,j\rangle$, with $i\in\{0,\ldots,N\}$ and $j\in\{0,1\}$. The nonzero matrix entries of this Hamiltonian are

for all $i$ and $j$. By construction, $|0,0\rangle$ is connected to either $|i,0\rangle$ or $|i,1\rangle$ for any $i$; it is connected to $|i,j\rangle$ if and only if $j=S_{0}\oplus S_{1}\oplus\ldots\oplus S_{i-1}$. Thus $|0,0\rangle$ is connected to either $|N,0\rangle$ or $|N,1\rangle$, and determining which is the case determines the parity of $S$. The graph of this Hamiltonian consists of two disjoint lines, one of which contains $|0,0\rangle$ and either $|N,0\rangle$ or $|N,1\rangle$ depending on the parity of $S$. Just as for $H_{1}$, starting with the state $|0,0\rangle$ and simulating $H_{2}$ for time $t=\pi N/2$ will give either $|N,0\rangle$ or $|N,1\rangle$, which determines the parity of $S$. Note that since $H_{2}$ is a permutation of $H_{1}\oplus H_{1}$, $\|{H_{2}}\|=\|{H_{1}}\|=1$.

Finally, we construct the dense Hamiltonian $H$ that has the properties stated in the theorem. As before, $H$ is generated from an $N$-bit string $S_{0}S_{1}\ldots S_{N-1}$. $H$ acts on vertices $|i,j,k\rangle$, with $i\in\{0,\ldots,N\}$, $j\in\{0,1\}$, and $k\in\{0,N-1\}$. The nonzero entries of $H$ are given by

for all $i$, $j$, $k$, and $k^{\prime}$. The graph of $H$ is similar to that of $H_{2}$, except that for each vertex in $H_{2}$, there are now $N$ copies of it in $H$. This Hamiltonian is dense because it has $\Theta(N^{2})$ vertices and each vertex is connected to all $N$ copies of its neighboring vertices, which gives at least $N$ nonzero entries in each row.

Now, just as before, the parity of $S$ can be determined by simulating $H$ for time $t=\pi N/2$. This gives the lower bound of $N/2$ queries.

A stronger limitation for dense Hamiltonians

The currently known dense Hamiltonian simulation algorithms rely on certain properties of the Hamiltonian that we call its structural properties. By this we mean the location of nonzero entries in $H$, which correspond to the location of edges in the graph of the Hamiltonian, and the magnitudes of the edge weights. (The remaining information about the Hamiltonian is the phase of each matrix entry $H_{ij}$.)

To show the lower bound, we need a black-box problem with an $\Omega(\sqrt{N/\log N})$ average-case lower bound, and a set of Hamiltonians whose simulation would solve this problem. We consider the problem of distinguishing strings $s\in\{-1,+1\}^{M}$ that have sum $-B$ or $+B$, given a black box for the entries of the string. When queried with an index $i\in\{1,2,\ldots,M\}$, the black box returns the value of $s_{i}\in\{-1,+1\}$, where $s=s_{1}s_{2}\ldots s_{M}$. The following lemma characterizes the query complexity of this problem.

Suppose we are given black-box access to a string $s\in\{-1,+1\}^{M}$, where $s$ is chosen uniformly at random from the set of strings with $\sum_{i}s_{i}\in\{-B,+B\}$. Then determining $\sum_{i}s_{i}$ has average-case quantum query complexity $\Theta(M/B)$.

Thus, determining whether the sum is $-\sqrt{M\log M}$ or $+\sqrt{M\log M}$, with the promise that one of these is the case, requires $\Omega(\sqrt{M/\log M})$ quantum queries on average. For each string $s$, we construct a Hamiltonian $H_{s}$ whose simulation for a particular time allows us to distinguish the two possible cases assuming $s$ satisfies the promise.

Let $H_{s}$ be a symmetric circulant matrix of size $N\times N$, where $N=2M+1$ is odd. (A circulant matrix is a matrix in which each row is rotated one element to the right relative to the preceding row.) In general, a circulant matrix is completely specified by its first row. However, since $H_{s}$ is a symmetric circulant matrix, it is completely specified by the first $M+1$ entries of the first row. Let the first entry of the first row be 0, and the next $M$ entries of the first row be $s_{1},s_{2},\ldots,s_{M}$. In other words, the first $M+1$ entries of the first row of $H_{s}$ are followed by the string $s$. Then the remaining entries of the first row are $s_{M},s_{M-1},\ldots,s_{1}$.

Given a black box for the entries of $s$, we can easily construct a black box for the entries of $H_{s}$. Indeed, one query to $H$ can be simulated with at most one query to the string $s$. Sometimes no query to $s$ is needed, since the diagonal entries of $H_{s}$ are always 0.

Since $H_{s}$ is a circulant matrix, it is diagonalized by the discrete Fourier transform. Its eigenvalues $\lambda_{0},\lambda_{1},\ldots,\lambda_{N-1}$ are

Thus the time evolution of $H_{s}$ can be used to learn whether $\sum_{j}s_{j}$ is $-\sqrt{M\log M}$ or $+\sqrt{M\log M}$. Since $\lambda_{0}=2\sum_{j}s_{j}$, the two cases can be distinguished by determining the sign of $\lambda_{0}$. Note that we know the eigenvector corresponding to $\lambda_{0}$: it is the first column of the discrete Fourier transform matrix, i.e., the uniform superposition over all computational basis states.

Using Lemmas 3 and 4, we can upper bound the probability that $\|{H_{s}}\|$ is large when $s$ is chosen uniformly at random from $\mathcal{P}$. If $X$ is the event that $\|{H_{s}}\|\geq 4d\sqrt{M\log M}$ and $Y$ is the event that $s\in\mathcal{P}$, then $\Pr(X)$ is given by Lemma 3 and $\Pr(Y)$ is given by Lemma 4. In these terms, we can compute an upper bound for $\Pr(X|Y)$ as follows:

Thus the average-case query complexity of the claimed algorithm is $O\left((\log M)^{c}\right)$, which violates the lower bound of $\Omega(\sqrt{M/\log M})$. ∎

The proof technique above can be extended to rule out algorithms with query complexity sub-exponential in $(\|{Ht}\|,\log N)$ as well, by changing the promised set (i.e., the value of $B$ used in Lemma 2) and choosing a larger value of $d$ in Lemma 3. Exponential functions of $(\|{Ht}\|,\log N)$ cannot be ruled out, of course, since any Hamiltonian can be simulated by making $O(N^{2})$ queries, which is exponential in $\log N$. On the other hand, if we insist that the query complexity of an algorithm depends only on $\|{Ht}\|$ (and not $\log N$), then the proof above can be modified to rule out algorithms whose time complexity is an arbitrary function of $\|{Ht}\|$. For example, there exists no Hamiltonian simulation algorithm that makes $\exp(\exp(\|{Ht}\|))$ queries.

Finally, we emphasize that even though the above proof involves average-case complexity and distributions over inputs, Theorem 4 is a statement about the worst-case complexity of simulating Hamiltonians.

Simulation complexity for structured Hamiltonians

The matrix $U$ is diagonal. To define $U_{ii}$, we arbitrarily fix some vertex as the root and consider the unique path from the root to vertex $i$. Let the path contain the vertices $i_{0},i_{1},\ldots,i_{p-1},i_{p},i$, where $i_{0}$ is the root and $i_{p}$ is the parent of $i$. For each nonzero entry of $H$, define \alpha_{ij}\mathrel{\mathchoice{\vbox{\hbox{\displaystyle:}}}{\vbox{\hbox{\textstyle:}}}{\vbox{\hbox{\scriptstyle:}}}{\vbox{\hbox{\scriptscriptstyle:}}}{=}}H_{ij}/|H_{ij}|. Then let U_{ii}\mathrel{\mathchoice{\vbox{\hbox{\displaystyle:}}}{\vbox{\hbox{\textstyle:}}}{\vbox{\hbox{\scriptstyle:}}}{\vbox{\hbox{\scriptscriptstyle:}}}{=}}1 if $i$ is the root and U_{ii}\mathrel{\mathchoice{\vbox{\hbox{\displaystyle:}}}{\vbox{\hbox{\textstyle:}}}{\vbox{\hbox{\scriptstyle:}}}{\vbox{\hbox{\scriptscriptstyle:}}}{=}}\alpha_{i_{0}i_{1}}\alpha_{i_{1}i_{2}}\cdots\alpha_{i_{p-1}i_{p}}\alpha_{i_{p}i} otherwise.

Since $U$ is diagonal, $(UHU^{\dagger})_{ij}=U_{ii}H_{ij}U_{jj}^{*}$. If $i$ and $j$ are not adjacent in the tree, then $(UHU^{\dagger})_{ij}=H_{ij}=0$ as required. Otherwise, suppose without loss of generality that $j$ is the parent of $i$. Then $U_{ii}U_{jj}^{*}=|\alpha_{i_{0}i_{1}}|^{2}|\alpha_{i_{1}i_{2}}|^{2}\cdots|\alpha_{i_{p-1}i_{p}}|^{2}\alpha_{ji}=H_{ji}/|H_{ij}|$, so $(UHU^{\dagger})_{ij}=|H_{ij}|$ as claimed. ∎

If $H$ can be expressed as the sum of a small number of Hamiltonians, each of whose graph is a forest, then $H$ can be efficiently simulated when $\|{H}\|$ is small. Recall that a graph is said to have arboricity $k$ if its adjacency matrix can be written as the sum of the adjacency matrices of $k$ forests, but not $k-1$ forests.

We begin by considering the case of a star graph. A star graph is a tree on $n$ vertices with one vertex having degree $n-1$ and the others having degree $1$ (i.e., the complete bipartite graph $K_{1,n-1}$). We show that if $S$ is a Hamiltonian whose graph is a star,

To show the result for graphs of arboricity $k$, we begin by showing how a rooted tree can be decomposed into the sum of two forests of stars. The first forest contains all the edges in which the parent vertex is at a even distance from the root. The second forest contains the rest of the edges. This decomposes a rooted tree into two forests of stars, and similarly decomposes a forest into two forests of stars. Since the Hamiltonian has arboricity $k$, it can be decomposed into $k$ forests, which can be decomposed into $2k$ forests of stars.

Open questions

Acknowledgments

We thank Aram Harrow for suggesting the use of Hoeffding’s inequality to simply the proof of Lemma 3. This work was supported by MITACS, NSERC, QuantumWorks, and the US ARO/DTO.

Appendix: Proofs of lemmas

In this appendix, we prove Lemmas 2, 3, and 4.

We show the lower bound by first showing the same lower bound for the worst-case problem using the quantum adversary method [Amb02] and then reducing the worst-case problem to the average-case problem.

For the worst-case lower bound, we use the notation of Theorem 2 of Ref. [Amb02]. We require two sets of inputs $X$ and $Y$ that have different outputs. Let $X$ be the set of all strings for which $\sum_{i}s_{i}=-B$, and $Y$ be the set for which $\sum_{i}s_{i}=+B$. We define a relation between the sets as follows. Let an element $x\in X$ be related to an element of $y\in Y$ if and only if $y$ can be reached from $x$ by changing exactly $B/2$ $-1$s to $+1$s in the string $x$. Note that a string in $X$ has exactly $(M/2+B/2)$ $-1$s and $(M/2-B/2)$ $+1$s.

Using these sets $X$ and $Y$, and the relation defined above, it is easy to see that

The adversary method now provides a lower bound of $\Omega\left(\sqrt{{mm^{\prime}}/{ll^{\prime}}}\right)=\Omega(M/B)$ for the worst-case query complexity of this problem.

The worst-case query complexity can now be reduced to the average-case query complexity under the uniform distribution over all input strings satisfying the promise. To do this, we first apply a uniformly random permutation to the input string, and then with probability $\frac{1}{2}$ multiply all the entries by $-1$ (and leave them unchanged with probability $\frac{1}{2}$). The resulting distribution is now uniform over all input strings satisfying the promise. If the string is not multiplied by $-1$, then the output of the permuted string is the same as the input string. If all the entries are multiplied by $-1$, then the output of the modified string is the opposite of that of the original input.

The lower bound is tight due to a matching upper bound provided by the algorithm for approximate quantum counting [BHMT02]. To distinguish the two types of inputs, we can approximately count the number of $+1$s to accuracy $\epsilon={B}/{2M}$, which requires $O(M/B)$ queries. ∎

The eigenvalues of $H_{s}$ are $\lambda_{r}=2\sum_{j=1}^{M}s_{j}\cos\frac{2\pi jr}{N}$, where $r\in\{0,1,\ldots,N-1\}$. We wish to bound the probability that $\lambda_{r}$ is large, so as to bound the probability of $\|{H_{s}}\|=\max_{r}|\lambda_{r}|$ being large. This is achieved by applying Hoeffding’s inequality [Hoe63, Theorem 2].

If $X_{1},X_{2},\ldots,X_{M}$ are independent and $a_{j}\leq X_{j}\leq b_{j}$ for all $1\leq j\leq M$, then for any $t>0$, we have

Since a similar inequality holds when $X_{j}$ is replaced by $-X_{j}$, we get

Finally, since $\|{H_{s}}\|=\max_{r}|\lambda_{r}|$, a union bound gives

Of the $2^{M}$ strings of length $M$, those with sum $-\sqrt{M\log M}$ or $+\sqrt{M\log M}$ have either $\frac{1}{2}(M+\sqrt{M\log M})$ $+1$s or $\frac{1}{2}(M+\sqrt{M\log M})$ $-1$s. Thus the total number of such strings is

We can asymptotically approximate this expression using a well-known approximation for the binomial coefficients (see for example equations 4.5 and 4.10 of Ref. [Odl95]), which states that

provided $|k-n/2|=o(n^{2/3})$. Applying this to (34), we get

References