Concepts Intermediate Free 44/53 in series 30 min

Quantum Walks and the Search Algorithm

How quantum walks on graphs achieve quadratic speedup for search problems, with the Szegedy framework explained and a Qiskit implementation of discrete-time quantum walk search.

What you'll learn

  • quantum walks
  • Szegedy walk
  • element distinctness
  • graph search
  • quadratic speedup

Prerequisites

  • Python proficiency
  • Beginner quantum computing concepts (superposition, entanglement)
  • Linear algebra basics

Classical Random Walks: Concrete Analysis

A classical random walk on a graph starts at some vertex and, at each step, moves to a uniformly random neighbor. Several quantities describe how the walk explores the graph. Understanding these quantities precisely is essential because quantum walks improve on each of them.

Hitting time. The hitting time h(u, v) is the expected number of steps to reach vertex v starting from vertex u. For a random walk on a line (path graph) with N vertices, the hitting time from vertex 0 to vertex N-1 is O(N^2). This quadratic scaling comes from the “drunkard’s walk” effect: the walker drifts back toward its origin as often as it moves forward.

Cover time. The cover time is the expected number of steps to visit every vertex at least once. For a line with N vertices, the cover time is O(N^2). For a 2D grid with N = n x n vertices, the cover time is O(N log N) = O(n^2 log(n^2)). The logarithmic factor in 2D arises because the last few unvisited vertices are hard to find by random exploration.

Return time. The expected return time to the starting vertex on a 1D line of N vertices with periodic boundaries is exactly N. On an N-vertex cycle, the stationary distribution assigns probability 1/N to each vertex, and by the Kac lemma the expected return time equals N.

Mixing time. The mixing time is the number of steps required for the walk’s distribution to converge (within total variation distance epsilon) to the stationary distribution. For a random walk on a graph with transition matrix P, the mixing time scales as 1/spectral_gap(P), where the spectral gap is 1 minus the second-largest eigenvalue magnitude of P.

For a line with N vertices (reflecting boundaries), the spectral gap of P is Theta(1/N^2), so the mixing time is Theta(N^2). For a 2D grid with N = n x n vertices, the spectral gap is Theta(1/n^2) = Theta(1/sqrt(N)^2), giving a mixing time of Theta(N).

Connection to linear algebra. The transition matrix P is a stochastic matrix with eigenvalues 1 = lambda_1 >= lambda_2 >= … >= lambda_N >= -1. The eigenvector for lambda_1 is the stationary distribution pi. The remaining eigenvectors describe the transient behavior: after t steps, the deviation from stationarity decays as |lambda_2|^t. The spectral gap delta = 1 - |lambda_2| controls how fast the walk forgets its starting position.

Random walks underlie PageRank, MCMC sampling, and classical search heuristics. A quantum walk replaces the random step with a unitary evolution, and interference between paths can dramatically speed up the spreading.

Discrete-Time Quantum Walks

In the discrete-time model, the walker has two registers: a coin register (representing which direction to move) and a position register (the current vertex). One step of the walk applies:

  1. A coin operator C on the coin register.
  2. A conditional shift operator S that moves the position according to the coin state.

The full step is U = S * (C tensor I_position). The coin operator creates superposition in the direction register, and the shift operator entangles the direction with position.

Coin choices

The Hadamard coin H = (1/sqrt(2)) * [[1,1],[1,-1]] is the standard choice for a walk on a line, where the coin has dimension 2 (left or right). But other coins produce different behavior.

Grover coin. For a d-dimensional coin (a walk on a graph where each vertex has degree d), the Grover coin is G = 2|s><s| - I, where |s> = (1/sqrt(d)) * (|0> + |1> + … + |d-1>). For d = 2, the Grover coin equals the Hadamard gate. For higher d, the Grover coin gives the fastest mixing on many regular graphs.

DFT coin. The quantum Fourier transform coin F_d, with entries F_{jk} = (1/sqrt(d)) * exp(2piijk/d), produces an asymmetric walk on the line. Unlike the Hadamard coin (which produces a symmetric distribution skewed slightly left), the DFT coin generates a distribution that moves predominantly in one direction.

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile
from qiskit_aer import AerSimulator
import numpy as np
import matplotlib.pyplot as plt

def quantum_walk_line(n_steps, n_position_qubits=6, coin_type="hadamard"):
    """
    Discrete-time quantum walk on a line (cyclic) with 2^n_position_qubits nodes.

    The position register uses periodic (cyclic) boundary conditions: incrementing
    past 2^n - 1 wraps to 0, and decrementing past 0 wraps to 2^n - 1.

    Parameters
    ----------
    n_steps : int
        Number of walk steps to apply.
    n_position_qubits : int
        Number of qubits for the position register. The line has 2^n_position_qubits vertices.
    coin_type : str
        "hadamard" for the standard Hadamard coin, "y" for the Y-rotation coin (balanced).
    """
    n = n_position_qubits
    coin = QuantumRegister(1, "coin")
    position = QuantumRegister(n, "pos")
    creg = ClassicalRegister(n, "c")
    qc = QuantumCircuit(coin, position, creg)

    # Start at the center of the line: position 2^(n-1)
    # This ensures the walk can spread in both directions without wrapping too early.
    qc.x(position[n - 1])

    # Initialize coin in |0> (no initial Hadamard).
    # For a symmetric walk, use the Y-rotation coin: R_y(pi/2).
    # The Hadamard coin produces a left-biased distribution when starting from |0>.

    for step in range(n_steps):
        # Apply coin operator
        if coin_type == "hadamard":
            qc.h(coin[0])
        elif coin_type == "y":
            # R_y(pi/2) coin produces a symmetric distribution
            qc.ry(np.pi / 2, coin[0])

        # Shift right when coin = |1>: controlled increment on position register.
        # Carry-ripple increment: flip bit i conditioned on coin=1 AND bits 0..i-1 all being 1.
        # Process from the most significant bit down to the least significant bit.
        for i in range(n - 1, 0, -1):
            # Multi-controlled X on position[i], controlled by coin[0] and position[0..i-1]
            controls = [coin[0]] + [position[j] for j in range(i)]
            qc.mcx(controls, position[i])
        qc.cx(coin[0], position[0])

        # Shift left when coin = |0>: controlled decrement on position register.
        # Decrement = flip coin polarity, increment, flip back.
        qc.x(coin[0])
        # Carry-ripple decrement: same structure as increment but with inverted coin.
        for i in range(n - 1, 0, -1):
            controls = [coin[0]] + [position[j] for j in range(i)]
            qc.mcx(controls, position[i])
        qc.cx(coin[0], position[0])
        qc.x(coin[0])

    qc.measure(position, creg)
    return qc

Why this shift operator works

The increment circuit is the trickiest part of the quantum walk implementation. It performs the mapping |x> -> |x+1 mod 2^n> on the position register, conditioned on the coin qubit being |1>.

Consider a 3-bit position register. To add 1 in binary:

  • Bit 0 (least significant) always flips.
  • Bit 1 flips if bit 0 was 1 (carry propagation).
  • Bit 2 flips if bits 0 and 1 were both 1 (carry propagation continues).

All these flips are conditioned on the coin qubit being |1>. So the gate sequence is:

  • MCX(coin, pos[0], pos[1]) controlled on coin and pos[0], targeting pos[2]
  • CX(coin, pos[0]) controlled on coin, targeting pos[1], but only if pos[0] = 1, so MCX(coin, pos[0], pos[1])
  • CX(coin, pos[0])

Wait, let us be precise. For the controlled increment on n position qubits, conditioned on coin = |1>:

  1. For bit n-1 (MSB): flip if coin=1 AND pos[0]=1 AND pos[1]=1 AND … AND pos[n-2]=1. This requires an (n)-controlled X gate.
  2. For bit n-2: flip if coin=1 AND pos[0]=1 AND … AND pos[n-3]=1. This requires an (n-1)-controlled X gate.
  3. Continue down to bit 0: flip if coin=1. This is a simple CX.

Processing from MSB to LSB is essential: if we flipped the LSB first, the carry conditions would read the already-modified bits.

Gate count per step. For n position qubits, the increment requires one CX gate and multi-controlled X gates with 2, 3, …, n controls. Each k-controlled X (Toffoli generalization) decomposes into O(k) Toffoli gates using ancillas, or O(k^2) Toffoli gates without ancillas. The decrement has the same cost. For a rough estimate without ancillas, each step costs about n*(n+1)/2 Toffoli-equivalent gates for the increment and the same for the decrement. For n = 5 position qubits, that is roughly 2 * 15 = 30 Toffoli-equivalent gates per step.

Ballistic vs Diffusive Spreading

The most striking difference between classical and quantum walks is how fast the walker spreads across the graph.

Classical random walk on a line: After t steps, the walker’s position follows a binomial distribution (approximately normal for large t). The standard deviation of the position grows as sqrt(t). This is diffusive spreading.

Quantum walk on a line: After t steps, the position distribution has a characteristic two-peaked shape, with peaks near positions +/- t/sqrt(2) from the origin. The standard deviation grows linearly as t/sqrt(2). This is ballistic spreading, quadratically faster than classical.

This difference is the root cause of quantum walk speedups. A quantum walker reaches distance d from its origin in O(d) steps, while a classical walker requires O(d^2) steps.

def compare_spreading():
    """
    Compare the spreading rate of classical vs quantum random walks on a line.
    Classical: standard deviation ~ sqrt(t).
    Quantum: standard deviation ~ t / sqrt(2).
    """
    sim = AerSimulator()
    steps_list = [5, 10, 15, 20, 25]
    n_pos = 7  # 128 positions, enough room for 25 steps from center
    shots = 8192

    quantum_stds = []
    classical_stds = []

    for t in steps_list:
        # Quantum walk
        qc = quantum_walk_line(t, n_position_qubits=n_pos, coin_type="hadamard")
        compiled = transpile(qc, sim)
        result = sim.run(compiled, shots=shots).result()
        counts = result.get_counts()

        center = 2 ** (n_pos - 1)  # starting position
        positions = []
        probs = []
        for bitstring, count in counts.items():
            pos = int(bitstring, 2)
            positions.append(pos - center)  # relative to start
            probs.append(count / shots)

        mean_q = sum(p * x for p, x in zip(probs, positions))
        var_q = sum(p * (x - mean_q) ** 2 for p, x in zip(probs, positions))
        quantum_stds.append(np.sqrt(var_q))

        # Classical random walk (exact computation)
        # After t steps with fair coin, position displacement follows
        # a distribution on {-t, -t+2, ..., t-2, t}.
        # Variance = t (each step contributes variance 1).
        classical_stds.append(np.sqrt(t))

    print("Steps | Quantum StdDev | Classical StdDev | Ratio (Q/C)")
    print("-" * 60)
    for i, t in enumerate(steps_list):
        ratio = quantum_stds[i] / classical_stds[i] if classical_stds[i] > 0 else 0
        print(f"  {t:3d} | {quantum_stds[i]:14.3f} | {classical_stds[i]:16.3f} | {ratio:.3f}")

    # The ratio should grow roughly as sqrt(t), confirming
    # quantum std ~ t while classical std ~ sqrt(t).
    print("\nQuantum std grows linearly with t (ballistic).")
    print("Classical std grows as sqrt(t) (diffusive).")

compare_spreading()

The output shows the quantum standard deviation growing roughly proportional to t, while the classical standard deviation grows as sqrt(t). The ratio quantum_std/classical_std increases with t, confirming the quadratic separation in spreading rate.

Common Mistakes in Quantum Walk Implementations

Before going further, here are pitfalls that catch many practitioners.

1. Boundary conditions matter. The shift operator above implements a cyclic shift: position wraps from 2^n - 1 back to 0 (and from 0 to 2^n - 1). This gives periodic boundary conditions. If you want reflective (hard-wall) boundaries instead, you must detect when the position is at 0 or 2^n - 1 and suppress the shift. Periodic and reflective boundaries produce very different distributions near the edges. For studying the bulk behavior of the walk, start the walker near the center of the line and use enough position qubits so the walker never reaches the boundary.

2. Coherent coin, not classical randomness. The Hadamard gate creates a coherent superposition of |0> and |1>. If you replaced it with a classical random bit (measure, then conditionally shift), you would get a classical random walk with diffusive spreading. The quantum speedup comes entirely from the coherent superposition creating interference.

3. Never measure between steps. Measuring the position after each step collapses the superposition and destroys all interference. The walker reverts to classical behavior. Run all t steps unitarily, then measure once at the end.

4. The Hadamard coin produces an asymmetric distribution. Starting with coin state |0> and applying the Hadamard coin gives a distribution skewed to the left (the |0> branch, which shifts left, has higher amplitude). For a symmetric distribution, either start the coin in (|0> - i|1>)/sqrt(2) or use the Y-rotation coin R_y(pi/2).

5. Grover search is not a quantum walk search. A common confusion: Grover’s algorithm operates on a flat, unstructured search space with no graph topology. Quantum walk search operates on a specific graph and exploits its structure. On some graphs (like the hypercube), the two coincide, but on others (like the 2D grid), quantum walk search achieves speedups that Grover’s algorithm cannot.

2D Quantum Walk on a Grid

A 2D quantum walk extends the 1D construction to a grid. The walker moves in four directions: up, down, left, right. This requires a 4-dimensional coin encoded in 2 qubits.

The position register splits into two parts: x-position and y-position. The shift operator moves:

  • Right (+x) when coin = |00>
  • Left (-x) when coin = |01>
  • Up (+y) when coin = |10>
  • Down (-y) when coin = |11>

For a 4x4 grid, each position coordinate uses 2 qubits (total: 2 coin + 2 x-position + 2 y-position = 6 qubits).

def quantum_walk_2d(n_steps, n_pos_qubits_per_dim=2):
    """
    Discrete-time quantum walk on a 2D grid with periodic boundary conditions.

    Grid size: 2^n x 2^n where n = n_pos_qubits_per_dim.
    Coin: 2 qubits encoding 4 directions.
    Coin operator: Grover diffusion on 2 qubits (= Hadamard on each, since for
    2 qubits the Grover coin G = 2|s><s| - I with |s> = |++> equals H tensor H).

    Actually, for a 4-dimensional coin, the Grover coin is:
    G = 2|s><s| - I where |s> = (|00> + |01> + |10> + |11>)/2.
    G has entries: G_ij = 2/4 - delta_ij = 1/2 - delta_ij.
    This is NOT the same as H tensor H.
    H tensor H has entries: (1/2)(-1)^(i.j), which is the Hadamard matrix.

    We use the Grover coin for the 2D walk, implemented via:
    G = H^{otimes 2} . diag(1,1,1,-1) . H^{otimes 2}
    Wait, that gives the wrong matrix. Let us use the direct construction.

    Actually, for a 2-qubit system:
    G = 2|++><++| - I = 2(H|0>)(H|0>)^dag tensor 2 ... this is simpler to
    implement as: H(q0), H(q1), then the 2-qubit Grover diffusion
    (X on both, CZ, X on both), then H(q0), H(q1).

    The 2-qubit Grover diffusion flips the phase of |00> and leaves others unchanged,
    then sandwiching with H gives the Grover coin.
    """
    n = n_pos_qubits_per_dim
    coin = QuantumRegister(2, "coin")
    pos_x = QuantumRegister(n, "x")
    pos_y = QuantumRegister(n, "y")
    creg_x = ClassicalRegister(n, "cx")
    creg_y = ClassicalRegister(n, "cy")
    qc = QuantumCircuit(coin, pos_x, pos_y, creg_x, creg_y)

    # Start at center of the grid
    if n > 1:
        qc.x(pos_x[n - 1])
        qc.x(pos_y[n - 1])

    for step in range(n_steps):
        # Grover coin on 2-qubit coin register
        # G = 2|s><s| - I where |s> = (|00>+|01>+|10>+|11>)/2
        # Implementation: H H, phase flip |00>, H H
        # Phase flip |00>: X on both, CZ, X on both. But CZ flips phase of |11>,
        # so X's convert |00> <-> |11>.
        qc.h(coin[0])
        qc.h(coin[1])
        qc.x(coin[0])
        qc.x(coin[1])
        qc.cz(coin[0], coin[1])
        qc.x(coin[0])
        qc.x(coin[1])
        qc.h(coin[0])
        qc.h(coin[1])

        # Shift in x-direction based on coin[0]
        # coin = |0.>: increment x (move right)
        # coin = |1.>: decrement x (move left)

        # Controlled increment on pos_x, conditioned on coin[0] = 0
        qc.x(coin[0])
        for i in range(n - 1, 0, -1):
            controls = [coin[0]] + [pos_x[j] for j in range(i)]
            qc.mcx(controls, pos_x[i])
        qc.cx(coin[0], pos_x[0])
        qc.x(coin[0])

        # Controlled decrement on pos_x, conditioned on coin[0] = 1
        for i in range(n - 1, 0, -1):
            controls = [coin[0]] + [pos_x[j] for j in range(i)]
            qc.mcx(controls, pos_x[i])
        qc.cx(coin[0], pos_x[0])

        # Shift in y-direction based on coin[1]
        # coin = |.0>: increment y (move up)
        # coin = |.1>: decrement y (move down)
        qc.x(coin[1])
        for i in range(n - 1, 0, -1):
            controls = [coin[1]] + [pos_y[j] for j in range(i)]
            qc.mcx(controls, pos_y[i])
        qc.cx(coin[1], pos_y[0])
        qc.x(coin[1])

        for i in range(n - 1, 0, -1):
            controls = [coin[1]] + [pos_y[j] for j in range(i)]
            qc.mcx(controls, pos_y[i])
        qc.cx(coin[1], pos_y[0])

    qc.measure(pos_x, creg_x)
    qc.measure(pos_y, creg_y)
    return qc

# Run a 4-step 2D quantum walk on a 4x4 grid
qc_2d = quantum_walk_2d(n_steps=4, n_pos_qubits_per_dim=2)
sim = AerSimulator()
result = sim.run(transpile(qc_2d, sim), shots=4096).result()
counts = result.get_counts()

print("2D quantum walk position distribution (x, y):")
grid = np.zeros((4, 4))
for bitstring, count in counts.items():
    # Qiskit returns bits in reverse register order
    # bitstring format: "cy1 cy0 cx1 cx0" (little-endian per register)
    bits = bitstring.replace(" ", "")
    # The classical registers are concatenated: cy bits then cx bits
    cx_bits = bits[2:]  # last n bits
    cy_bits = bits[:2]  # first n bits
    x = int(cx_bits, 2)
    y = int(cy_bits, 2)
    grid[y][x] += count

grid = grid / grid.sum()
for row in range(4):
    print("  " + "  ".join(f"{grid[row][col]:.3f}" for col in range(4)))

The 2D quantum walk distribution shows a characteristic diamond-shaped pattern, concentrated at positions roughly t steps from the origin along the diagonals. A classical 2D random walk would produce a circular Gaussian blob centered on the starting position, spreading as sqrt(t) in each dimension.

The Szegedy Walk Framework

The Szegedy walk generalizes discrete-time quantum walks to arbitrary graphs using a bipartite structure. This framework is the theoretical backbone of most quantum walk speedup results.

Setup. Start with a graph G on N vertices and its random walk transition matrix P, where P_{xy} = 1/deg(x) if (x, y) is an edge and 0 otherwise.

The Szegedy walk operates on a Hilbert space of dimension N^2, spanned by basis states |x>|y> for pairs of vertices. Define the superposition states:

|p_x> = sum_y sqrt(P_{xy}) |x>|y>

Each |p_x> is the “quantized” version of the classical random walk distribution from vertex x.

Reflections. The walk consists of two reflections:

R_A = 2 * (sum_x |p_x><p_x|) - I

R_B = SWAP * R_A * SWAP

Here SWAP exchanges the two registers: SWAP|x>|y> = |y>|x>. One step of the Szegedy walk is W = R_B * R_A.

Spectral connection. If the classical random walk has eigenvalues 1 = lambda_1 >= lambda_2 >= …, the Szegedy walk operator W has eigenvalues e^{+/- i*arccos(lambda_k)}. The eigenvalue gap of W is arccos(lambda_2), which is approximately sqrt(2 * (1 - lambda_2)) = sqrt(2 * delta) when the spectral gap delta is small. This square-root relationship is the origin of the quadratic speedup.

Worked example: Szegedy walk on a 4-vertex cycle

Consider the cycle graph C_4 with vertices {0, 1, 2, 3}. Each vertex connects to its two neighbors.

Transition matrix:

P = [[0,    1/2,  0,    1/2],
     [1/2,  0,    1/2,  0  ],
     [0,    1/2,  0,    1/2],
     [1/2,  0,    1/2,  0  ]]

Quantized states:

|p_0> = (1/sqrt(2))(|0>|1> + |0>|3>) |p_1> = (1/sqrt(2))(|1>|0> + |1>|2>) |p_2> = (1/sqrt(2))(|2>|1> + |2>|3>) |p_3> = (1/sqrt(2))(|3>|0> + |3>|2>)

Spectral analysis. The eigenvalues of P for C_4 are {1, 0, -1, 0}. The spectral gap is delta = 1 - |0| = 1. The Szegedy walk eigenvalues include e^{+/- iarccos(0)} = e^{+/- ipi/2} = {+i, -i} and e^{+/- iarccos(1)} = {1} and e^{+/- iarccos(-1)} = {-1}. The eigenvalue gap of W is pi/2.

def szegedy_walk_example():
    """
    Construct the Szegedy walk operator for C_4 (4-vertex cycle)
    and verify its eigenvalues.
    """
    # Transition matrix for C_4
    P = np.array([
        [0,   0.5, 0,   0.5],
        [0.5, 0,   0.5, 0  ],
        [0,   0.5, 0,   0.5],
        [0.5, 0,   0.5, 0  ]
    ])
    N = 4

    # Build the quantized states |p_x>
    # |p_x> = sum_y sqrt(P_xy) |x>|y>  in a N^2-dimensional space
    # Basis ordering: |0,0>, |0,1>, ..., |0,N-1>, |1,0>, ..., |N-1,N-1>
    p_states = np.zeros((N, N * N))
    for x in range(N):
        for y in range(N):
            if P[x, y] > 0:
                p_states[x, x * N + y] = np.sqrt(P[x, y])

    # R_A = 2 * sum_x |p_x><p_x| - I
    R_A = -np.eye(N * N)
    for x in range(N):
        px = p_states[x].reshape(-1, 1)
        R_A += 2 * (px @ px.T)

    # SWAP operator: SWAP|x,y> = |y,x>
    SWAP = np.zeros((N * N, N * N))
    for x in range(N):
        for y in range(N):
            SWAP[y * N + x, x * N + y] = 1

    # R_B = SWAP R_A SWAP
    R_B = SWAP @ R_A @ SWAP

    # Walk operator W = R_B R_A
    W = R_B @ R_A

    # Eigenvalues of W
    eigenvalues = np.linalg.eigvals(W)
    print("Szegedy walk eigenvalues for C_4:")
    for ev in sorted(eigenvalues, key=lambda z: -abs(z.real)):
        if abs(ev.imag) < 1e-10:
            print(f"  {ev.real:+.4f}")
        else:
            print(f"  {ev.real:+.4f} {ev.imag:+.4f}i")

    # Eigenvalues of P for comparison
    p_eigs = np.linalg.eigvals(P)
    print("\nClassical transition matrix eigenvalues:")
    for ev in sorted(p_eigs, key=lambda z: -z.real):
        print(f"  {ev.real:+.4f}")

    # Apply W to an initial state and observe the evolution
    # Start in |p_0>
    psi = p_states[0].copy()
    print("\nInitial state |p_0>:")
    for i in range(N * N):
        if abs(psi[i]) > 1e-10:
            x, y = divmod(i, N)
            print(f"  |{x},{y}>: {psi[i]:.4f}")

    psi = W @ psi
    print("\nAfter one step W|p_0>:")
    for i in range(N * N):
        if abs(psi[i]) > 1e-10:
            x, y = divmod(i, N)
            print(f"  |{x},{y}>: {psi[i]:.4f}")

szegedy_walk_example()

To search for a marked vertex m on a graph, the Szegedy walk is modified. The reflection R_A is replaced by R_A’ that omits the term for the marked vertex:

R_A’ = 2 * (sum_{x != m} |p_x><p_x|) - I

This is analogous to the oracle in Grover’s algorithm: the walk “avoids” the marked vertex, causing amplitude to accumulate there through interference.

Algorithm:

  1. Prepare the initial state |psi_0> = (1/sqrt(N)) * sum_x |p_x>.
  2. Repeat for O(sqrt(N/M)) steps: apply W’ = R_B * R_A’.
  3. Measure the first register. The marked vertex appears with high probability.

For a graph with N vertices and M marked vertices, this finds a marked vertex in O(sqrt(N/M)) steps, matching Grover’s optimal query complexity. The key insight is that the spectral gap of the modified walk operator controls the speed of amplitude concentration.

The geometric picture

The search algorithm has a clean geometric interpretation. Define two states:

|good> = (1/sqrt(M)) * sum_{x in marked} |p_x> |bad> = (1/sqrt(N - M)) * sum_{x not in marked} |p_x>

The initial state lies almost entirely in the |bad> direction. Each application of W’ rotates the state toward |good> by an angle proportional to sqrt(M/N). After O(sqrt(N/M)) rotations, the state aligns with |good>, and measurement yields a marked vertex.

Search on a hypercube

The hypercube {0,1}^n has N = 2^n vertices. Each vertex connects to n neighbors (obtained by flipping one bit). A quantum walk on the hypercube naturally coincides with Grover’s algorithm. The walk step (coin + shift on the hypercube) is equivalent to the Grover diffusion operator.

def quantum_walk_search_hypercube(n_qubits, marked_state, n_steps):
    """
    Quantum walk search on the hypercube {0,1}^n.

    On the hypercube, the quantum walk search reduces to Grover's algorithm.
    The coin register encodes which of the n neighbors to move to (which bit to flip).
    For the hypercube, the optimal coin is the Grover diffusion over n directions.

    For simplicity, this implementation uses the Grover oracle + diffusion formulation,
    which is equivalent to the quantum walk on the hypercube.

    Parameters
    ----------
    n_qubits : int
        Number of qubits (log2 of number of vertices).
    marked_state : str
        Binary string of the marked vertex, e.g. "1010".
    n_steps : int
        Number of walk steps (Grover iterations).
    """
    qc = QuantumCircuit(n_qubits, n_qubits)

    # Prepare uniform superposition (initial state of the walk)
    qc.h(range(n_qubits))

    for _ in range(n_steps):
        # Oracle: flip phase of the marked state
        for i, bit in enumerate(reversed(marked_state)):
            if bit == "0":
                qc.x(i)
        # Multi-controlled Z gate: flip phase when all qubits are |1>
        # Implemented as H on target, MCX, H on target
        qc.h(n_qubits - 1)
        qc.mcx(list(range(n_qubits - 1)), n_qubits - 1)
        qc.h(n_qubits - 1)
        for i, bit in enumerate(reversed(marked_state)):
            if bit == "0":
                qc.x(i)

        # Diffusion: inversion about the mean (= walk step on hypercube)
        qc.h(range(n_qubits))
        qc.x(range(n_qubits))
        qc.h(n_qubits - 1)
        qc.mcx(list(range(n_qubits - 1)), n_qubits - 1)
        qc.h(n_qubits - 1)
        qc.x(range(n_qubits))
        qc.h(range(n_qubits))

    qc.measure(range(n_qubits), range(n_qubits))
    return qc

# Search for marked state "1010" in a 4-qubit (16-vertex) hypercube
# Optimal number of steps: floor(pi/4 * sqrt(16)) = floor(pi) = 3
n_qubits = 4
optimal_steps = int(np.floor(np.pi / 4 * np.sqrt(2 ** n_qubits)))
print(f"Optimal steps for {2**n_qubits} vertices: {optimal_steps}")

qc = quantum_walk_search_hypercube(n_qubits, "1010", n_steps=optimal_steps)
sim = AerSimulator()
result = sim.run(transpile(qc, sim), shots=4096).result()
counts = result.get_counts()

# Display results
print("\nSearch results (top 5):")
for state, count in sorted(counts.items(), key=lambda x: -x[1])[:5]:
    print(f"  |{state}> = {int(state, 2):2d}: {count}/{4096} ({100*count/4096:.1f}%)")
print("Target '1010' should appear with high probability.")

Note the distinction: this works on the hypercube because the hypercube’s structure makes the walk equivalent to Grover’s algorithm. On other graphs, quantum walk search and Grover’s algorithm are genuinely different, and the walk can outperform Grover in structured settings.

Element Distinctness

The element distinctness problem illustrates quantum walk search on a non-trivial graph structure.

Problem. Given a list of N numbers x_1, …, x_N, determine whether any two are equal. You access the list through an oracle: querying index i returns x_i.

Classical complexity. O(N) queries using a hash table: insert each element and check for collisions.

Quantum complexity. Ambainis (2004) showed a quantum algorithm using O(N^{2/3}) queries. This is provably optimal (a matching lower bound exists).

The algorithm uses a quantum walk on the Johnson graph J(N, r):

  1. Vertices of J(N, r) are all r-element subsets of {1, …, N}. There are C(N, r) vertices.
  2. Edges connect two subsets that differ in exactly one element (swap one element for another).
  3. Choose r = N^{2/3}.

Algorithm steps:

  1. Setup: Pick a random r-subset S. Query x_i for all i in S. Store the results in quantum memory. This costs r = N^{2/3} queries.

  2. Check: Determine if any two elements in S are equal. If yes, output “collision found.” This costs O(r) time.

  3. Walk: If no collision in S, perform a quantum walk on J(N, r). Each walk step swaps one element of S for a new element from {1, …, N} \ S, queries the new element, and updates the quantum memory. Each step costs O(1) queries.

  4. Search: Modify the walk to search for subsets containing a collision. The walk runs for O(sqrt(C(N, r) / (number of marked vertices))) steps. A marked vertex is a subset S containing a pair (i, j) with x_i = x_j.

Why O(N^{2/3}) queries: The setup costs N^{2/3} queries. The number of walk steps is approximately sqrt(N/r) = sqrt(N^{1/3}) = N^{1/6}. Each step costs O(1) queries. The total is N^{2/3} + N^{1/6} = O(N^{2/3}).

Why this is optimal: Aaronson and Shi (2004) proved a lower bound of Omega(N^{2/3}) for element distinctness in the quantum query model, matching the upper bound.

The element distinctness algorithm extends to the k-element distinctness problem (find k equal elements), which requires O(N^{k/(k+1)}) quantum queries. For k=2 this gives N^{2/3}; for k=3 it gives N^{3/4}.

Continuous-Time Quantum Walks

An alternative formulation uses the adjacency matrix A of the graph as a Hamiltonian. The walker evolves continuously under:

U(t) = exp(-iAt)

There is no coin register; the spreading is driven directly by the graph structure. The walker starts at a single vertex and evolves as a wave packet across the graph.

Key properties:

  • The evolution is unitary: probability is conserved.
  • The eigenstates of A determine the dynamics. If A has eigendecomposition A = sum_k lambda_k |v_k><v_k|, then U(t) = sum_k exp(-ilambda_kt) |v_k><v_k|.
  • On a line graph, the continuous-time walk gives ballistic spreading, similar to the discrete-time walk.
from scipy.linalg import expm

def continuous_time_walk_path(n_vertices=8, times=None):
    """
    Simulate a continuous-time quantum walk on the path graph P_n.

    The Hamiltonian is the adjacency matrix of the path.
    Starting state: vertex 0 (leftmost).
    """
    if times is None:
        times = [0, 1.0, 2.0, 3.0, 5.0, 8.0]

    N = n_vertices
    # Adjacency matrix of the path graph P_N
    A = np.zeros((N, N))
    for i in range(N - 1):
        A[i, i + 1] = 1
        A[i + 1, i] = 1

    # Initial state: vertex 0
    psi0 = np.zeros(N)
    psi0[0] = 1.0

    print(f"Continuous-time quantum walk on P_{N}")
    print(f"Starting at vertex 0\n")

    for t in times:
        U = expm(-1j * A * t)
        psi_t = U @ psi0
        probs = np.abs(psi_t) ** 2

        prob_str = "  ".join(f"{p:.3f}" for p in probs)
        print(f"t = {t:5.1f}: [{prob_str}]")

    # Eigenvalues of A for the path graph
    eigenvalues = np.linalg.eigvalsh(A)
    print(f"\nEigenvalues of A: {np.round(eigenvalues, 4)}")
    print("The eigenvalue spread determines the walk speed.")
    print(f"Spectral range: {eigenvalues[-1] - eigenvalues[0]:.4f}")

continuous_time_walk_path(n_vertices=8)

The output shows the probability wave propagating from vertex 0 across the path graph. At specific times the probability concentrates on the far end (perfect state transfer), and at other times it spreads across the graph. This behavior contrasts with the classical continuous-time random walk (governed by the diffusion equation), which spreads monotonically and exponentially converges to the uniform distribution.

Comparison: continuous vs discrete time

Continuous-time and discrete-time quantum walks are closely related but not identical:

  • Universality. Both models are computationally universal (they can simulate any quantum computation).
  • Coin register. Discrete-time walks require a coin register; continuous-time walks do not. This makes continuous-time walks more natural for some applications.
  • Spectral connection. Discrete-time walk eigenvalues relate to arccos of the classical walk eigenvalues; continuous-time walk eigenvalues are the adjacency matrix eigenvalues directly.
  • Graph dependence. On regular graphs, the two models give similar speedups. On non-regular graphs, their behavior can differ.

Spatial Search on a 2D Grid

One of the most surprising results in quantum walk theory concerns spatial search on a 2D grid.

The problem. A walker lives on an n x n grid (N = n^2 vertices). One vertex is “marked.” Find the marked vertex.

Classical bound. A classical random walk on the 2D grid takes O(N log N) steps to cover all vertices (cover time), so classical walk-based search takes O(N log N).

Grover’s algorithm. Grover’s algorithm treats the N vertices as an unstructured set and finds the marked vertex in O(sqrt(N)) queries. But Grover ignores the grid structure entirely. The question is whether a quantum walk on the grid, which respects the locality of the graph, can match or beat this.

Naive quantum walk. The straightforward Szegedy walk on the 2D grid has a spectral gap of O(1/N), giving a search time of O(sqrt(N) * sqrt(1/delta)) = O(N). This is no better than classical!

The fix (Ambainis, Kempe, Rivosh, 2005). The 2D grid is in the “critical dimension” where the standard quantum walk fails. Ambainis, Kempe, and Rivosh showed that a modified algorithm achieves O(sqrt(N) * sqrt(log N)) steps. The modification adds an extra phase to the walk operator near the marked vertex. This remains the best known quantum walk search on a 2D grid.

Why 2D is special. The spectral gap of the random walk on a d-dimensional grid scales as O(1/n^2) = O(1/N^{2/d}). For d >= 3, the Szegedy walk search runs in O(sqrt(N)) steps, matching Grover. For d = 2, the gap is O(1/N), and the walk needs the logarithmic correction. For d = 1 (the line), the gap is O(1/N^2), and quantum walk search offers no speedup at all. The 2D grid sits exactly at the threshold between speedup and no-speedup.

This result shows that graph structure genuinely matters for quantum walk algorithms: the same algorithm that gives quadratic speedup on a 3D grid or hypercube gives only a sqrt(log N) advantage on a 2D grid and no advantage on a 1D line.

Circuit Resource Analysis

Understanding the gate cost of quantum walk circuits is essential for estimating feasibility on real hardware.

1D walk on a line with N = 2^n vertices, t steps.

Each step consists of:

  • Coin operation: 1 single-qubit gate (Hadamard).
  • Controlled increment (shift right): multi-controlled X gates with 1, 2, …, n controls. Each k-controlled X decomposes into O(k) Toffoli gates using one ancilla, or O(k^2) without ancillas.
  • Controlled decrement (shift left): same cost as the increment, plus 2 X gates on the coin.

Using the ancilla-free decomposition, the increment costs sum_{k=1}^{n} k = n(n+1)/2 Toffoli-equivalent gates. Each Toffoli decomposes into 6 CNOT + several single-qubit gates. The CNOT count for one step:

CNOT_per_step = 2 * n(n+1)/2 * 6 + overhead = 6n(n+1) + small corrections

Example: 10-step walk on 32 vertices (n = 5 position qubits).

  • Toffoli-equivalent gates per step: 2 * 5 * 6 / 2 = 30 (increment + decrement)
  • CNOT count per step: ~30 * 6 = 180 CNOTs
  • Total for 10 steps: ~1800 CNOTs + 10 single-qubit gates (coins)
  • Circuit depth (assuming limited parallelism): roughly 1000-1500 layers

Comparison with Grover’s algorithm on the same space.

Grover’s algorithm on 32 vertices (5 qubits) takes floor(pi/4 * sqrt(32)) approx 4 iterations. Each iteration requires one oracle (5-qubit multi-controlled Z, about 20 CNOTs) and one diffusion (about 20 CNOTs). Total: ~160 CNOTs.

Grover is much shallower for flat, unstructured search. Quantum walk search pays off when the graph structure encodes problem constraints that Grover cannot exploit (such as the Johnson graph for element distinctness).

When to prefer quantum walk search over Grover:

  • The search space has a natural graph structure that can be exploited.
  • The problem reduces to finding a marked vertex on a specific graph where the quantum walk speedup exceeds what Grover provides on the unstructured version.
  • The walk step can be implemented efficiently (O(polylog(N)) gates per step).

Qiskit Implementation: Quantum Walk Search on a Line

This implementation combines the quantum walk on a line with an oracle that marks a target position. The walk spreads the amplitude, and the oracle concentrates it on the target.

def quantum_walk_search_line(n_position_qubits, target_position, n_steps):
    """
    Quantum walk search on a cycle (periodic line) of 2^n vertices.

    This is a simplified demonstration. The oracle marks the target position
    by flipping its phase. The walk step spreads amplitude across the cycle.
    The combination of walk + oracle concentrates amplitude on the target.

    Note: this is NOT the same as Grover's algorithm. The walk operator
    respects the cycle topology, so the algorithm exploits the graph structure.
    On a cycle of N vertices, this takes O(sqrt(N)) steps to find the target.

    Parameters
    ----------
    n_position_qubits : int
        Number of position qubits. The cycle has 2^n vertices.
    target_position : int
        The marked vertex (integer in [0, 2^n - 1]).
    n_steps : int
        Number of walk-search iterations.
    """
    n = n_position_qubits
    N = 2 ** n
    coin = QuantumRegister(1, "coin")
    position = QuantumRegister(n, "pos")
    creg = ClassicalRegister(n, "meas")
    qc = QuantumCircuit(coin, position, creg)

    # Prepare uniform superposition over all positions
    qc.h(position)

    # Convert target to binary
    target_bits = format(target_position, f"0{n}b")

    for step in range(n_steps):
        # --- Oracle: flip phase of target position ---
        # Apply X gates to position qubits where the target bit is 0
        for i, bit in enumerate(reversed(target_bits)):
            if bit == "0":
                qc.x(position[i])
        # Multi-controlled Z on position register
        qc.h(position[n - 1])
        qc.mcx(list(position[:-1]), position[n - 1])
        qc.h(position[n - 1])
        # Undo X gates
        for i, bit in enumerate(reversed(target_bits)):
            if bit == "0":
                qc.x(position[i])

        # --- Quantum walk step ---
        # Coin flip
        qc.h(coin[0])

        # Shift right (controlled increment) when coin = |1>
        for i in range(n - 1, 0, -1):
            controls = [coin[0]] + [position[j] for j in range(i)]
            qc.mcx(controls, position[i])
        qc.cx(coin[0], position[0])

        # Shift left (controlled decrement) when coin = |0>
        qc.x(coin[0])
        for i in range(n - 1, 0, -1):
            controls = [coin[0]] + [position[j] for j in range(i)]
            qc.mcx(controls, position[i])
        qc.cx(coin[0], position[0])
        qc.x(coin[0])

    qc.measure(position, creg)
    return qc

# Search for position 5 on a 16-vertex cycle (4 position qubits)
n_pos = 4
target = 5
N = 2 ** n_pos
# Rough optimal steps for quantum walk search on a cycle
optimal_steps = int(np.ceil(np.sqrt(N)))
print(f"Searching for vertex {target} on a {N}-vertex cycle")
print(f"Using {optimal_steps} steps\n")

qc = quantum_walk_search_line(n_pos, target, optimal_steps)
sim = AerSimulator()
result = sim.run(transpile(qc, sim), shots=8192).result()
counts = result.get_counts()

print("Search results (top 5):")
for state, count in sorted(counts.items(), key=lambda x: -x[1])[:5]:
    pos = int(state, 2)
    marker = " <-- TARGET" if pos == target else ""
    print(f"  vertex {pos:2d}: {count}/{8192} ({100*count/8192:.1f}%){marker}")

Key Takeaways

Spreading. Classical random walks spread diffusively (standard deviation ~ sqrt(t)). Quantum walks spread ballistically (standard deviation ~ t). This quadratic speedup in spreading is the fundamental resource that quantum walk algorithms exploit.

Szegedy framework. For any graph with transition matrix P and spectral gap delta, the Szegedy walk has an eigenvalue gap of Theta(sqrt(delta)). This systematically converts classical random walk algorithms into quantum algorithms with a quadratic speedup, applicable to any problem modeled as a Markov chain.

Search complexity. Quantum walk search finds a marked vertex in O(sqrt(N/M)) steps on many graphs. The constant factors and feasibility depend on the graph structure: the hypercube gives clean Grover-like search, the 2D grid requires a logarithmic correction, and the 1D line offers no speedup.

Applications beyond search. Quantum walk algorithms currently achieve the best known quantum complexities for element distinctness (O(N^{2/3})), triangle finding (O(N^{5/4})), matrix product verification, and group commutativity testing.

Circuit costs. Quantum walk circuits are deeper than Grover’s algorithm for unstructured search. The walk pays off when the graph structure encodes problem constraints that reduce the effective search space. Understanding the gate count is critical for assessing near-term feasibility: a 10-step walk on 32 vertices requires roughly 1800 CNOT gates, well within the reach of current error-corrected architectures but challenging for near-term noisy devices.

Continuous vs discrete. Both continuous-time and discrete-time quantum walks are computationally universal. Continuous-time walks avoid the coin register and connect directly to Hamiltonian simulation. Discrete-time walks offer more flexibility in choosing the coin operator and are easier to implement on gate-based quantum computers.

Was this tutorial helpful?

Continue learning with these guides