Where does precisely the difficulty in exponentiating a Hamiltonian $H$ in the quantum simulation problem lay?

Question

I've read in the Nielsen's, Chuang's "Quantum Computation and Quantum Information":

Classical simulation begins with the realization that in solving a simple differential equation such as $dy/dt = f(y)$, to first order, it is known that $y(t + \Delta t) \approx y(t) + f (y)\Delta t$. Similarly, the quantum case is concerned with the solution of $id|\psi \rangle/dt = H|\psi \rangle$, which, for a time-independent $H$, is just $$|\psi(t)\rangle = e^{-iHt}|\psi(0)\rangle.\ \ \ \ \ \ (4.96)$$ Since H is usually extremely difficult to exponentiate (it may be sparse, but it is also exponentially large), a good beginning is the first order solution $|\psi(t + \Delta t)\rangle \approx (I − iH \Delta t)|\psi(t)\rangle$. This is tractable, because for many Hamiltonians $H$ it is straightforward to compose quantum gates to efficiently approximate $I − iH \Delta t$. However, such first order solutions are generally not very satisfactory.

Efficient approximation of the solution to Equation (4.96), to high order, is possible for many classes of Hamiltonian. For example, in most physical systems, the Hamiltonian can be written as a sum over many local interactions. Specifically, for a system of $n$ particles, $$ H = \sum_{k=1}^L H_k,\ \ \ \ \ \ \ (4.97)$$ where each $H_k$ acts on at most a constant c number of systems, and L is a polynomial in $n$. For example, the terms $H_k$ are often just two-body interactions such as $X_i X_j$ and one-body Hamiltonians such as $X_i$. [...] The important point is that although $e^{−iHt}$ is difficult to compute, $e^{−iH_kt}$ acts on a much smaller subsystem, and is straightforward to approximate, using quantum circuits.

This may be a silly question, but I'm stuck with this one. Does the difficulty of obtaining $e^{-iHt}$ lies only in its size? Both $e^{-iHt}$ and $e^{-iH_kt}$ can be seen as matrices (of course, the first one is immensely larger than the latter one) and a Taylor series can be used to approximate both of them. This in turn boils down to just making a number of multiplications of $H$ (with different coefficients standing by the consecutive matrices). So, it makes sense for a sparse matrix to be easier to obtain, because we just don't have to do a number of multiplications, which would at the end give 0.

There are two things that come to my mind. First of which is a divide-and-conquer approach, where obtainment of $e^{-iH_kt}$ is simple and all "small" results are combined to get a big one. In fact, I think that Trotterization is this kind of approach. The second thing is a guess, that maybe $e^{-iH_kt}$ can be computed in some different way, than using Taylor series (it's a really wild guess)?

Answer 1

TL;DR: Hamiltonian simulation does not just mean "exponentiating $H$". It means finding a quantum circuit $U$ that approximates the matrix exponentiation $e^{-iHt}$. More importantly, the size of the Hamiltonian matrix $H$ isn't the key concern here. The gate complexity (or query complexity, in case the Hamiltonian is described as an oracle) of matrix exponentiation is. Simulating arbitrary matrix exponentiations using quantum circuits is computationally very expensive. By imposing specific restrictions on the structure of the local interaction $H_k$ matrices (like, $H_k$ can act on at most constant $c$ number of systems), the complexity of the simulation can be reduced.

---

Definitions

The basic goal of the Hamiltonian simulation problem, given a Hamiltonian $H$ ($2^n\times 2^n$ Hermitian matrix acting on $n$ qubits), is to find an algorithm that approximates $U$ such that $||U - e^{-iHt}|| \leq \epsilon$, where $e^{-iHt}$ is the ideal evolution and $||.||$ is the operator norm (aka spectral norm)

$$||A|| := \mathrm{max}_{|\psi\rangle\neq 0}\frac{||A|\psi\rangle||}{|||\psi\rangle||},$$ where $||\psi\rangle|| = \sqrt{\langle \psi|\psi\rangle}$ is the usual Euclidean norm of $|\psi\rangle$.

The key that determines the hardness of this problem is not so much the size of $H$ but rather the gate complexity (gate complexity and time complexity are more or less proportional). In order for the Hamiltonian simulation to be efficient, we need $U$ to be approximable by a quantum circuit containing $\mathrm{poly}(n)$ gates. Most Hamiltonians $H$ do not satisfy this criterion and are thus not efficiently simulable (you may ask for the proof for this as a separate question if you're curious!).

Commuting Hamiltonians

Fortunately, most physically occurring Hamiltonians can indeed be simulated efficiently. A very special case is the class of $k$-local Hamiltonians - these Hamiltonians can be expressed as $H = \sum_{j=1}^{m} H_j$ where each $H_j$ acts non-trivially on at most $k$-qubits. We generally assume that $m \leq \binom{n}{k} = \mathcal{O}(n^k)$. Since $k$ is supposed to be a constant, $m$ is polynomial in $n$.

It can be proved from the Solovay-Kitaev theorem that each of the individual $H_j$ operators can be simulated efficiently (cf. Ashley Montanaro's lecture notes). Now if the $H_j$'s commute we have

$$\exp(-iHt) = \exp(-i(\sum_{j=1}^m H_j)) = \prod_{j=1}^{m}\exp(-iH_jt).$$

From here it can also be shown that for any $t$ there exists a quantum algorithm that approximates the operator $e^{-iHt}$ to within $\epsilon$ in time $\mathcal{O}(m \ \mathrm{polylog}(m/\epsilon))$.

Non-Commuting Hamiltonians

However, this technique does not work for non-commuting Hamiltonian $H_j$'s since the formula $e^{-i(A+B)t} = e^{-iAt}e^{-iBt}$ does not hold for non-commuting matrices $A$ and $B$. But the Lie-Trotter product formula comes to our rescue, which says

Let $A$ and $B$ be Hermitian matrices such that $||A|| \leq K$ and $||B|| \leq K$, for some real $K$. Then $e^{-iA}e^{-iB} = e^{-i(A+B)} + \mathcal{O}(k^2)$.

Applying this formula multiple times, for any Hermitian matrices $H_1, H_2, \ldots$ satisfying $||H_j|| \leq K \leq 1$ for all $j$.

$$e^{-iH_1}e^{-iH_2}\ldots e^{-iH_m} = e^{-i(H_1+\ldots + H_m)} + \mathcal{O}(m^3K^2).$$

Therefore, there is a universal constant $C$ such that if $r \geq Cm^3(Kt)^2/\epsilon$,

$$||e^{-iH_1t/r}e^{-iH_2t/r}\ldots e^{-iH_mt/r}|| \leq \epsilon/r.$$

Hence, for any such $n$

$$||(e^{-iH_1t/r}e^{-iH_2t/r}\ldots e^{-iH_mt/r})^r - e^{-i(H_1 + \ldots + H_m)}|| \leq \epsilon$$

follows from the lemma which states that if $(U_i), (V_i)$ are sequences of $m$ unitary operators satisfying $||U_i - V_i||$ for all $1\leq i \leq m$, then $||U_m\ldots U_1 - V_m\ldots V_1||\leq m\epsilon$.

Given this result, any $k$-local Hamiltonian can be simulated simply by simulating the evolution of each term for time $t/r$ to high enough accuracy and concatenating the individual simulations. Larger the $r$, more accurate the simulation. This can be formalized as

Let $H$ be a Hamiltonian which can be written as the sum of $m$ terms $H_j$, each acting non-trivially on $k = \mathcal{O}(1)$ qubits and satisfying $||H_j|| \leq K$ for some $K$. Then, for any $t$, there exists a quantum circuit which approximates the operator $e^{−iHt}$ to within $\epsilon$ in time $\mathcal{O}(m^3(Kt)^2/\epsilon)$, up to polylogarithmic factors.

The $t^2$ dependency can be further lowered, and the complexity can be improved to $\mathcal{O}(mkt)$. You're right that there are other techniques of Hamiltonian simulation like the Taylor series ($\mathcal{O}(\frac{t\log^2(t/\epsilon)}{\log \log \frac{t}{\epsilon}}$)) and quantum walk ($\mathcal{O}(\frac{t}{\sqrt{\epsilon}})$). With the quantum signal processing algorithm it is $\mathcal{O}(t + \log \frac{1}{\epsilon})$.

References

• Spring 2016: Advanced Quantum Information Theory by Ashley Montanaro (University of Bristol)

• Wikipedia: Hamiltonian Simulation

• Andrew Childs' lecture note and slides on quantum simulation

• Understanding Quantum Algorithms via Query Complexity (2017) by Andris Ambainis

• Quantum Computational Complexity (2008) by John Watrous