<jupyter_start><jupyter_text>Notebook 4: Quantum operations and distanceIn this notebook we will be taking a closer look at quantum operations, i.e. parts of a quantum circuit that are _not necessarily_ unitary.<jupyter_code>import numpy as np

# Import cirq, install it if it's not installed.
try:
    import cirq
except ImportError:
    print("installing cirq...")
    !pip install --quiet cirq
    print("installed cirq.")
    import cirq<jupyter_output><empty_output><jupyter_text>Working with density matricesTo work with quantum operations we need work with density matrices instead of pure states as we have been before. Let's first see how we can do simulations with a density matrix in a unitary quantum circuit using `DensityMatrixSimulator`.<jupyter_code>circuit = cirq.Circuit()
num_qubits = 2
qubits = cirq.LineQubit.range(num_qubits)

circuit.append([cirq.H(qubits[0])])
circuit.append([cirq.CNOT(qubits[0], qubits[1])])
print(circuit)

simulator = cirq.DensityMatrixSimulator()
result = simulator.simulate(circuit)
rho = result.final_density_matrix
rho<jupyter_output><empty_output><jupyter_text>The resulting density matrix is a hermitian positive semi-definite (PSD) matrix with trace equal to 1. Because the input is a pure state $|000\rangle$, and all the operations are unitary, the output should also be a pure state. Recall that a state $\rho$ is pure if and only if $\mathrm{tr}(\rho^2)=1$. Let's verify these properties.<jupyter_code>print("Trace =", np.real_if_close(np.trace(rho)))
print("Hermitian?: ", np.allclose(rho, rho.conjugate().T))
print("PSD?: ", np.all(np.linalg.eigvalsh(rho >= 0)))
# np.linalg.eigvals computes eigenvalues of a matrix
# np.linalg.eigvalsh computes eigenvalues of a hermitian matrix. The assumption that
# the matrix is hermitian allows for a faster more numerically stable computation.

print("Pure state?", np.trace(rho @ rho) > 1 - 1e-4)<jupyter_output><empty_output><jupyter_text>Noisy channelsUsually we need the quantum operator formalism because we want to model _noise_ in a quantum circuit. There are many types of noise that can occur in real-life quantum circuits. Perhaps the simplest type is the _bit-flip channel_, which flips the state (applies the X gate) of a single qubit with a certain probability.If $q$ is the probability of flipping the state, then the bit-flip channel acts as:$$    \rho \mapsto (1-q)\rho+qX\rho X^\dagger$$in the _operator sum formalism_, the _operation elements_ are therefore $\sqrt{1-q}I$ and $\sqrt{q}X$.Let's modify the circuit above to use the bit-flip channel on the first qubit, before the CNOT gate.<jupyter_code>circuit = cirq.Circuit()
num_qubits = 2
qubits = cirq.LineQubit.range(num_qubits)

circuit.append([cirq.H(qubits[0]), cirq.bit_flip(0.1)(qubits[1])])
circuit.append([cirq.CNOT(qubits[0], qubits[1])])
print(circuit)

simulator = cirq.DensityMatrixSimulator()
result = simulator.simulate(circuit)
rho = result.final_density_matrix
rho<jupyter_output><empty_output><jupyter_text>Unlike before, this circuit is _not_ unitary. And hence the output is not a pure state. Here's what happens if we compute $\mathrm{tr}(\rho^2)$:<jupyter_code>np.trace(rho @ rho)<jupyter_output><empty_output><jupyter_text>Exercise 1a> While $\rho$ is not a pure state, it is the mixture of two pure states: $\rho = 0.1|\psi\rangle\langle \psi| + 0.9|\varphi\rangle\langle \varphi|$. Use the eigenvalue decomposition `np.linalg.eigh` to find $|\psi\rangle $ and $|\varphi \rangle$. (Hint: look carefully at the eigenvalues `eigvals` to select the right eigenvectors). Use `cirq.qis.dirac_notation` to neatly format the resulting vector.<jupyter_code>eigvals, eigvects = np.linalg.eigh(rho)

# YOUR CODE HERE<jupyter_output><empty_output><jupyter_text>Exercise 1b> Denote the entire circuit above by $\mathcal E$, then we defined $\rho =\mathcal E(|0\rangle\langle 0|)$, and observed that $\mathrm{tr}(\rho^2)<1$. What happens if we iterate the circuit a few times? Use a for loop to show experimentally that $\mathrm{tr}(\mathcal E^n(\rho))$ converges to 0.5.- To apply the circuit multiple times, we can use the `initial_state=rho` keyword for the function `simulator.simulate`. This sets the initial state of the simulator to the density matrix `rho`. - If you use too many iterations, you might get this error:```py    ValueError: The density matrix is not hermitian.```This is because of accumulating numerical errors. To avoid this, simply use fewer iterations. The convergence should be pretty good after 10 iterations.<jupyter_code>rho = simulator.simulate(circuit).final_density_matrix

# Your code here<jupyter_output><empty_output><jupyter_text>Exercise 1c> We can get sates $\rho$ such that $\mathrm{tr}(\rho^2)$ is even smaller than 0.5. Modify the circuit by adding _a single_ `bit_flip(0.1)` gate to the circuit at the right place and repeat the experiment of Exercise 1a to converge to a state with  $\mathrm{tr}(\rho^2)\to 0.25$<jupyter_code>circuit = cirq.Circuit()
num_qubits = 2
qubits = cirq.LineQubit.range(num_qubits)

# YOUR CODE HERE<jupyter_output><empty_output><jupyter_text>Exercise 1d> The lowest value of $\mathrm{tr}(\rho^2)$ we can possibly achieve is when $\rho = I/d$, where d is the dimension of the system. Show that this state is a fixed point of the circuit $\mathcal E$; i.e. $\mathcal E(\rho) = \rho$.<jupyter_code>rho_worst = np.eye(4, dtype=np.complex64) / 4

# YOUR CODE HERE<jupyter_output><empty_output><jupyter_text>Trace distance and fidelityWe will investigate how different types of noise can affect the fidelity and trace distance between states. Your first job is to implement trace distance and fidelity. Exercise 2a> Recall that the trace distance is defined by $D(\rho,\sigma) = \mathrm{tr}|\rho-\sigma|$. Implement the trace distance in a function `trace_distance`. Here you can use the fact that for any hermitian matrix $A$ we have $\mathrm{tr}|A|=\sum_i \sigma_i(A)$ where $\sigma_i$ is the $i\!$ th _singular value_ of $A$. You can compute singular values using `scipy.linalg.svdvals`<jupyter_code>import scipy.linalg


def trace_distance(rho, sigma):
    # YOUR CODE HERE
    ...


# rho = |00><00|
rho = np.zeros((4, 4), dtype=np.complex64)
rho[0, 0] = 1

# sigma = E(rho)
sigma = simulator.simulate(circuit, initial_state=rho).final_density_matrix

print(trace_distance(rho, rho))  # should be 0
print(trace_distance(rho, sigma))  # should be around 1.33<jupyter_output><empty_output><jupyter_text>Exercise 2b> The fidelity is defined by $F(\rho,\sigma) = \mathrm{tr}\sqrt{\rho^{1/2}\sigma \rho^{1/2}}$. Implement the fidelity in a function `fidelity`. You can compute the square root of a matrix using `scipy.linalg.sqrtm`. Make sure that the result is a real number, possibly by using `np.abs( ... )` on the result.<jupyter_code>def fidelity(rho, sigma):
    # YOUR CODE HERE
    ...


print(fidelity(rho, rho))  # should be 1
print(fidelity(rho, sigma))  # should be around 0.67<jupyter_output><empty_output><jupyter_text>It is perhaps worth noting that while the definition  $F(\rho,\sigma) = \mathrm{tr}\sqrt{\rho^{1/2}\sigma \rho^{1/2}}$ is used in Nielsen-Chuang, the definition  $F(\rho,\sigma) = (\mathrm{tr}\sqrt{\rho^{1/2}\sigma \rho^{1/2}})^2$ is more common in contemporary literature. Since we follow the book, we will keep using the former definition.  Fidelity of a quantum operationWe often need to know how much a quantum operation (in particular noise) can distort a state. We can do this by computing the fidelity between a state and the result of applying the operation to the state. That is we consider $F(\rho,\mathcal E(\rho))$.This fidelity is going to be bigger for some states, and smaller for others. Therefore we are for example interested in the _minimum_ fidelity obtained among all states $\rho$. This is however not easy to compute in general. Instead we consider the _average_ fidelity$$\overline F(\mathcal E) := \int_{S^{n-1}}\! F(|\psi\rangle\langle \psi|,\mathcal E(|\psi\rangle\langle \psi|))\,\mathrm d\psi$$Here we took the average over all the pure states, but instead we could also take the average over all the mixed states. For now we can estimate this integral using Monte-Carlo integration. That is, we randomly sample over states $\rho$, compute $F(|\psi\rangle\langle \psi|,\mathcal E(|\psi\rangle\langle \psi|))$, and then average the result. Exercise 2c> Implement the function `average_fidelity(circuit, N)` that estimates the average fidelity of a quantum circuit `circuit` using `N` samples. You can use the function `random_pure_state(num_qubits)` to generate random pure states.<jupyter_code>circuit = cirq.Circuit()
qubit1 = cirq.LineQubit(0)
circuit.append(cirq.bit_flip(0.1)(qubit1))


def random_pure_state(num_qubits):
    n = 2**num_qubits

    # Vector of random normal complex numbers
    psi = np.random.normal(size=n) + 1j * np.random.normal(size=n)

    # Normalize
    psi = psi / np.linalg.norm(psi)
    psi = psi.astype(np.complex64)

    # Compute rank-1 matrix |psi><psi|
    state = np.outer(psi, psi.conj())

    return state


def average_fidelity(circuit, N):
    num_qubits = len(circuit.all_qubits())

    # YOUR CODE HERE

average_fidelity(circuit, 200)  # Should be around 0.967<jupyter_output><empty_output><jupyter_text>Now let's try to understand how the average fidelity changes if we use noise of different strength. Instead of the bit-flip channel we will be considering the _depolarizing channel_, which is defined as:$$\mathcal E(\rho) = \frac{pI}{2} + (1-p)\rho$$i.e. with probability $p$ we change the state into $I/2$ -- a completely random state, and with probability $1-p$ we leave the state unchanged. This channel can be implemented using the `depolarize(p)` gate. Exercise 2d> Using a `for` loop, define different circuits consisting of the depolarizing channel of strength $p$ on a single qubit. Then compute the average fidelity of this circuit, and add the result to the list `fidelities_list`. The result is then plotted for you. If done correctly, the two plotted lines should perfectly overlap.<jupyter_code>import matplotlib.pyplot as plt

fidelities_list = []
p_values = np.linspace(0, 1, 20)

for p in p_values:
    # YOUR CODE HERE
    ...

plt.plot(p_values, fidelities_list, "-o", label="Estimated")
plt.plot(p_values, np.sqrt(1 - 2 * p_values / 3), label="Theoretical")
plt.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10854058222479275, 'fraction_numerical': 0.011760169654906496, 'mean_word_length': 4.925185608223872, 'pattern_counts': {'":': 0, '<': 56, '<?xml version=': 0, '>': 65, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '5813', 'license': 'MIT', 'n_tokens_mistral': 3393, 'n_tokens_neox': 3144, 'n_words': 1202, 'path': 'notebook4.ipynb', 'repo_name': 'RikVoorhaar/seminar-applied-math-2022', 'size': '10374'}