|
| 1 | + |
| 2 | +Variational Quantum Eigensolver (VQE) Overview |
| 3 | + |
| 4 | + |
| 5 | +The Variational Quantum Eigensolver (VQE) is a hybrid quantum-classical algorithm that estimates the ground-state energy of a Hamiltonian by optimizing a parameterized quantum state . In VQE one (1) prepares a trial state $|\psi(\vec\theta)\rangle$ using a parameterized ansatz circuit, (2) computes the energy expectation $E(\vec\theta)=\langle\psi(\vec\theta)|H|\psi(\vec\theta)\rangle$, and (3) uses a classical optimizer to adjust the parameters $\vec\theta$ to minimize $E$ . By the variational principle $E(\vec\theta)\ge E_0$ (the true ground-state energy), so minimizing $E$ yields an estimate of the ground energy . In practice VQE proceeds by iterating these steps: |
| 6 | + |
| 7 | +Hamiltonian encoding: Map the problem Hamiltonian into a qubit Pauli form. |
| 8 | +Ansatz preparation: Choose a simple trial state, e.g. single-qubit rotations and entangling gates. |
| 9 | +Energy measurement: Compute $E(\vec\theta)=\langle\psi(\vec\theta)|H|\psi(\vec\theta)\rangle$ (classically by simulation here). |
| 10 | +Classical optimization: Use a routine like scipy.optimize.minimize to update $\vec\theta$ until convergence . |
| 11 | + |
| 12 | + |
| 13 | +Each of these steps is implemented below in Python using only NumPy and SciPy (no specialized quantum libraries). |
| 14 | + |
| 15 | + |
| 16 | +2D Ising and Heisenberg Models |
| 17 | + |
| 18 | + |
| 19 | +The 2D lattice Ising and Heisenberg models describe spins on a square lattice with nearest-neighbor interactions. For the transverse-field Ising model the Hamiltonian is |
| 20 | +H_{\rm Ising} = -J\sum_{\langle i,j\rangle} Z_iZ_j \;-\; h\sum_i X_i, |
| 21 | +where $Z_i,X_i$ are Pauli operators on site $i$, $J$ is the spin–spin coupling, $h$ is the transverse field, and $\langle i,j\rangle$ runs over nearest neighbors . In the isotropic Heisenberg model the Hamiltonian is |
| 22 | +H_{\rm Heis} = J\sum_{\langle i,j\rangle}(X_iX_j + Y_iY_j + Z_iZ_j), |
| 23 | +with $J>0$ for antiferromagnetic coupling . In this implementation we build the full $2^N\times2^N$ Hamiltonian matrix (for $N=L\times L$ spins) by tensor products of Pauli matrices. Specifically, we sum $Z_iZ_j$ or $X_iX_j+Y_iY_j+Z_iZ_j$ over all horizontal and vertical neighbor pairs on an $L\times L$ lattice, and add the single-spin $X_i$ terms for the Ising model. |
| 24 | +import numpy as np |
| 25 | + |
| 26 | +# Define 1-qubit Pauli matrices |
| 27 | +I = np.array([[1,0],[0,1]], dtype=complex) |
| 28 | +X = np.array([[0,1],[1,0]], dtype=complex) |
| 29 | +Y = np.array([[0,-1j],[1j,0]], dtype=complex) |
| 30 | +Z = np.array([[1,0],[0,-1]], dtype=complex) |
| 31 | + |
| 32 | +def get_neighbors(L): |
| 33 | + """List of nearest-neighbor pairs (i,j) for an LxL square lattice (open boundaries).""" |
| 34 | + neighbors = [] |
| 35 | + for i in range(L): |
| 36 | + for j in range(L): |
| 37 | + idx = i*L + j |
| 38 | + if j < L-1: |
| 39 | + neighbors.append((idx, idx+1)) |
| 40 | + if i < L-1: |
| 41 | + neighbors.append((idx, idx+L)) |
| 42 | + return neighbors |
| 43 | + |
| 44 | +def build_ising_hamiltonian(L, J, h): |
| 45 | + """Construct the 2^N x 2^N transverse-field Ising Hamiltonian for LxL spins.""" |
| 46 | + N = L*L |
| 47 | + dim = 2**N |
| 48 | + H = np.zeros((dim, dim), dtype=complex) |
| 49 | + for (i,j) in get_neighbors(L): |
| 50 | + # Build tensor-product for -J * Z_i Z_j |
| 51 | + ops = [Z if k==i or k==j else I for k in range(N)] |
| 52 | + term = ops[0] |
| 53 | + for op in ops[1:]: |
| 54 | + term = np.kron(term, op) |
| 55 | + H += -J * term |
| 56 | + for i in range(N): |
| 57 | + # -h * X_i term |
| 58 | + ops = [X if k==i else I for k in range(N)] |
| 59 | + term = ops[0] |
| 60 | + for op in ops[1:]: |
| 61 | + term = np.kron(term, op) |
| 62 | + H += -h * term |
| 63 | + return H |
| 64 | + |
| 65 | +def build_heisenberg_hamiltonian(L, J): |
| 66 | + """Construct the 2^N x 2^N isotropic Heisenberg Hamiltonian for LxL spins.""" |
| 67 | + N = L*L |
| 68 | + dim = 2**N |
| 69 | + H = np.zeros((dim, dim), dtype=complex) |
| 70 | + for (i,j) in get_neighbors(L): |
| 71 | + # Add J*(X_iX_j + Y_iY_j + Z_iZ_j) |
| 72 | + for P in (X, Y, Z): |
| 73 | + ops = [P if k==i or k==j else I for k in range(N)] |
| 74 | + term = ops[0] |
| 75 | + for op in ops[1:]: |
| 76 | + term = np.kron(term, op) |
| 77 | + H += J * term |
| 78 | + return H |
| 79 | +Each term above is built by a loop over qubit indices $k$, placing the appropriate Pauli matrix for sites $i,j$ or identity elsewhere, and taking the Kronecker product. This produces the exact Hamiltonian matrix for the small lattice. (For example, a $2\times2$ lattice has $N=4$ qubits and a $16\times16$ matrix.) |
| 80 | + |
| 81 | + |
| 82 | +Variational Ansatz and State Preparation |
| 83 | + |
| 84 | + |
| 85 | +We choose a simple hardware-efficient ansatz: start from the $|0\ldots0\rangle$ product state, apply parameterized single-qubit $R_y(\theta)$ rotations on each qubit, and then entangle neighboring qubits with CNOT gates. In vector form, an $R_y(\theta)$ rotation on qubit $q$ acts as |
| 86 | +R_y(\theta) = \begin{pmatrix} \cos(\theta/2) & -\sin(\theta/2)\\ \sin(\theta/2) & \cos(\theta/2) \end{pmatrix}, |
| 87 | +mixing the amplitudes of basis states differing on bit $q$. We implement these operations directly on the state vector. The functions below apply $R_y$ and CNOT to the $2^N$-dimensional state: |
| 88 | +def apply_RY(state, qubit, theta): |
| 89 | + """Apply a single-qubit RY(theta) rotation on the specified qubit of the state.""" |
| 90 | + N = int(np.log2(len(state))) |
| 91 | + new_state = np.zeros_like(state, dtype=complex) |
| 92 | + cos = np.cos(theta/2); sin = np.sin(theta/2) |
| 93 | + # Loop over basis states in binary; if qubit bit is 0, mix with its partner state |
| 94 | + for b in range(len(state)): |
| 95 | + if ((b >> qubit) & 1) == 0: |
| 96 | + partner = b | (1 << qubit) # flip that qubit bit |
| 97 | + new_state[b] += cos * state[b] - sin * state[partner] |
| 98 | + new_state[partner] += sin * state[b] + cos * state[partner] |
| 99 | + return new_state |
| 100 | + |
| 101 | +def apply_CNOT(state, control, target): |
| 102 | + """Apply a CNOT gate (control -> target) on the specified qubits of the state.""" |
| 103 | + new_state = np.zeros_like(state, dtype=complex) |
| 104 | + for b in range(len(state)): |
| 105 | + # If control bit is 1, flip the target bit |
| 106 | + if ((b >> control) & 1) == 1: |
| 107 | + b_new = b ^ (1 << target) |
| 108 | + else: |
| 109 | + b_new = b |
| 110 | + new_state[b_new] += state[b] |
| 111 | + return new_state |
| 112 | + |
| 113 | +def ansatz_state(params, L): |
| 114 | + """Construct the variational ansatz state for parameters `params` on an LxL lattice.""" |
| 115 | + N = L*L |
| 116 | + psi = np.zeros(2**N, dtype=complex) |
| 117 | + psi[0] = 1.0 # start in |00...0> |
| 118 | + # Apply RY on each qubit with corresponding parameter |
| 119 | + for i, theta in enumerate(params): |
| 120 | + psi = apply_RY(psi, i, theta) |
| 121 | + # Apply CNOT entanglers on each neighbor pair |
| 122 | + for (i,j) in get_neighbors(L): |
| 123 | + psi = apply_CNOT(psi, i, j) |
| 124 | + return psi |
| 125 | +Here params is an array of length $N=L^2$ (one parameter per qubit). After applying all $R_y$ rotations and CNOTs, ansatz_state(params,L) returns a state vector $|\psi(\vec\theta)\rangle$. This state incorporates both single-qubit rotations and two-qubit entanglement among neighbors. |
| 126 | + |
| 127 | + |
| 128 | +Energy Expectation and Optimization |
| 129 | + |
| 130 | + |
| 131 | +With the state prepared, the cost (energy) is simply $E(\vec\theta)=\langle\psi(\vec\theta)|H|\psi(\vec\theta)\rangle$. We implement this by a vector-matrix product using NumPy. A classical optimizer (e.g. SciPy’s minimize) then updates the parameters. For example: |
| 132 | +import scipy.optimize as opt |
| 133 | + |
| 134 | +def energy_expectation(params, H, L): |
| 135 | + """Compute expectation <psi(params)|H|psi(params)> for Hamiltonian H.""" |
| 136 | + psi = ansatz_state(params, L) |
| 137 | + # <psi|H|psi> using complex conjugate transpose |
| 138 | + return np.vdot(psi, H.dot(psi)).real |
| 139 | + |
| 140 | +# Example: find ground energy of 2x2 transverse Ising (J=1.0, h=0.5) |
| 141 | +L = 2 |
| 142 | +J = 1.0 |
| 143 | +h = 0.5 |
| 144 | +H_ising = build_ising_hamiltonian(L, J, h) |
| 145 | +# Initial random parameters |
| 146 | +np.random.seed(0) |
| 147 | +init_params = 0.1 * np.random.randn(L*L) |
| 148 | +# Optimize parameters to minimize energy |
| 149 | +result = opt.minimize(lambda th: energy_expectation(th, H_ising, L), |
| 150 | + init_params, method='COBYLA') |
| 151 | +estimated_energy = result.fun |
| 152 | +In this code, we pass the Hamiltonian matrix H_ising and lattice size L to the cost function. SciPy’s minimize routine (here using COBYLA) searches over the parameter vector th. The variational principle guarantees $E(\vec\theta)\ge E_0$, so the optimizer drives the estimate down toward the true ground energy . (One could also compute the exact ground-state energy by diagonalizing H_ising for comparison, since the matrix is small.) |
| 153 | + |
| 154 | + |
| 155 | +Results and Visualization |
| 156 | + |
| 157 | + |
| 158 | +Running the above VQE code for a small lattice (e.g. $2\times2$ or $3\times3$) yields an estimate of the ground-state energy as a function of the coupling parameters. For example, one can sweep the transverse field $h$ (for the Ising model) or coupling $J$ (for the Heisenberg model) and plot the converged energy. The curves typically follow the expected physical behavior (e.g. increasing field raises the ground energy). The code below sketches how to collect and plot these energies (not embedded here). |
| 159 | +import matplotlib.pyplot as plt |
| 160 | + |
| 161 | +# Example: Vary h for 2x2 Ising, J=1.0 |
| 162 | +hs = np.linspace(0.0, 2.0, 5) |
| 163 | +energies_vqe = [] |
| 164 | +for h_val in hs: |
| 165 | + H2 = build_ising_hamiltonian(2, 1.0, h_val) |
| 166 | + res = opt.minimize(lambda th: energy_expectation(th, H2, 2), |
| 167 | + init_params, method='COBYLA') |
| 168 | + energies_vqe.append(res.fun) |
| 169 | + |
| 170 | +plt.plot(hs, energies_vqe, marker='o') |
| 171 | +plt.xlabel('Transverse field h') |
| 172 | +plt.ylabel('Ground-state energy (VQE)') |
| 173 | +plt.title('2x2 Transverse-Field Ising (J=1)') |
| 174 | +plt.show() |
| 175 | +This would produce a plot of estimated ground-state energy versus $h$. A similar loop over $J$ can be done for the Heisenberg model. (In this example one can also compute the exact diagonalization energies for reference.) |
| 176 | + |
| 177 | +Modularity and Documentation: The code above is organized into functions for clarity: separate routines build the Hamiltonian, apply gates, form the ansatz state, and compute the cost. Each function is documented with comments. The overall procedure follows the standard VQE workflow and uses only NumPy/SciPy linear algebra, as required. |
| 178 | + |
| 179 | +Sources: The VQE method and variational principle are described in the quantum computing literature . The Hamiltonians for the transverse-field Ising and isotropic Heisenberg models are standard and given, e.g., in texts or reviews . SciPy’s optimization routine is documented in its manual . |
0 commit comments