Group operator terms: use the operator structure

A BosonOperator can tag each of its terms with a group index. The tags carry no meaning of their own: they only say which terms belong together. What the tags mean is decided by whatever consumes them. In this package the consumers are the time-evolution gate, which is split into one factor per group, and the randomized qDRIFT pass, which samples whole groups.

This guide shows how to assign groups, how to check the properties a consumer relies on, and how to lay the operator out so that picking a few groups out of a large one stays cheap.

Note

No automatic grouping strategy ships with this package. You assign the indices yourself, from what you know about your model.

Assign groups

The running example is a three-site Bose-Hubbard chain. Each bond contributes a hopping pair, and the on-site interactions form a diagonal piece. The simplest way to build a grouped operator is from_terms_with_groups(), which takes each term together with its group index as an (actions, coefficient, group) triple. You decide the group at the moment you write the term down, when you still know why it belongs there:

>>> from qiskit_bosons.operators import BosonOperator, ann, cre
>>>
>>> bonds = [(0, 1), (1, 2)]
>>> terms = []
>>> for group, (left, right) in enumerate(bonds):  # one group per bond
...     terms.append(([cre(left), ann(right)], -1.0, group))
...     terms.append(([cre(right), ann(left)], -1.0, group))
>>> for site in range(3):  # the on-site interactions share one group
...     terms.append(([cre(site), cre(site), ann(site), ann(site)], 1.0, len(bonds)))
>>> hamiltonian = BosonOperator.from_terms_with_groups(terms)
>>> len(hamiltonian), hamiltonian.num_groups()
(7, 3)

The terms are stored in the order you give them, so the group indices cannot end up attached to the wrong terms.

An operator that already exists has to be grouped after the fact, by assigning the groups attribute. That is where care is needed, because from_dict() does not specify the order in which it stores its terms. Writing the list out by position can then attach an index to the wrong term. Derive each index from the term itself instead:

>>> existing = BosonOperator.from_dict({tuple(actions): coeff for actions, coeff, _ in terms})
>>> existing.groups = [
...     len(bonds)
...     if len({mode for _, mode in term}) == 1
...     else bonds.index(tuple(sorted({mode for _, mode in term})))
...     for term, _ in existing.iter_terms()
... ]
>>> existing.num_groups(), existing.equiv(hamiltonian)
(3, True)

The number of indices has to match the number of terms. Anything else is rejected, because a too-short list would silently drop trailing terms wherever terms are iterated together with their groups:

>>> try:
...     hamiltonian.groups = [0, 1]
... except ValueError as error:
...     print(error)
expected one group index per term, but got 2 for 7 terms

Check what the consumer relies on

Three functions in qiskit_bosons.operators.terms answer one narrow question each about an existing grouping. None of them says a grouping is “correct”, because correctness is relative to the consumer.

Is each group Hermitian? This is the property that matters most. The time evolution of a group is only unitary if the group is Hermitian, and a lone hopping term \(b^\dagger_i b_j\) is not. Only the pair with its adjoint is Hermitian:

>>> from qiskit_bosons.operators.terms import groups_are_hermitian
>>> groups_are_hermitian(hamiltonian)
[True, True, True]

A grouping that splits a hopping pair is caught immediately:

>>> split = BosonOperator.from_terms_with_groups(
...     [(actions, coeff, index) for index, (actions, coeff, _) in enumerate(terms)]
... )  # one term per group
>>> [bool(flag) for flag in groups_are_hermitian(split)]
[False, False, False, False, True, True, True]

The four hopping terms are caught. The three True entries are the on-site terms \(b^\dagger_i b^\dagger_i b_i b_i\), which are Hermitian by themselves. A mode can legitimately repeat within a term, and nothing is dropped for doing so. The flags follow the order of terms, which is why each one can be read against the term it describes.

Are the coefficients uniform? A randomized product formula weighs a group by the mean magnitude of its coefficients, which is a faithful description of the group only if the magnitudes agree:

>>> from qiskit_bosons.operators.terms import group_coeff_means, groups_have_uniform_coeffs
>>> groups_have_uniform_coeffs(hamiltonian)
[True, True, True]
>>> group_coeff_means(hamiltonian)
[1.0, 1.0, 1.0]

Neither check implies the other: uniform coefficients do not make a group Hermitian, and a Hermitian group can mix magnitudes. By default the magnitudes are compared, because that is what group_coeff_means() averages. A Hermitian pair with complex coefficients has equal magnitudes but unequal coefficients, so it satisfies only the default form:

>>> conjugate = BosonOperator.from_terms_with_groups(
...     [([cre(0), ann(1)], 1.0j, 0), ([cre(1), ann(0)], -1.0j, 0)]
... )
>>> groups_have_uniform_coeffs(conjugate)
[True]
>>> groups_have_uniform_coeffs(conjugate, abs=False)
[False]

Note

groups_are_hermitian() is built on is_hermitian(), which compares an operator with its adjoint after normal ordering. That is an identity of the untruncated algebra, so a group that passes is Hermitian as a bosonic operator. See Design principles of operator representations for how products of ladder operators are ordered and compared.

Order the terms by group

split_out_groups() returns one operator per requested group. When the terms of a group are scattered through the operator, it has to scan every term to find them. When each group forms one contiguous run, it binary-searches the group boundaries instead, so the cost of a lookup follows the number of groups requested rather than the number of terms held.

group_order() produces that layout. It is a stable sort by group index: the value of the operator does not change, and terms within a group keep their relative order.

>>> from qiskit_bosons.operators.terms import group_order
>>> scattered = BosonOperator.from_terms_with_groups(
...     [(actions, coeff, group) for (actions, coeff, _), group in zip(terms, [1, 0, 1, 0, 2, 2, 2])]
... )
>>> ordered = group_order(scattered)
>>> ordered.groups
[0, 0, 1, 1, 2, 2, 2]
>>> ordered.equiv(scattered)
True

Both layouts return the same operators; only the cost differs:

>>> request = [2, 0]
>>> all(
...     a.equiv(b)
...     for a, b in zip(
...         scattered.split_out_groups(request), ordered.split_out_groups(request), strict=True
...     )
... )
True

Next steps

The grouping is chosen here, at the bosonic level, while the structure of the model is still visible.