|
| 1 | + |
| 2 | +Variational Quantum Eigensolver (VQE) with Jordan–Wigner |
| 3 | + |
| 4 | + |
| 5 | +The Variational Quantum Eigensolver (VQE) is a hybrid algorithm that uses a parameterized quantum circuit (ansatz) to estimate the ground-state energy of a Hamiltonian . In the molecular or fermionic context, one first maps the fermionic Hamiltonian (given in second quantization with creation/annihilation operators) to a qubit (spin) Hamiltonian, then optimizes circuit parameters to minimize the expectation value of the Hamiltonian. In code below we implement this fully in NumPy/SciPy: building the Hamiltonian via the Jordan–Wigner transform, preparing a trial state with a simple rotation-and-entangle ansatz, and using a classical optimizer to find the lowest energy. |
| 6 | + |
| 7 | + |
| 8 | +Jordan–Wigner Mapping of Fermionic Operators |
| 9 | + |
| 10 | + |
| 11 | +Fermionic operators obey anticommutation relations. The Jordan–Wigner (JW) transform maps fermionic creation/annihilation operators to tensor products of Pauli matrices by introducing a string of $Z$ operators on all preceding qubits . Concretely, in an $n$-qubit register (qubits labeled $0,\dots,n-1$), the annihilation operator $a_j$ on mode $j$ is represented as |
| 12 | +a_j \;\to\; (Z_0 Z_1 \cdots Z_{j-1}) \; \frac{X_j + iY_j}{2}, |
| 13 | +and similarly the creation operator $a_j^\dagger$ is represented by $(X_j - iY_j)/2$ with the same $Z$-chain . Here $X,Y,Z$ are the Pauli matrices: |
| 14 | +# Pauli matrices (2x2 identity and Pauli X,Y,Z) |
| 15 | +I = np.array([[1,0],[0,1]], dtype=complex) |
| 16 | +X = np.array([[0,1],[1,0]], dtype=complex) |
| 17 | +Y = np.array([[0,-1j],[1j,0]], dtype=complex) |
| 18 | +Z = np.array([[1,0],[0,-1]], dtype=complex) |
| 19 | +These definitions match the standard Pauli matrices . We can now build the full (dense) matrix for each fermionic operator. For example, the creation operator on qubit j in an n-qubit system is: |
| 20 | +def creation(n, j): |
| 21 | + """ |
| 22 | + Return the 2^n x 2^n matrix for the fermionic creation operator a_j^dagger |
| 23 | + on mode j, using the Jordan-Wigner transform. |
| 24 | + """ |
| 25 | + # Single-qubit raising operator (|1><0|) = (X - iY)/2 |
| 26 | + a_dag_j = (X - 1j * Y) / 2.0 |
| 27 | + op = 1 # start as 1x1 identity |
| 28 | + for k in range(n): |
| 29 | + if k < j: |
| 30 | + op = np.kron(op, Z) # apply Z on qubits 0..j-1 |
| 31 | + elif k == j: |
| 32 | + op = np.kron(op, a_dag_j) # apply (X - iY)/2 on qubit j |
| 33 | + else: |
| 34 | + op = np.kron(op, I) # apply identity on qubits j+1..n-1 |
| 35 | + return op |
| 36 | + |
| 37 | +def annihilation(n, j): |
| 38 | + """ |
| 39 | + Return the 2^n x 2^n matrix for the fermionic annihilation operator a_j on mode j. |
| 40 | + """ |
| 41 | + # Single-qubit lowering operator (|0><1|) = (X + iY)/2 |
| 42 | + a_j = (X + 1j * Y) / 2.0 |
| 43 | + op = 1 |
| 44 | + for k in range(n): |
| 45 | + if k < j: |
| 46 | + op = np.kron(op, Z) # Z on qubits 0..j-1 |
| 47 | + elif k == j: |
| 48 | + op = np.kron(op, a_j) # (X + iY)/2 on qubit j |
| 49 | + else: |
| 50 | + op = np.kron(op, I) # Identity on remaining qubits |
| 51 | + return op |
| 52 | +These routines construct the JW-mapped operators by tensoring $Z$ or $I$ on each qubit. One can verify they satisfy the required anticommutation relations (e.g. ${a_j,a_k^\dagger}=\delta_{jk}$) by inspection. |
| 53 | + |
| 54 | + |
| 55 | +Constructing the Qubit Hamiltonian |
| 56 | + |
| 57 | + |
| 58 | +Given a fermionic Hamiltonian specified by one-body coefficients $h_{pq}$ and two-body coefficients $h_{pqrs}$ (in second-quantized form $H = \sum_{pq}h_{pq} a_p^\dagger a_q + \frac12\sum_{pqrs}h_{pqrs}a_p^\dagger a_q^\dagger a_r a_s$), we build the full qubit Hamiltonian matrix in the computational basis. We sum over all terms, converting each product of fermionic operators into matrix form via the JW mapping: |
| 59 | +def build_qubit_hamiltonian(h1, h2): |
| 60 | + """ |
| 61 | + Build the full 2^n x 2^n qubit Hamiltonian matrix for a fermionic Hamiltonian |
| 62 | + with one-body terms h1[p,q] and two-body terms h2[p,q,r,s]. |
| 63 | + - h1 is an (n x n) matrix of coefficients for a_p^\dagger a_q. |
| 64 | + - h2 is an (n x n x n x n) tensor for a_p^\dagger a_q^\dagger a_r a_s. |
| 65 | + Returns: (H, n), where H is the Hamiltonian matrix and n is number of qubits. |
| 66 | + """ |
| 67 | + n = h1.shape[0] |
| 68 | + dim = 2**n |
| 69 | + H = np.zeros((dim, dim), dtype=complex) |
| 70 | + # One-body terms |
| 71 | + for p in range(n): |
| 72 | + for q in range(n): |
| 73 | + coeff = h1[p, q] |
| 74 | + if abs(coeff) > 1e-12: |
| 75 | + H += coeff * (creation(n,p) @ annihilation(n,q)) |
| 76 | + # Two-body terms |
| 77 | + for p in range(n): |
| 78 | + for q in range(n): |
| 79 | + for r in range(n): |
| 80 | + for s in range(n): |
| 81 | + coeff = h2[p, q, r, s] |
| 82 | + if abs(coeff) > 1e-12: |
| 83 | + term = (creation(n,p) @ creation(n,q) @ |
| 84 | + annihilation(n,r) @ annihilation(n,s)) |
| 85 | + H += coeff * term |
| 86 | + # Hermitian ensure (in case inputs weren't Hermitian) |
| 87 | + return (H + H.conj().T) / 2, n |
| 88 | +This function computes each matrix term by explicit matrix multiplication (@). The result H is a dense NumPy array of size $2^n\times 2^n$ acting on the $n$ qubits. (In practice this becomes large for $n>10$, but it meets the requirements.) |
| 89 | + |
| 90 | + |
| 91 | +Variational Ansatz (RY Rotations and Entanglement) |
| 92 | + |
| 93 | + |
| 94 | +We use a simple hardware-efficient ansatz: each qubit undergoes a parameterized $R_y(\theta)$ rotation, then a layer of CNOTs entangles qubit $j$ with $j+1$, and then a second layer of $R_y$ rotations on each qubit. (This kind of layered ansatz is common and can express a wide range of states.) The single-qubit rotation $R_y(\theta)$ has matrix |
| 95 | +R_y(\theta) = \begin{pmatrix} \cos(\theta/2) & -\sin(\theta/2) \\ \sin(\theta/2) & \cos(\theta/2) \end{pmatrix}, |
| 96 | +mixing $|0\rangle$ and $|1\rangle$ amplitudes . |
| 97 | + |
| 98 | +We implement this ansatz by manipulating the statevector directly. The helper functions below apply an $R_y$ gate or a CNOT gate to a state vector (length $2^n$). We index qubits so that qubit 0 is the most-significant bit of the basis index. |
| 99 | +import math |
| 100 | + |
| 101 | +def apply_RY(state, n, qubit, theta): |
| 102 | + """ |
| 103 | + Apply RY(theta) rotation on the specified qubit (0 = MSB) to the state vector. |
| 104 | + """ |
| 105 | + c = math.cos(theta/2) |
| 106 | + s = math.sin(theta/2) |
| 107 | + new_state = np.zeros_like(state, dtype=complex) |
| 108 | + for i in range(2**n): |
| 109 | + bit = (i >> (n-1-qubit)) & 1 |
| 110 | + j = i ^ (1 << (n-1-qubit)) # flip the target qubit bit |
| 111 | + if bit == 0: |
| 112 | + new_state[i] += c * state[i] - s * state[j] |
| 113 | + else: |
| 114 | + new_state[i] += s * state[j] + c * state[i] |
| 115 | + return new_state |
| 116 | + |
| 117 | +def apply_CNOT(state, n, control, target): |
| 118 | + """ |
| 119 | + Apply a CNOT with the given control and target qubit (0 = MSB) to the state vector. |
| 120 | + """ |
| 121 | + new_state = np.zeros_like(state, dtype=complex) |
| 122 | + for i in range(2**n): |
| 123 | + cbit = (i >> (n-1-control)) & 1 |
| 124 | + if cbit == 0: |
| 125 | + new_state[i] += state[i] # control=0: state unchanged |
| 126 | + else: |
| 127 | + # control=1: flip the target bit |
| 128 | + j = i ^ (1 << (n-1-target)) |
| 129 | + new_state[j] += state[i] |
| 130 | + return new_state |
| 131 | + |
| 132 | +def ansatz_state(theta, n): |
| 133 | + """ |
| 134 | + Prepare the ansatz state for n qubits given 2n parameters: |
| 135 | + First n parameters for RY rotations on each qubit, |
| 136 | + then entangling CNOTs (chain from qubit j to j+1), |
| 137 | + then another n RY rotations. |
| 138 | + Returns the 2^n statevector. |
| 139 | + """ |
| 140 | + assert len(theta) == 2*n |
| 141 | + # Start in the |00...0> state |
| 142 | + state = np.zeros(2**n, dtype=complex) |
| 143 | + state[0] = 1.0 |
| 144 | + # First layer of RY on all qubits |
| 145 | + for j in range(n): |
| 146 | + state = apply_RY(state, n, j, theta[j]) |
| 147 | + # Entangling layer: chain of CNOTs 0->1, 1->2, ..., n-2->n-1 |
| 148 | + for j in range(n-1): |
| 149 | + state = apply_CNOT(state, n, j, j+1) |
| 150 | + # Second layer of RY on all qubits |
| 151 | + for j in range(n): |
| 152 | + state = apply_RY(state, n, j, theta[n + j]) |
| 153 | + return state |
| 154 | +This ansatz uses $2n$ real parameters (each qubit has two $R_y$ angles, one before and one after entanglement). It is flexible enough to prepare various entangled states. |
| 155 | + |
| 156 | + |
| 157 | +Energy Expectation and Optimization |
| 158 | + |
| 159 | + |
| 160 | +To find the ground state energy, we compute the expectation value $\langle\psi(\theta)|H|\psi(\theta)\rangle$ of the Hamiltonian with respect to the ansatz state $|\psi(\theta)\rangle$, and then minimize it over the parameters $\theta$. Using the statevector, the expectation value is |
| 161 | +def energy_expectation(theta, H, n): |
| 162 | + """ |
| 163 | + Compute the expectation value <psi(theta)| H |psi(theta)> |
| 164 | + for the ansatz state on n qubits. |
| 165 | + """ |
| 166 | + psi = ansatz_state(theta, n) |
| 167 | + return np.real(np.vdot(psi, H @ psi)) |
| 168 | +We then call a classical optimizer from SciPy to minimize this expectation value. For example, using BFGS: |
| 169 | +from scipy.optimize import minimize |
| 170 | + |
| 171 | +# Example setup: define one- and two-body Hamiltonian coefficients h1, h2 (NumPy arrays) |
| 172 | +# h1 = ... (n x n), h2 = ... (n x n x n x n) |
| 173 | +H_qubit, n = build_qubit_hamiltonian(h1, h2) |
| 174 | + |
| 175 | +# Random initial guess for the 2n parameters |
| 176 | +initial_theta = np.random.rand(2*n) * np.pi |
| 177 | + |
| 178 | +result = minimize(lambda th: energy_expectation(th, H_qubit, n), |
| 179 | + initial_theta, method='BFGS') |
| 180 | +ground_energy = result.fun |
| 181 | +optimal_params = result.x |
| 182 | +print("Estimated ground-state energy:", ground_energy) |
| 183 | +This optimization finds (hopefully) the parameter set that minimizes the energy. In simple tests (small $n$), it should converge to the exact lowest eigenvalue of the Hamiltonian. |
| 184 | + |
| 185 | + |
| 186 | +Example Usage |
| 187 | + |
| 188 | + |
| 189 | +As a quick sanity check, consider a 2-qubit system ($n=2$) with a trivial Hamiltonian. For instance, let only orbital 0 have energy 1 and no two-body terms: $h_{00}=1, h_{11}=0, h_{pq}=0$ otherwise. The Hamiltonian is $a_0^\dagger a_0$, whose ground energy is 0 (vacuum state). Our code correctly finds this: |
| 190 | +n = 2 |
| 191 | +h1 = np.array([[1.0, 0.0], |
| 192 | + [0.0, 0.0]]) |
| 193 | +h2 = np.zeros((n,n,n,n)) |
| 194 | +H_qubit, _ = build_qubit_hamiltonian(h1, h2) |
| 195 | + |
| 196 | +# Optimize |
| 197 | +opt = minimize(lambda th: energy_expectation(th, H_qubit, n), |
| 198 | + np.zeros(2*n), method='BFGS') |
| 199 | +print("Computed ground energy ≈", opt.fun) |
| 200 | +# Output: Computed ground energy ≈ 0.0 |
| 201 | +This demonstrates the VQE routine finds (within numerical tolerance) the true ground energy. More complex Hamiltonians (including nonzero two-body terms) can be handled similarly by providing the appropriate h1 and h2 arrays. |
| 202 | + |
| 203 | +References: The Jordan–Wigner mapping and VQE approach are standard in quantum computation. See Nielsen & Chuang or Nielsen’s notes on JW transform and general VQE descriptions . Our code follows these established formulas to build and optimize the variational energy. |
0 commit comments