In HHL algorithm, there is a step in Quantum Phase Estimation where we have to apply powers of $e^{iAt}$ to the register (see pic). I am not able to understand how to find the quantum circuit corresponding to this unitary transformation i.e. if I know $A$, then how to find the decomposition of this operator in terms of quantum gates.

As an exapmle, Cao et al. v2 , Cao et al. v3 - in these papers the authors have decomposed the operator in terms of basic gates, but I am not able to understand, how they have done that. So, I am not able to understand how to decompose the operator $e^{iAt}$ in terms of basic gates for any Hermitian Matrix $A$ .

enter image description here

  • 4
    $\begingroup$ Converting a given Hamiltonian $A$ into a sequence of gates $e^{iAt}$ is called Hamiltonian simulation, which is a topic unto itself. This depends very much on the nature of the Hamiltonian - e.g. how local or how sparse it is, how easy it is to determine entries of the matrix, etc. HHL left the Hamiltonian simulation of $A$ as a black box, but there are lots of algorithms now known to efficiently perform Hamiltonian simulation. Wikipedia lists a couple - many of these were optimized after HHL's paper. $\endgroup$ Commented Mar 22, 2021 at 15:04
  • 1
    $\begingroup$ +1 to what @MarkS said. In addition, the circuit implementing $e^{iAt}$ depends on the choice of basic gates. In fact, the circuit ranges from trivial (when $U=e^{iAt}$ is one of the basic gates) to potentially impossible (when the set of basic gates is not universal). $\endgroup$ Commented Mar 22, 2021 at 17:15
  • 1
    $\begingroup$ I recommend you ask another question specific to the paper(s). But briefly at the bottom of FIG. 4 of version 2 of the paper the authors provide the unitary which is equivalent to their $e^{iAt}$ for their example matrix $A$ given in equation $1$ of their paper. They appear to rely on Group Leader Optimization for their Hamiltonian simulation, which I've not heard of but is discussed on this site here. Any further details you'd likely have to be more specific. $\endgroup$ Commented Mar 22, 2021 at 19:57
  • 1
    $\begingroup$ @quankid Note that you can decompose $A$ as follow: $A = 3.75II + 1.25XZ + 0.75YY + 2.25 ZX $. Now you can use trotterization... but maybe there is a clever way to do this for this specific problem. In general, given a Hermitian matrix $A$, you always can decomposed $A$ in term of pauli string, then use trotterization to approximate $e^{iAt}$. If you are using qiskit then this can be done very quickly. $\endgroup$
    – KAJ226
    Commented Mar 22, 2021 at 20:13
  • 1
    $\begingroup$ @KAJ226 May I know how did you find the coefficients of Pauli matrices. Also, can you explain in detail how can this be done using qiskit. I would be thankful if you could write a detailed answer below. $\endgroup$
    – quankid
    Commented Mar 22, 2021 at 20:20

1 Answer 1


As requested through the comment by the OP.

Given a Hermitian matrix $H$, we can always write it as linear combination of Pauli strings. That is,

$$ H = \sum_i \alpha_i P_i \hspace{1 cm} P_i \in \{I,X,Y,Z\}^{\otimes n} $$

note this linear combination can have up to $4^n$ terms. Thus, for a general Hamiltonian splitting it into linear combination of Pauli terms can be hard and inefficient. However, many interesting Hamiltonian have very efficient decomposition (the linear combination scales polynomial). For instance, the electronic structure Hamiltonian for molecule.

Now, the coefficients $\alpha_i$ can be determined through the formula $$ \alpha_i = \dfrac{1}{2^n}Tr(P_i H)$$

So here, we are given that $H = A =\dfrac{1}{4} \begin{pmatrix} 15 & 9 & 5 & -3\\ 9 & 15 & 3 & -5\\ 5 & 3 & 15 & -9 \\ -3 & -5 & -9 & 15 \end{pmatrix} $.

so for instance, $$ \dfrac{1}{2^2}Tr(II \cdot H)=\dfrac{1}{2^2} Tr\bigg(\begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \dfrac{1}{4} \begin{pmatrix} 15 & 9 & 5 & -3\\ 9 & 15 & 3 & -5\\ 5 & 3 & 15 & -9 \\ -3 & -5 & -9 & 15 \end{pmatrix} \bigg)= \dfrac{1}{2^2} Tr\bigg(\begin{pmatrix} 3.75 & 2.25 & 1.25 & -0.75\\ 2.25 & 3.75 & 0.75 & -1.25\\ 1.25 & 0.75 & 3.75 & -2.25 \\ -0.75 & -1.25 & -2.25 & 3.75\end{pmatrix} \bigg) = 3.75 $$

Similarly, $\dfrac{1}{4}Tr( XZ \cdot H) = 1.25$, $\dfrac{1}{4}Tr( ZX \cdot H) = 2.25$, and $\dfrac{1}{4}Tr( YY \cdot H ) = 0.75$.

The rest of the coefficients are $0$. For instance, $\dfrac{1}{4}Tr( IX) = 0$, and $\dfrac{1}{4}Tr( YX) = 0$, etc.

Once you have the Pauli string decomposition, you can input this into quiskit to create a circuit that approximate (if the pauli terms don't commute with each other) the operation $e^{iAt}$ through Trotterization (look at page 207 here to know more about this technique) pretty quickly using their WeightedPauliOperator as follow:

from qiskit.aqua.operators import WeightedPauliOperator
pauli_dict = {'paulis': [{"coeff": {"imag": 0, "real": 3.75 }, "label": "II"},
                         {"coeff": {"imag": 0, "real": 1.25 }, "label": "XZ"},
                         {"coeff": {"imag": 0, "real": 0.75 }, "label": "YY"},
                         {"coeff": {"imag": 0, "real": 2.25 }, "label": "ZX"},
operator = WeightedPauliOperator.from_dict(pauli_dict)
approx_circuit = operator.evolve(evo_time= 1, num_time_slices=1).decompose()

Note that you have option to change the evolution time and number of time slices in the evolve method. But note that all the Pauli terms here commute with each other so this is not an approximation but exact representation of the circuit. The reason is if $A$ and $B$ are two commuting matrix, $[A,B] = AB - BA = 0$, then $e^{A+B} = e^Ae^B $. Here $II$ commute with everything, and $[XZ,ZX] = 0$, $[ZX,YY] = 0$, and $[XZ,YY] = 0$ .

Which would output something like:

enter image description here

and if you want to transpiled to a specific gate set that you want to use, you can use the transpile function:

from qiskit.compiler import transpile
transpiled_circuit = transpile(approx_circuit, basis_gates=['cx', 'u3'], optimization_level = 3)

enter image description here

Then now, if you want to create the controlled version of $Controlled-e^{iAt}$ then you can create a controlled gate to all the gates you have in the circuit. But then you would still need to decompose those controlled operations to native gates that matches the device's native gate set. You can also do this automatically with qiskit too, actually. That is, you can just create a controlled version of the entire circuit above. Let's suppose the circuit to simulate $e^{iAt}$ is described as the transpiled_circuit we have above, then what I can do is:

from qiskit import QuantumRegister, QuantumCircuit
xs_gate = transpiled_circuit.to_gate()
cxs_gate = xs_gate.control()
qreg_q = QuantumRegister(3, 'q')
circuit = QuantumCircuit(qreg_q)
circuit.append(cxs_gate, [0,1,2])

This will create the circuit:

enter image description here

where the circuit231 represents the entire transpiled_circuit we had above. You can then again* decompose this circuit to see what it would look like in term of more elementary gates:

decomposed_circuit = circuit.decompose().decompose() 

which can give you something ridiculously long:

enter image description here

What I described is for a general case. For specific problem, you can use different tricks to make this neater. For instance, see this amazing answer by Adam.

  • $\begingroup$ KAJ226, I think there is a typo: it should be $\frac{1}{4}Tr(XZ \cdot H)$ instead of $\frac{1}{4}Tr(XZ)$ and similarly, for the other terms, am I right? $\endgroup$ Commented Mar 23, 2021 at 15:07
  • 1
    $\begingroup$ @DavitKhachatryan Yes! Thank you! :) $\endgroup$
    – KAJ226
    Commented Mar 23, 2021 at 15:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.