核心量子算法¶
本节介绍教科书级量子算法：搜索（Grover）、相位估计（QPE）、振幅估计， 以及变分量子本征求解器（VQD）。这些是理解更复杂算法的前置基础。
Grover’s search algorithm — complete example.¶
Source: examples/2_advanced/algorithms/grover.py
Status: pass
Demonstrates:
Oracle construction for an n-qubit Grover search
Diffusion operator (amplitude amplification)
Running the algorithm with UnifiedQuantum simulators
Using the measurement module for result analysis
Usage: python grover.py [–n-qubits N] [–marked-state STATE] [–shots N]
References: Grover, L. K. (1996). “A fast quantum mechanical algorithm for database search.” STOC ‘96. https://arxiv.org/abs/quant-ph/9605043
Source code
#!/usr/bin/env python
"""Grover's search algorithm — complete example.
Demonstrates:
* Oracle construction for an n-qubit Grover search
* Diffusion operator (amplitude amplification)
* Running the algorithm with UnifiedQuantum simulators
* Using the measurement module for result analysis
Usage:
python grover.py [--n-qubits N] [--marked-state STATE] [--shots N]
References:
Grover, L. K. (1996). "A fast quantum mechanical algorithm
for database search." STOC '96.
https://arxiv.org/abs/quant-ph/9605043
[doc-require: ]
"""
import argparse
import sys
import math
import numpy as np
# Add parent directory to path so we can import uniqc when running as a script
sys.path.insert(0, str(__file__.rsplit("/", 2)[0]))
from uniqc import Circuit
from uniqc.simulator import Simulator
from uniqc import pauli_expectation
def build_oracle(n_qubits: int, marked_state: int) -> tuple[Circuit, list[int]]:
"""Build a Grover oracle that flips the phase of the marked basis state.
The oracle applies a Z gate on an ancilla qubit (last qubit) conditioned
on all data qubits being in the marked state. This implements the standard
phase-kickback oracle:
U_φ|x⟩|−⟩ = (−1)^{[x==marked]}|x⟩|−⟩
where |−⟩ = (|0⟩−|1⟩)/√2 is the ancilla in the magic basis.
Args:
n_qubits: Number of data qubits.
marked_state: Integer encoding the marked basis state (0 to 2^n−1).
Returns:
A tuple of (circuit, ancilla_qubits) where ancilla_qubits is the
list containing the ancilla qubit index.
"""
total_qubits = n_qubits + 1
c = Circuit()
# Initialise ancilla to |−⟩ = X·H·|0⟩
c.x(n_qubits)
c.h(n_qubits)
# Apply multi-controlled Z (MCZ) targeting the ancilla qubit
# This flips the ancilla phase iff all data qubits are in |1⟩
# We achieve this with a CCZ (Toffoli equivalent) followed by some X gates
# to encode the marked-state condition.
#
# Strategy: flip data qubits that are |0⟩ in the marked state,
# then apply MCZ from all data qubits to ancilla, then flip back.
marked_bits = [(marked_state >> (n_qubits - 1 - i)) & 1 for i in range(n_qubits)]
# Flip qubits that should be |0⟩ in the marked state
for i, bit in enumerate(marked_bits):
if bit == 0:
c.x(i)
# Apply multi-controlled Z gate to ancilla
# We use a sequence of Toffoli (CCX) gates to implement MCZ
# For n data qubits, we need n-1 ancillas for an efficient MCZ
# Here we use a simpler approach: apply CCZ if n_qubits=2, or chain Toffolis
if n_qubits == 1:
c.cz(i=n_qubits, target=n_qubits)
elif n_qubits == 2:
c.toffoli(0, 1, n_qubits)
c.z(n_qubits) # CCZ = CCX with final Z
c.toffoli(0, 1, n_qubits)
else:
# General multi-controlled Z using multiple Toffoli stages
# MCZ: apply CNOT cascade then H on target then CNOT cascade
# Using the standard linear-depth circuit
_apply_mcz(c, list(range(n_qubits)), n_qubits)
# Flip back
for i, bit in enumerate(marked_bits):
if bit == 0:
c.x(i)
# Return circuit with measurements on data qubits only
c.measure(*list(range(n_qubits)))
return c, list(range(n_qubits))
def _apply_mcz(circuit: Circuit, controls: list[int], target: int) -> None:
"""Apply multi-controlled Z gate using the standard linear-depth circuit.
Circuit:
(H on target)
CNOT cascade up
CNOT cascade down
(H on target)
This is equivalent to applying Z on the target iff all controls are |1⟩.
"""
n = len(controls)
if n == 0:
circuit.z(target)
return
if n == 1:
circuit.cz(controls[0], target)
return
# H on target
circuit.h(target)
# Cascade CNOTs from controls to target
circuit.cx(controls[0], controls[1])
for i in range(1, n - 1):
circuit.cx(controls[i], controls[i + 1])
circuit.cx(controls[-1], target)
# Cascade back
for i in range(n - 2, 0, -1):
circuit.cx(controls[i], controls[i + 1])
circuit.cx(controls[0], controls[1])
# H on target
circuit.h(target)
def build_diffusion(n_qubits: int) -> Circuit:
"""Build the Grover diffusion operator for amplitude amplification.
D = H^{⊗n} · Z · H^{⊗n} · A_0
where A_0 = H^{⊗n}·Z^{⊗n} is the initialisation (uniform superposition).
D = 2|s⟩⟨s| − I where |s⟩ = H^{⊗n}|0⟩^{⊗n}.
Args:
n_qubits: Number of qubits.
Returns:
A Circuit implementing the diffusion operator.
"""
c = Circuit()
# Apply Hadamard to all
for i in range(n_qubits):
c.h(i)
# Apply X to all
for i in range(n_qubits):
c.x(i)
# Multi-controlled Z (phase flip)
if n_qubits == 1:
c.z(0)
elif n_qubits == 2:
c.toffoli(0, 1, n_qubits)
c.z(n_qubits)
c.toffoli(0, 1, n_qubits)
else:
# Use ancilla at position n_qubits as target for MCZ
_apply_mcz(c, list(range(n_qubits)), n_qubits)
# Flip back X
for i in range(n_qubits):
c.x(i)
# Flip back H
for i in range(n_qubits):
c.h(i)
return c
def run_grover(
n_qubits: int,
marked_state: int,
shots: int = 4096,
) -> dict[str, float]:
"""Run Grover's search for the given marked state.
Args:
n_qubits: Number of data qubits (search space = 2^n qubits).
marked_state: The integer index of the target state.
shots: Number of measurement shots.
Returns:
Dictionary of measurement outcome frequencies.
"""
total_qubits = n_qubits + 1
# Initialise all data qubits to |+⟩ (uniform superposition)
c = Circuit()
for i in range(n_qubits):
c.h(i)
# Optimal number of Grover iterations
n_iterations = int(math.pi / 4 * math.sqrt(2**n_qubits))
n_iterations = max(1, min(n_iterations, 10)) # cap for small n
for _ in range(n_iterations):
# Oracle
oracle_c = Circuit()
for i in range(n_qubits):
oracle_c.cx(i, n_qubits) # CCNOT-style phase flip
# Marked-state flip: flip the marked state's amplitude
# We encode the marked state as a sequence of X gates
marked_bits = [(marked_state >> (n_qubits - 1 - i)) & 1 for i in range(n_qubits)]
for i, bit in enumerate(marked_bits):
if bit == 0:
oracle_c.x(i)
oracle_c.toffoli(0, 1, n_qubits)
oracle_c.z(n_qubits)
oracle_c.toffoli(0, 1, n_qubits)
for i, bit in enumerate(marked_bits):
if bit == 0:
oracle_c.x(i)
# (End oracle)
# Diffusion
for i in range(n_qubits):
oracle_c.h(i)
for i in range(n_qubits):
oracle_c.x(i)
if n_qubits == 2:
oracle_c.toffoli(0, 1, n_qubits)
oracle_c.z(n_qubits)
oracle_c.toffoli(0, 1, n_qubits)
else:
_apply_mcz(oracle_c, list(range(n_qubits)), n_qubits)
for i in range(n_qubits):
oracle_c.x(i)
for i in range(n_qubits):
oracle_c.h(i)
c.add_circuit(oracle_c)
c.measure(*list(range(n_qubits)))
# Simulate
sim = Simulator(least_qubit_remapping=False)
result = sim.simulate_shots(c.qasm, shots=shots)
total = sum(result.values())
return {f"{k:0{n_qubits}b}": v / total for k, v in result.items()}
def main():
parser = argparse.ArgumentParser(description="Grover's Search Algorithm")
parser.add_argument("--n-qubits", type=int, default=3, help="Number of data qubits")
parser.add_argument(
"--marked-state",
type=int,
default=5,
help="Integer index of the marked state (0 to 2^n−1)",
)
parser.add_argument("--shots", type=int, default=4096, help="Number of shots")
args = parser.parse_args()
n = args.n_qubits
marked = args.marked_state
if marked >= 2**n:
print(f"Error: marked_state {marked} is out of range for {n} qubits (max {2**n-1})")
sys.exit(1)
print(f" Grover's Search — {n} data qubits")
print(f" Marked state: {marked} ({marked:0{n}b})")
print(f" Search space: {2**n} states")
print()
result = run_grover(n, marked, shots=args.shots)
marked_str = f"{marked:0{n}b}"
marked_prob = result.get(marked_str, 0.0)
print(f" Results (top 5 most probable states):")
sorted_results = sorted(result.items(), key=lambda x: x[1], reverse=True)
for state, prob in sorted_results[:5]:
marker = " ← TARGET" if state == marked_str else ""
print(f" |{state}⟩ {prob*100:5.1f}%{marker}")
print()
print(f" Target probability: {marked_prob*100:.1f}%")
print(f" Expected (ideal): ~{95.0:.1f}% (after optimal iterations)")
print()
print(f" ✓ Run complete.")
if __name__ == "__main__":
main()
Stdout
Grover's Search — 3 data qubits
Marked state: 5 (101)
Search space: 8 states
Results (top 5 most probable states):
|110⟩ 26.6%
|101⟩ 25.8% ← TARGET
|011⟩ 25.5%
|000⟩ 22.1%
Target probability: 25.8%
Expected (ideal): ~95.0% (after optimal iterations)
✓ Run complete.
Quantum Phase Estimation (QPE) — complete example.¶
Source: examples/2_advanced/algorithms/qpe.py
Status: pass
Demonstrates:
QPE circuit construction with phase register + eigenstate register
Inverse Quantum Fourier Transform (QFTdagger)
Running QPE with UnifiedQuantum simulators
Using the measurement module to extract phase bits
Connecting the estimated phase to the eigenvalue
Usage: python qpe.py [–n-precision N] [–unitary TYPE] [–shots N]
References: Nielsen & Chuang, “Quantum Computation and Quantum Information”, Chapter 5. Cleve et al. (1998), “Efficient Discrete Random Unitary Circuits for Approximating the Quantum Fourier Transform.” https://arxiv.org/abs/quant-ph/9904026
Source code
#!/usr/bin/env python
"""Quantum Phase Estimation (QPE) — complete example.
Demonstrates:
* QPE circuit construction with phase register + eigenstate register
* Inverse Quantum Fourier Transform (QFTdagger)
* Running QPE with UnifiedQuantum simulators
* Using the measurement module to extract phase bits
* Connecting the estimated phase to the eigenvalue
Usage:
python qpe.py [--n-precision N] [--unitary TYPE] [--shots N]
References:
Nielsen & Chuang, "Quantum Computation and Quantum Information", Chapter 5.
Cleve et al. (1998), "Efficient Discrete Random Unitary Circuits for Approximating
the Quantum Fourier Transform." https://arxiv.org/abs/quant-ph/9904026
[doc-require: ]
"""
import argparse
import sys
import cmath
import numpy as np
sys.path.insert(0, str(__file__.rsplit("/", 2)[0]))
from uniqc import Circuit
from uniqc.simulator import Simulator
from uniqc import basis_rotation_measurement
def apply_cu1(circuit: Circuit, control: int, target: int, theta: float) -> None:
"""Apply controlled-U1 gate, where U1(θ) = diag(1, e^{iθ}).
In QASM: cu1(theta) q[control], q[target]
Decomposition: U1(θ) = Rz(θ) = [[1,0],[0,e^{iθ}]]
We implement this as:
u1(theta/2) on target
cx from control to target
u1(-theta/2) on target
cx from control to target
which gives: CU1(θ) = [[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,e^{iθ}]]
"""
circuit.u1(target, theta / 2)
circuit.cx(control, target)
circuit.u1(target, -theta / 2)
circuit.cx(control, target)
def apply_qft(circuit: Circuit, qubits: list[int]) -> None:
"""Apply the Quantum Fourier Transform to a register of qubits.
QFT on |x⟩ = (1/√2^n) ∑_{k=0}^{2^n-1} e^{2πikx/2^n} |k⟩
Circuit for n qubits (leftmost = most significant):
for i = 0 to n-1:
H on q[i]
for j = i+1 to n-1:
CU1(π/2^{j-i}) on (q[j], q[i])
"""
n = len(qubits)
for i in range(n):
circuit.h(qubits[i])
for j in range(i + 1, n):
theta = np.pi / (2 ** (j - i))
apply_cu1(circuit, qubits[j], qubits[i], theta)
def apply_qft_dagger(circuit: Circuit, qubits: list[int]) -> None:
"""Apply QFTdagger (inverse QFT) — the final step of QPE.
QFTdagger = apply_qft in reverse order with adjoint rotations:
for i = n-1 down to 0:
for j = i+1 to n-1:
CU1(-π/2^{j-i}) on (q[j], q[i])
H on q[i]
"""
n = len(qubits)
for i in range(n - 1, -1, -1):
for j in range(n - 1, i, -1):
theta = -np.pi / (2 ** (j - i))
apply_cu1(circuit, qubits[j], qubits[i], theta)
circuit.h(qubits[i])
def build_qpe_circuit(
n_precision: int,
unitary_matrix: np.ndarray,
eigenstate: list[int],
) -> tuple[Circuit, list[int]]:
"""Build the complete QPE circuit.
Args:
n_precision: Number of precision (ancilla) qubits.
The phase is estimated to n_precision bits.
unitary_matrix: 2^n × 2^n unitary matrix representing U.
U must be diagonal in the computational basis for
clean phase estimation; otherwise the result is a
superposition of phases (equivalent to eigenstate decomposition).
eigenstate: The n-qubit basis state |u⟩ that is an eigenstate of U.
Given as a list of bits, e.g. [1, 0] for |10⟩.
Returns:
A tuple (circuit, precision_qubits).
"""
n_system = len(eigenstate)
total_qubits = n_precision + n_system
q_precision = list(range(n_precision))
q_system = list(range(n_precision, total_qubits))
c = Circuit()
# Step 1: Initialise eigenstate on system qubits
for i, bit in enumerate(eigenstate):
if bit == 1:
c.x(q_system[i])
# Step 2: Apply Hadamard on precision qubits (uniform superposition)
for q in q_precision:
c.h(q)
# Step 3: Controlled powers of U
# For k = 0 to n_precision-1, apply U^{2^k} controlled by precision qubit k
for k, q_ctrl in enumerate(q_precision):
power = 2**k
# Compute U^{2^k} by repeated matrix multiplication
U_power = unitary_matrix
for _ in range(power - 1):
U_power = U_power @ unitary_matrix
# Apply controlled-U^{2^k} via basis decomposition
# Since U is diagonal, controlled-U is just CU1 gates
_apply_controlled_unitary(c, q_ctrl, q_system, U_power)
# Step 4: Inverse QFT on precision qubits
apply_qft_dagger(c, q_precision)
# Measure precision qubits
c.measure(*q_precision)
return c, q_precision
def _apply_controlled_unitary(
circuit: Circuit,
control: int,
system_qubits: list[int],
U: np.ndarray,
) -> None:
"""Apply controlled-U where U is diagonal in the computational basis.
For a diagonal U = diag(e^{iθ_0}, e^{iθ_1}, ...):
CU |c⟩|s⟩ = |c⟩ ⊗ U|s⟩ if c=1, else |s⟩
= exp(iθ_s) |c⟩|s⟩ if c=1
This is implemented as a sequence of CU1 gates encoding the phases.
"""
n = len(system_qubits)
dim = 2**n
# Extract diagonal elements
diags = np.diag(U)
phases = [cmath.phase(d) for d in diags]
# Apply CU1 gates encoding the phases
# U|s⟩ = exp(iθ_s)|s⟩
# For a 1-qubit system: CU1(θ) on (control, sys) implements e^{iθ}|1⟩|s⟩
# For multi-qubit: decompose into CU1s using binary phase encoding
if n == 1:
# CU1 on (control, system_qubits[0]) encodes phase for |1⟩ on control
theta = phases[1] - phases[0] # relative phase
apply_cu1(circuit, control, system_qubits[0], theta)
# Also need to apply the global phase e^{iθ_0} — handled by QPE automatically
else:
# Multi-qubit controlled-U: decompose using Toffoli cascade
# U|s_0 s_1 ...⟩ = e^{iθ_s}|s⟩
# We encode each basis state's phase as a sum of CU1 angles
# Using Gray code decomposition
_apply_cu_cascade(circuit, control, system_qubits, phases)
def _apply_cu_cascade(
circuit: Circuit,
control: int,
system_qubits: list[int],
phases: list[float],
) -> None:
"""Apply multi-controlled phase gate using a cascade of CU1 gates.
For n system qubits with 2^n basis states, each having phase θ_k,
we decompose this into n CU1 gates using the following trick:
The phase on state |k⟩ is encoded in binary, and each qubit
contributes a phase conditioned on higher-order bits.
"""
n = len(system_qubits)
# For simplicity with n=2 qubits:
if n == 2:
# U|ab⟩ = e^{iθ_ab}|ab⟩
# CU1(θ_01) on (ctrl, q1) gives e^{iθ_ab} when ctrl=1 and q1=b=1
# This depends on q1 only, not q0.
# We use the decomposition:
# e^{iθ_00}|00⟩ |e^{iθ_01}|01⟩ |e^{iθ_10}|10⟩ |e^{iθ_11}|11⟩
# = e^{iθ_00}|00⟩ ⊗ (e^{i(θ_01-θ_00)} for |01⟩) ⊗ (e^{i(θ_10-θ_00)} for |10⟩) ⊗
# (e^{i(θ_11-θ_01+θ_00)} for |11⟩)
#
# CU1(ctrl, q1) only applies phase when ctrl=1 AND q1=1
# CU1(ctrl, q0) only applies phase when ctrl=1 AND q0=1
#
# We decompose using the standard method:
# θ_0 = θ_00 (base phase, absorbed into global)
# θ_1 = θ_01 - θ_00 (q1 phase relative to q0)
# θ_2 = θ_10 - θ_00 (q0 phase relative to q1, for |10⟩)
# θ_3 = θ_11 + θ_00 - θ_01 - θ_10 (correction for |11⟩)
base = phases[0]
d1 = phases[1] - phases[0] # relative to |01⟩
d2 = phases[2] - phases[0] # relative to |10⟩
d3 = phases[3] + phases[0] - phases[1] - phases[2] # correction
apply_cu1(circuit, control, system_qubits[0], d3) # q1 is MSB for |11⟩
apply_cu1(circuit, control, system_qubits[1], d2) # q0
else:
# For n=1 case
apply_cu1(circuit, control, system_qubits[0], phases[1] - phases[0])
def run_qpe(
n_precision: int,
unitary: str = "t",
shots: int = 4096,
) -> tuple[float, float, np.ndarray]:
"""Run QPE to estimate the phase of a quantum gate.
Args:
n_precision: Number of precision qubits (bits of phase precision).
unitary: Type of unitary to estimate:
"t" → T = diag(1, e^{iπ/4}), phase = π/8 (T-gate)
"z" → Z = diag(1, -1), phase = 0.5
"s" → S = diag(1, i), phase = 0.25
"rz" → Rz(π/3), phase = π/6 ≈ 0.5236
shots: Number of measurement shots.
Returns:
Tuple of (estimated_phase, true_phase, phase_counts dict).
"""
# Define the unitary matrices
if unitary == "t":
# T gate: |0⟩ → phase 0 (eigenvalue 1), |1⟩ → phase π/8 (eigenvalue e^{iπ/4})
U = np.diag([1, np.exp(1j * np.pi / 4)])
true_phase = 0 # phase of |0⟩ eigenstate
eigenstate = [0] # |0⟩ is eigenstate with eigenvalue 1 (phase 0)
elif unitary == "z":
# Z gate: phase = 0.5 (since eigenvalue = -1 = e^{iπ})
U = np.diag([1, -1])
true_phase = 0.5 # eigenvalue = e^{2πi*0.5} = -1
eigenstate = [0] # |0⟩ → phase 0
elif unitary == "s":
# S gate: phase = 0.25 (eigenvalue = e^{iπ/2} = i)
U = np.diag([1, 1j])
true_phase = 0.25
eigenstate = [0]
elif unitary == "rz":
# Rz(π/3) = diag(e^{-iπ/6}, e^{iπ/6}), phase = π/6 / 2π = 1/12 ≈ 0.0833
U = np.diag([np.exp(-1j * np.pi / 6), np.exp(1j * np.pi / 6)])
true_phase = 1 / 12
eigenstate = [0]
else:
raise ValueError(f"Unknown unitary: {unitary}")
# Build QPE circuit
c, q_precision = build_qpe_circuit(n_precision, U, eigenstate)
# Simulate
sim = Simulator(least_qubit_remapping=False)
counts = sim.simulate_shots(c.qasm, shots=shots)
total = sum(counts.values())
# Extract phase from measurement results
# Precision qubits are ordered with q[0] = most significant (leftmost in binary)
# QPE gives: measured integer m → phase ≈ m / 2^n_precision
phase_counts = {f"{k:0{n_precision}b}": v / total for k, v in counts.items()}
# Find most likely outcome
m_estimated = max(counts, key=counts.get)
m_int = m_estimated if isinstance(m_estimated, int) else int(m_estimated, 2)
estimated_phase = m_int / (2**n_precision)
return estimated_phase, true_phase, phase_counts
def main():
parser = argparse.ArgumentParser(description="Quantum Phase Estimation")
parser.add_argument(
"--n-precision", type=int, default=4,
help="Number of precision qubits (default 4)"
)
parser.add_argument(
"--unitary",
type=str,
default="t",
choices=["t", "z", "s", "rz"],
help="Unitary to estimate (default: t gate)",
)
parser.add_argument(
"--shots", type=int, default=4096,
help="Number of measurement shots"
)
args = parser.parse_args()
print(f" Quantum Phase Estimation")
print(f" Precision qubits: {args.n_precision}")
print(f" Phase precision: 1/{2**args.n_precision} = {1/2**args.n_precision:.4f}")
print(f" Unitary: {args.unitary}")
print()
est, true, counts = run_qpe(args.n_precision, args.unitary, args.shots)
print(f" Measurement results:")
for state, prob in sorted(counts.items(), key=lambda x: x[1], reverse=True):
m_int = int(state, 2)
phase = m_int / (2**args.n_precision)
marker = " ← most likely" if prob == max(counts.values()) else ""
print(f" |{state}⟩ prob={prob*100:5.1f}% phase={phase:.4f}{marker}")
print()
print(f" Estimated phase: {est:.4f}")
print(f" True phase: {true:.4f}")
print(f" Absolute error: {abs(est - true):.4f}")
print(f" ✓ QPE complete.")
if __name__ == "__main__":
main()
Stdout
Quantum Phase Estimation
Precision qubits: 4
Phase precision: 1/16 = 0.0625
Unitary: t
Measurement results:
|0111⟩ prob= 41.2% phase=0.4375 ← most likely
|0000⟩ prob= 40.9% phase=0.0000
|0011⟩ prob= 5.1% phase=0.1875
|0100⟩ prob= 5.0% phase=0.2500
|0010⟩ prob= 2.4% phase=0.1250
|0101⟩ prob= 2.2% phase=0.3125
|0110⟩ prob= 1.6% phase=0.3750
|0001⟩ prob= 1.5% phase=0.0625
Estimated phase: 0.4375
True phase: 0.0000
Absolute error: 0.4375
✓ QPE complete.
Quantum Amplitude Estimation (QAE) — complete example.¶
Source: examples/2_advanced/circuits/amplitude_estimation.py
Status: pass
Demonstrates:
Building a simple oracle for amplitude estimation
Running QAE to estimate the probability of “good” states
Using the amplitude_estimation_result function to extract the estimate
Usage: python amplitude_estimation.py [–n-qubits N] [–n-eval-qubits M] [–shots N]
References: Brassard, G., Høyer, P., Mosca, M. & Tapp, A. (2002). “Quantum Amplitude Amplification and Estimation.” AMS Contemporary Mathematics, 305, 53–74.
Source code
#!/usr/bin/env python
"""Quantum Amplitude Estimation (QAE) — complete example.
Demonstrates:
* Building a simple oracle for amplitude estimation
* Running QAE to estimate the probability of "good" states
* Using the amplitude_estimation_result function to extract the estimate
Usage:
python amplitude_estimation.py [--n-qubits N] [--n-eval-qubits M] [--shots N]
References:
Brassard, G., Høyer, P., Mosca, M. & Tapp, A. (2002).
"Quantum Amplitude Amplification and Estimation."
AMS Contemporary Mathematics, 305, 53–74.
[doc-require: ]
[doc-warning-ignore: DeprecationWarning]
"""
import argparse
import sys
import math
# Add parent directory to path so we can import uniqc when running as a script
sys.path.insert(0, str(__file__.rsplit("/", 2)[0]))
from uniqc import Circuit
from uniqc.simulator import Simulator
from uniqc import (
amplitude_estimation_circuit,
amplitude_estimation_result,
)
def build_simple_oracle(qubits: list, marked_state: int) -> Circuit:
"""Build a phase-flip oracle that marks a specific computational basis state.
Flips the phase of |marked_state⟩ by:
1. X gates on qubits where marked_state has 0
2. Multi-controlled Z on all qubits
3. Undo X gates
Args:
qubits: Indices of the search-register qubits the oracle should act on.
marked_state: Integer representing the marked basis state.
Returns:
Oracle circuit.
"""
oracle = Circuit()
n_qubits = len(qubits)
bits = format(marked_state, f"0{n_qubits}b")
# Flip qubits where bit is 0
for i, bit in enumerate(bits):
if bit == "0":
oracle.x(qubits[i])
# Apply multi-controlled Z
if n_qubits == 1:
oracle.z(qubits[0])
elif n_qubits == 2:
oracle.cz(qubits[0], qubits[1])
else:
# For n > 2, use toffoli cascade
oracle.h(qubits[-1])
oracle.toffoli(qubits[0], qubits[1], qubits[-1])
oracle.h(qubits[-1])
# Undo X flips
for i, bit in enumerate(bits):
if bit == "0":
oracle.x(qubits[i])
return oracle
def run_qae(n_qubits: int, n_eval_qubits: int, marked_state: int, shots: int = 4096) -> dict:
"""Run QAE and return results.
Args:
n_qubits: Number of search register qubits.
n_eval_qubits: Number of evaluation qubits (precision).
marked_state: The marked state for the oracle.
shots: Number of measurement shots.
Returns:
Dictionary with estimated probability and measurement counts.
"""
total_qubits = n_qubits + n_eval_qubits
search_qubits = list(range(n_eval_qubits, total_qubits))
eval_qubits = list(range(n_eval_qubits))
oracle = build_simple_oracle(search_qubits, marked_state)
c = amplitude_estimation_circuit(
oracle,
qubits=search_qubits,
eval_qubits=eval_qubits,
)
# Simulate via OriginIR (its toffoli decomposition is wider; QASM2 export
# rejects nested controlled-CCZ produced by the QAE wrapper).
from uniqc.simulator import Simulator
sim = Simulator()
result = sim.simulate_shots(c.originir, shots=shots)
# Convert to bit-strings
counts = {f"{int(k):0{n_eval_qubits}b}": v for k, v in result.items()}
# Estimate amplitude
estimated_a = amplitude_estimation_result(counts, n_eval_qubits)
# True probability (for uniform superposition with one marked state)
true_a = 1.0 / (2 ** n_qubits)
return {
"estimated": estimated_a,
"true": true_a,
"counts": counts,
}
def main():
parser = argparse.ArgumentParser(description="Quantum Amplitude Estimation (QAE)")
parser.add_argument("--n-qubits", type=int, default=2, help="Number of search qubits")
parser.add_argument(
"--n-eval-qubits", type=int, default=3, help="Number of evaluation qubits"
)
parser.add_argument("--shots", type=int, default=4096, help="Number of shots")
args = parser.parse_args()
n = args.n_qubits
m = args.n_eval_qubits
marked = 0 # Mark |0...0⟩
print(f" Quantum Amplitude Estimation")
print(f" Search qubits: {n}, Eval qubits: {m}")
print(f" Marked state: |{marked:0{n}b}⟩")
print()
result = run_qae(n, m, marked, args.shots)
print(f" Estimated probability: {result['estimated']:.6f}")
print(f" True probability: {result['true']:.6f}")
print(f" Error: {abs(result['estimated'] - result['true']):.6f}")
print()
print(" Measurement counts (eval register):")
for outcome, count in sorted(result["counts"].items()):
print(f" |{outcome}⟩: {count}")
if __name__ == "__main__":
main()
Stdout
Quantum Amplitude Estimation
Search qubits: 2, Eval qubits: 3
Marked state: |00⟩
Estimated probability: 0.853553
True probability: 0.250000
Error: 0.603553
Measurement counts (eval register):
|000⟩: 2
|001⟩: 316
|010⟩: 141
|011⟩: 515
|100⟩: 420
|101⟩: 558
|110⟩: 1282
|111⟩: 862
Variational Quantum Deflation (VQD) — complete example.¶
Source: examples/2_advanced/circuits/vqd.py
Status: pass
Demonstrates:
Finding the ground state of H = Z₀ + Z₁ using VQE (layer 0)
Finding the first excited state using VQD with overlap penalty
Using scipy.optimize.minimize as the classical optimiser
Usage: python vqd.py [–n-qubits N] [–n-layers L] [–penalty B]
References: Higgott, O., Wang, D. & Brierley, S. (2019). “Variational Quantum Computation of Excited States.” Quantum 3, 156.
Source code
#!/usr/bin/env python
"""Variational Quantum Deflation (VQD) — complete example.
Demonstrates:
* Finding the ground state of H = Z₀ + Z₁ using VQE (layer 0)
* Finding the first excited state using VQD with overlap penalty
* Using scipy.optimize.minimize as the classical optimiser
Usage:
python vqd.py [--n-qubits N] [--n-layers L] [--penalty B]
References:
Higgott, O., Wang, D. & Brierley, S. (2019).
"Variational Quantum Computation of Excited States."
Quantum 3, 156.
[doc-require: ]
[doc-warning-ignore: DeprecationWarning]
"""
import argparse
import sys
import numpy as np
# Add parent directory to path so we can import uniqc when running as a script
sys.path.insert(0, str(__file__.rsplit("/", 2)[0]))
from uniqc import Circuit
from uniqc.simulator import Simulator
from uniqc import (
vqd_circuit,
vqd_overlap_circuit,
)
def _make_hamiltonian_matrix(n_qubits: int) -> np.ndarray:
"""Build the matrix for H = sum_i Z_i on *n_qubits*."""
dim = 2**n_qubits
H = np.zeros((dim, dim), dtype=complex)
for i in range(n_qubits):
# Z_i = I⊗...⊗Z⊗...⊗I
mat = np.array([1.0])
for j in range(n_qubits):
if j == i:
mat = np.kron(mat, np.diag([1.0, -1.0]))
else:
mat = np.kron(mat, np.eye(2))
H += mat
return H
def _state_from_hea(ansatz_params: np.ndarray, n_qubits: int, n_layers: int) -> np.ndarray:
"""Compute the state vector produced by the HEA ansatz with given params."""
dim = 2**n_qubits
state = np.zeros(dim, dtype=complex)
state[0] = 1.0
# Apply HEA gates to state vector
idx = 0
for _ in range(n_layers):
for q in range(n_qubits):
state = _apply_ry(state, ansatz_params[idx], q, n_qubits)
idx += 1
for q in range(n_qubits - 1):
state = _apply_cnot(state, q, q + 1, n_qubits)
return state
def _apply_ry(state: np.ndarray, theta: float, qubit: int, n_qubits: int) -> np.ndarray:
"""Apply Ry(theta) on qubit to state vector."""
dim = 2**n_qubits
new_state = state.copy()
cos_t = np.cos(theta / 2)
sin_t = np.sin(theta / 2)
for i in range(dim):
if (i >> (n_qubits - 1 - qubit)) & 1 == 0:
j = i | (1 << (n_qubits - 1 - qubit))
new_state[i] = cos_t * state[i] - sin_t * state[j]
new_state[j] = sin_t * state[i] + cos_t * state[j]
return new_state
def _apply_cnot(state: np.ndarray, control: int, target: int, n_qubits: int) -> np.ndarray:
"""Apply CNOT(control, target) to state vector."""
dim = 2**n_qubits
new_state = state.copy()
for i in range(dim):
c_bit = (i >> (n_qubits - 1 - control)) & 1
if c_bit == 1:
j = i ^ (1 << (n_qubits - 1 - target))
new_state[i] = state[j]
new_state[j] = state[i]
return new_state
def vqe_objective(params: np.ndarray, n_qubits: int, n_layers: int, H: np.ndarray) -> float:
"""VQE cost function: ⟨ψ(θ)|H|ψ(θ)⟩."""
state = _state_from_hea(params, n_qubits, n_layers)
return float(np.real(state.conj() @ H @ state))
def vqd_objective(
params: np.ndarray,
n_qubits: int,
n_layers: int,
H: np.ndarray,
prev_states: list,
penalty: float,
) -> float:
"""VQD cost function: ⟨ψ(θ)|H|ψ(θ)⟩ + Σᵢ β|⟨ψ(θ)|φᵢ⟩|²."""
state = _state_from_hea(params, n_qubits, n_layers)
energy = float(np.real(state.conj() @ H @ state))
overlap_penalty = 0.0
for ps in prev_states:
overlap_penalty += penalty * abs(np.dot(state.conj(), ps))**2
return energy + overlap_penalty
def run_vqd(n_qubits: int, n_layers: int, penalty: float) -> None:
"""Run full VQD example: VQE ground state → VQD first excited state."""
from scipy.optimize import minimize
H = _make_hamiltonian_matrix(n_qubits)
n_params = n_qubits * n_layers
# Exact eigenvalues for reference
eigenvalues = np.linalg.eigvalsh(H)
# Step 1: VQE for ground state
print(f"=== VQD Example: H = Z₀ + ... + Z_{{n-1}} on {n_qubits} qubits ===")
print(f"n_layers={n_layers}, penalty={penalty}, n_params={n_params}")
print()
print(f"Exact eigenvalues: {eigenvalues}")
print()
print("--- Step 1: VQE (ground state) ---")
best_energy = float("inf")
best_params = None
for trial in range(5):
x0 = np.random.uniform(0, 2 * np.pi, n_params)
res = minimize(vqe_objective, x0, args=(n_qubits, n_layers, H), method="COBYLA")
if res.fun < best_energy:
best_energy = res.fun
best_params = res.x
gs_state = _state_from_hea(best_params, n_qubits, n_layers)
print(f"VQE ground state energy: {best_energy:.6f} (exact: {eigenvalues[0]:.6f})")
print()
# Step 2: VQD for first excited state
print("--- Step 2: VQD (first excited state) ---")
best_vqd_energy = float("inf")
best_vqd_params = None
for trial in range(5):
x0 = np.random.uniform(0, 2 * np.pi, n_params)
res = minimize(
vqd_objective,
x0,
args=(n_qubits, n_layers, H, [gs_state], penalty),
method="COBYLA",
)
if res.fun < best_vqd_energy:
best_vqd_energy = res.fun
best_vqd_params = res.x
es_state = _state_from_hea(best_vqd_params, n_qubits, n_layers)
es_energy = float(np.real(es_state.conj() @ H @ es_state))
overlap_with_gs = abs(np.dot(es_state.conj(), gs_state))**2
print(f"VQD first excited state energy: {es_energy:.6f} (exact: {eigenvalues[1]:.6f})")
print(f"Overlap with ground state: {overlap_with_gs:.6f}")
print()
print("✓ Run complete.")
def main():
parser = argparse.ArgumentParser(description="Variational Quantum Deflation (VQD)")
parser.add_argument("--n-qubits", type=int, default=2, help="Number of qubits (default: 2)")
parser.add_argument("--n-layers", type=int, default=2, help="Number of HEA layers (default: 2)")
parser.add_argument("--penalty", type=float, default=10.0, help="Overlap penalty β (default: 10.0)")
args = parser.parse_args()
np.random.seed(42)
run_vqd(args.n_qubits, args.n_layers, args.penalty)
if __name__ == "__main__":
main()
Stdout
=== VQD Example: H = Z₀ + ... + Z_{n-1} on 2 qubits ===
n_layers=2, penalty=10.0, n_params=4
Exact eigenvalues: [-2. 0. 0. 2.]
--- Step 1: VQE (ground state) ---
VQE ground state energy: -2.000000 (exact: -2.000000)
--- Step 2: VQD (first excited state) ---
VQD first excited state energy: -0.000000 (exact: 0.000000)
Overlap with ground state: 0.000000
✓ Run complete.