diff --git a/analysis/tests/test_pauli_string_lcu.py b/analysis/tests/test_pauli_string_lcu.py new file mode 100644 index 00000000..6062e766 --- /dev/null +++ b/analysis/tests/test_pauli_string_lcu.py @@ -0,0 +1,143 @@ +""" +Test construction of Pauli String LCU block encodings. +""" + +import numpy as np +import pytest + +from pyLIQTR.ProblemInstances.ProblemInstance import ProblemInstance + +from qhat.analysis.config_types import ( + GeneralConfiguration, + GeneralConfigurationUser, +) +from qhat.analysis.hamiltonian import ( + Hamiltonian, + LinearCombinationOfPauliStrings, +) +from qhat.analysis.unitary import PauliStringLCU, PyLIQTRPauliStringLCU + + +def test_pauli_string_lcu1(): + """Test converting a small, simple Hamiltonian to Pauli String LCU + block encoding.""" + # Create a simple 1-qubit Hamiltonian + pauli_data = { + 'X': 1.0, + 'Z': 1.0, + } + lcps = LinearCombinationOfPauliStrings(num_qubits=1, dense=pauli_data) + H = Hamiltonian(lcps) + + # Convert Hamiltonian to matrix + H_matrix = H.to_matrix(memory_threshold_gb=1.0) + + # Create unitary operator and convert to matrix + unitaryop = PauliStringLCU(H, 'AS', probability_eps=0.5) + unitarymx = unitaryop.tensor_contract() + + # Verify that the upper corner of our unitary is equal to the + # original Hamiltonian matrix, scaled + alpha = np.sum([v for v in H.get_all_pauli_strings().values()]) + scale = 1. / alpha + np.testing.assert_array_almost_equal(unitarymx[0:2,0:2], scale * H_matrix) + + +def test_pauli_string_lcu2(): + """Test converting a small, slightly more complex Hamiltonian to + Pauli String LCU block encoding.""" + # Create a simple 2-qubit Hamiltonian + pauli_data = { + 'XX': 1.0, + 'YZ': 1.0, + 'ZY': 2.0, + } + lcps = LinearCombinationOfPauliStrings(num_qubits=2, dense=pauli_data) + H = Hamiltonian(lcps) + + # Convert Hamiltonian to matrix + H_matrix = H.to_matrix(memory_threshold_gb=1.0) + + # Create unitary operator and convert to matrix + unitaryop = PauliStringLCU(H, 'AS', probability_eps=0.1) + unitarymx = unitaryop.tensor_contract() + + # Verify that the upper corner of our unitary is equal to the + # original Hamiltonian matrix, scaled + alpha = np.sum([v for v in H.get_all_pauli_strings().values()]) + scale = 1. / alpha + np.testing.assert_array_almost_equal(unitarymx[0:4,0:4], scale * H_matrix) + + +# Define a wrapper class to make a QHAT PauliString Hamiltonian behave +# like a PyLIQTR ProblemInstance +class PauliStringInstance(Hamiltonian, ProblemInstance): + def __init__(self, hamiltonian): + super().__init__(hamiltonian) + + def __str__(self): + return str(self.get_all_pauli_strings(return_as='strings')) + + def n_terms(self, **kwargs): + return len(self.get_all_pauli_strings()) + + def n_qubits(self): + return self.num_qubits() + + def get_alpha(self): + return np.sum([v for v in self.get_all_pauli_strings().values()]) + + def yield_PauliLCU_Info(self, do_pad=0, return_as='strings'): + for t in self.get_all_pauli_strings(return_as=return_as).items(): + yield t + + +def test_pauli_string_lcu_pyliqtr1(): + """Test converting a small, simple Hamiltonian to Pauli String LCU + block encoding.""" + # Create a simple 1-qubit Hamiltonian + pauli_data = { + 'X': 1.0, + 'Z': 1.0, + } + lcps = LinearCombinationOfPauliStrings(num_qubits=1, dense=pauli_data) + inst = PauliStringInstance(lcps) + + # Convert Hamiltonian to matrix + H_matrix = inst.to_matrix(memory_threshold_gb=1.0) + + # Create unitary operator and convert to matrix + unitaryop = PyLIQTRPauliStringLCU(inst, 'AS', probability_eps=0.5) + unitarymx = unitaryop.tensor_contract() + + # Verify that the upper corner of our unitary is equal to the + # original Hamiltonian matrix, scaled + scale = 1. / inst.get_alpha() + np.testing.assert_array_almost_equal(unitarymx[0:2,0:2], scale * H_matrix) + + +def test_pauli_string_lcu_pyliqtr2(): + """Test converting a small, slightly more complex Hamiltonian to + Pauli String LCU block encoding.""" + # Create a simple 2-qubit Hamiltonian + pauli_data = { + 'XX': 1.0, + 'YZ': 1.0, + 'ZY': 2.0, + } + lcps = LinearCombinationOfPauliStrings(num_qubits=2, dense=pauli_data) + inst = PauliStringInstance(lcps) + + # Convert Hamiltonian to matrix + H_matrix = inst.to_matrix(memory_threshold_gb=1.0) + + # Create unitary operator and convert to matrix + unitaryop = PyLIQTRPauliStringLCU(inst, 'AS', probability_eps=0.1) + unitarymx = unitaryop.tensor_contract() + + # Verify that the upper corner of our unitary is equal to the + # original Hamiltonian matrix, scaled + scale = 1. / inst.get_alpha() + np.testing.assert_array_almost_equal(unitarymx[0:4,0:4], scale * H_matrix) + + diff --git a/analysis/unitary.py b/analysis/unitary.py index 9b132862..817adfb0 100644 --- a/analysis/unitary.py +++ b/analysis/unitary.py @@ -1,15 +1,22 @@ import logging import math +from typing import Dict import cirq import numpy as np +from qualtran import ( + Bloq, + BloqBuilder, + SoquetT, +) +from qualtran._infra.registers import Signature from qualtran.bloqs.block_encoding import LCUBlockEncoding from qualtran.bloqs.multiplexers.select_pauli_lcu import SelectPauliLCU from qualtran.bloqs.state_preparation import StatePreparationAliasSampling from pyLIQTR.BlockEncodings.DoubleFactorized import DoubleFactorized from pyLIQTR.BlockEncodings.LinearT import Fermionic_LinearT -from pyLIQTR.BlockEncodings.PauliStringLCU import PauliStringLCU as PyLIQTRPauliStringLCU +from pyLIQTR.BlockEncodings.PauliStringLCU import PauliStringLCU as PyLIQTRPauliStringLCU_orig from pyLIQTR.ProblemInstances.ChemicalHamiltonian import ChemicalHamiltonian from qhat.analysis.config_types import UnitaryConfiguration @@ -36,10 +43,8 @@ def __init__(self, hamiltonian, prepare_type=None, probability_eps=0.002, **kwar n_tot = 2**(int(np.ceil(np.log2(n_terms)))) n_pad = n_tot - n_terms - alphas = [np.sqrt(np.abs(t.coefficient)) for t in pauli_terms] - alpha = np.sum([a**2 for a in alphas]) - alphas_scaled = [a/np.sqrt(alpha) for a in alphas] - alphas_scaled.extend([0.0 for i in range(n_pad)]) + weights = [np.abs(t.coefficient) for t in pauli_terms] + alpha = np.sum(weights) selection_bitsize = int(np.ceil(np.log2(n_tot))) @@ -54,23 +59,13 @@ def __init__(self, hamiltonian, prepare_type=None, probability_eps=0.002, **kwar # see https://github.com/quantumlib/Qualtran/issues/1045 #prepare = StatePreparationViaRotations( # phase_bitsize = 4, - # state_coefficients = alphas_scaled, + # state_coefficients = weights, # ) raise NotImplementedError("PauliStringLCU") elif prepare_type=='AS': - for t in pauli_terms: - coeff = np.real(t.coefficient) - if coeff < 0: - logger.warning("Alias sampling preparation with negative coefficients is not " - "supported yet. Circuits and estimates will assume positive " - "coefficients.") - break - prepare = StatePreparationAliasSampling.from_lcu_probs( - lcu_probabilities=[ - np.abs(np.real(t.coefficient)) for t in pauli_terms - ], + lcu_probabilities=weights, probability_epsilon=probability_eps, ) @@ -88,6 +83,26 @@ def _select_gate(self): def _prepare_gate(self): return self.prepare + # Backport corrected LCUBlockEncoding method from v0.5.0 + def build_composite_bloq(self, bb: 'BloqBuilder', **soqs: SoquetT) -> Dict[str, 'SoquetT']: + def _extract_soqs(bloq: Bloq) -> Dict[str, 'SoquetT']: + return {reg.name: soqs.pop(reg.name) for reg in bloq.signature.lefts()} + + soqs |= bb.add_d(self.prepare, **_extract_soqs(self.prepare)) + soqs |= bb.add_d(self.select, **_extract_soqs(self.select)) + soqs |= bb.add_d(self.prepare.adjoint(), **_extract_soqs(self.prepare.adjoint())) + return soqs + +# ------------------------------------------------------------------------------------------------- + +# add bugfix: correct ordering of Signature registers +class PyLIQTRPauliStringLCU(PyLIQTRPauliStringLCU_orig): + @property + def signature(self): + return Signature( + [*self.control_registers, *self.selection_registers, + *self.junk_registers, *self.target_registers] ) + # ------------------------------------------------------------------------------------------------- def encode_linear_t(