Chapter 3: The Mathematics of Quantum Computing: State Vectors, Bra-Ket Notation, Unitary Operators, and the Linear Algebra You Need

Learning Objectives

By the end of this chapter, you will be able to:

• Manipulate quantum states using bra-ket notation fluently
• Compute inner products, outer products, and tensor products by hand and in NumPy
• Understand the properties of unitary and Hermitian operators and why they matter for quantum mechanics
• Diagonalize operators via spectral decomposition
• Represent the Pauli matrices and identity as the fundamental basis of single-qubit operations
• Recognize that quantum computing is, at its core, applied linear algebra over the complex numbers
• Compute measurement probabilities and expectation values using projectors
• Understand the polar decomposition and its role in quantum gate synthesis
• Derive key identities involving Pauli matrices from first principles
• Implement all mathematical operations in Python/NumPy and verify results

3.1 Introduction: Quantum Mechanics Is Linear Algebra

If you take one idea from this chapter, let it be this: quantum computing is linear algebra with complex numbers, dressed in physical interpretation. Every qubit is a vector. Every gate is a matrix. Every measurement is a projection. The notation looks exotic — $| \psi \rangle$, $\langle \phi | \psi \rangle$, $\hat{H}$ — but underneath it all, you are doing exactly what you did in your first linear algebra course: multiplying matrices by vectors, finding eigenvalues, and taking dot products.

This chapter builds the mathematical foundation for everything that follows. We will move methodically from the Dirac notation through the operator theory you need to read quantum computing papers and implement algorithms. Along the way, we will ground every abstract concept in Python/NumPy code so you can compute with these ideas, not just admire them.

Recurring Theme: Quantum mechanics is not magic. It is a well-defined mathematical framework. The "weirdness" arises from the fact that the vector spaces are over complex numbers and the measurement postulate is probabilistic — but the math itself is crisp, consistent, and computable.

3.1.1 Why This Mathematics Matters

Before we dive in, let's motivate why you need to master this material. Here's what each concept enables:

Every single one of these concepts will be used repeatedly throughout the rest of this book. If you master the math now, the quantum algorithms that follow will make sense. If you skip it, everything will seem like magic—and we've established that quantum computing is not magic.

3.2 Bra-Ket Notation: The Language of Quantum States

3.2.1 Kets

A ket $| \psi \rangle$ is a column vector in a complex Hilbert space $\mathcal{H}$. For an $n$-dimensional system, we write:

$$| \psi \rangle = \begin{pmatrix} \alpha_1 \\ \alpha_2 \\ \vdots \\ \alpha_n \end{pmatrix}, \quad \alpha_i \in \mathbb{C}$$

The computational basis states for a single qubit are:

$$|0\rangle = \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \qquad |1\rangle = \begin{pmatrix} 0 \\ 1 \end{pmatrix}$$

An arbitrary single-qubit state is a superposition:

$$| \psi \rangle = \alpha |0\rangle + \beta |1\rangle = \begin{pmatrix} \alpha \\ \beta \end{pmatrix}, \quad |\alpha|^2 + |\beta|^2 = 1$$

The normalization condition $|\alpha|^2 + |\beta|^2 = 1$ ensures that when we measure, probabilities sum to 1.

For two qubits, the computational basis states are:

$$|00\rangle = \begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \end{pmatrix}, \quad |01\rangle = \begin{pmatrix} 0 \\ 1 \\ 0 \\ 0 \end{pmatrix}, \quad |10\rangle = \begin{pmatrix} 0 \\ 0 \\ 1 \\ 0 \end{pmatrix}, \quad |11\rangle = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 1 \end{pmatrix}$$

And a general two-qubit state is:

$$| \psi \rangle = \alpha_{00}|00\rangle + \alpha_{01}|01\rangle + \alpha_{10}|10\rangle + \alpha_{11}|11\rangle = \begin{pmatrix} \alpha_{00} \\ \alpha_{01} \\ \alpha_{10} \\ \alpha_{11} \end{pmatrix}$$

3.2.2 Bras

A bra $\langle \psi |$ is the conjugate transpose (Hermitian conjugate) of the corresponding ket:

$$\langle \psi | = (| \psi \rangle)^\dagger = \begin{pmatrix} \alpha_1^* & \alpha_2^* & \cdots & \alpha_n^* \end{pmatrix}$$

For our single qubit: $\langle \psi | = \begin{pmatrix} \alpha^* & \beta^* \end{pmatrix}$.

Why conjugate and transpose? In standard linear algebra, the inner product on $\mathbb{C}^n$ is defined as $\langle u, v \rangle = \sum_i u_i^* v_i$, which ensures positive-definiteness: $\langle v, v \rangle = \sum_i |v_i|^2 \geq 0$ with equality iff $v = 0$. The bra-ket notation makes this natural: $\langle \psi | \phi \rangle$ is computed by taking the bra (conjugate transpose of the ket) and multiplying by the ket (column vector), which gives exactly the standard complex inner product.

Common Misconception: "The bra is just the transpose of the ket." No—it's the conjugate transpose. If $|\psi\rangle = \begin{pmatrix} 1+i \\ 2-3i \end{pmatrix}$, then $\langle\psi| = \begin{pmatrix} 1-i & 2+3i \end{pmatrix}$, not $\begin{pmatrix} 1+i & 2-3i \end{pmatrix}$. The conjugation is essential for the inner product to satisfy positive-definiteness.

3.2.3 The Inner Product

The inner product (or "bra-ket") $\langle \phi | \psi \rangle$ is the complex dot product:

$$\langle \phi | \psi \rangle = \sum_i \phi_i^* \psi_i$$

This yields a complex scalar. The squared magnitude $|\langle \phi | \psi \rangle|^2$ is the probability of measuring state $| \psi \rangle$ as $| \phi \rangle$ (if $| \phi \rangle$ is a basis state). Orthonormal basis states satisfy $\langle i | j \rangle = \delta_{ij}$.

Key properties of the inner product:

1. Conjugate symmetry: $\langle \phi | \psi \rangle = \langle \psi | \phi \rangle^*$
2. Linearity in the second argument: $\langle \phi | a\psi_1 + b\psi_2 \rangle = a\langle \phi | \psi_1 \rangle + b\langle \phi | \psi_2 \rangle$
3. Positive-definiteness: $\langle \psi | \psi \rangle \geq 0$ with equality iff $|\psi\rangle = 0$
4. Cauchy-Schwarz inequality: $|\langle \phi | \psi \rangle|^2 \leq \langle \phi | \phi \rangle \langle \psi | \psi \rangle$

Worked Example: Compute $\langle+|\psi\rangle$ where $|\psi\rangle = \frac{3}{5}|0\rangle + \frac{4i}{5}|1\rangle$ and $|+\rangle = \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)$.

$$\langle+|\psi\rangle = \frac{1}{\sqrt{2}}(\langle 0| + \langle 1|)\left(\frac{3}{5}|0\rangle + \frac{4i}{5}|1\rangle\right) = \frac{1}{\sqrt{2}}\left(\frac{3}{5} + \frac{4i}{5}\right) = \frac{3+4i}{5\sqrt{2}}$$

$$|\langle+|\psi\rangle|^2 = \frac{9+16}{50} = \frac{25}{50} = \frac{1}{2}$$

import numpy as np

# Define basis kets
ket0 = np.array([[1], [0]], dtype=complex)
ket1 = np.array([[0], [1]], dtype=complex)

# Define an arbitrary normalized state
alpha = (1 + 1j) / np.sqrt(3)
beta  = 1 / np.sqrt(3)
psi   = np.array([[alpha], [beta]], dtype=complex)

# Verify normalization
norm_sq = np.vdot(psi, psi)  # conjugate transpose dot product
print(f"⟨ψ|ψ⟩ = {norm_sq:.4f}")  # Should be 1.0

# Inner product with basis states
print(f"⟨0|ψ⟩ = {np.vdot(ket0, psi)[0,0]:.4f}")
print(f"⟨1|ψ⟩ = {np.vdot(ket1, psi)[0,0]:.4f}")

# Probability of measuring |0⟩
prob0 = np.abs(np.vdot(ket0, psi))**2
print(f"P(|0⟩) = {prob0:.4f}")

# Compute ⟨+|ψ⟩
ket_plus = np.array([[1], [1]], dtype=complex) / np.sqrt(2)
overlap = np.vdot(ket_plus, psi)
print(f"⟨+|ψ⟩ = {overlap[0,0]:.4f}")
print(f"|⟨+|ψ⟩|² = {np.abs(overlap[0,0])**2:.4f}")

3.2.4 The Outer Product

The outer product $| \psi \rangle \langle \phi |$ produces a matrix (an operator):

$$| \psi \rangle \langle \phi | = \begin{pmatrix} \alpha_1 \\ \vdots \\ \alpha_n \end{pmatrix} \begin{pmatrix} \beta_1^* & \cdots & \beta_n^* \end{pmatrix} = \begin{pmatrix} \alpha_1\beta_1^* & \cdots & \alpha_1\beta_n^* \\ \vdots & \ddots & \vdots \\ \alpha_n\beta_1^* & \cdots & \alpha_n\beta_n^* \end{pmatrix}$$

Outer products are the building blocks of projectors and density matrices. The projector onto $| \psi \rangle$ is $P_\psi = | \psi \rangle \langle \psi |$.

Properties of projectors:

1. $P_\psi^2 = P_\psi$ (idempotent)
2. $P_\psi^\dagger = P_\psi$ (Hermitian)
3. $\text{Tr}(P_\psi) = 1$ (for a rank-1 projector)
4. $P_\psi |\phi\rangle = \langle\psi|\phi\rangle |\psi\rangle$ (projects onto $|\psi\rangle$)

Worked Example: Compute $|+\rangle\langle+|$ and verify it acts as a projector.

$$|+\rangle = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 \\ 1 \end{pmatrix}, \quad \langle+| = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \end{pmatrix}$$

$$|+\rangle\langle+| = \frac{1}{2}\begin{pmatrix} 1 \\ 1 \end{pmatrix}\begin{pmatrix} 1 & 1 \end{pmatrix} = \frac{1}{2}\begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}$$

Verify idempotency: $(|+\rangle\langle+|)^2 = \frac{1}{4}\begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}\begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix} = \frac{1}{4}\begin{pmatrix} 2 & 2 \\ 2 & 2 \end{pmatrix} = \frac{1}{2}\begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix} = |+\rangle\langle+|$ ✓

Apply to $|0\rangle$: $|+\rangle\langle+|0\rangle = \frac{1}{2}\begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}\begin{pmatrix} 1 \\ 0 \end{pmatrix} = \frac{1}{2}\begin{pmatrix} 1 \\ 1 \end{pmatrix} = \frac{1}{\sqrt{2}}|+\rangle$. And indeed, $\langle+|0\rangle = 1/\sqrt{2}$, so $|+\rangle\langle+|0\rangle = \langle+|0\rangle |+\rangle = \frac{1}{\sqrt{2}}|+\rangle$ ✓

# Outer product |0⟩⟨0|
outer_00 = ket0 @ ket0.conj().T
print("|0⟩⟨0| =\n", outer_00)

# Projector onto |+⟩ = (|0⟩ + |1⟩)/√2
ket_plus = np.array([[1], [1]], dtype=complex) / np.sqrt(2)
proj_plus = ket_plus @ ket_plus.conj().T
print("|+⟩⟨+| =\n", proj_plus)

# Verify projector property: P^2 = P
print("P² = P?", np.allclose(proj_plus @ proj_plus, proj_plus))

# Apply projector to |0⟩
projected = proj_plus @ ket0
print("P₊|0⟩ =", projected.flatten())
print("⟨+|0⟩|+⟩ =", (1/np.sqrt(2)) * ket_plus.flatten())

3.2.5 Completeness Relation

For any orthonormal basis $\{|i\rangle\}$, the completeness relation holds:

$$\sum_i |i\rangle \langle i| = I$$

This is the identity operator expressed as a sum of projectors. It is the workhorse of many quantum derivations — whenever you see an identity inserted into an expression, this is what's happening.

Proof: For any state $|\psi\rangle = \sum_i c_i |i\rangle$ where $c_i = \langle i|\psi\rangle$:

$$|\psi\rangle = \sum_i c_i |i\rangle = \sum_i |i\rangle\langle i|\psi\rangle = \left(\sum_i |i\rangle\langle i|\right)|\psi\rangle$$

Since this holds for all $|\psi\rangle$, we must have $\sum_i |i\rangle\langle i| = I$.

Application: Resolving the identity. The completeness relation lets us insert identity into any expression. For example:

$$\langle\phi|\psi\rangle = \langle\phi|I|\psi\rangle = \langle\phi|\left(\sum_i |i\rangle\langle i|\right)|\psi\rangle = \sum_i \langle\phi|i\rangle\langle i|\psi\rangle$$

This is just the statement that the inner product can be computed in any basis. We'll use this technique extensively when computing measurement probabilities and deriving gate decompositions.

# Completeness for single qubit
I_from_basis = ket0 @ ket0.conj().T + ket1 @ ket1.conj().T
print("∑|i⟩⟨i| =\n", I_from_basis)
print("Equals identity?", np.allclose(I_from_basis, np.eye(2)))

3.2.6 Worked Example: Computing in Bra-Ket Notation

Let's work through a complete calculation using all the bra-ket operations we've learned.

Problem: Given $|\psi\rangle = \frac{1}{\sqrt{5}}|0\rangle + \frac{2}{\sqrt{5}}|1\rangle$ and $|\phi\rangle = \frac{1}{\sqrt{2}}(|0\rangle + i|1\rangle)$, compute:

(a) $\langle\phi|\psi\rangle$

(b) $|\langle\phi|\psi\rangle|^2$

(c) The outer product $|\psi\rangle\langle\phi|$

(d) The projector $|\psi\rangle\langle\psi|$ and its eigenvalues

Solution (a):

$$\langle\phi|\psi\rangle = \frac{1}{\sqrt{2}}(\langle 0| - i\langle 1|)\left(\frac{1}{\sqrt{5}}|0\rangle + \frac{2}{\sqrt{5}}|1\rangle\right) = \frac{1}{\sqrt{10}}(1 + 2i \cdot (-i)) = \frac{1}{\sqrt{10}}(1 + 2) = \frac{3}{\sqrt{10}}$$

Wait, let me redo this more carefully:

$$\langle\phi| = \frac{1}{\sqrt{2}}(\langle 0| + (-i)\langle 1|) = \frac{1}{\sqrt{2}}(\langle 0| - i\langle 1|)$$

Actually, if $|\phi\rangle = \frac{1}{\sqrt{2}}(|0\rangle + i|1\rangle)$, then $\langle\phi| = \frac{1}{\sqrt{2}}(\langle 0| - i\langle 1|)$.

$$\langle\phi|\psi\rangle = \frac{1}{\sqrt{2}} \cdot \frac{1}{\sqrt{5}} \cdot 1 + \frac{1}{\sqrt{2}} \cdot (-i) \cdot \frac{2}{\sqrt{5}} = \frac{1}{\sqrt{10}} - \frac{2i}{\sqrt{10}} = \frac{1-2i}{\sqrt{10}}$$

Solution (b):

$$|\langle\phi|\psi\rangle|^2 = \frac{|1-2i|^2}{10} = \frac{1+4}{10} = \frac{1}{2}$$

Solution (c):

$$|\psi\rangle\langle\phi| = \begin{pmatrix} 1/\sqrt{5} \\ 2/\sqrt{5} \end{pmatrix}\begin{pmatrix} 1/\sqrt{2} & -i/\sqrt{2} \end{pmatrix} = \frac{1}{\sqrt{10}}\begin{pmatrix} 1 & -i \\ 2 & -2i \end{pmatrix}$$

Solution (d):

$$|\psi\rangle\langle\psi| = \frac{1}{5}\begin{pmatrix} 1 \\ 2 \end{pmatrix}\begin{pmatrix} 1 & 2 \end{pmatrix} = \frac{1}{5}\begin{pmatrix} 1 & 2 \\ 2 & 4 \end{pmatrix}$$

The eigenvalues of a rank-1 projector are 1 and 0. Let's verify: $\text{Tr}(\rho) = 1/5 + 4/5 = 1$ ✓, and $\det(\rho) = 4/25 - 4/25 = 0$ ✓ (confirming one eigenvalue is 0). The other eigenvalue must be $\text{Tr} - 0 = 1$.

3.3 Hilbert Spaces

A Hilbert space $\mathcal{H}$ is a complete inner product space over the complex numbers. For quantum computing, we primarily work with finite-dimensional Hilbert spaces:

• A single qubit lives in $\mathcal{H} \cong \mathbb{C}^2$
• $n$ qubits live in $\mathcal{H}^{\otimes n} \cong \mathbb{C}^{2^n}$

The key properties are:

1. Linearity: $| \psi \rangle, | \phi \rangle \in \mathcal{H} \implies a| \psi \rangle + b| \phi \rangle \in \mathcal{H}$ for $a,b \in \mathbb{C}$
2. Inner product: $\langle \cdot | \cdot \rangle : \mathcal{H} \times \mathcal{H} \to \mathbb{C}$, conjugate-symmetric, linear in second argument, positive-definite
3. Completeness: Every Cauchy sequence converges within the space (guaranteed for finite dimensions)

The dimension of the state space grows exponentially with the number of qubits — this is the source of both quantum computing's power and the difficulty of classical simulation.

3.3.1 The Dimension Explosion

Let's make the exponential growth concrete with a table:

This table illustrates why simulating more than about 50 qubits is infeasible on classical computers, and why quantum computers with even modest numbers of high-quality qubits can perform classically intractable computations.

Recurring Theme: Quantum advantage is problem-specific. The exponential growth of the state space is necessary but not sufficient for quantum speedup. You also need the ability to manipulate amplitudes in a way that produces useful interference patterns. Not all problems in exponentially large spaces benefit from quantum computation.

3.4 Operators on Hilbert Space

3.4.1 Linear Operators

A linear operator $\hat{A}: \mathcal{H} \to \mathcal{H}$ satisfies:

$$\hat{A}(a| \psi \rangle + b| \phi \rangle) = a\hat{A}| \psi \rangle + b\hat{A}| \phi \rangle$$

In a given basis, $\hat{A}$ is represented by a matrix $A$ with entries $A_{ij} = \langle i | \hat{A} | j \rangle$.

Worked Example: Compute the matrix elements of the operator $\hat{A} = |0\rangle\langle+|$ in the computational basis.

$$A_{00} = \langle 0|A|0\rangle = \langle 0|0\rangle\langle+|0\rangle = 1 \cdot \frac{1}{\sqrt{2}} = \frac{1}{\sqrt{2}}$$ $$A_{01} = \langle 0|A|1\rangle = \langle 0|0\rangle\langle+|1\rangle = 1 \cdot \frac{1}{\sqrt{2}} = \frac{1}{\sqrt{2}}$$ $$A_{10} = \langle 1|A|0\rangle = \langle 1|0\rangle\langle+|0\rangle = 0 \cdot \frac{1}{\sqrt{2}} = 0$$ $$A_{11} = \langle 1|A|1\rangle = \langle 1|0\rangle\langle+|1\rangle = 0 \cdot \frac{1}{\sqrt{2}} = 0$$

So $A = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 0 & 0 \end{pmatrix}$.

3.4.2 The Adjoint (Hermitian Conjugate)

The adjoint $\hat{A}^\dagger$ is defined by:

$$\langle \psi | \hat{A}^\dagger | \phi \rangle = \langle \phi | \hat{A} | \psi \rangle^*$$

In matrix terms: conjugate transpose. $(A^\dagger)_{ij} = A_{ji}^*$.

Key properties of the adjoint:

1. $(AB)^\dagger = B^\dagger A^\dagger$ (note the order reversal!)
2. $(aA + bB)^\dagger = a^*A^\dagger + b^*B^\dagger$ (scalars get conjugated)
3. $(A^\dagger)^\dagger = A$ (double adjoint returns the original)
4. $\text{Tr}(A^\dagger) = (\text{Tr}(A))^*$

Worked Example: Let $A = \begin{pmatrix} 1+2i & 3 \\ 0 & 4-1i \end{pmatrix}$. Compute $A^\dagger$.

$$A^\dagger = \begin{pmatrix} (1+2i)^* & 0^* \\ 3^* & (4-1i)^* \end{pmatrix} = \begin{pmatrix} 1-2i & 0 \\ 3 & 4+i \end{pmatrix}$$

A = np.array([[1+2j, 3], [0, 4-1j]], dtype=complex)
A_dag = A.conj().T
print("A† =\n", A_dag)
print("(A†)† = A?", np.allclose(A_dag.conj().T, A))

3.4.3 Hermitian Operators

An operator is Hermitian (self-adjoint) if $\hat{A} = \hat{A}^\dagger$. Hermitian operators have real eigenvalues and orthogonal eigenvectors — they represent observables in quantum mechanics (things you can measure).

Properties of Hermitian operators:

1. All eigenvalues are real: $\lambda_i \in \mathbb{R}$
2. Eigenvectors with different eigenvalues are orthogonal: $\langle\lambda_i|\lambda_j\rangle = \delta_{ij}$
3. The operator can be diagonalized by a unitary transformation: $A = U\Lambda U^\dagger$
4. $\langle\psi|A|\psi\rangle \in \mathbb{R}$ for all $|\psi\rangle$ (expectation values are real)

Worked Example: Show that the Pauli Z matrix is Hermitian and find its eigenvalues and eigenvectors.

$$Z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}$$

$Z^\dagger = Z^{*T} = Z$ ✓ (all entries are real and the matrix is symmetric... actually, Z is diagonal, so trivially $Z = Z^\dagger$.)

Eigenvalues: $\det(Z - \lambda I) = (1-\lambda)(-1-\lambda) = 0$, so $\lambda_1 = 1, \lambda_2 = -1$.

Eigenvectors: For $\lambda_1 = 1$: $(Z-I)v = 0 \Rightarrow v \propto |0\rangle$. For $\lambda_2 = -1$: $(Z+I)v = 0 \Rightarrow v \propto |1\rangle$.

This confirms that the computational basis $\{|0\rangle, |1\rangle\}$ is the eigenbasis of Z, with eigenvalues +1 and -1 respectively. When we "measure a qubit in the computational basis," we are really measuring the observable $\hat{Z}$.

3.4.4 Unitary Operators

An operator $\hat{U}$ is unitary if:

$$\hat{U}^\dagger \hat{U} = \hat{U} \hat{U}^\dagger = I$$

Equivalently, $\hat{U}^{-1} = \hat{U}^\dagger$. Unitary operators preserve inner products:

$$\langle \hat{U}\psi | \hat{U}\phi \rangle = \langle \psi | \hat{U}^\dagger \hat{U} | \phi \rangle = \langle \psi | \phi \rangle$$

This is why quantum gates must be unitary: they must preserve the normalization of the state vector. Unitary evolution is the only kind of evolution allowed for a closed quantum system (the Schrödinger equation yields unitary time evolution).

Key properties of unitary matrices:

1. Rows and columns form orthonormal bases
2. Eigenvalues lie on the unit circle: $|\lambda| = 1$ (i.e., $\lambda = e^{i\theta}$ for some $\theta$)
3. The product of two unitaries is unitary
4. The set of $n \times n$ unitary matrices forms the unitary group $U(n)$
5. $\det(U) = e^{i\theta}$ (determinant has unit modulus)
6. Unitary operators preserve probabilities: if $|\psi\rangle$ is normalized, so is $U|\psi\rangle$

Proof that unitary operators preserve normalization:

If $\langle\psi|\psi\rangle = 1$, then $\langle U\psi|U\psi\rangle = \langle\psi|U^\dagger U|\psi\rangle = \langle\psi|\psi\rangle = 1$ ✓

Worked Example: Verify that the Hadamard gate is unitary.

$$H = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}$$

$$H^\dagger = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} = H$$

(since H is real and symmetric, $H^\dagger = H^T = H$.)

$$H^\dagger H = H^2 = \frac{1}{2}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} = \frac{1}{2}\begin{pmatrix} 2 & 0 \\ 0 & 2 \end{pmatrix} = I$$

So $H$ is both unitary and self-inverse: $H^{-1} = H^\dagger = H$.

Common Misconception: "Unitary means the matrix has determinant 1." Not necessarily. A unitary matrix has determinant $|det| = 1$, meaning $\det(U) = e^{i\theta}$ for some angle $\theta$. The special unitary group $SU(n)$ consists of unitary matrices with $\det = 1$. The Hadamard gate has $\det(H) = \frac{1}{2}((-1)(1) - (1)(1)) = -1$, so $H \in U(2)$ but $H \notin SU(2)$.

# Verify a Hadamard gate is unitary
H = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)
print("H†H =\n", H.conj().T @ H)
print("Is unitary?", np.allclose(H.conj().T @ H, np.eye(2)))

# Eigenvalues of a unitary
eigvals = np.linalg.eigvals(H)
print("Eigenvalues:", eigvals)
print("|λ| =", np.abs(eigvals))  # All should be 1

3.4.5 The Polar Decomposition

Every operator $A$ on a Hilbert space can be written as:

$$A = U|A|$$

where $U$ is unitary and $|A| = \sqrt{A^\dagger A}$ is positive semidefinite. This is the polar decomposition, analogous to writing a complex number as $z = e^{i\theta}|z|$.

For quantum computing, this means every invertible matrix can be decomposed into a rotation (unitary) and a scaling (positive). Since quantum gates must be unitary (to preserve normalization), only the unitary part is physically realizable as a gate. This is why we focus on unitary operators.

3.5 Eigenvalues, Eigenvectors, and Spectral Decomposition

3.5.1 The Eigenvalue Equation

For an operator $\hat{A}$, if there exists a non-zero vector $|v\rangle$ and scalar $\lambda$ such that:

$$\hat{A} |v\rangle = \lambda |v\rangle$$

then $|v\rangle$ is an eigenvector and $\lambda$ is the corresponding eigenvalue.

Finding eigenvalues: The eigenvalues of an $n \times n$ matrix $A$ are the roots of the characteristic polynomial:

$$\det(A - \lambda I) = 0$$

For a $2 \times 2$ matrix $A = \begin{pmatrix} a & b \\ c & d \end{pmatrix}$:

$$\det(A - \lambda I) = (a-\lambda)(d-\lambda) - bc = \lambda^2 - (a+d)\lambda + (ad-bc) = 0$$

The trace $a + d$ is the sum of eigenvalues, and the determinant $ad - bc$ is the product.

3.5.2 Spectral Decomposition

Any normal operator ($\hat{A}\hat{A}^\dagger = \hat{A}^\dagger\hat{A}$) can be diagonalized. For a Hermitian operator with eigenvalues $\lambda_i$ and orthonormal eigenvectors $|v_i\rangle$:

$$\hat{A} = \sum_i \lambda_i |v_i\rangle \langle v_i|$$

This is the spectral theorem — the operator is a weighted sum of projectors onto its eigenstates. This decomposition is central to understanding measurement in quantum mechanics: measuring an observable projects the state onto one of the operator's eigenstates, yielding the corresponding eigenvalue as the measurement outcome.

Worked Example: Spectral decomposition of the Hadamard gate.

The Hadamard gate is Hermitian ($H = H^\dagger$), so it admits a spectral decomposition.

Eigenvalues: $\det(H - \lambda I) = \frac{1}{2}(1-\lambda)(-1-\lambda) - \frac{1}{2} = 0 \Rightarrow \lambda^2 - 1 = 0 \Rightarrow \lambda = \pm 1$.

For $\lambda = +1$: $(H - I)v = 0 \Rightarrow v \propto \begin{pmatrix} 1+\sqrt{2} \\ 1 \end{pmatrix} \propto |+\rangle = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 \\ 1 \end{pmatrix}$... wait, let me compute this more carefully.

Actually, $H|+\rangle = \frac{1}{\sqrt{2}}H(|0\rangle + |1\rangle) = \frac{1}{\sqrt{2}}(|+\rangle + |-\rangle) = \frac{1}{\sqrt{2}} \cdot \frac{1}{\sqrt{2}}((|0\rangle + |1\rangle) + (|0\rangle - |1\rangle)) = |0\rangle$. Hmm, that gives $H|+\rangle = |0\rangle$, which is not an eigenvalue equation.

Let me try directly. $Hv = \lambda v$ means:

$$\frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\begin{pmatrix} v_1 \\ v_2 \end{pmatrix} = \lambda\begin{pmatrix} v_1 \\ v_2 \end{pmatrix}$$

For $\lambda = 1$: $v_1 + v_2 = \sqrt{2} v_1$ and $v_1 - v_2 = \sqrt{2} v_2$. From the first: $v_2 = (\sqrt{2}-1)v_1$. The normalized eigenvector is $|v_+\rangle \propto \begin{pmatrix} 1 \\ \sqrt{2}-1 \end{pmatrix}$.

Actually, let me just compute this numerically:

H = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)
eigvals, eigvecs = np.linalg.eigh(H)
print("Eigenvalues:", eigvals)
print("Eigenvectors:\n", eigvecs)

# Reconstruct H from spectral decomposition
H_reconstructed = sum(eigvals[i] * np.outer(eigvecs[:,i], eigvecs[:,i].conj()) for i in range(2))
print("H reconstructed:\n", H_reconstructed)
print("Matches original?", np.allclose(H, H_reconstructed))

3.5.3 Spectral Decomposition and Measurement

The spectral decomposition is intimately connected to quantum measurement. When we measure an observable $\hat{A}$ with spectral decomposition $A = \sum_i \lambda_i |v_i\rangle\langle v_i|$, the possible outcomes are the eigenvalues $\lambda_i$, and the probability of obtaining $\lambda_i$ when the system is in state $|\psi\rangle$ is:

$$p(\lambda_i) = |\langle v_i|\psi\rangle|^2$$

The expectation value (average value over many measurements) is:

$$\langle A \rangle = \langle\psi|A|\psi\rangle = \sum_i \lambda_i |\langle v_i|\psi\rangle|^2$$

And the variance is:

$$(\Delta A)^2 = \langle A^2 \rangle - \langle A \rangle^2 = \langle\psi|A^2|\psi\rangle - (\langle\psi|A|\psi\rangle)^2$$

Worked Example: Compute the expectation value and variance of measuring $\hat{Z}$ on state $|+\rangle$.

$$\langle Z \rangle = \langle+|Z|+\rangle = \frac{1}{2}(\langle 0| + \langle 1|)\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}(|0\rangle + |1\rangle) = \frac{1}{2}(1 - 1) = 0$$

$$\langle Z^2 \rangle = \langle+|Z^2|+\rangle = \langle+|I|+\rangle = 1$$

$$(\Delta Z)^2 = \langle Z^2 \rangle - \langle Z \rangle^2 = 1 - 0 = 1$$

So measuring Z on $|+\rangle$ gives expectation 0 and standard deviation 1, consistent with equal probability of outcomes +1 and -1.

3.5.4 Functions of Operators

Using spectral decomposition, we can define functions of operators:

$$f(\hat{A}) = \sum_i f(\lambda_i) |v_i\rangle \langle v_i|$$

This is how we compute matrix exponentials $e^{i\theta A}$, square roots $\sqrt{A}$, and logarithms — all essential for quantum simulation and gate synthesis.

Worked Example: Compute $e^{i\theta Z}$.

The eigenvalues of Z are $\lambda_1 = 1$ and $\lambda_2 = -1$, with eigenvectors $|0\rangle$ and $|1\rangle$.

$$e^{i\theta Z} = e^{i\theta \cdot 1}|0\rangle\langle 0| + e^{i\theta \cdot (-1)}|1\rangle\langle 1| = e^{i\theta}|0\rangle\langle 0| + e^{-i\theta}|1\rangle\langle 1|$$

$$= \begin{pmatrix} e^{i\theta} & 0 \\ 0 & e^{-i\theta} \end{pmatrix}$$

This is the $R_z(\theta)$ gate (up to a global phase)!

from scipy.linalg import expm, sqrtm

# Matrix exponential e^{iθZ} (a common phase gate form)
theta = np.pi / 4
Z = np.array([[1, 0], [0, -1]], dtype=complex)
U = expm(1j * theta * Z)
print("e^{iπZ/4} =\n", U)
print("Is unitary?", np.allclose(U.conj().T @ U, np.eye(2)))
print("Compare with R_z formula:")
Rz = np.array([[np.exp(1j*theta), 0], [0, np.exp(-1j*theta)]], dtype=complex)
print("R_z(π/4) =\n", Rz)

3.5.5 The Operator-Square-Root and Unitary Gate Synthesis

A crucial application of functions of operators is gate synthesis: expressing any single-qubit unitary as a product of elementary gates.

Theorem (Single-qubit decomposition): Any single-qubit unitary $U$ can be written as:

$$U = e^{i\alpha} R_z(\beta) R_y(\gamma) R_z(\delta)$$

for some angles $\alpha, \beta, \gamma, \delta$. This means any single-qubit quantum gate can be decomposed into rotations around the Z and Y axes, plus a global phase. This is the basis for how quantum computers implement arbitrary single-qubit operations using a small set of physically realizable gates.

3.6 Tensor Products: Building Multi-Qubit Systems

The tensor product (Kronecker product) $\otimes$ is how we combine independent quantum systems. If system A is in state $| \psi \rangle_A$ and system B is in state $| \phi \rangle_B$, the joint state is:

$$| \psi \rangle_A \otimes | \phi \rangle_B$$

For matrices, the Kronecker product of $A$ ($m \times n$) and $B$ ($p \times q$) is the $mp \times nq$ block matrix:

$$A \otimes B = \begin{pmatrix} a_{11}B & a_{12}B & \cdots & a_{1n}B \\ a_{21}B & a_{22}B & \cdots & a_{2n}B \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1}B & a_{m2}B & \cdots & a_{mn}B \end{pmatrix}$$

Key properties:

1. $(A \otimes B)(C \otimes D) = AC \otimes BD$ (mixed-product property)
2. $(A \otimes B)^\dagger = A^\dagger \otimes B^\dagger$
3. $\text{Tr}(A \otimes B) = \text{Tr}(A) \cdot \text{Tr}(B)$
4. If $A$ and $B$ are unitary, $A \otimes B$ is unitary
5. $\det(A \otimes B) = \det(A)^n \det(B)^m$ where $A$ is $m \times m$ and $B$ is $n \times n$

Worked Example: Compute $H \otimes X$ where $H = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}$ and $X = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}$.

$$H \otimes X = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 \cdot X & 1 \cdot X \\ 1 \cdot X & -1 \cdot X \end{pmatrix} = \frac{1}{\sqrt{2}}\begin{pmatrix} 0 & 1 & 0 & 1 \\ 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & -1 \\ 1 & 0 & -1 & 0 \end{pmatrix}$$

Now apply this to $|00\rangle$:

$$H \otimes X |00\rangle = \frac{1}{\sqrt{2}}\begin{pmatrix} 0 & 1 & 0 & 1 \\ 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & -1 \\ 1 & 0 & -1 & 0 \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \end{pmatrix} = \frac{1}{\sqrt{2}}\begin{pmatrix} 0 \\ 1 \\ 0 \\ 1 \end{pmatrix} = \frac{1}{\sqrt{2}}(|01\rangle + |11\rangle)$$

This makes sense: the first qubit goes through $H$ (becoming $|+\rangle$) and the second goes through $X$ (flipping $|0\rangle$ to $|1\rangle$), so the joint state is $|+\rangle \otimes |1\rangle = \frac{1}{\sqrt{2}}(|01\rangle + |11\rangle)$.

# Two-qubit computational basis states
ket00 = np.kron(ket0, ket0)  # |00⟩
ket01 = np.kron(ket0, ket1)  # |01⟩
ket10 = np.kron(ket1, ket0)  # |10⟩
ket11 = np.kron(ket1, ket1)  # |11⟩

print("|00⟩ =\n", ket00)

# A Bell state: (|00⟩ + |11⟩)/√2
bell = (ket00 + ket11) / np.sqrt(2)
print("|Φ+⟩ =\n", bell)

# Applying a CNOT gate: CNOT = |0⟩⟨0| ⊗ I + |1⟩⟨1| ⊗ X
X = np.array([[0, 1], [1, 0]], dtype=complex)
I = np.eye(2, dtype=complex)
proj0 = ket0 @ ket0.conj().T
proj1 = ket1 @ ket1.conj().T
CNOT = np.kron(proj0, I) + np.kron(proj1, X)
print("CNOT =\n", CNOT)

# Apply CNOT to |10⟩ → |11⟩
result = CNOT @ ket10
print("CNOT|10⟩ =\n", result)
print("Equals |11⟩?", np.allclose(result, ket11))

# Create a Bell state: CNOT (H ⊗ I) |00⟩
H = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)
H_I = np.kron(H, I)
bell_from_circuit = CNOT @ H_I @ ket00
print("\nBell state from circuit: CNOT (H⊗I) |00⟩ =")
print(bell_from_circuit)
print("Equals (|00⟩ + |11⟩)/√2?", np.allclose(bell_from_circuit, bell))

3.6.1 The Exponential Dimension

With $n$ qubits, the state vector has $2^n$ complex amplitudes. A 50-qubit state vector would require $2^{50} \approx 1.13 \times 10^{15}$ complex numbers — over 16 petabytes of memory. This is why we cannot classically simulate large quantum computers, and it is precisely this exponential space that quantum computers exploit.

3.6.2 Entanglement and Non-Separable States

A multi-qubit state is separable (or product) if it can be written as a tensor product of single-qubit states:

$$|\psi\rangle = |\phi_1\rangle \otimes |\phi_2\rangle \otimes \cdots \otimes |\phi_n\rangle$$

A state that is not separable is entangled. The Bell state $|\Phi^+\rangle = \frac{1}{\sqrt{2}}(|00\rangle + |11\rangle)$ is entangled—there is no way to write it as $|\psi\rangle_A \otimes |\phi\rangle_B$.

Proof of entanglement: Suppose for contradiction that $|\Phi^+\rangle = (a|0\rangle + b|1\rangle) \otimes (c|0\rangle + d|1\rangle)$. Expanding:

$$ac|00\rangle + ad|01\rangle + bc|10\rangle + bd|11\rangle = \frac{1}{\sqrt{2}}|00\rangle + \frac{1}{\sqrt{2}}|11\rangle$$

Comparing coefficients: $ac = 1/\sqrt{2}$, $ad = 0$, $bc = 0$, $bd = 1/\sqrt{2}$.

From $ad = 0$: either $a = 0$ or $d = 0$. - If $a = 0$: then $ac = 0 \neq 1/\sqrt{2}$. Contradiction. - If $d = 0$: then $bd = 0 \neq 1/\sqrt{2}$. Contradiction.

So no such factorization exists. The state is entangled. ✓

Recurring Theme: Entanglement is the quintessential quantum resource. It has no classical analogue, and it's what makes quantum computing fundamentally different from classical computing. Without entanglement, a quantum computer can be efficiently simulated classically (via the Gottesman-Knill theorem for stabilizer states, and more generally because product states can be described efficiently).

3.7 The Pauli Matrices and the Identity

The four matrices $\{I, X, Y, Z\}$ form a basis for all $2 \times 2$ complex matrices. Every single-qubit operator can be written as a linear combination:

$$\hat{A} = a_0 I + a_1 X + a_2 Y + a_3 Z, \quad a_i \in \mathbb{C}$$

3.7.1 The Matrices

$$I = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \qquad X = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \qquad Y = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} \qquad Z = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}$$

3.7.2 Algebraic Properties

The Pauli matrices satisfy:

1. Hermiticity: $\sigma_i = \sigma_i^\dagger$ for $i \in \{X, Y, Z\}$
2. Unitarity: $\sigma_i^\dagger \sigma_i = I$
3. Involutory: $\sigma_i^2 = I$
4. Commutation relations: $[\sigma_i, \sigma_j] = 2i\varepsilon_{ijk}\sigma_k$
5. Anti-commutation relations: $\{\sigma_i, \sigma_j\} = 2\delta_{ij}I$

Derivation of the commutation relation $[X, Y] = 2iZ$:

$$XY = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}\begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} = \begin{pmatrix} i & 0 \\ 0 & -i \end{pmatrix} = iZ$$

$$YX = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} = \begin{pmatrix} -i & 0 \\ 0 & i \end{pmatrix} = -iZ$$

$$[X, Y] = XY - YX = iZ - (-iZ) = 2iZ \checkmark$$

The commutation relation $[X, Y] = 2iZ$ is the mathematical origin of the uncertainty principle. Non-commuting observables cannot be simultaneously measured to arbitrary precision.

Derivation of the anti-commutation relation $\{X, Z\} = 0$:

$$XZ = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}$$

$$ZX = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}$$

$$\{X, Z\} = XZ + ZX = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} + \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix} = \begin{pmatrix} 0 & 0 \\ 0 & 0 \end{pmatrix} \checkmark$$

3.7.3 The Pauli Group and Its Significance

The Pauli group on $n$ qubits $\mathcal{P}_n$ consists of all $n$-fold tensor products of Pauli matrices with phases $\{\pm 1, \pm i\}$:

$$\mathcal{P}_n = \{\pm 1, \pm i\} \times \{I, X, Y, Z\}^{\otimes n}$$

This group is fundamental to quantum error correction and stabilizer formalism. For $n = 1$, the Pauli group has 8 elements: $\{\pm I, \pm iI, \pm X, \pm iX, \pm Y, \pm iY, \pm Z, \pm iZ\}$.

Key fact: Any two elements $P, Q \in \mathcal{P}_n$ either commute ($[P, Q] = 0$) or anti-commute ($\{P, Q\} = 0$). This binary commutation structure is the foundation of the stabilizer formalism, which we'll encounter in later chapters on quantum error correction.

3.7.4 Decomposing the Hadamard Gate

As an application of the Pauli basis decomposition, let's decompose the Hadamard gate:

$$H = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}$$

Using the decomposition $H = a_0 I + a_1 X + a_2 Y + a_3 Z$:

$$a_0 = \frac{1}{2}\text{Tr}(H \cdot I) = \frac{1}{2}\text{Tr}(H) = \frac{1}{2} \cdot 0 = 0$$

Wait, $\text{Tr}(H) = \frac{1}{\sqrt{2}}(1 + (-1)) = 0$. Let me recalculate using the correct formula for the coefficients:

$$a_i = \frac{1}{2}\text{Tr}(H \cdot \sigma_i)$$

where $\sigma_0 = I$:

$$a_0 = \frac{1}{2}\text{Tr}(H) = \frac{1}{2} \cdot 0 = 0$$

$$a_1 = \frac{1}{2}\text{Tr}(HX) = \frac{1}{2}\text{Tr}\left(\frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}\right) = \frac{1}{2}\text{Tr}\left(\frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ -1 & 1 \end{pmatrix}\right) = \frac{1}{2} \cdot \frac{2}{\sqrt{2}} = \frac{1}{\sqrt{2}}$$

$$a_2 = \frac{1}{2}\text{Tr}(HY) = \frac{1}{2}\text{Tr}\left(\frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix}\right) = \frac{1}{2}\text{Tr}\left(\frac{1}{\sqrt{2}}\begin{pmatrix} i & -i \\ -i & -i \end{pmatrix}\right) = \frac{1}{2} \cdot 0 = 0$$

$$a_3 = \frac{1}{2}\text{Tr}(HZ) = \frac{1}{2}\text{Tr}\left(\frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}\right) = \frac{1}{2}\text{Tr}\left(\frac{1}{\sqrt{2}}\begin{pmatrix} 1 & -1 \\ 1 & 1 \end{pmatrix}\right) = \frac{1}{2} \cdot \frac{2}{\sqrt{2}} = \frac{1}{\sqrt{2}}$$

So $H = \frac{1}{\sqrt{2}}(X + Z)$. Let's verify:

$$\frac{1}{\sqrt{2}}\left(\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} + \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}\right) = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} = H \checkmark$$

This decomposition tells us that the Hadamard gate is an equal superposition of X and Z rotations—it's a $\pi$ rotation around the axis $(\hat{x} + \hat{z})/\sqrt{2}$ on the Bloch sphere.

X = np.array([[0, 1], [1, 0]], dtype=complex)
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
I = np.eye(2, dtype=complex)

# Verify properties
print("X² = I?", np.allclose(X @ X, I))
print("Y² = I?", np.allclose(Y @ Y, I))
print("Z² = I?", np.allclose(Z @ Z, I))

# Commutator [X, Y] = XY - YX
comm_XY = X @ Y - Y @ X
print("[X, Y] =\n", comm_XY)
print("Equals 2iZ?", np.allclose(comm_XY, 2j * Z))

# Decompose Hadamard in Pauli basis
H = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)
a0 = 0.5 * np.trace(H @ I)
a1 = 0.5 * np.trace(H @ X)
a2 = 0.5 * np.trace(H @ Y)
a3 = 0.5 * np.trace(H @ Z)
print(f"\nH = {a0:.4f}·I + {a1:.4f}·X + {a2:.4f}·Y + {a3:.4f}·Z")
print(f"H = (X + Z)/√2")
H_reconstructed = a0*I + a1*X + a2*Y + a3*Z
print("Reconstruction correct?", np.allclose(H, H_reconstructed))

3.8 The Bloch Sphere Representation

Any single-qubit pure state can be written as:

$$| \psi \rangle = \cos\frac{\theta}{2} |0\rangle + e^{i\phi} \sin\frac{\theta}{2} |1\rangle$$

with $\theta \in [0, \pi]$ and $\phi \in [0, 2\pi)$. These are spherical coordinates on the Bloch sphere:

                     |0⟩ (north pole)
                      *
                     /|\
                    / | \
                   /  |  \     θ = polar angle (0 to π)
                  /   |   \    φ = azimuthal angle (0 to 2π)
                 /    |    \   r = (sinθcosφ, sinθsinφ, cosθ)
                /     |     \
               /      |      * ← state |ψ⟩
              /       |     / \
             /     r  |    /   \
            /        |   /     \
           /         |  /       \
          *----------+-/---------*-----> x-axis (equator: |+⟩)
                     / \
                    /   \
                   /     \
                  /       \
                 /         \
                *-----------*
               |1⟩ (south pole)

The Pauli matrices generate rotations on the Bloch sphere:

$$R_x(\theta) = e^{-i\theta X/2} = \cos\frac{\theta}{2}I - i\sin\frac{\theta}{2}X$$ $$R_y(\theta) = e^{-i\theta Y/2} = \cos\frac{\theta}{2}I - i\sin\frac{\theta}{2}Y$$ $$R_z(\theta) = e^{-i\theta Z/2} = \cos\frac{\theta}{2}I - i\sin\frac{\theta}{2}Z$$

Derivation of $R_z(\theta)$: Using the spectral decomposition, $Z = |0\rangle\langle 0| - |1\rangle\langle 1|$:

$$e^{-i\theta Z/2} = e^{-i\theta/2}|0\rangle\langle 0| + e^{i\theta/2}|1\rangle\langle 1| = \begin{pmatrix} e^{-i\theta/2} & 0 \\ 0 & e^{i\theta/2} \end{pmatrix}$$

This is a rotation around the Z-axis of the Bloch sphere. The factor of 1/2 in the exponent accounts for the fact that a full $2\pi$ rotation on the Bloch sphere corresponds to a $4\pi$ rotation in physical space (a deep fact about spin-1/2 particles).

Worked Example: Show that $R_x(\pi) = -iX$, and explain why this is a $\pi$-pulse (a bit flip).

$$R_x(\pi) = \cos\frac{\pi}{2}I - i\sin\frac{\pi}{2}X = 0 \cdot I - i \cdot X = -iX$$

Since global phase is irrelevant, $R_x(\pi) \sim X$ (a bit flip). Up to a global phase, this is the X gate.

def Rx(theta):
    return np.cos(theta/2)*I - 1j*np.sin(theta/2)*X

def Ry(theta):
    return np.cos(theta/2)*I - 1j*np.sin(theta/2)*Y

def Rz(theta):
    return np.cos(theta/2)*I - 1j*np.sin(theta/2)*Z

# Verify unitarity
for gate in [Rx(0.7), Ry(1.2), Rz(2.1)]:
    print("Unitary?", np.allclose(gate.conj().T @ gate, np.eye(2)))

# Verify Rx(π) = -iX
print("\nRx(π) =\n", Rx(np.pi))
print("-iX =\n", -1j*X)
print("Equal (up to global phase)?", np.allclose(Rx(np.pi), -1j*X))

# Any single-qubit unitary can be decomposed as Rz(α)Ry(β)Rz(γ)
# This is the Z-Y decomposition theorem
# Let's verify for the Hadamard gate
print("\nHadamard = Rz(π) Ry(π/2) Rz(0):")
U_decomposed = Rz(np.pi) @ Ry(np.pi/2) @ Rz(0)
print("U_decomposed =\n", U_decomposed)
print("Matches H (up to global phase)?", np.allclose(np.abs(U_decomposed), np.abs(H)))

3.9 The Trace and Partial Trace

The trace of an operator is the sum of its diagonal elements:

$$\text{Tr}(\hat{A}) = \sum_i \langle i | \hat{A} | i \rangle$$

Key properties of the trace:

1. Cyclic: $\text{Tr}(AB) = \text{Tr}(BA)$
2. Linear: $\text{Tr}(aA + bB) = a\text{Tr}(A) + b\text{Tr}(B)$
3. Basis-independent: The trace has the same value in any orthonormal basis
4. Of a density matrix: $\text{Tr}(\rho) = 1$
5. Of a projector: $\text{Tr}(|\psi\rangle\langle\psi|) = \langle\psi|\psi\rangle = 1$
6. Of a tensor product: $\text{Tr}(A \otimes B) = \text{Tr}(A) \cdot \text{Tr}(B)$

Worked Example: Show that $\text{Tr}(|\psi\rangle\langle\phi|) = \langle\phi|\psi\rangle$.

$$\text{Tr}(|\psi\rangle\langle\phi|) = \sum_i \langle i|\psi\rangle\langle\phi|i\rangle = \sum_i \langle\phi|i\rangle\langle i|\psi\rangle = \langle\phi|\left(\sum_i |i\rangle\langle i|\right)|\psi\rangle = \langle\phi|\psi\rangle$$

where we used the completeness relation $\sum_i |i\rangle\langle i| = I$.

3.9.1 The Partial Trace

The partial trace is how we describe subsystems of a composite quantum system. For a bipartite state $\rho_{AB}$, the reduced state of system A is:

$$\rho_A = \text{Tr}_B(\rho_{AB}) = \sum_i (I_A \otimes \langle i_B|) \rho_{AB} (I_A \otimes |i_B\rangle)$$

This is essential for understanding entanglement and decoherence. The partial trace "traces out" system B, leaving us with a description of system A alone.

Worked Example: Compute the reduced density matrix of qubit A for the Bell state $|\Phi^+\rangle = \frac{1}{\sqrt{2}}(|00\rangle + |11\rangle)$.

$$\rho_{AB} = |\Phi^+\rangle\langle\Phi^+| = \frac{1}{2}(|00\rangle + |11\rangle)(\langle 00| + \langle 11|) = \frac{1}{2}(|00\rangle\langle 00| + |00\rangle\langle 11| + |11\rangle\langle 00| + |11\rangle\langle 11|)$$

$$\rho_A = \text{Tr}_B(\rho_{AB}) = \sum_{i=0}^{1} (\langle i|_B \otimes I_A)\rho_{AB}(|i\rangle_B \otimes I_A)$$

For $i = 0$: $$(\langle 0|_B \otimes I_A)\rho_{AB}(|0\rangle_B \otimes I_A) = \frac{1}{2}(|0\rangle_A\langle 0|_A \cdot 1 + 0) = \frac{1}{2}|0\rangle\langle 0|$$

For $i = 1$: $$(\langle 1|_B \otimes I_A)\rho_{AB}(|1\rangle_B \otimes I_A) = \frac{1}{2}(0 + |1\rangle_A\langle 1|_A \cdot 1) = \frac{1}{2}|1\rangle\langle 1|$$

$$\rho_A = \frac{1}{2}|0\rangle\langle 0| + \frac{1}{2}|1\rangle\langle 1| = \frac{1}{2}I$$

The reduced state is the maximally mixed state! This confirms that the Bell state is maximally entangled—each qubit individually contains no information, but together they are perfectly correlated.

# Partial trace: given a 2-qubit density matrix, trace out qubit B
rho_AB = bell @ bell.conj().T  # density matrix of Bell state

def partial_trace(rho, dims, keep):
    rho_tensor = rho.reshape(dims + dims)
    trace_over = [i for i in range(len(dims)) if i not in keep]
    for axis in sorted(trace_over, reverse=True):
        rho_tensor = np.trace(rho_tensor, axis1=axis, axis2=axis + len(dims))
    return rho_tensor

rho_A = partial_trace(rho_AB, [2, 2], [0])
print("Reduced state of qubit A:\n", rho_A)
print("Is maximally mixed?", np.allclose(rho_A, I/2))

# Compute purity
purity_A = np.trace(rho_A @ rho_A).real
print(f"Purity of reduced state: {purity_A:.4f}")
print(f"Purity of Bell state (full): {np.trace(rho_AB @ rho_AB).real:.4f}")

The purity of the reduced state is 0.5 (mixed), while the purity of the full Bell state is 1.0 (pure). This is the hallmark of entanglement: a pure global state with mixed local states.

3.10 Important Identities and Theorems

Let's collect some important identities that will be used throughout the book.

3.10.1 The Hadamard Conjugation Identity

For any Pauli matrix $\sigma \in \{X, Y, Z\}$:

$$H\sigma H = \begin{cases} Z & \text{if } \sigma = X \\ -Y & \text{if } \sigma = Y \\ X & \text{if } \sigma = Z \end{cases}$$

This shows that $H$ swaps X and Z (and flips Y). This is why the Hadamard gate converts between the computational basis (Z eigenbasis) and the Hadamard basis (X eigenbasis).

3.10.2 The Pauli Twirl Identity

For any $2 \times 2$ matrix $A$:

$$\frac{1}{4}\sum_{\sigma \in \{I,X,Y,Z\}} \sigma A \sigma = \frac{\text{Tr}(A)}{2} I$$

This identity is crucial in quantum error correction: it shows that averaging over all Pauli operations "twirls" any operator into a multiple of the identity.

3.10.3 The Cauchy-Schwarz Inequality

For any two states $|\psi\rangle$ and $|\phi\rangle$:

$$|\langle\psi|\phi\rangle|^2 \leq \langle\psi|\psi\rangle \cdot \langle\phi|\phi\rangle$$

Equality holds if and only if $|\psi\rangle$ and $|\phi\rangle$ are linearly dependent. This inequality underlies many proofs in quantum information, including the no-cloning theorem.

3.11 Qiskit: Implementing the Math

Let's put all these mathematical concepts together in a comprehensive Qiskit example:

from qiskit import QuantumCircuit, transpile
from qiskit_aer import AerSimulator
from qiskit.quantum_info import Statevector, Operator, DensityMatrix
import numpy as np

# Example 1: Verify all Pauli algebraic properties
print("=== Pauli Matrix Properties ===")
X = np.array([[0, 1], [1, 0]], dtype=complex)
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
I = np.eye(2, dtype=complex)

# Commutation relations
print(f"[X,Y] = 2iZ? {np.allclose(X@Y - Y@X, 2j*Z)}")
print(f"[Y,Z] = 2iX? {np.allclose(Y@Z - Z@Y, 2j*X)}")
print(f"[Z,X] = 2iY? {np.allclose(Z@X - X@Z, 2j*Y)}")

# Anti-commutation relations
print(f"{{X,Y}} = 0? {np.allclose(X@Y + Y@X, np.zeros((2,2)))}")
print(f"{{Y,Z}} = 0? {np.allclose(Y@Z + Z@Y, np.zeros((2,2)))}")
print(f"{{Z,X}} = 0? {np.allclose(Z@X + X@Z, np.zeros((2,2)))}")

# Example 2: Create and verify a Bell state
print("\n=== Bell State Verification ===")
bell_qc = QuantumCircuit(2)
bell_qc.h(0)
bell_qc.cx(0, 1)
bell_state = Statevector.from_instruction(bell_qc)
print(f"Bell state: {bell_state.data}")

# Verify entanglement: partial trace should give maximally mixed state
bell_dm = DensityMatrix(bell_state)
print(f"Full density matrix purity: {np.trace(bell_dm.data @ bell_dm.data).real:.4f}")

# Example 3: Spectral decomposition of Z
print("\n=== Spectral Decomposition of Z ===")
eigvals, eigvecs = np.linalg.eigh(Z)
Z_reconstructed = sum(eigvals[i] * np.outer(eigvecs[:,i], eigvecs[:,i].conj()) for i in range(2))
print(f"Eigenvalues: {eigvals}")
print(f"Reconstruction matches? {np.allclose(Z, Z_reconstructed)}")

# Example 3: Demonstrate measurement and projection
from qiskit import QuantumCircuit, transpile
from qiskit_aer import AerSimulator
from qiskit.quantum_info import Statevector
import numpy as np

# Prepare |ψ⟩ = (3/5)|0⟩ + (4/5)|1⟩
theta = 2 * np.arcsin(4/5)  # sin(θ/2) = 4/5
qc = QuantumCircuit(1)
qc.ry(theta, 0)

state = Statevector.from_instruction(qc)
print(f"State: {state.data}")
print(f"Probability of |0⟩: {abs(state.data[0])**2:.4f}")
print(f"Probability of |1⟩: {abs(state.data[1])**2:.4f}")

# Measure
qc.measure_all()
simulator = AerSimulator()
result = simulator.run(transpile(qc, simulator), shots=10000).result()
counts = result.get_counts()
print(f"Measurement results: {counts}")

# Example 4: Verify unitarity of rotation gates
from scipy.linalg import expm
print("\n=== Rotation Gate Unitarity ===")
for angle in [0.5, 1.0, np.pi/4]:
    for label, matrix in [("X", X), ("Y", Y), ("Z", Z)]:
        U = expm(-1j * angle/2 * matrix)
        is_unitary = np.allclose(U.conj().T @ U, np.eye(2))
        print(f"R_{label}({angle:.2f}) unitary? {is_unitary}")

3.12 The Schmidt Decomposition

The Schmidt decomposition is a powerful tool for analyzing bipartite (two-party) quantum states. It tells us that any pure state of a bipartite system can be written in a canonical form:

Theorem (Schmidt decomposition): For any pure state $|\psi\rangle_{AB}$ of a bipartite system with dimensions $d_A$ and $d_B$, there exist orthonormal states $\{|i\rangle_A\}$ for system A and $\{|i\rangle_B\}$ for system B such that:

$$|\psi\rangle_{AB} = \sum_{i=1}^{r} \sqrt{\lambda_i} |i\rangle_A |i\rangle_B$$

where $r = \min(d_A, d_B)$ is the Schmidt rank, and $\lambda_i \geq 0$ with $\sum_i \lambda_i = 1$ are the Schmidt coefficients.

The number of non-zero Schmidt coefficients is the Schmidt number, which quantifies the entanglement of the state:

• Schmidt number = 1: the state is separable (not entangled)
• Schmidt number > 1: the state is entangled
• Schmidt number = $\min(d_A, d_B)$: the state is maximally entangled

Worked Example: Find the Schmidt decomposition of $|\psi\rangle = \frac{1}{2}|00\rangle + \frac{1}{2}|01\rangle + \frac{1}{2}|10\rangle + \frac{1}{2}|11\rangle = |+\rangle \otimes |+\rangle$.

Since this is already a product state, the Schmidt decomposition has Schmidt number 1:

$$|\psi\rangle = 1 \cdot |+\rangle_A \otimes |+\rangle_B$$

The Schmidt coefficients are just $\{1\}$, confirming no entanglement.

Worked Example: Find the Schmidt decomposition of $|\psi\rangle = \frac{3}{5}|00\rangle + \frac{4}{5}|11\rangle$.

This is already in Schmidt form with orthonormal states $|0\rangle_A, |0\rangle_B$ and $|1\rangle_A, |1\rangle_B$. The Schmidt coefficients are $(3/5)^2 = 9/25$ and $(4/5)^2 = 16/25$... wait, let me be more careful.

The state is $|\psi\rangle = \frac{3}{5}|00\rangle + \frac{4}{5}|11\rangle$. In matrix form, the density matrix of the full system is $\rho_{AB} = |\psi\rangle\langle\psi|$.

To find the Schmidt decomposition, we compute the reduced density matrix of system A:

$$\rho_A = \text{Tr}_B(\rho_{AB}) = \left(\frac{3}{5}\right)^2|0\rangle\langle 0| + \left(\frac{4}{5}\right)^2|1\rangle\langle 1| = \frac{9}{25}|0\rangle\langle 0| + \frac{16}{25}|1\rangle\langle 1|$$

The eigenvalues of $\rho_A$ are $9/25$ and $16/25$, which are the Schmidt coefficients squared. So the Schmidt coefficients are $\lambda_1 = 3/5$ and $\lambda_2 = 4/5$, and the Schmidt decomposition is:

$$|\psi\rangle = \frac{3}{5}|0\rangle_A|0\rangle_B + \frac{4}{5}|1\rangle_A|1\rangle_B$$

which is what we started with (already in Schmidt form).

The entanglement entropy is defined as $S = -\sum_i \lambda_i^2 \log_2(\lambda_i^2)$, which for this state is:

$$S = -\frac{9}{25}\log_2\frac{9}{25} - \frac{16}{25}\log_2\frac{16}{25} \approx 0.94 \text{ bits}$$

For a maximally entangled state like $|\Phi^+\rangle$, the entanglement entropy is $\log_2(2) = 1$ bit.

import numpy as np

# Schmidt decomposition of (3/5)|00⟩ + (4/5)|11⟩
psi = np.array([3/5, 0, 0, 4/5], dtype=complex).reshape(2, 2)

# SVD: psi = U @ diag(s) @ V†
U, s, Vh = np.linalg.svd(psi)
print("Schmidt coefficients:", s)
print("Schmidt number:", np.sum(s > 1e-10))

# Entanglement entropy
probs = s**2
entropy = -np.sum(probs * np.log2(probs + 1e-15))
print(f"Entanglement entropy: {entropy:.4f} bits")

# For Bell state |Φ+⟩
bell = np.array([1, 0, 0, 1], dtype=complex).reshape(2, 2) / np.sqrt(2)
U, s, Vh = np.linalg.svd(bell)
probs_bell = s**2
entropy_bell = -np.sum(probs_bell * np.log2(probs_bell + 1e-15))
print(f"\nBell state Schmidt coefficients: {s}")
print(f"Bell state entanglement entropy: {entropy_bell:.4f} bits")

3.13 The Density Matrix Formalism (Extended)

3.13.1 Pure States vs. Mixed States

A pure state is a state that can be written as a ket $|\psi\rangle$. Its density matrix is $\rho = |\psi\rangle\langle\psi|$, a rank-1 projector.

A mixed state is a statistical ensemble of pure states: $\rho = \sum_i p_i |\psi_i\rangle\langle\psi_i|$ where $p_i \geq 0$ and $\sum_i p_i = 1$. The ensemble $\{(p_i, |\psi_i\rangle)\}$ represents: "with probability $p_i$, the system is in state $|\psi_i\rangle$."

Key distinction: The pure state $|+\rangle$ and the mixed state $\frac{1}{2}|0\rangle\langle 0| + \frac{1}{2}|1\rangle\langle 1|$ both give 50/50 measurement probabilities in the computational basis, but they are physically different:

3.13.2 Properties of Density Matrices

A valid density matrix must satisfy:

1. Hermiticity: $\rho = \rho^\dagger$
2. Positive semidefiniteness: $\langle\psi|\rho|\psi\rangle \geq 0$ for all $|\psi\rangle$ (all eigenvalues $\geq 0$)
3. Unit trace: $\text{Tr}(\rho) = 1$

Theorem: A state is pure if and only if $\text{Tr}(\rho^2) = 1$. A state is mixed if and only if $\text{Tr}(\rho^2) < 1$.

Proof: If $\rho = |\psi\rangle\langle\psi|$, then $\rho^2 = |\psi\rangle\langle\psi|\psi\rangle\langle\psi| = |\psi\rangle\langle\psi| = \rho$, so $\text{Tr}(\rho^2) = \text{Tr}(\rho) = 1$.

If $\rho = \sum_i p_i |\psi_i\rangle\langle\psi_i|$ with more than one non-zero $p_i$, then $\text{Tr}(\rho^2) = \sum_{i,j} p_i p_j |\langle\psi_i|\psi_j\rangle|^2 \leq \sum_i p_i^2 < (\sum_i p_i)^2 = 1$ (by the Cauchy-Schwarz inequality applied to the probability distribution).

3.13.3 The Von Neumann Entropy

The von Neumann entropy of a density matrix is:

$$S(\rho) = -\text{Tr}(\rho \log_2 \rho) = -\sum_i \lambda_i \log_2 \lambda_i$$

where $\lambda_i$ are the eigenvalues of $\rho$.

Properties: - $S(\rho) = 0$ for pure states (only one non-zero eigenvalue, which is 1) - $S(\rho) = \log_2 d$ for the maximally mixed state $\rho = I/d$ - $S(\rho) \leq \log_2 d$ for any state in a $d$-dimensional space

For a single qubit, the maximum entropy is $\log_2 2 = 1$ bit, achieved by the maximally mixed state $I/2$.

# Compare purity and entropy for different states
import numpy as np

def von_neumann_entropy(rho):
    """Compute von Neumann entropy S = -Tr(ρ log₂ ρ)"""
    eigvals = np.linalg.eigvalsh(rho)
    # Filter out zero eigenvalues
    eigvals = eigvals[eigvals > 1e-15]
    return -np.sum(eigvals * np.log2(eigvals))

# Pure state: |+⟩
rho_plus = np.array([[0.5, 0.5], [0.5, 0.5]], dtype=complex)
print(f"|+⟩ purity: {np.trace(rho_plus @ rho_plus).real:.4f}")
print(f"|+⟩ entropy: {von_neumann_entropy(rho_plus):.4f} bits")

# Mixed state: I/2
rho_mixed = np.array([[0.5, 0.0], [0.0, 0.5]], dtype=complex)
print(f"\nI/2 purity: {np.trace(rho_mixed @ rho_mixed).real:.4f}")
print(f"I/2 entropy: {von_neumann_entropy(rho_mixed):.4f} bits")

# Partially mixed state: 3/4|0⟩⟨0| + 1/4|1⟩⟨1|
rho_partial = np.array([[0.75, 0.0], [0.0, 0.25]], dtype=complex)
print(f"\nρ_partial purity: {np.trace(rho_partial @ rho_partial).real:.4f}")
print(f"ρ_partial entropy: {von_neumann_entropy(rho_partial):.4f} bits")

3.14 The Heisenberg Uncertainty Principle for Qubits

The most famous consequence of non-commuting observables is the Heisenberg uncertainty principle. For two observables $A$ and $B$:

$$\Delta A \cdot \Delta B \geq \frac{1}{2}|\langle[A, B]\rangle|$$

where $\Delta A = \sqrt{\langle A^2 \rangle - \langle A \rangle^2}$ is the standard deviation of $A$ in the state $|\psi\rangle$.

Worked Example: For $A = X$ and $B = Z$, we have $[X, Z] = -2iY$, so:

$$\Delta X \cdot \Delta Z \geq \frac{1}{2}|\langle -2iY \rangle| = |\langle Y \rangle|$$

For the state $|+\rangle$: - $\langle X \rangle = \langle+|X|+\rangle = 1$ - $\langle Z \rangle = \langle+|Z|+\rangle = 0$ - $\langle X^2 \rangle = \langle+|I|+\rangle = 1$ - $\langle Z^2 \rangle = \langle+|I|+\rangle = 1$ - $\Delta X = \sqrt{1 - 1} = 0$ - $\Delta Z = \sqrt{1 - 0} = 1$ - $\Delta X \cdot \Delta Z = 0$ - $|\langle Y \rangle| = |\langle+|Y|+\rangle| = |\frac{1}{2}(\langle 0| + \langle 1|)Y(|0\rangle + |1\rangle)| = |\frac{1}{2}(\langle 0|Y|1\rangle + \langle 1|Y|0\rangle)| = |\frac{1}{2}(-i + i)| = 0$

So $\Delta X \cdot \Delta Z = 0 \geq 0 = |\langle Y \rangle|$. The inequality is satisfied with equality because $|+\rangle$ is an eigenstate of $X$, so $\Delta X = 0$ and $\Delta Z$ is maximal.

This illustrates the complementarity of $X$ and $Z$: if you know $X$ precisely ($\Delta X = 0$), you know nothing about $Z$ ($\Delta Z = 1$), and vice versa.

Common Misconception: "The uncertainty principle is about measurement disturbance—you can't measure position without disturbing momentum." While this is one interpretation, the uncertainty principle is more fundamental: it's a mathematical property of the operators themselves. Even without any measurement, a quantum state cannot be simultaneously an eigenstate of non-commuting operators. The uncertainty is inherent in the state, not just in the measurement process.

3.15 Change of Basis

A fundamental operation in quantum computing is changing the basis in which we express a state or operator. Given two orthonormal bases $\{|i\rangle\}$ and $\{|j'\rangle\}$ related by a unitary transformation $U$ (where $|j'\rangle = \sum_i U_{ji}|i\rangle$), we can transform:

States: $|\psi'\rangle = U|\psi\rangle$

Operators: $A' = UAU^\dagger$

Worked Example: Transform the Z operator from the computational basis to the Hadamard basis.

In the computational basis: $Z = |0\rangle\langle 0| - |1\rangle\langle 1|$

In the Hadamard basis: $Z' = HZH$

$$Z' = HZH = \frac{1}{2}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} = \frac{1}{2}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\begin{pmatrix} 1 & 1 \\ -1 & 1 \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} = X$$

So in the Hadamard basis, the Z operator becomes the X operator! This makes perfect sense: Z measures "computational basis spin" while X measures "Hadamard basis spin," and changing basis swaps the two.

Similarly: $X' = HXH = Z$ and $Y' = HYH = -Y$.

This is why the Hadamard gate converts between Z-measurements and X-measurements—it changes the measurement basis.

import numpy as np

H = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)
X = np.array([[0, 1], [1, 0]], dtype=complex)
Y = np.array([[0, -1j], [1j, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)

# Change of basis: Z in Hadamard basis = HZH
Z_prime = H @ Z @ H
print(f"Z in Hadamard basis = HZH =\n{Z_prime}")
print(f"This equals X? {np.allclose(Z_prime, X)}")

# Change of basis: X in Hadamard basis = HXH
X_prime = H @ X @ H
print(f"\nX in Hadamard basis = HXH =\n{X_prime}")
print(f"This equals Z? {np.allclose(X_prime, Z)}")

# Change of basis: Y in Hadamard basis = HYH
Y_prime = H @ Y @ H
print(f"\nY in Hadamard basis = HYH =\n{Y_prime}")
print(f"This equals -Y? {np.allclose(Y_prime, -Y)}")

3.16 Qiskit: Putting It All Together

Let's create a comprehensive example that ties together all the mathematical concepts from this chapter:

from qiskit import QuantumCircuit, transpile
from qiskit_aer import AerSimulator
from qiskit.quantum_info import Statevector, Operator, DensityMatrix, partial_trace
import numpy as np

print("=" * 60)
print("COMPREHENSIVE LINEAR ALGEBRA REVIEW FOR QUANTUM COMPUTING")
print("=" * 60)

# 1. State vectors and inner products
print("\n--- 1. State Vectors and Inner Products ---")
ket_psi = np.array([[3/5], [4/5]], dtype=complex)  # NOT normalized!
norm = np.linalg.norm(ket_psi)
ket_psi_normalized = ket_psi / norm
print(f"|ψ⟩ = {ket_psi_normalized.flatten()}")
print(f"⟨ψ|ψ⟩ = {np.vdot(ket_psi_normalized, ket_psi_normalized):.4f}")

# 2. Outer products and projectors
print("\n--- 2. Outer Products and Projectors ---")
proj = ket_psi_normalized @ ket_psi_normalized.conj().T
print(f"|ψ⟩⟨ψ| =\n{proj}")
print(f"P² = P? {np.allclose(proj @ proj, proj)}")
print(f"Tr(P) = {np.trace(proj).real:.4f}")

# 3. Unitary operations
print("\n--- 3. Unitary Operations ---")
theta = np.pi / 3
Rz_gate = np.array([[np.exp(-1j*theta/2), 0], [0, np.exp(1j*theta/2)]], dtype=complex)
print(f"R_z(π/3) =\n{Rz_gate}")
print(f"Unitary? {np.allclose(Rz_gate.conj().T @ Rz_gate, np.eye(2))}")
result = Rz_gate @ ket_psi_normalized
print(f"R_z(π/3)|ψ⟩ = {result.flatten()}")

# 4. Spectral decomposition
print("\n--- 4. Spectral Decomposition ---")
A = np.array([[3, 1+1j], [1-1j, 2]], dtype=complex)
print(f"A =\n{A}")
print(f"Hermitian? {np.allclose(A, A.conj().T)}")
eigvals, eigvecs = np.linalg.eigh(A)
print(f"Eigenvalues: {eigvals}")
A_reconstructed = sum(eigvals[i] * np.outer(eigvecs[:,i], eigvecs[:,i].conj()) for i in range(2))
print(f"Reconstruction matches? {np.allclose(A, A_reconstructed)}")

# 5. Tensor products and entanglement
print("\n--- 5. Tensor Products and Entanglement ---")
qc_bell = QuantumCircuit(2)
qc_bell.h(0)
qc_bell.cx(0, 1)
bell_state = Statevector.from_instruction(qc_bell)
print(f"Bell state: {bell_state.data}")
bell_dm = DensityMatrix(bell_state)
rho_A = partial_trace(bell_dm, [1], dims=[2, 2])
print(f"Reduced density matrix of qubit A:\n{rho_A.data}")
print(f"Purity of Bell state: {np.trace(bell_dm.data @ bell_dm.data).real:.4f}")
print(f"Purity of reduced state: {np.trace(rho_A.data @ rho_A.data).real:.4f}")
print(f"Reduced state is maximally mixed? {np.allclose(rho_A.data, np.eye(2)/2)}")

This comprehensive example demonstrates all the key mathematical concepts: state vectors, inner products, outer products, projectors, unitary operations, spectral decomposition, tensor products, entanglement, and partial traces. Every concept in this chapter has a concrete numerical implementation.

3.17 Common Pitfalls and How to Avoid Them

As you work with the mathematics of quantum computing, be aware of these common mistakes:

Pitfall 1: Confusing Bra and Ket Order

The inner product $\langle\phi|\psi\rangle$ is a complex number (scalar). The outer product $|\psi\rangle\langle\phi|$ is a matrix (operator). Getting the order wrong changes the result completely.

Rule: If the bra is on the left, it's an inner product (number). If the ket is on the left, it's an outer product (matrix).

Pitfall 2: Forgetting Complex Conjugation

$\langle\phi|\psi\rangle \neq \langle\psi|\phi\rangle$. They are complex conjugates: $\langle\phi|\psi\rangle = \langle\psi|\phi\rangle^*$. This matters when computing probabilities: $|\langle\phi|\psi\rangle|^2 = \langle\phi|\psi\rangle \cdot \langle\phi|\psi\rangle^* = \langle\phi|\psi\rangle\langle\psi|\phi\rangle$.

Pitfall 3: Treating Global Phase as Physical

Two states that differ by a global phase $|\psi\rangle$ and $e^{i\gamma}|\psi\rangle$ are physically indistinguishable. No experiment can tell them apart. But relative phase (between components of a superposition) is physical and has measurable consequences.

Pitfall 4: Confusing Tensor Product Order

$|a\rangle \otimes |b\rangle \neq |b\rangle \otimes |a\rangle$ in general. The qubit ordering convention matters: in Qiskit, qubit 0 is the least significant (rightmost) bit.

Pitfall 5: Assuming All States Can Be Factored

Not all multi-qubit states can be written as $|\psi\rangle_A \otimes |\phi\rangle_B$. Entangled states like $|\Phi^+\rangle = \frac{1}{\sqrt{2}}(|00\rangle + |11\rangle)$ cannot be factored—this is what makes them entangled.

Pitfall 6: Treating the Density Matrix as Just Any Matrix

A density matrix must be Hermitian, positive semidefinite, and have trace 1. Not every $2 \times 2$ matrix satisfies these conditions. Always verify these properties.

# Quick check: is a matrix a valid density matrix?
def is_valid_density_matrix(rho):
    is_hermitian = np.allclose(rho, rho.conj().T)
    is_positive = all(eigval >= -1e-10 for eigval in np.linalg.eigvalsh(rho))
    is_unit_trace = np.allclose(np.trace(rho), 1)
    return is_hermitian and is_positive and is_unit_trace

# Test
rho_pure = np.array([[0.5, 0.5], [0.5, 0.5]], dtype=complex)  # |+⟩⟨+|
rho_mixed = np.array([[0.5, 0.0], [0.0, 0.5]], dtype=complex)  # I/2
rho_invalid = np.array([[0.5, 0.5], [0.5, 0.5]], dtype=complex) * 2  # trace = 2, not 1

print(f"|+⟩⟨+| valid? {is_valid_density_matrix(rho_pure)}")
print(f"I/2 valid? {is_valid_density_matrix(rho_mixed)}")
print(f"Invalid matrix valid? {is_valid_density_matrix(rho_invalid)}")