Advantage of simulating sparse Hamiltonians

Question

In @DaftWullie's answer to this question he showed how to represent in terms of quantum gates the matrix used as example in this article. However, I believe it to be unlikely to have such well structured matrices in real life examples, therefore I was trying to look at other methods to simulate an Hamiltonian. I have found in several articles a reference to this one by Aharonov and Ta-Shma in which, among other things they state that it is possible to have some advantage in simulating sparse hamiltonians. After reading the article, however, I haven't understood how the simulation of sparse hamiltonians could be performed. The problem is usually presented as one of graph coloring, however also looking at the presentation that @Nelimee suggested to read to study matrix exponentiation, this all falls down the silmulation through product formula.

To make an example, let's take a random matrix like:

$$ A = \left[\begin{matrix} 2 & 0 & 0 & 0\\ 8 & 5 & 0 & 6\\ 0 & 0 & 7 & 0\\ 0 & 5 & 3 & 4 \end{matrix}\right]; $$ this is not hermitian, but using the suggestion from Harrow,Hassidim and Lloyd we can construct an hermitian matrix starting from it:

$$ C = \left[ \begin{matrix} 0 & A\\ A^{\dagger} & 0 \end{matrix} \right] = \left[\begin{matrix} 0 & 0 & 0 & 0 & 2 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 8 & 5 & 0 & 6\\ 0 & 0 & 0 & 0 & 0 & 0 & 7 & 0\\ 0 & 0 & 0 & 0 & 0 & 5 & 3 & 4\\ 2 & 8 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 5 & 0 & 5 & 0 & 0 & 0 & 0\\ 0 & 0 & 7 & 3 & 0 & 0 & 0 & 0 \\ 0 & 6 & 0 & 4 & 0 & 0 & 0 & 0 \\ \end{matrix}\right]. $$

Now that I have an 8x8, 2-sparse hermitian matrix:

• Can I simulate its evolution in other ways than the product formula method?
• Even if I use the product formula, how do I exploit the fact that it is sparse? Is it just because there are less non-zero entries and therefore it should be easier to find the product of basic gates?

Answer 1

The insight that suggests that sparse matrices are useful goes along the lines of: for any $H$, we can decompose it in terms of a set of $H_i$ whose individual components all commute (making diagonalisation straightforward), $$ H=\sum_{i=1}^mH_i. $$ If the matrix is sparse, then you shouldn't need too many distinct $H_i$. Then you can simulate the Hamiltonian evolution $$ e^{-iHt}=\prod_{j=1}^Ne^{-iH_m\delta t}e^{-iH_{m-1}\delta t}\ldots e^{-iH_{1}\delta t}, $$ where $t=N\delta t$. For example, in your case, you can have $$ H_1=\frac14 X\otimes(18\mathbb{I}-6Z\otimes Z-4Z\otimes\mathbb{I}) \\ H_2=\frac14(X\otimes(11\mathbb{I}+5Z)\otimes X+Y\otimes(11\mathbb{I}+5Z)\otimes Y)\\ H_3=\frac14(11X\otimes X-Y\otimes Y)\otimes(\mathbb{I}-Z) $$ (the 3 terms corresponding to the fact that it's a 3-sparse Hamiltonian). I believe there's a strategy here: you go through all the non-zero matrix elements of your Hamiltonian and group them so that if I write their coordinates as $(i,j)$ (and I always include their complex conjugate pair), I continue adding other elements to my set $(k,l)$ provided neither $k$ nor $l$ equal $i$ or $j$.. This would mean for an $m$-sparse Hamiltonian, you have $m$ different $H_i$.

The problem is this doesn't necessarily work this straightforwardly in practice. For one thing, there's still exponentially many matrix elements that you have to go through, but that's always going to be the case with the way you're setting it up.

The way that people get around this is they set up an oracle. One possible oracle is essentially a function $f(j,l)$ which returns the position and value of the $l^{th}$ non-zero entry on the $j^{th}$ row. This can be built into a full on quantum algorithm. There are a few papers on this topic (none of which I've completely understood yet). For example, here and here. Let me try to give a crude description of the way they work.

The first step is to decompose the Hamiltonian as a set of unitaries, multiplied by positive scale factors $\alpha_i$: $$ H=\sum_i\alpha_iU_i $$ For simplicity, let's assume $H=U_1+\alpha U_2$. It might be assumed that you're given this decomposition. One then defines an operation (constructed out of controlled-$U_1$ and controlled-$U_2$) that implements $V=|0\rangle\langle 0|\otimes U_1+|1\rangle\langle 1|\otimes U_2$. If we input a particular state $|0\rangle+\sqrt{\alpha}|1\rangle$ (up to normalisation) on the control qubit, apply $V$, then measure the control qubit, post-selecting on it being in the state $|0\rangle+\sqrt{\alpha}|1\rangle$, then if the post-selection succeeds, we have implemented $U_1+\alpha U_2$, which happens with a probability at least $(1-\alpha)^2/(1+\alpha)^2$. You can do exactly the same with multiple terms, and indeed with exponentials of Hamiltonians (think about the series expansion), although in practice some better series expansions are used based on Bessel functions.