Add unit tests for PauliStringLCU block encoding

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)
+
+

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(