|
| 1 | +""" |
| 2 | +Reproduction of "A density-matrix renormalization group algorithm for simulating |
| 3 | +quantum circuits with a finite fidelity" |
| 4 | +Link: https://arxiv.org/abs/2207.05612 |
| 5 | +
|
| 6 | +Description: |
| 7 | +This script reproduces Figure 2(a) from the paper. |
| 8 | +It simulates a Sycamore-like random quantum circuit using both exact state vector simulation |
| 9 | +and an MPS-based simulator (DMRG-like algorithm) with varying bond dimensions. |
| 10 | +The script plots the infidelity (1 - Fidelity) as a function of the bond dimension. |
| 11 | +""" |
| 12 | + |
| 13 | +import time |
| 14 | +import logging |
| 15 | +import numpy as np |
| 16 | +import matplotlib.pyplot as plt |
| 17 | +import tensorcircuit as tc |
| 18 | + |
| 19 | +# Configure logging |
| 20 | +logging.basicConfig(level=logging.INFO) |
| 21 | +logger = logging.getLogger(__name__) |
| 22 | + |
| 23 | +# Use numpy backend for broad compatibility |
| 24 | +K = tc.set_backend("numpy") |
| 25 | + |
| 26 | + |
| 27 | +def generate_sycamore_like_circuit(rows, cols, depth, seed=42): |
| 28 | + """ |
| 29 | + Generates a random quantum circuit with a structure similar to Sycamore circuits. |
| 30 | + Since real Sycamore gates are specific, we use a simplified model: |
| 31 | + - 2D Grid connectivity |
| 32 | + - Layers of random single-qubit gates |
| 33 | + - Layers of two-qubit gates (CZ or similar) in a specific pattern |
| 34 | + """ |
| 35 | + np.random.seed(seed) |
| 36 | + n_qubits = rows * cols |
| 37 | + c = tc.Circuit(n_qubits) |
| 38 | + |
| 39 | + def q(r, col): |
| 40 | + return r * cols + col |
| 41 | + |
| 42 | + for d in range(depth): |
| 43 | + # Single qubit gates |
| 44 | + for i in range(n_qubits): |
| 45 | + # Random single qubit gate (e.g., Rx, Ry, Rz) |
| 46 | + # Simplification: Use a general random unitary or specific rotations |
| 47 | + # Sycamore uses sqrt(X), sqrt(Y), sqrt(W) |
| 48 | + # We'll use random rotations for generic RQC behavior |
| 49 | + theta = np.random.uniform(0, 2 * np.pi) |
| 50 | + phi = np.random.uniform(0, 2 * np.pi) |
| 51 | + lam = np.random.uniform(0, 2 * np.pi) |
| 52 | + c.rz(i, theta=phi) |
| 53 | + c.ry(i, theta=theta) |
| 54 | + c.rz(i, theta=lam) |
| 55 | + |
| 56 | + # Two-qubit gates |
| 57 | + # Sycamore uses a specific pattern (ABCD...) |
| 58 | + # We will use a simple alternating pattern |
| 59 | + # Layer A: horizontal even |
| 60 | + # Layer B: horizontal odd |
| 61 | + # Layer C: vertical even |
| 62 | + # Layer D: vertical odd |
| 63 | + # We cycle through these patterns based on depth |
| 64 | + |
| 65 | + layer_type = d % 4 |
| 66 | + |
| 67 | + if layer_type == 0: # Horizontal (col, col+1) for even cols |
| 68 | + for r in range(rows): |
| 69 | + for col in range(0, cols - 1, 2): |
| 70 | + c.cz(q(r, col), q(r, col + 1)) |
| 71 | + elif layer_type == 1: # Horizontal (col, col+1) for odd cols |
| 72 | + for r in range(rows): |
| 73 | + for col in range(1, cols - 1, 2): |
| 74 | + c.cz(q(r, col), q(r, col + 1)) |
| 75 | + elif layer_type == 2: # Vertical (row, row+1) for even rows |
| 76 | + for col in range(cols): |
| 77 | + for r in range(0, rows - 1, 2): |
| 78 | + c.cz(q(r, col), q(r + 1, col)) |
| 79 | + elif layer_type == 3: # Vertical (row, row+1) for odd rows |
| 80 | + for col in range(cols): |
| 81 | + for r in range(1, rows - 1, 2): |
| 82 | + c.cz(q(r, col), q(r + 1, col)) |
| 83 | + |
| 84 | + return c |
| 85 | + |
| 86 | + |
| 87 | +def run_mps_simulation(c, bond_dim): |
| 88 | + """ |
| 89 | + Runs the simulation using MPSCircuit with a maximum bond dimension. |
| 90 | + """ |
| 91 | + # Create MPSCircuit from the circuit operations |
| 92 | + # We need to re-apply the gates to an MPSCircuit |
| 93 | + # Or we can just build the MPSCircuit directly in the generation function. |
| 94 | + # But to ensure exactly the same circuit, we can iterate over the qir of the tc.Circuit |
| 95 | + n = c._nqubits |
| 96 | + mps = tc.MPSCircuit(n) |
| 97 | + |
| 98 | + # Set truncation rules |
| 99 | + # We use `max_singular_values` to control the bond dimension (chi) |
| 100 | + mps.set_split_rules({"max_singular_values": bond_dim}) |
| 101 | + |
| 102 | + for gate in c._qir: |
| 103 | + index = gate["index"] |
| 104 | + params = gate.get("parameters", {}) |
| 105 | + |
| 106 | + # Construct the gate on MPS |
| 107 | + # tc.Circuit stores 'gatef' which returns a Gate object when called with params |
| 108 | + g_obj = gate["gatef"](**params) |
| 109 | + |
| 110 | + # Apply to MPS |
| 111 | + mps.apply(g_obj, *index) |
| 112 | + |
| 113 | + return mps |
| 114 | + |
| 115 | + |
| 116 | +def calculate_fidelity(exact_c, mps_c): |
| 117 | + """ |
| 118 | + Calculates the fidelity between the exact state and the MPS state. |
| 119 | + F = |<psi_exact | psi_mps>|^2 |
| 120 | + """ |
| 121 | + # Get exact state vector |
| 122 | + psi_exact = exact_c.state() |
| 123 | + |
| 124 | + # Get MPS state vector (converted to full tensor) |
| 125 | + # Note: For large N, this will OOM. We should keep N small (e.g. <= 20). |
| 126 | + psi_mps = mps_c.wavefunction() |
| 127 | + |
| 128 | + # Compute overlap |
| 129 | + # exact state is (2^N,) or (1, 2^N) |
| 130 | + # mps state is (1, 2^N) usually |
| 131 | + |
| 132 | + psi_exact = K.reshape(psi_exact, (-1,)) |
| 133 | + psi_mps = K.reshape(psi_mps, (-1,)) |
| 134 | + |
| 135 | + overlap = K.tensordot(K.conj(psi_exact), psi_mps, axes=1) |
| 136 | + fidelity = np.abs(overlap) ** 2 |
| 137 | + return float(fidelity) |
| 138 | + |
| 139 | + |
| 140 | +def main(): |
| 141 | + # Parameters |
| 142 | + ROWS = 3 |
| 143 | + COLS = 4 # 12 qubits |
| 144 | + DEPTH = 8 |
| 145 | + # Bond dimensions to sweep |
| 146 | + # For 12 qubits, full bond dimension is 2^6=64. |
| 147 | + # So we should see perfect fidelity at 64, and error below it. |
| 148 | + BOND_DIMS = [2, 4, 8, 16, 32, 64] |
| 149 | + |
| 150 | + logger.info(f"Generating random circuit: {ROWS}x{COLS} grid, Depth {DEPTH}") |
| 151 | + circuit = generate_sycamore_like_circuit(ROWS, COLS, DEPTH, seed=42) |
| 152 | + |
| 153 | + # 1. Exact Simulation |
| 154 | + logger.info("Running Exact Simulation...") |
| 155 | + start_time = time.time() |
| 156 | + # Force state calculation |
| 157 | + _ = circuit.state() |
| 158 | + logger.info(f"Exact simulation done in {time.time() - start_time:.4f}s") |
| 159 | + |
| 160 | + infidelities = [] |
| 161 | + |
| 162 | + # 2. MPS Simulation with varying bond dimension |
| 163 | + logger.info("Running MPS Simulations...") |
| 164 | + for chi in BOND_DIMS: |
| 165 | + start_time = time.time() |
| 166 | + mps = run_mps_simulation(circuit, chi) |
| 167 | + |
| 168 | + # Calculate Fidelity |
| 169 | + fid = calculate_fidelity(circuit, mps) |
| 170 | + infidelity = 1.0 - fid |
| 171 | + # Avoid log(0) |
| 172 | + if infidelity < 1e-15: |
| 173 | + infidelity = 1e-15 |
| 174 | + |
| 175 | + infidelities.append(infidelity) |
| 176 | + |
| 177 | + logger.info( |
| 178 | + f"Bond Dim: {chi}, Fidelity: {fid:.6f}, Infidelity: {infidelity:.4e}, Time: {time.time() - start_time:.4f}s" |
| 179 | + ) |
| 180 | + |
| 181 | + # 3. Plotting |
| 182 | + plt.figure(figsize=(8, 6)) |
| 183 | + plt.loglog(BOND_DIMS, infidelities, "o-", label="Total Infidelity (1-F)") |
| 184 | + |
| 185 | + plt.xlabel("Bond Dimension (chi)") |
| 186 | + plt.ylabel("Infidelity (1 - F)") |
| 187 | + plt.title( |
| 188 | + f"MPS Simulation Accuracy vs Bond Dimension\n{ROWS}x{COLS} Circuit, Depth {DEPTH}" |
| 189 | + ) |
| 190 | + plt.grid(True, which="both", ls="--") |
| 191 | + plt.legend() |
| 192 | + |
| 193 | + output_path = ( |
| 194 | + "examples/reproduce_papers/2022_dmrg_circuit_simulation/outputs/result.png" |
| 195 | + ) |
| 196 | + plt.savefig(output_path) |
| 197 | + logger.info(f"Plot saved to {output_path}") |
| 198 | + |
| 199 | + |
| 200 | +if __name__ == "__main__": |
| 201 | + main() |
0 commit comments