Chapter 4

Here I document the approach via Sympy. This is not really a feasable approach since Matrix Multiplication in Sympy is extremely slow. Of course this is not entirely unexpected since sympy is not meant to be used in a brute force environment.

In fact, finding the circuit of Figure 4.8 via brute force is out of reach for this approach. It can't be done within "reasonable" time limits. Going through all combinations of just two gates (just 25=5*5 cases for the five admissible gates) already takes one second. For three gates the number raises to almost nine seconds. So several minutes are to be expected in case of five gates. I count this as "unreasonable" since this is still a very small problem.

Conclusion: Do not use sympy to calculate lots of matrix products.

IMPORTANT: Complex numbers in numpy are always implemented via floats. To avoid the typical floating point arithmetic issues we assume that all numbers are either integrals or more generally multiples of some 2**(-n). In that case floating point arithmetic is exact (up to overflow).

The eigenvalues are \(\pm1\) for all Pauli matrices. The eigenvectors (tip: Z.eigenvects()) are:

Recall the correspondence between the state vector and the parameterization of the Bloch Sphere:

and (using \(\sin\theta=2sc\), \(\cos\theta=c^2-s^2\))

From this it is not hard to see that the \(+1\) eigenvectors of X, Y, Z correspond to \(\hat{x}\), \(\hat{y}\), \(\hat{z}\). The \(-1\) eigenvectors correspond to \(-\hat{x}\), \(-\hat{y}\), \(-\hat{z}\).

The equation follows directly from the polynomial series of \(\exp\), \(\sin\), and \(\cos\).

It is easy to see that \(T = \exp(i\pi/8) \cdot R_z(\pi/4)\). We let sympy do the job:

rz = Rz.subs(theta, pi/4)
exp(I*pi/8) * rz - T  # should be zero

Matrix([[0, 0], [0, 0]])

QED.

Clearly a mere product of two factors can't produce the Hadamard gate. Hence we try to find \(\alpha,\beta,\gamma\) such that \(R_z(\alpha)\cdot R_x(\gamma) \cdot R_z(\beta)\) is the Hadamard up to the phase factor. Due to the sqrt of 2 factor appearing in \(H\) we deduce that \(\gamma=\pm\pi/2\) is a good bet. Knowing what comes we choose \(\gamma=\pi/2\).

a, b = sp.symbols('\\alpha \\beta')

ra = Rz.subs(theta, a)
rb = Rz.subs(theta, b)
rx = Rx.subs(theta, pi/2)

h0 = ra * rx * rb

A short calculation leads to (use the code above to do it):

\[ R_z(\alpha)\cdot R_x(\pi/2) \cdot R_z(\beta) = \frac{1}{\sqrt{2}} \cdot \left[\begin{matrix}e^{0.5 i \left(- \alpha - \beta\right)} & - 1.0 i e^{- 0.5 i \left(\alpha - \beta\right)}\\- 1.0 i e^{0.5 i \left(\alpha - \beta\right)} & e^{0.5 i \left(\alpha + \beta\right)}\end{matrix}\right] \]

Hence setting \(\alpha = \beta = \pi/2\) and \(\phi = \pi/2\) we get

\[ H = e^{i\pi/2} \cdot R_z(\pi/2)\cdot R_x(\pi/2) \cdot R_z(\pi/2) \]

This is easy to see from the following calculation:

nx, ny, nz = sp.symbols('n_x n_y n_z')
r = nx*X + ny*Y + nz*Z

sp.simplify(r*r)

Matrix([[n_x**2 + n_y**2 + n_z**2, 0], [0, n_x**2 + n_y**2 + n_z**2]])

By assumption we have \(n_x^2 + n_y^2 + n_z^2 = 1\), which implies the claim.

Alternatively one can avoid calculating with matrices by using the well known (anti) commutator relations between the Pauli Operators like \(XY=iZ=-YX\). This shows that in the expanded form of \((\hat{n}\cdot\sigma)^2\) only the squares of the Pauli Operators survive. Now use \(X^2=Y^2=Z^2=I\) to obtain the claim.

We divide the proof into several steps.

Let \(|\psi\rangle = \cos\left(\theta /2\right) |0 \rangle \, + \, e^{i\varphi} \sin\left(\theta /2\right) |1\rangle\) and recall the parameterization \(p = (\sin\theta \cos\varphi,\; \sin\theta \sin\varphi,\; \cos\theta)\) of the sphere.

PROOF: Clearly every unitary Matrix has the form:

\[ U = e^{i\gamma} \begin{pmatrix} e & -\overline{f} \\ f & \overline{e} \end{pmatrix} \text{ where } |e|^2 + |f|^2 = 1, \]

since the columns have to be orthogonal. In the following we show that the rotation operators are characterized as those unitary operators which look like the RHS without the phase factor.

We have:

\[ R_{\hat{n}}(\theta) = e^{-i\theta N/2} = \left(\begin{matrix} \cos(\theta/2) - iz \sin(\theta/2) & -i\overline{a} \sin(\theta/2) \\ -ia \sin(\theta/2) & \cos(\theta/2) + iz \sin(\theta/2) \end{matrix}\right) \]

where \(\hat{n} = (x, y, z)\) and \(a = x + iy\). It is sufficient to show that the First column of the rotation opterator can be made equal to \((e, f)\), since the second column of \(U\) is fixed by unitarity - up to a phase factor (this is where the \(\alpha\) kicks in).

Abbreviating \(\zeta = \cos(\theta/2) - iz \sin(\theta/2)\) and decomposing \(a = e^{i\varphi}|a|\) we see that the first columnt of the rotation is:

\[ (\zeta, -i e^{i \varphi} \sqrt{1 - |\zeta|^2}) \]

Clearly \(\zeta\) can be chosen to be any complex number with modulus at most \(1\). Once this is fixed, the second component can be made any number with modulus \(\sqrt{1-|\zeta|^2}\). Hence the system

has a (unique) solution, which proves the claim. QED.

We follow the proof of Part 1 and first represent \(H\) in this special form:

\[ H = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} = \frac{i}{\sqrt{2}} \begin{pmatrix} -i & -i \\ -i & i \end{pmatrix} \]

Hence \(\zeta=-i/\sqrt{2}\) and \(e^{i\varphi}=1\). This in turn leads to \(x=z=1/\sqrt{2}\), \(z=0\) and \(\theta=\pi\). In other words:

\[ H = i R_{(\hat{x}+\hat{z})/\sqrt{2}}(\pi) = \frac{1}{\sqrt{2}} \left( X + Z \right). \]

Recall \(Z = i R_z(\pi)\), hence:

\[ S = e^{i\pi/4} R_z(\pi/2). \]

By Theorem 4.1 we find the following representation of \(HUH\):

\[ HUH = e^{i\alpha} R_z(\beta)\cdot R_y(\gamma)\cdot R_z(\delta) \]

Conjugating this again with \(H\) we get:

\[ U = e^{i\alpha} R_x(\beta)\cdot R_y(-\gamma)\cdot R_x(\delta) \]

QED.

Since we consider only products of Gates we may identify Gates resp. Qubits which are equivalent up to a phase factor. That is we work on the quotient space \(\CC^2/\CC\) for the Qubits and \(\CC^{2\times2}/\CC\) for the Gates.

We have to proof that:

\[ U \equiv R_{\hat{n}}(\beta) R_{\hat{m}}(\gamma) R_{\hat{n}}(\delta) \]

In exercise 4.8 we proved that \(U\equiv R_{\hat{k}}(\varepsilon)\) for some axis and angle. Let us denote by \(r\) the (ordinary) rotations in the Bloch Sphere. By exercise 4.6 we know that the quantum roations are isomorphic to the corresponding ordinary rotations. That is, we have to show:

\[ r_{\hat{k}}(\varepsilon) = r_{\hat{n}}(\beta) r_{\hat{m}}(\gamma) r_{\hat{n}}(\delta) . \]

On the other hand this is a well known fact about the group of rotations (real orthogonal matrices with determinant equal to 1) in three dimensions. We prove it here anyway:

The above equation is equivalent to

\[ r_{\hat{n}}(-\beta) = r_{\hat{m}}(\gamma) r_{\hat{n}}(\delta) r_{\hat{k}}(-\varepsilon) . \]

Since \(\beta\) is a free parameter it suffices to show that the RHS maps \(\hat{n}\) to itself, for appropriate choices of \(\gamma\) and \(\delta\). Therefore let \(\hat{t} = r_{\hat{k}}(-\varepsilon)\hat{n}\) and consider:

\[ \hat{n} = r_{\hat{m}}(\gamma) r_{\hat{n}}(\delta) \hat{t} . \]

Consider the Meridian \(M_1\) around \(\hat{n}\) which contains \(\hat{t}\) and the Meridian \(M_2\) around \(\hat{m}\) which contains \(\hat{n}\). Since the two axes are not parallel the two Meridians have two intersection points \(a\) and \(b\) (they can be equal in the trivial case \(\hat{t}=\hat{n}\)). Clearly we can choose \(\delta\) in such a way that \(\hat{t}\) moves to one of the intersection points, say \(a\). Now, by definition of \(M_2\) we can choose \(\gamma\) in sich a way that \(a\) moves to \(\hat{n}\). QED.

First let us find the parameters in

\[ H = e^{i\alpha} R_z(\beta) R_y(\gamma) R_z(\delta) \]

The RHS is given explicitly in (4.12) from where it is immediat that \(\gamma=\pi/2\). The other three parameters satisfy a linear system of equations, which can be solved easily: \(\alpha=\pi/2\), \(\beta=0\), and \(\delta=\pi\).

Hence:

The first two identities follow from the fact that \(H\) is (unitary and) self-adjoint and contains eigenvectors of X as columns (and rows). The third identity follows from

\[ H = \frac{1}{\sqrt{2}} \left( X + Z \right). \]

and the commutator properties of the Paulis, e.g. \(XY=iZ\) and its cyclic variations. Of course the first two identities could also be proved like that.

This follows from \(T\equiv R_z(\pi/4)\) (up to phase factor) and \(HZH=X\).

We only proof the first assertion since the second one follows trivially from the first.

We use the previously established fact that the 3D-rotations are isomorphic to the Pauli rotation operators (exercise 4.6). Along that way we essentially establish the relation between Quaternions and 3D rotations (and Pauli rotations).

For a vector \(\hat{n}\) let us define \(N_{\hat{n}}=-i(n_1X+n_2Y+n_3Z)\). The composition of the two rotations is given by:

\[ (c_{12} + s_{12} N_{\hat{n}_{12}}) = (c_2 + s_2 N_{\hat{n}_2}) \cdot (c_1 + s_1 N_{\hat{n}_1}) . \]

To spare us from the tedious work to expand the RHS by hand we use sympy for that. We are only interested in algebraic operations and in particular we don't want to see any cumbersome matrix expression. To help sympy we use the fact that the Quaternions can be modelled by the Pauli Matrices. Just use the Identity matrix as the unit (the real number 1) and \(-iX\), \(-iY\), \(-iZ\) as the three imaginary units. Keep in mind that we only use the Quaternions as a trick to tell sympy that we are not interested in Matrices. Essentially we still use the Pauli Rotations to prove the claim.

x1, y1, z1 = sp.symbols('x1 y1 z1')
x2, y2, z2 = sp.symbols('x2 y2 z2')
c1, c2, s1, s2 = sp.symbols('c1 c2 s1 s2')

N1 = Quaternion(0, x1, y1, z1)
N2 = Quaternion(0, x2, y2, z2)

R1 = c1 + s1*N1
R2 = c2 + s2*N2

result = R2 * R1

# Formatting for readability:
str(result).replace(' + (', '\n+ (')

(c1*c2 - s1*s2*x1*x2 - s1*s2*y1*y2 - s1*s2*z1*z2)
+ (c1*s2*x2 + c2*s1*x1 - s1*s2*y1*z2 + s1*s2*y2*z1)*i
+ (c1*s2*y2 + c2*s1*y1 + s1*s2*x1*z2 - s1*s2*x2*z1)*j
+ (c1*s2*z2 + c2*s1*z1 - s1*s2*x1*y2 + s1*s2*x2*y1)*k

From here we can directly read out the claim! QED.

According the convention of the book, these two circuits are represented by \(I\otimes H\), \(H\otimes I\) (in that order). We note here that Qiskit uses the reverse order. The Kronecker Product yields the matrices (since the basis vectors are ordered in a canonical way: 00, 01, 10, 11):

\[ I\otimes H = \frac{1}{\sqrt{2}} \,\begin{bmatrix} 1 & 1 & 0 & 0 \\ 1 & -1 & 0 & 0 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 1 & -1 \,\end{bmatrix} \]

\[ H\otimes I = \frac{1}{\sqrt{2}} \,\begin{bmatrix} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ 1 & 0 & -1 & 0 \\ 0 & 1 & 0 & -1 \,\end{bmatrix} \]

The general form of a CU-gate is (\(P_j\) being projections):

\[ CU = P_1 \otimes I + P_2 \otimes U \]

Hence by \(HZH=X\):

\[ CX = I\otimes H \cdot CZ \cdot I\otimes H \]

This follows from

\[ CZ = P_1 \otimes I + P_2 \otimes Z = I \otimes P_1 + Z \otimes P_2 . \]

The density matrix for two qubits in the computational basis looks as follows:

\[ \rho = \sum_{i,j=0}^1 p_{ij} |ij\rangle . \]

CX maps \(|0j\rangle\) to itself, and it swaps \(|10\rangle\) with \(|11\rangle\). Hence, as a mapping on the density matrix it acts as follows on the density matrix:

\[ p_{0i} \mapsto p_{0i},\quad p_{10} \mapsto p_{11},\quad p_{11} \mapsto p_{10} . \]

The statement after "Moreover" is direct consequence of the first statement. Therefore we only prove the first one.

We have to show that \(H\otimes H \cdot CX \cdot H\otimes H\) is the CX gate with the qubits reversed.

First, let us recall the projections on the eigenbasis of X:

Consider

\[ H\otimes H \cdot CX \cdot H\otimes H = P_+ \otimes I + P_- \otimes Z . \]

We have to show that, after swapping the two qubits the RHS equals CX. In fact:

\[ I \otimes P_+ + Z \otimes P_- = \frac{1}{2} \,\begin{bmatrix} 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 1 & 1 \,\end{bmatrix} \,+ \frac{1}{2} \,\begin{bmatrix} 1 & -1 & 0 & 0 \\ -1 & 1 & 0 & 0 \\ 0 & 0 & -1 & 1 \\ 0 & 0 & 1 & -1 \,\end{bmatrix} = \,\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \,\end{bmatrix} = C(X) . \]

The easiest solution might be to just try out all four possibilities to plugin the computational basis into the first two wires.

Another approach is to just calculate the corresponding algebraic expression:

Using \(XP_0X=P_1\) and \(XP_1X=P_0\) we get, that this is:

\[ P_0 \otimes P_0 \otimes I + P_1 \otimes P_0 \otimes I + P_0 \otimes P_1 \otimes I + P_1 \otimes P_1 \otimes V^2 . \]

Which is indeed \(C^2(U)\).

If we combine the construction of \(C^2(U)\) from figure 4.8 (the \(U=V^2\) thing) with the construction of \(C(U)\) from figure 4.6 (the \(AXBXC=U\) thing), we immediately see that we can construct a circuit with six CX and ten one-qubit gates (two of the three phase gates cancel each other). However, we can reduce the number of one qubit gates by two, since two pairs can actually be merged into a single one-qubit operator.

For the y-rotation this is easy:

\[ C^1(R_y(\theta)) = C(X) \cdot R_{y,2}(-\theta/2) \cdot C(X) \cdot R_{y,2}(\theta/2) . \]

q_0: ──────────────■────────────────■──
     ┌──────────┐┌─┴─┐┌──────────┐┌─┴─┐
q_1: ┤ Ry(+π/2) ├┤ X ├┤ Ry(-π/2) ├┤ X ├
     └──────────┘└───┘└──────────┘└───┘

We get the same result if we use the proof of Corollary 4.2. For the x-rotation we first deduce from exercise 4.6 (interpretation as 3D rotations):

\[ R_x(\theta) = R_z(-\pi/2) \cdot R_y(\theta) \cdot R_z(\pi/2) . \]

Now the proof of Corollary 4.2 shows:

\[ C^1(R_x(\theta)) = R_{z,2}(-\pi/2) C(X) \cdot R_{y,2}(-\theta/2) \cdot C(X) \cdot R_{z,2}(-\pi/2) R_{y,2}(\theta/2) . \]

q_0: ─────────────────────────■────────────────■─────────────
     ┌──────────┐┌─────────┐┌─┴─┐┌──────────┐┌─┴─┐┌─────────┐
q_1: ┤ Ry(+π/2) ├┤ Rz(2Δt) ├┤ X ├┤ Ry(-π/2) ├┤ X ├┤ Rz(2Δt) ├
     └──────────┘└─────────┘└───┘└──────────┘└───┘└─────────┘

Preliminaries: Recall that \(T=e^{i\pi/8}R_z(\pi/4)\) and \(XZX=-Z\). It suffices to prove the claim for states of the form \(|ij\rangle\otimes|\psi\rangle\) (the basis states in the first two wires), by linearity.

The first wire always acts as \(T\). This does nothing to \(|0\rangle\) and applies a phase shift to \(|1\rangle\).

Now we consider the first two wires (the upper wires) states of the form \(|01\rangle\) (that is, the control qubit is turned off). Then, the second wire acts as \(ST^\dagger T^\dagger=I\). Hence \(|01\rangle\) is mapped to itself. Now consider \(|1i\rangle\). The second wire now acts as

\[ SX T\dagger XT\dagger = e^{-i\pi/4} S XR_z(\pi/4)XR_z(\pi/4) = e^{-i\pi/4} S . \]

Hence \(|1i\rangle\) is mapped to \(|1\rangle\otimes S|i\rangle\). We summarize: The first two wires act as:

\[ \,|ij\rangle \mapsto \begin{cases} i|11\rangle & \text{if } i=j=1 \\ \,|ij\rangle & \text{else.} \,\end{cases} \]

Now consider the third wire. It is not hard to see that if one of the first two qubits is \(|0\rangle\), then all operators cancel each other (equalize to \(I\)) due to \(X^2=-1\) and \(T^\dagger T=I\).

So the only interesting case is if the first to qubits are both on (equal state 1). In that case we have to show that the operator in the third row is equal to \(-iX\). We have:

\[ T X T^\dagger X = R_z(\pi/4) X R_z(-\pi/4) X = R_z(\pi/2) . \]

Hence the operator in the third row is:

\[ H R_z(\pi/2) R_z(\pi/2) H = R_x(\pi) = -iX. \]

QED.

Taking the swap-gate construction from Figure 1.7 we may deduce that the following circuit might do what we want:

q_0: ──■────■────■──
       │  ┌─┴─┐  │
q_1: ──■──┤ X ├──■──
     ┌─┴─┐└─┬─┘┌─┴─┐
q_2: ┤ X ├──■──┤ X ├
     └───┘     └───┘

In fact, a short calculation yields the desired result.

T01 = make_CnU(3, [0, 1], 2, X)
T02 = make_CnU(3, [0, 2], 1, X)

assert T01 * T02 * T01 == Fred, "Exercise 4.25, Solution to 1."
"PASSED"

PASSED

The following circuit equals the Fredkin Gate too:

q_0: ───────■───────
          ┌─┴─┐
q_1: ──■──┤ X ├──■──
     ┌─┴─┐└─┬─┘┌─┴─┐
q_2: ┤ X ├──■──┤ X ├
     └───┘     └───┘

In fact:

C12 = make_CU(3, 1, 2, X)
T02 = make_CnU(3, [0, 2], 1, X)

assert C12 * T02 * C12 == Fred, "Exercise 4.25, Solution to 2."
"PASSED"

PASSED

Apriori this new circuit uses 8 = 6 + 2 two-qubit gates. However the first two of these can be merged into one (together with the Hadamard in between). Observe that the a similar fact is true for the last two two-qubit gates of the construction from Figure 4.9. Hence we only need 6 such gates.

So far I have no solution. My first guess was to use the construction from Firgure 4.8 (the thing with \(V^2=X\)). With the same trick as in Solution 3 (merging adjacent gates) one gets a solution with 6 two-qubit gates. I tried to move around gates at the right end (via commutator relations) hoping the find an equivalent circuit with on gate less - so far without success.

The following code saves us from the tedious calculations. It shows that \(?=4\) (everywhere) is a possible choice and that

\[ \theta(c_1, c_2, t) = \begin{cases} \pi & \text{if } c_1=1, \text{ and } c_2=t=0 \\ 0 & \text{else.} \,\end{cases} \]

def get_circuit(n: int):
    """The n is the divisor missing in the Figure of Exercise 4.26 (the '?')."""
    R = tprod(Id, Id, Ry.subs(theta, pi / n))
    Ri = tprod(Id, Id, Ry.subs(theta, -pi / n))

    CX12 = make_CU(3, 1, 2, X)
    CX02 = make_CU(3, 0, 2, X)

    return sp.simplify(R * CX12 * R * CX02 * Ri * CX12 * Ri)


C = get_circuit(4)
Phase = tprod(Id, Id, Id) - 2*tprod(P1, P0, P0)  # Multiplies |100> by -1

assert Phase * C == Toff, "Exercise 4.26: Circuit should be Toffoli up to Phase Factor"
"PASSED"

PASSED

Since the exercise already speaks about permutations lets dive into the permutation group formalism. This works since all involved matrices (CNOT, Toffoli, and the matrix from the exercise text) are all permutation matrices. These form a group which is isomorphic to a group of permuations. The procedure below also works analogously if we would directly use the matrix representations but I think it is nice to think it throught in a different formalism.

Keep in mind that the isomorphism reverses the order of multiplication, that is: if Matrices \(A\) and \(B\) correspond to permutations \(a\) and \(b\) then \(AB\) corresponds to \(ba\) (meaning apply \(b\) first and then \(a\)). This is not just a sympy quirk but a normal convention in the mathematics of permutations. Maybe the rationale is that the action of a permutation is sometimes written in superscript notation: \(x^{ba}=(x^b)^a\) where \(x\) is some set which gets permutated (see also right-group-action on wikipedia, which does not use superscript notation though).

Recall the Matrix representation for the CNOT gate with control at 0 and target at 1 (on three wires):

\[ C_{01}(X) = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \,\end{bmatrix} \]

If we enumerate the basis vectors from \(0\) (for \(|000\rangle\)) to \(7\) (for \(|111\rangle\)), then this corresponds to the permuation \((4\;6)(5\;7)\) (in cycle-notation). The other variations of CX (with different control/target) can be found by noting that permuting the three wires induces a permutation on the basis vectors which in turn translates to an action on the permutation group itself. See the code block below for the results (or use make_CU to make life easier).

The "standard" Toffoli Gate is the one with target at 2:

\[ T_2 = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \,\end{bmatrix} , \]

which corresponds to the permuation \((6\;7)\). The other three Toffolis (with different targets) are listed in the code block below.

# cx3ij means: control i, target j
cx301 = Perm(4, 6)(5, 7)
cx302 = Perm(4, 5)(6, 7)
cx312 = Perm(2, 3)(6, 7)
cx310 = Perm(2, 6)(3, 7)
cx320 = Perm(1, 5)(3, 7)
cx321 = Perm(1, 3)(5, 7)

# t3i means: target at qubit i
t30 = Perm(3, 7)
t31 = Perm(5, 7)
t32 = Perm(6, 7)

The following code defines an exhaustive search for solutions with a fixed number of gates.

def search_ex4_27(n: int):
    """Find all solutions to exercise 4.27 with n gates."""
    perms = [cx301, cx302, cx312, cx310, cx320, cx321, t30, t31, t32]
    names = ["cx301", "cx302", "cx312", "cx310", "cx320", "cx321", "t30", "t31", "t32"]
    indices = list(range(len(perms)))

    for idxs in product(*([indices]*n)):
        ps = [perms[i] for i in idxs]
        prod = reduce((lambda x, y: x*y), ps)
        if prod == Perm(1, 2, 3, 4, 5, 6, 7):
            print(" * ".join([names[i] for i in idxs]))

Trying out some values of \(n\) we see that \(n=5\) is the minimal value which yields a solution. There are six solutions in fact but they are similar in the sense that the last three permutations (gates) actually commute. The solutions are:

search_ex4_27(5)

t30 * cx321 * cx302 * cx312 * t32
t30 * cx321 * cx302 * t32 * cx312
t30 * cx321 * cx312 * cx302 * t32
t30 * cx321 * cx312 * t32 * cx302
t30 * cx321 * t32 * cx302 * cx312
t30 * cx321 * t32 * cx312 * cx302

One benefit of the "reversed" multiplication of the permutations is that it is the same ordering like in typical drawings of circuits. Let us just draw the third one:

The idea for the proof is surprisingly simple - watch out for the determinant of the circuit! The credits for the idea go to this stackoverflow answer. The overall discussion contains also a reference to the pseudo-solution mentioned in the remark above.

Consider \(n+1\) wires and recall that with all controls being on the first \(n\) wires we have:

\[ C^n(U) = \mathrm{diag}(I,I,\ldots,I,U) \text{ with } 2^n - 1 \text{ copies of } I=\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}. \]

That is, \(C^n(U)\) is a diagonal matrix up to a tiny 2 by 2 sub-matrix \(U\) in the lower right corner. This follows from the ordering of the basis vectors (counting in binary, like 000, 001, 010, …, 111). Hence:

\[ \det(C^n(U)) = \det(U). \]

Note that this formula also holds if the controls are elsewhere. This follows from the fact that the determinant is invariant with respect to a change of basis.

Now consider the more general case of \(C^m(W)\) (for a unitary \(W\)) with \(m\leq n\) on \(n+1\) wires. Assume for the moment that \(W\) sits on the last wire and the controls being immediately before. Again by the ordering of the basis vectors we see that:

\[ C^m(W) = \mathrm{diag}(W_m, W_m,\ldots, W_m) \text{ with } 2^{n-m} \text{ copies of } W_m, \]

where \(W_m\) looks like in the special case \(m=n\) above. Since the determinant of \(W_m\) is \(\det(W)\), we deduce:

\[ \det(C^m(W)) = \det(W)^{2^{n-m}} . \]

Again this still holds if the controls and the target are elsewhere on the \(n\) wires.

Now suppose the construction was possible:

\[ C^n(U) = \prod_i G_i, \]

where each \(G_i\) is one of the allowed gates. By well known properties of the determinant we have:

\[ \det(U) = \prod_i \det(G_i) . \]

The analysis above, together with \(n-m\geq 2\), shows that \(\det(G_i)=\det(W)^{4k}\) for some \(k\) and \(W\) one of \(V\), \(V^\dagger\), \(U\) or \(X\). Since \(\det(X)=-1\) we have:

\[ \det(U) = \left(\det(V)^{k_V} \cdot \det(V^\dagger)^{k_{V^\dagger}} \cdot \det(U)^{k_U}\right)^4. \]

This can't be satisfied for arbitrary \(U\). To be specific consider \(U=X\). In that case the LHS is \(-1\) and the RHS is \(+1\) because \(\det(V),\det(V^\dagger)\in\{+i,-i\}\). QED.

Preparations: Let \(N\) stand for X, Y, or Z. Recall that \(P_0\) and \(P_1\) are the projections onto the computational basis, and that:

In particular: \(X_1 X_2=X\otimes X\).

Identity (1) follows by considering only states of the form \(|0j\rangle\), \(|1j\rangle\), noting that \(X_1\) swaps the computational basis.

Identitiy (2) follows similarly if one observes that \(Y_1\) acts like \(X_1\) up to different phases, which do not play a role here.

Identity (3) follows from the fact that \(Z_1\) commutes with \(P_0\) and \(P_1\) (and hence with \(C\)), and that \(C^2=I\).

Identity (4) follows similarly from the fact that \(X_2\) also commutes with \(C\).

Identity (5) follows from \(YX=-iZ=-XY\), and \(P_0P_1=P_1P_0=0\), and \(X^2=I\), and \(Z=P_0-P_1\).

Identity (6) follows the same way as (5).

Identity (7) follows from the fact that Z, and hence \(R_z\), commutes with \(P_0\) and \(P_1\).

Identitiy (8) follows from the fact X, and hence \(R_x\), commutes with X.

First of all observe that the formulation of the exercise clearly identifies \(P_i\) with \(I_1\otimes P_i\). We keep this identification in the following.

If the interpret the measurement postulate "the right way" we see that it immediately implies that a measurement with respect to \(P\) leaves the system in the post measurement state

\[ \left\{ \left(\trace{P_0\rho}, P_0 \rho P_0\right), \left(\trace{P_1\rho}, P_1 \rho P_1\right) \right\} , \]

which is an ensemble. I suppose it only "collapses" into one of its members if "some classical information flows". I suppose that "an observer who did not learn the measurement result" means that this flow of information did not happen, hence the state of the system is still in the ensemble. The corresponding density operator is:

\[ P_0 \rho P_0 + P_1 \rho P_1 = \rho' . \]

This proves the first claim. To prove the second claim we may assume wlog that \(\rho=\alpha\otimes\beta\) - by linearity of the trace. Indeed:

\[ \sum_i \ptrace{2}{P_i\,\rho\,P_i} = \sum_i \ptrace{2}{I\otimes P_i \, \alpha\otimes\beta \, I\otimes P_i} = \alpha \sum_i \trace{P_i \beta P_i} = \alpha = \ptrace{2}{\rho} , \]

which concludes the solution of the exercise. QED.

Let \(x,y\in\{0,1\}\). Observe

\[ CX \cdot H\otimes I \; |xy\rangle = \frac{1}{\sqrt{2}} (|0y\rangle + (-1)^x |1\overline{y}) = |\beta_{xy}\rangle . \]

Hence the circuit first maps \(|\beta_{xy}\rangle\) to \(|xy\rangle\) and then measures in the computational basis. This already shows that the circuit indeed "measures in the Bell Basis" in some sense.

The measurement operators can be read directly from the circuit:

\[ Q_{ij} = P_{ij} \cdot H\otimes I \cdot CX =: P_{ij} \cdot U , \]

where the \(P_{ij}\) are the standard projections (for the measurement in the computational basis after the barrier). Observe that these are no projections (\(U\) and hence the \(P_{ij}U\) map the Bell Basis to the standard basis as implied by the introductory calculations). This is in line with the observation that two times applying the circuit results in a potentially different outcome.

The corresponding POVMs are indeed the projections on the Bell Basis:

\[ Q_{ij}^\dagger Q_{ij} = U^\dagger P_{ij} U , \]

since \(U\) maps the Bell Basis to the standard basis. QED.

Observe that the projections on the eigenvalues of \(U\) are given by:

\[ Q_{\pm} = \frac{1}{2} (I \pm U) . \]

This and a short calculation shows that the result of the circuit at the position of the barrier (before measurement) is:

\[ H_1 \cdot C(U) \cdot H_1 \; |0\psi\rangle = |0\rangle \otimes Q_+ |\psi\rangle + |1\rangle \otimes Q_- |\psi\rangle . \]

Now we see that finally applying the measurement in the first qubit "reveals" the corresponding projection of \(|\psi\rangle\) in the second wire. QED.

Let us consider a generic state \(|\psi\rangle=\sum_{ij}\alpha_{ij}|ij\rangle\).

Case 1 (Measure first): Directly after measuring, the state of the system is an ensemble (of the possible measurement outcomes):

\[ \{ (p_0, |\psi_0\rangle), (p_1, |\psi_1\rangle) \} , \]

where

\[ p_i = \sqrt{|\alpha_{i0}|^2 + |\alpha_{i1}|^2}, \quad \,|\psi_i\rangle = \frac{1}{p_i} (\alpha_{i0} |i0\rangle + \alpha_{i1} |i1\rangle) . \]

The conditional application of \(U\) leaves the system in state

\[ \{ (p_0, |\psi_0\rangle), (p_1, I\otimes U \; |\psi_1\rangle) \} . \]

Case 2 (CU first): After applying the controlled \(U\) the state is:

\[ C(U) |\psi\rangle = \sum_j (\alpha_{0j} |0j\rangle + \alpha_{1j} |1\rangle \otimes U|j\rangle) . \]

Measuring this, leaves us in the same final state as in case 1:

\[ \{ (p_0, |\psi_0\rangle), (p_1, I\otimes U \; |\psi_1\rangle) \} . \]

QED.

To find a possible solution let us write the operation in a slightly more explicit bit-wise form:

\[ \, |x_0,x_1,y_0,y_1\rangle \mapsto |x_0,x_1,x_0\oplus y_0, x_1\oplus c \oplus y_1\rangle , \]

where \(c=x_0\cdot y_0\) is the carry-bit. Recall that \(C(X)\) implements binary addition modulo 2. Hence the following circuit, which is inspired by the above equation, does what we want:

q_0: ──■────■───────
       │    │
q_1: ──┼────┼────■──
       │  ┌─┴─┐  │
q_2: ──■──┤ X ├──┼──
     ┌─┴─┐└───┘┌─┴─┐
q_3: ┤ X ├─────┤ X ├
     └───┘     └───┘

We semi-automate the task with sympy. First we have to set up the matrix:

F = Matrix([[I**(x*y) for x in range(4)] for y in range(4)]) / 2

Using the the utility function make_twolevel it is relatively easy to follow the procedure from the book (note that in general we also had to use make_onelevel, but here we don't need it)):

U1 = sp.simplify(make_twolevel(4, [1,0], [1, -1]))
U2 = sp.simplify(make_twolevel(4, [2,0], [1, -sqrt(2)]))
U3 = sp.simplify(make_twolevel(4, [3,0], [1, -sqrt(3)]))
U4 = sp.simplify(make_twolevel(4, [2,1], [sqrt(3) * (3 + I), -3 * (1 - I)]))
U5 = sp.simplify(make_twolevel(4, [3,1], [I, -sqrt(2)]))

F0 = make_CU(2, 0, 1, Matrix([[1, I], [-1, I]]) / sqrt(2))
assert sp.simplify(U5*U4*U3*U2*U1*F) == F0
"PASSED"

PASSED

From this we can easily readout the decomposition:

\[ \left[\begin{matrix}\frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2} & 0 & 0\\\frac{\sqrt{2}}{2} & - \frac{\sqrt{2}}{2} & 0 & 0\\0 & 0 & 1 & 0\\0 & 0 & 0 & 1\end{matrix}\right] \cdot \left[\begin{matrix}\frac{\sqrt{6}}{3} & 0 & \frac{\sqrt{3}}{3} & 0\\0 & 1 & 0 & 0\\\frac{\sqrt{3}}{3} & 0 & - \frac{\sqrt{6}}{3} & 0\\0 & 0 & 0 & 1\end{matrix}\right] \cdot \left[\begin{matrix}\frac{\sqrt{3}}{2} & 0 & 0 & \frac{1}{2}\\0 & 1 & 0 & 0\\0 & 0 & 1 & 0\\\frac{1}{2} & 0 & 0 & - \frac{\sqrt{3}}{2}\end{matrix}\right] \cdot \left[\begin{matrix}1 & 0 & 0 & 0\\0 & - \frac{\sqrt{3} i \left(1 + i\right)}{4} & \frac{3}{4} - \frac{i}{4} & 0\\0 & \frac{3}{4} + \frac{i}{4} & - \frac{\sqrt{3} \cdot \left(1 + i\right)}{4} & 0\\0 & 0 & 0 & 1\end{matrix}\right] \cdot \left[\begin{matrix}1 & 0 & 0 & 0\\0 & \frac{\sqrt{6}}{3} & 0 & - \frac{\sqrt{3} i}{3}\\0 & 0 & 1 & 0\\0 & \frac{\sqrt{3} i}{3} & 0 & - \frac{\sqrt{6}}{3}\end{matrix}\right] \cdot \left[\begin{matrix}1 & 0 & 0 & 0\\0 & 1 & 0 & 0\\0 & 0 & \frac{\sqrt{2}}{2} & \frac{\sqrt{2} i}{2}\\0 & 0 & - \frac{\sqrt{2}}{2} & \frac{\sqrt{2} i}{2}\end{matrix}\right] \]

The following code was used to produce the above latex expression:

def get_latex_for_exercise_4_37():
    Vs = [U1.H, U2.H, U3.H, U4.H, U5.H, F0]
    return " \\cdot ".join([sp.latex(V) for V in Vs])

For each unitary matrix \(U\) let us define an equivalence relation \(\sim_U\) on \(\{0,\ldots,d-1\}\) as follows. First of all define:

\[ i \sim^0_U j :\Leftrightarrow \langle i|U|j\rangle \neq 0 \]

Now let \(\sim_U\) be the reflexive and transitive closure \(\mathrm{closure}(\sim^0_U)\) of \(\sim^0_U\). We remark that we actually only care about \(i,j\) with \(i\neq j\). But it is convenient to have an equivalence relation, hence we require reflexivity.

Let us call equivalence classes with at least two elements non-trivial. In a sense, the non-trivial equivalence classes correspond to those subspaces, spanned by base vectors, on which the matrix acts non-trivially. Note that this is a slightly stricter notion then in the book. Merely multiplying a base vector with a complex number of modulus one does not count. Instead we require some "interaction" between the basis vectors.

Note that two-Level matrices have at most one non-trivial equivalence class \(A\), and this \(A\) has exactly two elements.

Let us define the score of an equivalence relation \(\sim\) by:

\[ \mathrm{score}(\sim) = \sum_A (|A| - 1), \]

where \(A\) ranges over all (non-trivial) equivalence classes of \(\sim\). Observe that the score of a two-level matrix is at most \(1\). Let us write \(\mathrm{score}(U):=\mathrm{score}(\sim_U)\). The score has the following easy to prove properties (we only sketch their proofs)

monotonicity

\(\sim_1\;\subset\;\sim_2\quad\Rightarrow\quad\mathrm{score}(\sim_1)\leq\mathrm{score}(\sim_2)\)

sub-additivity

\(\mathrm{score}(\mathrm{closure}(\sim_1\cup\sim_2))\leq\mathrm{score}(\sim_1)+\mathrm{score}(\sim_2)\).

The key to prove monotonicity is to observe that the equivalence classes of \(\sim_1\) are subsets of equivalence classes of \(\sim_2\). Note that the claim would even be true without the \(-1\) "punishment" for each equivalence class.

For sub-additivity first observe that the equivalence classes of the combined relation is "essentially" the set of the "old" equivalence classes. In the extreme case where this is "realy" true equality holds. In general, there can be old equivalence classes which overlap. This means that we have to merge them. This "merger" of two old equivalence classes decreases the score by \(m-1\), where \(m\) is the number of common elements of the two old equivalence classes. Of course the merging has to be repeated until the equivalence classes are disjoint.

To say something about \(\sim_{UV}\) recall:

\[ \langle i|UV|j\rangle = \sum_k \langle i|U|k\rangle\langle k|V|j\rangle . \]

This implies that

\[ \sim^0_{UV} \; \subseteq \; \sim_U \circ \sim_V \; \subseteq \; \mathrm{closure}(\sim_U \cup \sim_V) \; =: \; \sim_{U,V} . \]

In particular \(\sim_{UV}\;\subseteq\;\sim_{U,V}\) (note the comma between \(U\) and \(V\) in the RHS). Hence, by monotonicity and sub-additivity:

\[ \mathrm{score}(UV) \leq \mathrm{score}(\sim_{U,V}) \leq \mathrm{score}(U) + \mathrm{score}(V) . \]

In particular the score of a matrix which can be decomposed into \(k\) two-level matrices is at most \(k\). On the other hand there are matrices whose score is \(d-1\) (the maximal possible value). Take for example the cyclic permuation which maps \(|i\rangle\) to \(|i+1\mod d\rangle\). Such a matrix cannot be expressend by fewer than \(d-1\) matrices. QED.

Let us call the matrix \(M\). Clearly \(M\) represents a two-level operation at \(s=010\), \(t=111\). The following gray code connects the two bit-patterns:

Let \(P\) be the permutation matrix corresponding to the permuation \((7)(2,3)\) in cycle notation (it swaps \(|g_1\rangle\) and \(|g_2\rangle\)). It can be implemented by a controlled X. \(PMP\) is a controlled \(\tilde{U}\) gate with target \(0\) and controls at all remaining wires. Hence the following circuit implements \(M\):

          ┌───┐
q_0: ──o──┤ Ũ ├──o──
       │  └─┬─┘  │
q_1: ──■────■────■──
     ┌─┴─┐  │  ┌─┴─┐
q_2: ┤ X ├──■──┤ X ├
     └───┘     └───┘

Figure 4.9 shows how to implement the controlled X from fundamental gates (take also into account the "trick" from Figure 4.11). Use Figure 4.8 and Figure 4.6 to implement the controlled \(\tilde{U}\).

This directly follows from the fact that the spectral radius of a normal operator is equal to its operator norm (for non-normal operators it might be strictly smaller). For diagonal matrices this is easy to see. For the general case, recall that normal operators are those operators which can be diagonalized by unitary operators.

Let us return to the particular exercise. By the spectral theorem, and the fact that the rotation operators are unitary we have:

\[ \norm{R_{\hat{n}}(\alpha) - R_{\hat{n}}(\alpha + \beta)} = \norm{ R_{\hat{n}}(\alpha) (I-R_{\hat{n}}(\beta)) } = \norm{I-R_{\hat{n}}(\beta)} . \]

Again by the spectral theorem the eigenvalues of \(I-R_{\hat{n}}(\beta)\) are \(1-\exp(\mp i\beta/2)\) (because the eigenvalues of \(\hat{n}\cdot(X,Y,Z)\) are \(\pm1\)). Hence the spectral radius is \(|1-\exp(i\beta/2)|\). Putting everything together we proved the claim.

Formula (4.76) follows directly from this and \(R_{\hat{n}}(\theta)^n=R_{\hat{n}}(n\theta)\) (again spectral theorem) for an appropriate choice of \(n\) (so that \(\alpha\) and \(n\theta\) are close enough (as elements of the unit circle)). QED.

The initial state is \(|00\psi\rangle\), however let us consider the more general situation where the initial state is \(|ab\psi\rangle\) for \(a,b\in\{0,1\}\). Observe that

\[ H\otimes H |ab\rangle = \frac{1}{2} \sum_{ij} (-1)^{ai+bj} |ij\rangle . \]

It is easy to calculate the state directly after the second CCX. It is:

\[ \, |\zeta_0\rangle = \frac{1}{2} \sum_{ij} \, (-1)^{ai+bj} \, |ij\rangle \, X^{ij} S X^{ij} \, |\psi\rangle . \]

Before applying the last two Hadamard gates it comes in handy to represent this in a slightly different way:

\[ \, |\zeta_0\rangle = \frac{1}{2} \sum_{ij} \, (-1)^{ai+bj} \, |ij\rangle \, S \, |\psi\rangle \, + \, \frac{1}{2} \, (-1)^{a+b} \, |11\rangle \, (XSX - S) \, |\psi\rangle . \]

Applying the last two Hadamards leads to

Note that the thus defined linear operator \(A\) is not unitary but invertible. It has the nice property that \(\norm{A|ab\rangle}=\sqrt{3}/4\) (independent of \(a,b\)) and that \(A|ab\rangle\) is orthogonal to \(|ab\rangle\). On the other hand these properties do not generalize to arbitrary states.

Clearly

\[ XSX - S = (i-1) Z . \]

The other operator which is applied to the first occurrence of \(|\psi\rangle\) is a bit more tricky to make sense of. We make use of the functional calculus for Z (more specifically: we use the Pauli-Rotation formalism). First of all recall that \(S=e^{i\pi/4}R_z(\pi/2)\) and \(R_z(-\pi)=iZ\). Using this, with \(XZX=-Z\) we get:

Observe that \(\sin\theta_0=2\sin(\theta_0/2)\cos(\theta_0/2)=-3/5\). Let \(\theta=\theta_0+\pi/2\). Hence \(\cos(\theta)=3/5\) and

\[ \frac{3}{4}S+\frac{1}{4}XSX = \frac{\sqrt{10}}{4} e^{i\pi/4} R_z(\theta) . \]

Let us summarize what we found out. The final state of the qubits is:

\[ \, |\zeta\rangle = |ab\rangle \, \frac{\sqrt{10}}{4} e^{i\pi/4} R_z(\theta) |\psi\rangle \,+\, A|ab\rangle (i-1) Z |\psi\rangle \]

Conclusion: If the initial state is \(|ab\rangle|\psi\rangle\) then we measure \(|ab\rangle\) with probability \(5/8\). In that case the final state of the "payload" is \(R_z(\theta)|\psi\rangle\) (up to a global phase). Otherwise the payload ends up in state \(Z|\psi\rangle\) (up to a global phase). The angle satisfies \(\cos\theta=3/5\). This proves the first part of the exercise.

Let us sketch a possibility to improve the probability of success. One can introduce two more ancillary qubits (initialized to \(|0\rangle\)) and put another such circuit on them (connected to the "main" qubit as the first one). To cancel the Z from the first circuit, in case it "fails", one can use three CCZ gates (for the three "failing" measurement outcomes). To prevent a double rotation (if both circuits succeed) one can put a controls on the \(S\) gate (or simpler: put a controlled \(S^3\) next to it) of the second circuit making it dependent on success of the first one. My the principle of deferred measurement all measurements can be done at the end of the circuit. To improve the probability even further, just repeat the procedure.

The exercise essentially contains the solution. If \(\theta\) was rational, then there would be integers \(k,m\) such that \(\theta=2k\pi/m\). Hence

\[ (3 + 4i)^m = 5^m \cdot e^{2k\pi i} = 5^m . \]

On the other hand

\[ (3+4i)^2 = -7 + 24i \equiv 3 + 4i \mod 5 . \]

Clearly the same is true for higher powers. Hence \((3+4i)^m\) cannot be real - which contradicts a previous statement. Hence \(\theta\) can not be rational. QED.

We already know that if we add the \(T\) gate then we have a universal set of gates. Hence it suffices to be able to approximate the \(T\) gate. Note that the \(T\) gate is \(R_z(\pi/4)\) up to a global phase. Exercises 4.41 and 4.42 show that any \(R_z(\theta)\) can be approximated arbitrarily well. QED.

Of course, contrary to the construction from the main-text, this construction with the Toffoli gate needs some ancillary qubits. So we proved a weaker kind of universality then what is proved for the "original" set of universal gates!

We show universality for the situation where one has access to a reservoir of ancillary qubits in state \(|0\rangle\) or \(|1\rangle\).

By irrationality of \(\alpha\) one can approximate \(C^2(R_x(\theta))\) for any angle \(\theta\) (we drop the factor \(i\) because we ignore the global phase here and in the following). By the use of ancillary qubits in state \(|1\rangle\) one can thus implement \(R_x(\theta)\) for any \(\theta\).

Moreover the Toffoli CCX can be implemented. By the use of ancillary qubits in state \(|1\rangle\) it is easy to implement CX and X from it.

To proof universality it suffices to implement \(R_z(\theta)\) for arbitrary \(\theta\). To see this, first recall that according to exercise 4.11 any single qubit gate has the form \(R_z(\beta)R_x(\gamma)R_z(\delta)\) (up to a global phase). On the other hand, the book shows that single qubit gates together with CX is a universal set.

Let \(|\pm\rangle\) be eigenstates of X corresponding to the eigenvalues \(\pm1\). By the spectral theorem we have:

\[ C(R_x(\theta)) |k-\rangle = \,\begin{cases} \, |0-\rangle & \text{for } k=0 \\ e^{i\theta/2} \, |1-\rangle & \text{for } k=1 \,\end{cases} = (e^{i\theta/4} R_z(\theta/2) |k\rangle) \otimes |-\rangle . \]

That is, if we had access to an ancillary qubit in state \(|-\rangle\) we could implement all \(R_z(\theta)\) gates. To see that, in a sense, we have access to such a qubit, consider:

\[ C = R_x(\pi/2)_1 \cdot C(X) \cdot R_x(\pi/2)_1 . \]

As usual the index \(1\) means that we apply the rotation gate to the first qubit. Observe that

\[ C |00\rangle = \frac{1}{\sqrt{2}} (|0\rangle|-\rangle - i|1\rangle|+\rangle) . \]

That is, the circuit \(C\), feeded with the initial state \(|00\rangle\) and measured at the first qubit produces \(|-\rangle\) with a probability of \(1/2\) (in that case \(0\) is the measurement outcome). One can make the probability as close to one as needed by running sufficiently many of these circuits, but with the same "output" qubit. Of course, we have to prevent "multiple successes". To accomplish this, we add additional controls to the \(C(X)\) gate making it dependent on the measurement outcomes of the "previous" runs of \(C\). Of course measurements can be deferred to the end.

For the argument above it is of course sufficient to be able to generate \(|-\rangle\) with merely high probability. This proves the claim. QED.

Observe that except for the Hadamard gate all mentioned gates have complex integer coefficients (and Kronecker products of these gates with the identity obviously does not change this property). Moreover we have

\[ H = 2^{-1/2} \tilde{H} , \]

for a matrix \(\tilde{H}\) with (real - but this does not matter) integer coefficients. Hence we can choose \(k\) to be the number of Hadamard gates and \(M\) to be the product of the matrices where \(H\) is replaced by \(\tilde{H}\).

Now consider the \(T\) gate. Observe that

\[ T = \frac{1}{\sqrt{2}} \begin{bmatrix} \sqrt{2} & 0 \\ 0 & i+1 \end{bmatrix} . \]

By a similar reasoning as above (technically by mathematical induction) we can easily show that any gate composed of the mentioned base gates has the form

\[ 2^{-k/2} (M_1 + \sqrt{2} M_2) , \]

where \(M_1\) and \(M_2\) have only complex integer coefficients and \(k\) is the number of \(T\) gates involved. QED.

We know that \(\rho\in\CC^{2n\times2n}\). Hence, a priori we need \(2\cdot4^n\) real numbers to describe it. Since \(\rho\) is (positive) hermitian, the diagonal entries are real and the off-diagonal entries above the diagonl determine the off-diagonal entries below the diagonal completely. Hence we are left with \(4^n\) reals. This again reduces by one to \(4^n-1\) by the trace condition.

First of all let us observe that it is sufficient to prove the statement for the special case \(L=2\) since the general case directly follows from it. Let us restate a "cleaned up" version of the statement:

Let \(A\), \(B\) two commuting bounded operators (on the same Hilbert Space). Then \(e^{A+B}=e^Ae^B\). We will prove this statement.

The proof has two key ingredients. The first one is the Binomial Formula

\[ (A + B)^k = \sum_{i=0}^k \binom{k}{i} A^{k-i} B^i . \]

This formula can be easily proved by mathematical induction on \(k\). The commutativity of \(A\) and \(B\) is needed in the proof and in fact it allows us to take the standard proof, where \(A\) and \(B\) are real (or complex) numbers, as a template.

The second ingredient is the Cauchy Product Formula:

\[ \left( \sum_{i=0}^\infty a_i \right) \cdot \left( \sum_{i=0}^\infty b_j \right) = \sum_{k=0}^\infty \sum_{i=0}^k a_{k-i} b_i . \]

This formula holds under rather general conditions. One sufficient condition is that the two series on the left are both absolutely convergent series (in a space of bounded operators over a Hilbert Space). In fact, no commutativity is required and even associativity is unnecessary (the results holds in abstract Banach Algebras).

Now let us come back to our problem. From the binomial formula we deduce:

\[ e^{A+B} = \sum_{k=0}^\infty \frac{1}{k!} (A+B)^k = \sum_{k=0}^\infty \frac{1}{k!} \sum_{i=0}^k \binom{k}{i} A^{k-i} B^i = \sum_{k=0}^\infty \sum_{i=0}^k \frac{1}{(k-i)!} A^{k-i} \cdot \frac{1}{i!} B^i . \]

The RHS has the form of the RHS of the Cauchy Product formula. Hence

\[ e^{A+B} = e^A e^B \]

as desired. QED.

I suppose that \(n\) is the number of qubits. We assume that any particle is described by at most \(m\) qubits. This implies that each Hamiltonian involves at most \(mc\) qubits.

Moreover we assume that there are only a finite number \(K\) of types of interactions. As an example consider a special case of the (quantum) Ising Model (aka Heisenberg Model) with periodic boundary conditions:

\[ H = -\frac{J}{2} \sum_{i=0}^n X_i X_{i+1} - h \sum_{j=1}^n Z_j . \]

Here we have \(c=2\), \(m=1\) and \(K=2\). More generally the two particle interaction could also involve \(Y_iY_{i+1}\) and \(Z_iZ_{i+1}\) contributions, in which case we had \(K=4\).

In the extreme case that each combination of \(c\) particles gets a Hamiltonian of each type we get the maximum value for \(L\). So we get the upper bound

\[ K \binom{n/m}{c} = O(Km^{-c} n^c) \]

if each particle needs exactly \(m\) qubits. In general at least some particles could need less qubits. So the save upper bound is:

\[ L \leq K \binom{n}{c} = O(K n^c) . \]

Note that not each type of Hamiltonian needs to act on all \(c\) particles. In that case the inequality is already strict, since then the given binomial coefficient overestimates the possible combinations. QED.

Let us first prove formula (1) - because it is the simplest. In fact, it directly follows from first order approximating all involved exponential series:

\[ e^{(A+B)t} = I + (A+B)t + O(t^2) = e^{At} e^{Bt} + O(t^2) . \]

Remark: From this we easily obtain the Trotter Formula (Theorem 4.3) by replacing \(t\,\leftarrow\,t/n\) and taking the \(n\)​-th power. We even get a uniform error estimate \(O(1/n)\) if \(t\) is restricted to a fixed interval \([0,t_0]\). Of course the original proof shows the same (including the error estimate), but it contains an error: It is not true that \(\binom{n}{k}/n^k=(1+O(1/n))/k!\). Set \(k=n\) to see this.

The other two formulas are also showed by comparing power series. This rather boring task can be done by sympy. Hence let us define three symbols. Two non-commuting symbols for the operators and a commuting symbol for the time:

A, B = sp.symbols('A B', commutative=False)
t = sp.symbols('t', commutative=True)  # commutes with A and B

For a warm up let us return to formula (1):

E1 = exp((A+B)*t) - exp(A*t)*exp(B*t)
E1.series(t, 0, 3)

t**2*(-A*B/2 + B*A/2) + O(t**3)

This of course shows the claim again. But it is slightly more precise. The \(t^2\) term as the coefficient \(-[A,B]/2\). In a sense this already proves the BCH formula (since the exponential's leading term is \(I\)). But let us check it anyway. We increase the order of the series expansion to \(4\) since we are curious what the first non-vanishing term looks like:

E2 = exp((A+B)*t) - exp(A*t)*exp(B*t)*exp(-(A*B-B*A)*t**2/2)
E2.series(t, 0, 4)

t**3*(-A*B*A/3 - A*B**2/3 + A**2*B/6 + 2*B*A*B/3 + B*A**2/6 - B**2*A/3) + O(t**4)

Now let us check formula (2):

E3 = exp((A+B)*t) - exp(B*t/2)*exp(A*t)*exp(B*t/2)
E3.series(t, 0, 4)

t**3*(A*B*A/6 + A*B**2/24 - A**2*B/12 - B*A*B/12 - B*A**2/12 + B**2*A/24) + O(t**4)

Comparing this with the third order term of the BCH formula it looks like this formula has a slightly better error estimate (the coefficients are smaller). But this is speculation at this point!

QED.

The first part is a direct consequence of one of the formulas proved in Exercise 4.49:

\[ e^{(A + B)t} = e^{Bt/2} e^{At} e^{Bt/2} + O(t^3) . \]

The rest of the proof is just mathematical induction on \(L\). The base case \(L=0\) is trivial (holds even with a zero-error). The induction step uses the above formula.

The case \(m=1\) of the second part is essentially a reformulation of the first part. In fact the first part shows that there is a \(t_0>0\) such that the inequality holds for all \(t\in[-t_0,t_0]\). But for \(t>t_0\) we can use that the LHS is bounded by \(2\). Hence the inequality also holds in that range after possibly adjusting \(\alpha\). The general case follows from the case \(m=1\) because all involved operators are unitary and hence the error estimates just add up according to Box 4.1. QED.

Let

\[ H = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix}, \quad K = \frac{1}{\sqrt{2}} \begin{bmatrix} i & -1 \\ -1 & i \end{bmatrix}. \]

The first matrix is the Hadamard gate. Recall \(HXH=Z\). The second matrix is chosen so that \(K^\dagger YK=Z\). To construct it we just calculated the eigenvectors of Y and used them as columns in \(K\). Hence

\[ \mathcal{H} = H\otimes K^\dagger\otimes I \cdot Z_1\otimes Z_2\otimes Z_3 \cdot H\otimes K\otimes I . \]

Observe that \(K=iR_x(\pi/2)\). Note that the factor \(i\) in the rotation operator representation of \(K\) and \(K^\dagger\) cancels out when diagonalizing: \(R_x(\pi/2)YR_x(-\pi/2)=Z\). Hence

\[ \mathcal{H} = H\otimes R_x(\pi/2)\otimes I \cdot Z_1\otimes Z_2\otimes Z_3 \cdot H\otimes R_x(-\pi/2)\otimes I , \]

and thus:

\[ e^{-i\mathcal{H} \Delta t} = H\otimes R_x(\pi/2)\otimes I \cdot e^{-i Z_1\otimes Z_2\otimes Z_3 \Delta t} \cdot H\otimes R_x(-\pi/2)\otimes I . \]

From the book we already know how to construct the operator in the middle (Figure 4.19). The complete circuit is:

                       ░    ┌───┐                                                ┌───┐
x_0: ──────────────────░────┤ H ├──────■───────────────────────────────────■─────┤ H ├────
                       ░ ┌──┴───┴───┐  │                                   │  ┌──┴───┴───┐
x_1: ──────────────────░─┤ Rx(+π/2) ├──┼────■─────────────────────────■────┼──┤ Rx(-π/2) ├
                       ░ └──────────┘  │    │                         │    │  └──────────┘
x_2: ──────────────────░───────────────┼────┼────■───────────────■────┼────┼──────────────
     ┌───────────────┐ ░             ┌─┴─┐┌─┴─┐┌─┴─┐┌─────────┐┌─┴─┐┌─┴─┐┌─┴─┐
anc: ┤ Initialize(0) ├─░─────────────┤ X ├┤ X ├┤ X ├┤ Rz(2Δt) ├┤ X ├┤ X ├┤ X ├────────────
     └───────────────┘ ░             └───┘└───┘└───┘└─────────┘└───┘└───┘└───┘