QubitizationWalkOperator add attribute: sum of LCU coefficients (#890)

* `QubitizationWalkOperator` attribute: `sum_of_lcu_coefficients`

defaults to `None`

* add property to `PrepareOracle` instead

* ising: return true LCU coeffs

* fix symbolic ham sim example

* walk op docstrings

* fix docstring nits

Author: anurudhp

qualtran/bloqs/chemistry/ising/hamiltonian.py

@@ -90,5 +90,4 @@ def get_1d_ising_lcu_coeffs(
     spins = cirq.LineQubit.range(n_spins)
     ham = get_1d_ising_hamiltonian(spins, j_zz_strength, gamma_x_strength)
     coeffs = np.array([term.coefficient.real for term in ham])
-    lcu_coeffs = coeffs / np.sum(abs(coeffs))
-    return lcu_coeffs
+    return coeffs

qualtran/bloqs/for_testing/qubitization_walk_test.py

@@ -12,18 +12,19 @@
 #  See the License for the specific language governing permissions and
 #  limitations under the License.
 from functools import cached_property
-from typing import Tuple
+from typing import Optional, Tuple
 
 import attrs
 import cirq
 import scipy
-from numpy._typing import NDArray
+from numpy.typing import NDArray
 
 from qualtran import Signature
 from qualtran.bloqs.multiplexers.select_pauli_lcu import SelectPauliLCU
 from qualtran.bloqs.qubitization_walk_operator import QubitizationWalkOperator
 from qualtran.bloqs.select_and_prepare import PrepareOracle
 from qualtran.bloqs.state_preparation import PrepareUniformSuperposition
+from qualtran.resource_counting.symbolic_counting_utils import SymbolicFloat
 
 
 @attrs.frozen
@@ -32,6 +33,7 @@ class PrepareUniformSuperpositionTest(PrepareOracle):
     cvs: Tuple[int, ...] = attrs.field(
         converter=lambda v: (v,) if isinstance(v, int) else tuple(v), default=()
     )
+    qlambda: Optional[float] = None
 
     @cached_property
     def selection_registers(self) -> Signature:
@@ -41,6 +43,10 @@ def selection_registers(self) -> Signature:
     def junk_registers(self) -> Signature:
         return Signature.build()
 
+    @cached_property
+    def l1_norm_of_coeffs(self) -> Optional['SymbolicFloat']:
+        return self.qlambda
+
     def decompose_from_registers(
         self, *, context: cirq.DecompositionContext, **quregs: NDArray[cirq.Qid]  # type: ignore[type-var]
     ) -> cirq.OP_TREE:
@@ -57,7 +63,7 @@ def get_uniform_pauli_qubitized_walk(target_bitsize: int):
     ham_dps = [ps.dense(q) for ps in ham]
 
     assert scipy.linalg.ishermitian(ham.matrix())
-    prepare = PrepareUniformSuperpositionTest(len(ham_coeff))
+    prepare = PrepareUniformSuperpositionTest(len(ham_coeff), qlambda=sum(ham_coeff))
     select = SelectPauliLCU(
         (len(ham_coeff) - 1).bit_length(), select_unitaries=ham_dps, target_bitsize=target_bitsize
     )

qualtran/bloqs/for_testing/random_select_and_prepare.py

@@ -56,6 +56,10 @@ def __attrs_post_init__(self):
     def selection_registers(self) -> tuple[Register, ...]:
         return (Register('selection', BoundedQUInt(bitsize=self.U.bitsize)),)
 
+    @property
+    def l1_norm_of_coeffs(self) -> float:
+        return 1.0
+
     @classmethod
     def random(cls, bitsize: int, *, random_state: np.random.RandomState):
         """Generate a random unitary s.t. the first column has all real amplitudes"""

qualtran/bloqs/hamiltonian_simulation/hamiltonian_simulation_by_gqsp.ipynb

@@ -85,7 +85,6 @@
     "#### Parameters\n",
     " - `walk_operator`: qubitization walk operator of $H$ constructed from SELECT and PREPARE oracles.\n",
     " - `t`: time to simulate the Hamiltonian, i.e. $e^{-iHt}$\n",
-    " - `alpha`: the $1$-norm of the coefficients of the unitaries comprising the Hamiltonian $H$.\n",
     " - `precision`: the precision $\\epsilon$ to approximate $e^{it\\cos\\theta}$ to a polynomial.\n"
    ]
   },
@@ -123,9 +122,7 @@
     "from qualtran.bloqs.hubbard_model import get_walk_operator_for_hubbard_model\n",
     "\n",
     "walk_op = get_walk_operator_for_hubbard_model(2, 2, 1, 1)\n",
-    "hubbard_time_evolution_by_gqsp = HamiltonianSimulationByGQSP(\n",
-    "    walk_op, t=5, alpha=1, precision=1e-7\n",
-    ")"
+    "hubbard_time_evolution_by_gqsp = HamiltonianSimulationByGQSP(walk_op, t=5, precision=1e-7)"
    ]
   },
   {
@@ -142,10 +139,9 @@
     "from qualtran.bloqs.hubbard_model import get_walk_operator_for_hubbard_model\n",
     "\n",
     "walk_op = get_walk_operator_for_hubbard_model(2, 2, 1, 1)\n",
-    "t, alpha, inv_eps = sympy.symbols(\"t alpha N\")\n",
-    "symbolic_hamsim_by_gqsp = HamiltonianSimulationByGQSP(\n",
-    "    walk_op, t=t, alpha=alpha, precision=1 / inv_eps\n",
-    ")"
+    "\n",
+    "t, inv_eps = sympy.symbols(\"t N\")\n",
+    "symbolic_hamsim_by_gqsp = HamiltonianSimulationByGQSP(walk_op, t=t, precision=1 / inv_eps)"
    ]
   },
   {
@@ -191,8 +187,8 @@
    },
    "outputs": [],
    "source": [
-    "from qualtran.resource_counting.generalizers import generalize_rotation_angle, ignore_split_join, ignore_alloc_free\n",
-    "hubbard_time_evolution_by_gqsp_g, hubbard_time_evolution_by_gqsp_sigma = hubbard_time_evolution_by_gqsp.call_graph(generalizer=[ignore_split_join, ignore_alloc_free, generalize_rotation_angle])\n",
+    "from qualtran.resource_counting.generalizers import ignore_split_join\n",
+    "hubbard_time_evolution_by_gqsp_g, hubbard_time_evolution_by_gqsp_sigma = hubbard_time_evolution_by_gqsp.call_graph(max_depth=1, generalizer=ignore_split_join)\n",
     "show_call_graph(hubbard_time_evolution_by_gqsp_g)\n",
     "show_counts_sigma(hubbard_time_evolution_by_gqsp_sigma)"
    ]

qualtran/bloqs/hamiltonian_simulation/hamiltonian_simulation_by_gqsp.py

@@ -88,18 +88,26 @@ class HamiltonianSimulationByGQSP(GateWithRegisters):
     Args:
         walk_operator: qubitization walk operator of $H$ constructed from SELECT and PREPARE oracles.
         t: time to simulate the Hamiltonian, i.e. $e^{-iHt}$
-        alpha: the $1$-norm of the coefficients of the unitaries comprising the Hamiltonian $H$.
         precision: the precision $\epsilon$ to approximate $e^{it\cos\theta}$ to a polynomial.
     """
 
     walk_operator: QubitizationWalkOperator
     t: SymbolicFloat = field(kw_only=True)
-    alpha: SymbolicFloat = field(kw_only=True)
     precision: SymbolicFloat = field(kw_only=True)
 
+    def __attrs_post_init__(self):
+        if self.walk_operator.sum_of_lcu_coefficients is None:
+            raise ValueError(
+                f"Missing attribute `sum_of_ham_coeffs` for {self.walk_operator}, cannot implement Hamiltonian Simulation"
+            )
+
     def is_symbolic(self):
         return is_symbolic(self.t, self.alpha, self.precision)
 
+    @property
+    def alpha(self):
+        return self.walk_operator.sum_of_lcu_coefficients
+
     @cached_property
     def degree(self) -> SymbolicInt:
         r"""degree of the polynomial to approximate the function e^{it\cos(\theta)}"""
@@ -179,9 +187,7 @@ def _hubbard_time_evolution_by_gqsp() -> HamiltonianSimulationByGQSP:
     from qualtran.bloqs.hubbard_model import get_walk_operator_for_hubbard_model
 
     walk_op = get_walk_operator_for_hubbard_model(2, 2, 1, 1)
-    hubbard_time_evolution_by_gqsp = HamiltonianSimulationByGQSP(
-        walk_op, t=5, alpha=1, precision=1e-7
-    )
+    hubbard_time_evolution_by_gqsp = HamiltonianSimulationByGQSP(walk_op, t=5, precision=1e-7)
     return hubbard_time_evolution_by_gqsp
 
 
@@ -192,10 +198,9 @@ def _symbolic_hamsim_by_gqsp() -> HamiltonianSimulationByGQSP:
     from qualtran.bloqs.hubbard_model import get_walk_operator_for_hubbard_model
 
     walk_op = get_walk_operator_for_hubbard_model(2, 2, 1, 1)
-    t, alpha, inv_eps = sympy.symbols("t alpha N")
-    symbolic_hamsim_by_gqsp = HamiltonianSimulationByGQSP(
-        walk_op, t=t, alpha=alpha, precision=1 / inv_eps
-    )
+
+    t, inv_eps = sympy.symbols("t N")
+    symbolic_hamsim_by_gqsp = HamiltonianSimulationByGQSP(walk_op, t=t, precision=1 / inv_eps)
     return symbolic_hamsim_by_gqsp
 
 

qualtran/bloqs/hamiltonian_simulation/hamiltonian_simulation_by_gqsp_test.py

@@ -25,7 +25,7 @@
     check_polynomial_pair_on_random_points_on_unit_circle,
     verify_generalized_qsp,
 )
-from qualtran.bloqs.qubitization_walk_operator import _walk_op, QubitizationWalkOperator
+from qualtran.bloqs.qubitization_walk_operator import QubitizationWalkOperator
 
 from .hamiltonian_simulation_by_gqsp import (
     _hubbard_time_evolution_by_gqsp,
@@ -52,10 +52,11 @@ def test_generalized_qsp_with_exp_cos_approx_on_random_unitaries(
     bitsize: int, t: float, precision: float
 ):
     random_state = np.random.RandomState(42)
+    W, _ = random_qubitization_walk_operator(1, 1, random_state=random_state)
 
     for _ in range(5):
         U = MatrixGate.random(bitsize, random_state=random_state)
-        gqsp = HamiltonianSimulationByGQSP(_walk_op, t=t, alpha=1, precision=precision).gqsp
+        gqsp = HamiltonianSimulationByGQSP(W, t=t, precision=precision).gqsp
         P, Q = gqsp.P, gqsp.Q
 
         check_polynomial_pair_on_random_points_on_unit_circle(
@@ -65,16 +66,11 @@ def test_generalized_qsp_with_exp_cos_approx_on_random_unitaries(
 
 
 def verify_hamiltonian_simulation_by_gqsp(
-    W: QubitizationWalkOperator,
-    H: NDArray[np.complex_],
-    *,
-    t: float,
-    alpha: float,
-    precision: float,
+    W: QubitizationWalkOperator, H: NDArray[np.complex_], *, t: float, precision: float
 ):
     N = H.shape[0]
 
-    W_e_iHt = HamiltonianSimulationByGQSP(W, t=t, alpha=alpha, precision=precision)
+    W_e_iHt = HamiltonianSimulationByGQSP(W, t=t, precision=precision)
     result_unitary = cirq.unitary(W_e_iHt)
 
     expected_top_left = scipy.linalg.expm(-1j * H * t)
@@ -96,4 +92,4 @@ def test_hamiltonian_simulation_by_gqsp(
         W, H = random_qubitization_walk_operator(
             select_bitsize, target_bitsize, random_state=random_state
         )
-        verify_hamiltonian_simulation_by_gqsp(W, H.matrix(), t=t, alpha=1, precision=precision)
+        verify_hamiltonian_simulation_by_gqsp(W, H.matrix(), t=t, precision=precision)

qualtran/bloqs/hubbard_model.py

@@ -47,7 +47,7 @@
 See the documentation for `PrepareHubbard` and `SelectHubbard` for details.
 """
 from functools import cached_property
-from typing import Collection, Optional, Sequence, Tuple, Union
+from typing import Collection, Optional, Sequence, Tuple, TYPE_CHECKING, Union
 
 import attrs
 import cirq
@@ -67,6 +67,9 @@
     PrepareUniformSuperposition,
 )
 
+if TYPE_CHECKING:
+    from qualtran.resource_counting.symbolic_counting_utils import SymbolicFloat
+
 
 @attrs.frozen
 class SelectHubbard(SelectOracle):
@@ -309,6 +312,13 @@ def selection_registers(self) -> Tuple[Register, ...]:
     def junk_registers(self) -> Tuple[Register, ...]:
         return (Register('temp', QAny(2)),)
 
+    @cached_property
+    def l1_norm_of_coeffs(self) -> 'SymbolicFloat':
+        # https://arxiv.org/abs/1805.03662v2 equation 60
+        N = self.x_dim * self.y_dim * 2
+        qlambda = 2 * N * self.t + (N * self.mu) // 2
+        return qlambda
+
     @cached_property
     def signature(self) -> Signature:
         return Signature([*self.selection_registers, *self.junk_registers])
@@ -321,8 +331,7 @@ def decompose_from_registers(
         temp = quregs['temp']
 
         N = self.x_dim * self.y_dim * 2
-        qlambda = 2 * N * self.t + (N * self.mu) // 2
-        yield cirq.Ry(rads=2 * np.arccos(np.sqrt(self.t * N / qlambda))).on(*V)
+        yield cirq.Ry(rads=2 * np.arccos(np.sqrt(self.t * N / self.l1_norm_of_coeffs))).on(*V)
         yield cirq.Ry(rads=2 * np.arccos(np.sqrt(1 / 5))).on(*U).controlled_by(*V)
         yield PrepareUniformSuperposition(self.x_dim).on_registers(controls=[], target=p_x)
         yield PrepareUniformSuperposition(self.y_dim).on_registers(controls=[], target=p_y)
@@ -352,4 +361,5 @@ def get_walk_operator_for_hubbard_model(
 ) -> 'QubitizationWalkOperator':
     select = SelectHubbard(x_dim, y_dim)
     prepare = PrepareHubbard(x_dim, y_dim, t, mu)
+
     return QubitizationWalkOperator(select=select, prepare=prepare)

qualtran/bloqs/multiplexers/select_pauli_lcu_test.py

@@ -172,6 +172,7 @@ def test_select_application_to_eigenstates():
     # right now we only handle positive interaction term values
     ham = get_1d_ising_hamiltonian(target, 1, 1)
     dense_pauli_string_hamiltonian = [tt.dense(target) for tt in ham]
+    qubitization_lambda = sum(xx.coefficient.real for xx in dense_pauli_string_hamiltonian)
     # built select with unary iteration gate
     bloq = SelectPauliLCU(
         selection_bitsize=selection_bitsize,
@@ -187,11 +188,10 @@ def test_select_application_to_eigenstates():
     all_qubits = select_circuit.all_qubits()
 
     coeffs = get_1d_ising_lcu_coeffs(num_sites, 1, 1)
-    prep_circuit = _fake_prepare(np.sqrt(coeffs), selection)
+    prep_circuit = _fake_prepare(np.sqrt(coeffs / qubitization_lambda), selection)
     turn_on_control = cirq.Circuit(cirq.X.on(control))
 
     ising_eigs, ising_wfns = np.linalg.eigh(ham.matrix())
-    qubitization_lambda = sum(xx.coefficient.real for xx in dense_pauli_string_hamiltonian)
     for iw_idx, ie in enumerate(ising_eigs):
         eigenstate_prep = cirq.Circuit()
         eigenstate_prep.append(

qualtran/bloqs/qubitization_walk_operator.ipynb

@@ -38,14 +38,16 @@
     "## `QubitizationWalkOperator`\n",
     "Constructs a Szegedy Quantum Walk operator using LCU oracles SELECT and PREPARE.\n",
     "\n",
+    "For a Hamiltonian $H = \\sum_l w_l H_l$ (s.t. $w_l > 0$ and $H_l$ are unitaries),\n",
     "Constructs a Szegedy quantum walk operator $W = R_{L} . SELECT$, which is a product of\n",
     "two reflections $R_{L} = (2|L><L| - I)$ and $SELECT=\\sum_{l}|l><l|H_{l}$.\n",
     "\n",
     "The action of $W$ partitions the Hilbert space into a direct sum of two-dimensional irreducible\n",
     "vector spaces. For an arbitrary eigenstate $|k>$ of $H$ with eigenvalue $E_k$, $|\\ell>|k>$ and\n",
     "an orthogonal state $\\phi_{k}$ span the irreducible two-dimensional space that $|\\ell>|k>$ is\n",
     "in under the action of $W$. In this space, $W$ implements a Pauli-Y rotation by an angle of\n",
-    "$-2arccos(E_{k} / \\lambda)$ s.t. $W = e^{i arccos(E_k / \\lambda) Y}$.\n",
+    "$-2arccos(E_{k} / \\lambda)$ s.t. $W = e^{i arccos(E_k / \\lambda) Y}$,\n",
+    "where $\\lambda = \\sum_l w_l$.\n",
     "\n",
     "Thus, the walk operator $W$ encodes the spectrum of $H$ as a function of eigenphases of $W$\n",
     "s.t. $spectrum(H) = \\lambda cos(arg(spectrum(W)))$ where $arg(e^{i\\phi}) = \\phi$.\n",

qualtran/bloqs/qubitization_walk_operator.py

@@ -35,14 +35,16 @@
 class QubitizationWalkOperator(GateWithRegisters):
     r"""Constructs a Szegedy Quantum Walk operator using LCU oracles SELECT and PREPARE.
 
+    For a Hamiltonian $H = \sum_l w_l H_l$ (s.t. $w_l > 0$ and $H_l$ are unitaries),
     Constructs a Szegedy quantum walk operator $W = R_{L} . SELECT$, which is a product of
     two reflections $R_{L} = (2|L><L| - I)$ and $SELECT=\sum_{l}|l><l|H_{l}$.
 
     The action of $W$ partitions the Hilbert space into a direct sum of two-dimensional irreducible
     vector spaces. For an arbitrary eigenstate $|k>$ of $H$ with eigenvalue $E_k$, $|\ell>|k>$ and
     an orthogonal state $\phi_{k}$ span the irreducible two-dimensional space that $|\ell>|k>$ is
     in under the action of $W$. In this space, $W$ implements a Pauli-Y rotation by an angle of
-    $-2arccos(E_{k} / \lambda)$ s.t. $W = e^{i arccos(E_k / \lambda) Y}$.
+    $-2arccos(E_{k} / \lambda)$ s.t. $W = e^{i arccos(E_k / \lambda) Y}$,
+    where $\lambda = \sum_l w_l$.
 
     Thus, the walk operator $W$ encodes the spectrum of $H$ as a function of eigenphases of $W$
     s.t. $spectrum(H) = \lambda cos(arg(spectrum(W)))$ where $arg(e^{i\phi}) = \phi$.
@@ -92,6 +94,11 @@ def signature(self) -> Signature:
     def reflect(self) -> ReflectionUsingPrepare:
         return ReflectionUsingPrepare(self.prepare, control_val=self.control_val, global_phase=-1)
 
+    @cached_property
+    def sum_of_lcu_coefficients(self) -> Optional[float]:
+        r"""value of $\lambda$, i.e. sum of absolute values of coefficients $w_l$."""
+        return self.prepare.l1_norm_of_coeffs
+
     def decompose_from_registers(
         self,
         context: cirq.DecompositionContext,
@@ -132,17 +139,13 @@ def controlled(
         ):
             c_select = self.select.controlled(control_values=control_values)
             assert isinstance(c_select, SelectOracle)
-            return QubitizationWalkOperator(
-                c_select, self.prepare, control_val=control_values[0], power=self.power
-            )
+            return attrs.evolve(self, select=c_select, control_val=control_values[0])
         raise NotImplementedError(
             f'Cannot create a controlled version of {self} with control_values={control_values}.'
         )
 
     def with_power(self, new_power: int) -> 'QubitizationWalkOperator':
-        return QubitizationWalkOperator(
-            self.select, self.prepare, control_val=self.control_val, power=new_power
-        )
+        return attrs.evolve(self, power=new_power)
 
     def __pow__(self, power: int):
         return self.with_power(self.power * power)