diff --git a/README.md b/README.md index 96b7613..0b5471d 100644 --- a/README.md +++ b/README.md @@ -69,7 +69,7 @@ Python Program"] - **Bernstein-Vazirani**: Determining hidden linear functions - **Deutsch-Jozsa**: Distinguishing constant from balanced functions - **Quantum Random Number Generator (QRNG)**: True quantum randomness -- **Iterative Amplitude Estimation**: Maximum-likelihood amplitude estimation for Monte Carlo-style workflows (canonical QPE-based variant planned) +- **Amplitude Estimation**: Canonical QPE-based amplitude estimation, plus the iterative (maximum-likelihood) variant for Monte Carlo-style workflows ### Quantum Information & States diff --git a/qpiai_quantum/algorithms/README.md b/qpiai_quantum/algorithms/README.md index b408d23..b338ad2 100644 --- a/qpiai_quantum/algorithms/README.md +++ b/qpiai_quantum/algorithms/README.md @@ -165,9 +165,9 @@ print(f"10 random values: {values}") ``` ### 9. Amplitude Estimation -Quantum speedup for Monte Carlo estimation. The SDK provides the **iterative -(maximum-likelihood) variant**; the canonical QPE-based variant is planned but -not yet implemented. +Quantum speedup for Monte Carlo estimation. The SDK provides both the +**canonical QPE-based variant** (`AmplitudeEstimation`) and the **iterative +(maximum-likelihood) variant** (`IterativeAmplitudeEstimation`). ```python from qpiai_quantum.algorithms.amplitude_estimation import ( @@ -193,10 +193,45 @@ amplitude = iae.estimate(problem, shots=2000) print(f"Estimated amplitude: {amplitude:.4f}") ``` -> **Note:** The canonical QPE-based `AmplitudeEstimation` class is not yet -> available and raises `NotImplementedError`. Use -> `IterativeAmplitudeEstimation` for currently supported amplitude estimation -> workflows. +The canonical QPE-based variant estimates the amplitude from the phase +measured on an evaluation register: + +```python +from qpiai_quantum.algorithms.amplitude_estimation import ( + AmplitudeEstimation, + EstimationProblem, +) +from qpiai_quantum.circuit import Circuit + +A = Circuit(1) +A.ry(0, 0.8) +problem = EstimationProblem(state_preparation=A, objective_qubits=[0]) + +# More evaluation qubits gives finer phase resolution; 5 or more is recommended. +qae = AmplitudeEstimation(num_evaluation_qubits=6) +amplitude = qae.estimate(problem, shots=1000, device_name="QpiAI-QSV-Local") +print(f"Estimated amplitude: {amplitude:.4f}") # true value: sin^2(0.4) = 0.1516 +``` + +`AmplitudeEstimation` honours a custom `is_good_state` predicate: it is +enumerated over the state register and synthesized into an explicit oracle. +That enumeration is supported for up to `AmplitudeEstimation.MAX_ORACLE_QUBITS` +state qubits; problems using the default `objective_qubits` marking have no +such limit. `IterativeAmplitudeEstimation` does not yet build its Grover +operator from the predicate, so custom predicates are not supported there: + +```python +# Bitstrings are MSB first, so bitstring[-1] is qubit 0. +B = Circuit(2) +B.ry(0, 0.8) +B.ry(1, 1.1) + +problem = EstimationProblem( + state_preparation=B, + objective_qubits=[0], + is_good_state=lambda bitstring: bitstring.count("1") == 1, # exactly one qubit set +) +``` ## Common Operations @@ -359,7 +394,7 @@ for n in [2, 3, 4, 5]: - **Complexity**: O(1/ε) queries vs O(1/ε²) classical Monte Carlo - **Use Cases**: Option pricing, risk analysis, counting problems - **Key Feature**: Quadratic speedup for estimating expectation values -- **Note**: Only the iterative (maximum-likelihood) variant is currently implemented. The canonical QPE-based variant is planned. +- **Note**: Both the canonical QPE-based variant (`AmplitudeEstimation`) and the iterative (maximum-likelihood) variant (`IterativeAmplitudeEstimation`) are implemented. Custom `is_good_state` predicates are currently supported by the canonical variant only. ## Performance Tips diff --git a/qpiai_quantum/algorithms/amplitude_estimation.py b/qpiai_quantum/algorithms/amplitude_estimation.py index 12fde2c..f6a2870 100644 --- a/qpiai_quantum/algorithms/amplitude_estimation.py +++ b/qpiai_quantum/algorithms/amplitude_estimation.py @@ -52,6 +52,7 @@ def __init__( self.state_preparation = state_preparation self.objective_qubits = objective_qubits self.is_good_state: Callable[[str], bool] + self._default_is_good_state: Callable[[str], bool] | None if is_good_state is None: def default_is_good_state(bitstring: str) -> bool: @@ -65,43 +66,255 @@ def default_is_good_state(bitstring: str) -> bool: return True self.is_good_state = default_is_good_state + self._default_is_good_state = default_is_good_state else: self.is_good_state = is_good_state + self._default_is_good_state = None @property def num_qubits(self) -> int: return self.state_preparation.num_qubits + @property + def has_custom_good_state(self) -> bool: + """ + True when is_good_state is anything other than the default marking. -class AmplitudeEstimation(QuantumAlgorithm): - """ - Canonical Quantum Amplitude Estimation (QAE) — **not yet implemented**. + Compared by identity rather than recorded at construction time, so + reassigning is_good_state afterwards is still honoured by the circuit + rather than only by the result counting. + """ + return self.is_good_state is not self._default_is_good_state - The canonical QPE-based variant is not yet available. Canonical QAE - requires constructing controlled Grover operators, which in turn need - a general controlled-subcircuit compiler not yet present in the SDK. - Use :class:`IterativeAmplitudeEstimation` for currently supported - amplitude estimation workflows. +class AmplitudeEstimation(QuantumAlgorithm): + """ + Canonical Quantum Amplitude Estimation (QAE) using Quantum Phase Estimation (QPE). + + QAE uses an evaluation register of m qubits initialized in superposition, + applies controlled powers of the Grover operator C-Q^(2^(m-1-j)) on + evaluation qubit j, and then applies Inverse QFT on the evaluation register + to estimate the amplitude a = sin^2(theta_a). + + Custom ``is_good_state`` predicates are supported: the predicate is + enumerated over the state register and synthesized into an explicit oracle + when it differs from the default all-objective-qubits-are-1 marking. That + enumeration is bounded by ``MAX_ORACLE_QUBITS``; problems relying on the + default marking are unaffected by the bound. """ + # Upper bound on the state register for which an explicit oracle can be + # synthesized from a custom is_good_state predicate (2^n enumeration). + # Only reached when a caller supplies their own predicate. + MAX_ORACLE_QUBITS = 16 + def __init__(self, num_evaluation_qubits: int): self.num_evaluation_qubits = num_evaluation_qubits super().__init__(num_qubits=0, name="Amplitude Estimation") self.description = "Canonical Amplitude Estimation using QPE" - def build_circuit(self, problem: EstimationProblem) -> Circuit: - raise NotImplementedError( - "Canonical QAE requires native implementation of controlled circuits." - ) + @staticmethod + def _apply_multi_controlled_z(circuit: Circuit, controls: list[int], target: int): + """Apply a phase flip on the all-ones state of controls + target.""" + if not controls: + circuit.z(target) + elif len(controls) == 1: + circuit.cz(controls[0], target) + elif len(controls) == 2: + circuit.h(target) + circuit.ccx(controls[0], controls[1], target) + circuit.h(target) + else: + circuit.h(target) + circuit.mcx(controls, target) + circuit.h(target) + + def _resolve_good_states(self, problem: EstimationProblem) -> list[int] | None: + """ + Resolve ``problem.is_good_state`` into the basis states it marks. + + Returns None when the predicate marks exactly the states covered by the + default "every objective qubit is 1" rule, so the compact + objective-qubit encoding can be used. Otherwise returns the integer + encoding of every marked basis state, where bit q is qubit q. + + Only a caller-supplied predicate has to be enumerated, so the state + register size is bounded on that path alone; the default marking is + emitted directly on the objective qubits at any register size. + """ + if not problem.has_custom_good_state: + return None + + n = problem.num_qubits + if n > self.MAX_ORACLE_QUBITS: + raise ValueError( + f"A custom is_good_state predicate can only be synthesized into " + f"an oracle for at most {self.MAX_ORACLE_QUBITS} state qubits, " + f"got {n}. Express the marked set through objective_qubits " + f"instead, which has no such limit." + ) + + objective_mask = 0 + for q in problem.objective_qubits: + objective_mask |= 1 << q + + good: list[int] = [] + default: list[int] = [] + for value in range(2**n): + # Same layout as measurement counts: MSB first, qubit q at index n-1-q. + if problem.is_good_state(format(value, f"0{n}b")): + good.append(value) + if value & objective_mask == objective_mask: + default.append(value) + + return None if good == default else good + + def _apply_controlled_s_chi( + self, + circuit: Circuit, + control_qubit: int, + state_offset: int, + problem: EstimationProblem, + good_states: list[int] | None = None, + ): + """ + Apply controlled phase-flip S_chi on good states. + + With ``good_states`` None the default "all objective qubits are 1" + marking is emitted directly on the objective qubits. Otherwise one + phase flip is emitted per marked basis state, so custom + ``is_good_state`` predicates are honoured rather than silently ignored. + """ + if good_states is None: + obj_qubits = [state_offset + q for q in problem.objective_qubits] + self._apply_multi_controlled_z( + circuit, [control_qubit] + obj_qubits[:-1], obj_qubits[-1] + ) + return + + n_state = problem.num_qubits + state_qubits = [state_offset + q for q in range(n_state)] + controls = [control_qubit] + state_qubits[:-1] + target = state_qubits[-1] + + for value in good_states: + zeros = [state_qubits[q] for q in range(n_state) if not (value >> q) & 1] + for q in zeros: + circuit.x(q) + self._apply_multi_controlled_z(circuit, controls, target) + for q in zeros: + circuit.x(q) + + def _apply_controlled_s_0( + self, + circuit: Circuit, + control_qubit: int, + state_offset: int, + n_state_qubits: int, + ): + """Apply controlled phase-flip S_0 on zero state |0...0>.""" + state_qubits = [state_offset + q for q in range(n_state_qubits)] - def estimate(self, problem: EstimationProblem) -> float: - raise NotImplementedError( - "Canonical QAE requires native implementation of controlled sub-circuits, " - "which are not yet fully supported by the basic Circuit builder. " - "Please use IterativeAmplitudeEstimation instead." + for q in state_qubits: + circuit.x(q) + + self._apply_multi_controlled_z( + circuit, [control_qubit] + state_qubits[:-1], state_qubits[-1] ) + for q in state_qubits: + circuit.x(q) + + def _apply_controlled_grover( + self, + circuit: Circuit, + control_qubit: int, + state_offset: int, + problem: EstimationProblem, + power: int, + good_states: list[int] | None = None, + ): + r""" + Apply controlled Grover operator Q^power controlled by control_qubit. + + Controlled-Q = A · (Controlled-S_0) · A^\dagger · (Controlled-S_chi) + """ + n_state = problem.num_qubits + state_qubits = [state_offset + q for q in range(n_state)] + A = problem.state_preparation + A_inv = A.inverse() + + for _ in range(power): + circuit.z(control_qubit) + self._apply_controlled_s_chi( + circuit, control_qubit, state_offset, problem, good_states + ) + circuit.compose(A_inv, state_qubits) + self._apply_controlled_s_0(circuit, control_qubit, state_offset, n_state) + circuit.compose(A, state_qubits) + + def build_circuit(self, problem: EstimationProblem) -> Circuit: + m = self.num_evaluation_qubits + n = problem.num_qubits + total_qubits = m + n + + self.num_qubits = total_qubits + self.circuit = Circuit(total_qubits, m) + + # 1. Initialize evaluation register to |+> + for i in range(m): + self.circuit.h(i) + + # 2. Initialize state register with A + state_qubits = [m + q for q in range(n)] + self.circuit.compose(problem.state_preparation, state_qubits) + + # 3. Apply controlled Grover powers C-Q^(2^(m-1-j)). + # Evaluation qubit j must accumulate the phase weight the inverse QFT + # expects at that position, matching QuantumPhaseEstimation. + good_states = self._resolve_good_states(problem) + for j in range(m): + power = 2 ** (m - 1 - j) + self._apply_controlled_grover( + self.circuit, j, m, problem, power, good_states + ) + + # 4. Apply Inverse QFT on evaluation register + QFT.apply_inverse_qft_to_circuit(self.circuit, start=0, n=m) + + # 5. Measure evaluation register + for i in range(m): + self.circuit.measure(i, i) + + return self.circuit + + def estimate( + self, + problem: EstimationProblem, + shots: int = 1024, + device_name: str = "QpiAI-QSV-Local", + ) -> float: + """ + Estimate amplitude a = sin^2(theta_a) using Canonical QAE. + """ + m = self.num_evaluation_qubits + circ = self.build_circuit(problem) + result = circ.run(shots=shots, device_name=device_name) + counts = result.get()["counts"] + + if not counts: + return 0.0 + + total_shots = sum(counts.values()) + weighted_amplitude = 0.0 + + for bitstring, count in counts.items(): + y = int(bitstring, 2) + theta = math.pi * y / (2**m) + weighted_amplitude += (count / total_shots) * (math.sin(theta) ** 2) + + return float(weighted_amplitude) + class IterativeAmplitudeEstimation(QuantumAlgorithm): """ diff --git a/tests/algorithms/test_amplitude_estimation.py b/tests/algorithms/test_amplitude_estimation.py index a462878..5f97cfb 100644 --- a/tests/algorithms/test_amplitude_estimation.py +++ b/tests/algorithms/test_amplitude_estimation.py @@ -58,13 +58,31 @@ def make_mock_result(counts): assert abs(estimated_prob - expected_prob) < 0.05 -def test_canonical_qae_not_implemented(): +def test_canonical_qae_build_circuit(): circuit = Circuit(1) + circuit.ry(0, 0.6) problem = EstimationProblem(state_preparation=circuit, objective_qubits=[0]) qae = AmplitudeEstimation(num_evaluation_qubits=3) - with pytest.raises(NotImplementedError): - qae.estimate(problem) + built_circuit = qae.build_circuit(problem) + # Total qubits = 3 evaluation + 1 state qubit = 4 + assert built_circuit.num_qubits == 4 + assert built_circuit.num_clbits == 3 + + +def test_canonical_qae_estimation_local(): + theta = 0.8 + expected_prob = math.sin(theta / 2) ** 2 + + circuit = Circuit(1) + circuit.ry(0, theta) + problem = EstimationProblem(state_preparation=circuit, objective_qubits=[0]) + + qae = AmplitudeEstimation(num_evaluation_qubits=6) + estimated_prob = qae.estimate(problem, shots=1000, device_name="QpiAI-QSV-Local") + + # Allow 0.05 tolerance with 6 evaluation qubits + assert abs(estimated_prob - expected_prob) < 0.05 @pytest.mark.skipif( @@ -82,3 +100,171 @@ def test_amplitude_estimation_correctness(): iae = IterativeAmplitudeEstimation(epsilon_target=0.08, alpha=0.05) estimated_prob = iae.estimate(problem, shots=200, device_name="QpiAI-QSV-Local") assert abs(estimated_prob - expected_prob) < 0.15 + + +# Regression coverage for the canonical QAE controlled-power ordering. +# Evaluation qubit j must receive Q^(2^(m-1-j)) to match the inverse-QFT +# convention used by QuantumPhaseEstimation. With the reversed assignment +# theta=0.8 still happens to pass, so the sweep below is what actually pins +# the convention down: theta=0.2 returns ~0.378 instead of ~0.010. +@pytest.mark.parametrize("theta", [0.2, 0.5, 0.8, 1.5, 2.2, 2.8]) +def test_canonical_qae_amplitude_sweep(theta): + expected_prob = math.sin(theta / 2) ** 2 + + circuit = Circuit(1) + circuit.ry(0, theta) + problem = EstimationProblem(state_preparation=circuit, objective_qubits=[0]) + + qae = AmplitudeEstimation(num_evaluation_qubits=6) + estimated_prob = qae.estimate(problem, shots=4000, device_name="QpiAI-QSV-Local") + + assert abs(estimated_prob - expected_prob) < 0.05 + + +@pytest.mark.parametrize("theta,expected_prob", [(0.0, 0.0), (math.pi, 1.0)]) +def test_canonical_qae_boundary_amplitudes(theta, expected_prob): + circuit = Circuit(1) + circuit.ry(0, theta) + problem = EstimationProblem(state_preparation=circuit, objective_qubits=[0]) + + qae = AmplitudeEstimation(num_evaluation_qubits=6) + estimated_prob = qae.estimate(problem, shots=4000, device_name="QpiAI-QSV-Local") + + assert abs(estimated_prob - expected_prob) < 0.05 + + +def _two_qubit_problem(is_good_state): + """RY(0.8) on qubit 0, RY(1.1) on qubit 1, with a custom predicate.""" + circuit = Circuit(2) + circuit.ry(0, 0.8) + circuit.ry(1, 1.1) + return EstimationProblem( + state_preparation=circuit, + objective_qubits=[0], + is_good_state=is_good_state, + ) + + +def test_canonical_qae_honours_custom_good_state(): + # Counts layout is MSB first, so bitstring[0] is qubit 1 and bitstring[1] + # is qubit 0. "exactly one qubit is 1" cannot be expressed as the default + # all-objective-qubits-are-1 marking, so this fails outright if + # is_good_state is ignored (it would estimate P(qubit 0 = 1) = 0.1516). + p0 = math.sin(0.8 / 2) ** 2 + p1 = math.sin(1.1 / 2) ** 2 + expected_prob = p0 * (1 - p1) + (1 - p0) * p1 + + problem = _two_qubit_problem(lambda bitstring: bitstring.count("1") == 1) + qae = AmplitudeEstimation(num_evaluation_qubits=6) + estimated_prob = qae.estimate(problem, shots=4000, device_name="QpiAI-QSV-Local") + + assert abs(estimated_prob - expected_prob) < 0.05 + + +def test_canonical_qae_custom_good_state_marking_zeros(): + p0 = math.sin(0.8 / 2) ** 2 + p1 = math.sin(1.1 / 2) ** 2 + expected_prob = (1 - p0) * (1 - p1) + + problem = _two_qubit_problem(lambda bitstring: bitstring == "00") + qae = AmplitudeEstimation(num_evaluation_qubits=6) + estimated_prob = qae.estimate(problem, shots=4000, device_name="QpiAI-QSV-Local") + + assert abs(estimated_prob - expected_prob) < 0.05 + + +def test_resolve_good_states_uses_compact_marking_by_default(): + qae = AmplitudeEstimation(num_evaluation_qubits=3) + circuit = Circuit(2) + circuit.ry(0, 0.8) + circuit.ry(1, 1.1) + + default = EstimationProblem(state_preparation=circuit, objective_qubits=[0, 1]) + assert qae._resolve_good_states(default) is None + + # A custom predicate that happens to match the default must not pay for an + # explicitly synthesised oracle either. + equivalent = EstimationProblem( + state_preparation=circuit, + objective_qubits=[0], + is_good_state=lambda bitstring: bitstring[-1] == "1", + ) + assert qae._resolve_good_states(equivalent) is None + + +def test_resolve_good_states_enumerates_custom_predicate(): + qae = AmplitudeEstimation(num_evaluation_qubits=3) + circuit = Circuit(2) + + problem = EstimationProblem( + state_preparation=circuit, + objective_qubits=[0], + is_good_state=lambda bitstring: bitstring == "10", + ) + # "10" is qubit 1 = 1 and qubit 0 = 0, i.e. integer 0b10. + assert qae._resolve_good_states(problem) == [0b10] + + +class _OversizedStatePreparation: + num_qubits = AmplitudeEstimation.MAX_ORACLE_QUBITS + 1 + + +def test_canonical_qae_rejects_oversized_custom_predicate(): + """Only a custom predicate needs enumerating, so only it is size limited.""" + problem = EstimationProblem( + state_preparation=_OversizedStatePreparation(), + objective_qubits=[0], + is_good_state=lambda bitstring: bitstring.count("1") == 1, + ) + qae = AmplitudeEstimation(num_evaluation_qubits=3) + + with pytest.raises(ValueError, match="at most"): + qae._resolve_good_states(problem) + + +def test_canonical_qae_allows_oversized_default_marking(): + """The default objective-qubit marking is emitted directly, so a large + state register must not be rejected and must not be enumerated.""" + problem = EstimationProblem( + state_preparation=_OversizedStatePreparation(), objective_qubits=[0] + ) + qae = AmplitudeEstimation(num_evaluation_qubits=3) + + assert qae._resolve_good_states(problem) is None + + +@pytest.mark.parametrize( + "is_good_state,expected_prob", + [ + (lambda bitstring: False, 0.0), + (lambda bitstring: True, 1.0), + ], +) +def test_canonical_qae_degenerate_good_state(is_good_state, expected_prob): + """A predicate marking no states or every state must still be exact.""" + problem = _two_qubit_problem(is_good_state) + qae = AmplitudeEstimation(num_evaluation_qubits=6) + estimated_prob = qae.estimate(problem, shots=4000, device_name="QpiAI-QSV-Local") + + assert abs(estimated_prob - expected_prob) < 0.05 + + +def test_good_state_reassignment_is_picked_up_by_the_circuit(): + """Reassigning is_good_state must reach the oracle, not just the counting. + + Otherwise the circuit and the final calculation would use different + definitions of a good state. + """ + circuit = Circuit(2) + circuit.ry(0, 0.8) + circuit.ry(1, 1.1) + problem = EstimationProblem(state_preparation=circuit, objective_qubits=[0]) + qae = AmplitudeEstimation(num_evaluation_qubits=3) + + assert problem.has_custom_good_state is False + assert qae._resolve_good_states(problem) is None + + problem.is_good_state = lambda bitstring: bitstring == "10" + + assert problem.has_custom_good_state is True + assert qae._resolve_good_states(problem) == [0b10]