Compute PauliStringLCU tolerances from energy error

analysis/algorithm.py

@@ -136,9 +136,15 @@ def build_qpe_qualtran_qubitized(
     logger.verbose(
             "Build a QPE algorithm with Qualtran's \"QubitizationQPE\" method.")
 
+    P = config_algorithm.num_phase_qubits
+    if P is None:
+        dE = config_algorithm.energy_error
+        alpha = unitary.alpha
+        P = int(math.ceil(math.log2(math.pi * alpha / (2 * dE))))
+
     return NewQubitizationQPE(QubitizationWalkOperator(unitary._select_gate,
                                                        unitary._prepare_gate),
-                              config_algorithm.num_phase_qubits)
+                              P)
 
 # -------------------------------------------------------------------------------------------------
 
@@ -149,9 +155,15 @@ def build_qpe_pyliqtr_qubitized(
     logger.verbose(
             "Build a QPE algorithm with pyLIQTR's \"QubitizedPhaseEstimation\" method.")
 
+    P = config_algorithm.num_phase_qubits
+    if P is None:
+        dE = config_algorithm.energy_error
+        alpha = unitary.alpha
+        P = int(math.ceil(math.log2(math.pi * alpha / (2 * dE))))
+
     # TODO: The name and signature suggest that this may _only_ be valid for block-encoded
     #       unitaries.  Is that true?
-    return QubitizedPhaseEstimation(block_encoding=unitary, prec=config_algorithm.num_phase_qubits)
+    return QubitizedPhaseEstimation(block_encoding=unitary, prec=P)
 
 # -------------------------------------------------------------------------------------------------
 

analysis/config.py

@@ -73,7 +73,7 @@
 elif my_method == "pauli-lcu":
     algorithm.method = "QPE: pyLIQTR qubitized"
 #    algorithm.method = "QPE: qualtran qubitization"
-    algorithm.num_phase_qubits = 12
+    algorithm.energy_error = 0.5 * energy_error
 elif my_method == "double-factorization":
     algorithm.method = "QPE: pyLIQTR qubitized"
     algorithm.num_phase_qubits = 12

analysis/tests/test_pauli_string_lcu.py

@@ -33,7 +33,8 @@ def test_pauli_string_lcu1():
     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)
+    # energy error = 5.0 => probability_eps ~= 0.44
+    unitaryop = PauliStringLCU(H, 'AS', energy_error=5.0)
     unitarymx = unitaryop.tensor_contract()
 
     # Verify that the upper corner of our unitary is equal to the
@@ -59,7 +60,8 @@ def test_pauli_string_lcu2():
     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)
+    # energy error = 5.0 => probability_eps ~= 0.10
+    unitaryop = PauliStringLCU(H, 'AS', energy_error=5.0)
     unitarymx = unitaryop.tensor_contract()
 
     # Verify that the upper corner of our unitary is equal to the
@@ -78,7 +80,7 @@ def __init__(self, hamiltonian):
     def __str__(self):
         return str(self.get_all_pauli_strings(return_as='strings'))
 
-    def n_terms(self, **kwargs):
+    def n_terms(self, encoding, **kwargs):
         return len(self.get_all_pauli_strings())
 
     def n_qubits(self):
@@ -107,7 +109,8 @@ def test_pauli_string_lcu_pyliqtr1():
     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)
+    # energy error = 5.0 => probability_eps ~= 0.44
+    unitaryop = PyLIQTRPauliStringLCU(inst, 'AS', energy_error=5.0)
     unitarymx = unitaryop.tensor_contract()
 
     # Verify that the upper corner of our unitary is equal to the
@@ -132,7 +135,8 @@ def test_pauli_string_lcu_pyliqtr2():
     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)
+    # energy error = 5.0 => probability_eps ~= 0.10
+    unitaryop = PyLIQTRPauliStringLCU(inst, 'AS', energy_error=5.0)
     unitarymx = unitaryop.tensor_contract()
 
     # Verify that the upper corner of our unitary is equal to the

analysis/unitary.py

@@ -14,6 +14,7 @@
 from qualtran.bloqs.multiplexers.select_pauli_lcu import SelectPauliLCU
 from qualtran.bloqs.state_preparation import StatePreparationAliasSampling
 
+from pyLIQTR.BlockEncodings import VALID_ENCODINGS
 from pyLIQTR.BlockEncodings.DoubleFactorized import DoubleFactorized
 from pyLIQTR.BlockEncodings.LinearT import Fermionic_LinearT
 from pyLIQTR.BlockEncodings.PauliStringLCU import PauliStringLCU as PyLIQTRPauliStringLCU_orig
@@ -29,7 +30,7 @@
 
 class PauliStringLCU(LCUBlockEncoding):
 
-    def __init__(self, hamiltonian, prepare_type=None, probability_eps=0.002, **kwargs):
+    def __init__(self, hamiltonian, prepare_type=None, energy_error=None, **kwargs):
 
         pauli_terms = [cirq.DensePauliString(term[0],coefficient=term[1])
                           for term in hamiltonian.get_all_pauli_strings(return_as='strings').items()
@@ -48,6 +49,8 @@ def __init__(self, hamiltonian, prepare_type=None, probability_eps=0.002, **kwar
 
         selection_bitsize = int(np.ceil(np.log2(n_tot)))
 
+        eps = energy_error / (2 * np.sqrt(2) * alpha * n_terms)
+
         select = SelectPauliLCU(
                      selection_bitsize=selection_bitsize,
                      target_bitsize=hamiltonian.num_qubits(),
@@ -66,11 +69,11 @@ def __init__(self, hamiltonian, prepare_type=None, probability_eps=0.002, **kwar
         elif prepare_type=='AS':
             prepare = StatePreparationAliasSampling.from_lcu_probs(
                           lcu_probabilities=weights,
-                          probability_epsilon=probability_eps,
+                          probability_epsilon=eps,
                       )
 
         super().__init__(
-            alpha=alpha, epsilon=probability_eps,
+            alpha=alpha, epsilon=eps,
             select=select, prepare=prepare
         )
 
@@ -95,8 +98,16 @@ def _extract_soqs(bloq: Bloq) -> Dict[str, 'SoquetT']:
 
 # -------------------------------------------------------------------------------------------------
 
-# add bugfix:  correct ordering of Signature registers
+# wrapper to add bugfixes to PyLIQTR's class
 class PyLIQTRPauliStringLCU(PyLIQTRPauliStringLCU_orig):
+    # use energy_error to determine probability_eps
+    def __init__(self,ProblemInstance, prepare_type=None, energy_error=None, **kwargs):
+        alpha = ProblemInstance.get_alpha()
+        n_terms = ProblemInstance.n_terms(VALID_ENCODINGS.PauliLCU)
+        eps = energy_error / (2 * np.sqrt(2) * alpha * n_terms)
+        super().__init__(ProblemInstance, prepare_type, probability_eps=eps, **kwargs)
+
+    # correct ordering of Signature registers
     @property
     def signature(self):
         return Signature(
@@ -160,7 +171,7 @@ def encode_pauli_lcu_pyliqtr(
             ProblemInstance=problem_instance,
             # TODO: Setting `prepare_type='AS'` might give the alias sampling Obenland mentioned?
             prepare_type='AS',
-            # TODO: What to do with energy_error?
+            energy_error=config_unitary.energy_error
             )
 
 # -------------------------------------------------------------------------------------------------
@@ -175,7 +186,7 @@ def encode_pauli_lcu_qualtran(
     return PauliStringLCU(
             hamiltonian=hamiltonian,
             prepare_type='AS',
-            # TODO: What to do with energy_error?
+            energy_error=config_unitary.energy_error
             )
 
 # -------------------------------------------------------------------------------------------------