Hamiltonian Unitary Builder

The HamiltonianUnitaryBuilder algorithm in QDK/Chemistry constructs a unitary based on the Hamiltonian, such as time simulation unitary \(U(t) = e^{-iHt}\) or block-encoded unitary \(U = \frac{H}{\|H\|}\). Following QDK/Chemistry’s algorithm design principles, it takes a QubitOperator and produces a UnitaryRepresentation as output.

Overview

Building unitary from Hamiltonian — such as the Hamiltonian simulation unitary \(U(t) = e^{-iHt}\) or block encoding unitary \(U = \frac{H}{\|H\|}\) — is a central subroutine in many quantum algorithms. The HamiltonianUnitaryBuilder provides a unified interface for methods that construct this operator from a QubitOperator.

QDK/Chemistry currently provides two families of implementations for this task: Trotter-Suzuki product formulas and block encoding.

The resulting UnitaryRepresentation objects wrap either a PauliProductFormulaContainer (Trotter) or an LCUContainer (block encoding).

Using the HamiltonianUnitaryBuilder

Note

This algorithm is currently available only in the Python API.

This section demonstrates how to create, configure, and run a Hamiltonian unitary builder. The run method returns a UnitaryRepresentation object that can be used by any algorithm that requires a Hamiltonian simulation unitary (e.g., PhaseEstimation).

Input requirements

The HamiltonianUnitaryBuilder requires the following inputs:

QubitOperator

A QubitOperator containing the Pauli-string representation of the Hamiltonian. This can be obtained from the QubitMapper algorithm, constructed from a model Hamiltonian, or built directly.

Creating a builder

Python API

from qdk_chemistry.algorithms import create

# Create a time evolution builder
trotter = create("hamiltonian_unitary_builder", "trotter")

Configuring settings

Settings vary by implementation. See Available implementations below for implementation-specific options.

Python API

# Configure a second-order Trotter builder with automatic step count
trotter = create("hamiltonian_unitary_builder", "trotter")
trotter.settings().set("order", 2)
trotter.settings().set("target_accuracy", 1e-3)
trotter.settings().set("error_bound", "commutator")

Running the builder

Python API

import numpy as np
from qdk_chemistry.algorithms import create
from qdk_chemistry.data import Structure

# 1. Setup molecule
coords = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 1.4]])
symbols = ["H", "H"]
structure = Structure(coords, symbols=symbols)

# 2. SCF
scf_solver = create("scf_solver")
E_scf, wfn_scf = scf_solver.run(
    structure, charge=0, spin_multiplicity=1, basis_or_guess="sto-3g"
)

# 3. Hamiltonian construction
hamiltonian_constructor = create("hamiltonian_constructor")
hamiltonian = hamiltonian_constructor.run(wfn_scf.get_orbitals())

# 4. Qubit mapping
from qdk_chemistry.data import MajoranaMapping

n_spin_orbitals = 2 * hamiltonian.get_orbitals().get_num_molecular_orbitals()
qubit_mapper = create("qubit_mapper")
qubit_ham = qubit_mapper.run(
    hamiltonian, MajoranaMapping.jordan_wigner(n_spin_orbitals)
)

# 5. Build time evolution unitary
trotter = create("hamiltonian_unitary_builder", "trotter", order=2, time=0.1)
evolution = trotter.run(qubit_ham)

print(f"Container type: {evolution.get_container_type()}")
print(f"Number of qubits: {evolution.get_num_qubits()}")
print(evolution.get_summary())

Available implementations

QDK/Chemistry’s HamiltonianUnitaryBuilder provides a unified interface for Hamiltonian simulation methods. You can discover available implementations programmatically:

Python API

from qdk_chemistry.algorithms import registry

# List all registered time evolution builder implementations
implementations = registry.available("hamiltonian_unitary_builder")
print(implementations)  # e.g. ['trotter', 'qdrift', 'partially_randomized', 'lcu']

Trotter-Suzuki product formulas

Factory name: "trotter"

The Trotter-Suzuki decomposition approximates the time-evolution operator as a product of individual Pauli exponentials. For a Hamiltonian \(H = \sum_j \alpha_j P_j\):

First-order Trotter (\(p = 1\)):

\[e^{-iHt} \approx \left[\prod_j e^{-i\alpha_j P_j t / N}\right]^N\]

Second-order Trotter (Strang splitting, \(p = 2\)):

\[e^{-iHt} \approx \left[\prod_{j=1}^{L-1} e^{-i\alpha_j P_j t / 2N} \cdot e^{-i\alpha_L P_L t / N} \cdot \prod_{j=L-1}^{1} e^{-i\alpha_j P_j t / 2N}\right]^N\]

Higher even orders are constructed via the recursive Suzuki composition [Suz92]:

\[S_{2k}(t) = S_{2k-2}(u_k t)^2 \, S_{2k-2}\!\bigl((1-4u_k) t\bigr) \, S_{2k-2}(u_k t)^2\]

where \(u_k = 1/(4 - 4^{1/(2k-1)})\).

The number of Trotter steps \(N\) can be specified directly (num_divisions) or computed automatically from a target_accuracy using one of two error bounds:

Commutator bound (default)

A tighter bound [CST+21] based on nested commutators: \(N = \lceil \frac{t^{2}}{2\epsilon} \sum_{j<k}\lVert[\alpha_jP_j,\alpha_kP_k]\rVert \rceil\)

Naive bound

A looser triangle-inequality bound: \(N = \lceil (\sum_j|\alpha_j|)^{2}t^{2}/\epsilon \rceil\)

When both num_divisions and target_accuracy are specified, the builder uses whichever requires more Trotter steps.

Settings

Setting	Type	Description
order	int	Trotter-Suzuki order (1 for first-order, 2+ for higher even orders). Default is 1.
target_accuracy	float	Target approximation error \(\epsilon\). When set to 0.0 (default), automatic step-count estimation is disabled.
num_divisions	int	Explicit number of Trotter steps \(N\). When set to 0 (default), determined from target_accuracy.
error_bound	str	Error bound strategy: "commutator" (default, tighter) or "naive" (simpler).
weight_threshold	float	Coefficient threshold below which Pauli terms are discarded. Default is 1e-12.

Zassenhaus product formulas

Factory name: "zassenhaus"

The Zassenhaus decomposition approximates the time-evolution operator by recursively generating commutator correction terms. For a Hamiltonian \(H = \sum_j \alpha_j P_j\), the first-order approximation is the Trotter product, and higher-order approximations introduce explicit product-formula factors for lower-order commutator corrections:

First-order Zassenhaus (\(p = 1\)): Matches the first-order Trotter product formula.

Second-order Zassenhaus (\(p = 2\)):

\[e^{-iHt} \approx \left[ \prod_{j > k} e^{-\frac{1}{2} [H_j, H_k] t^2} \cdot \prod_j e^{-i H_j t} \right]^N\]

Higher orders (up to \(p = 4\)) recursively construct explicit nested-commutator terms. Unlike Trotter-Suzuki formulas where step count is used to reduce commutator errors, the Zassenhaus builder computes low-order corrections explicitly.

The number of divisions can be specified directly (num_divisions) or computed automatically from a target_accuracy using one of two bounds:

Commutator bound (default)

A tighter bound based on nested commutators of the first omitted exponent: \(N = \lceil \frac{\|C_{p+1}(-iH_1, \dots, -iH_m)\|^{1/p} t^{1+1/p}}{\epsilon^{1/p}} \rceil\)

Naive bound

A looser bound using the absolute coefficient sum of the first omitted exponent and the 1-norm of the Hamiltonian: \(N = \lceil \frac{(\kappa_{p+1} 2^p \|H\|_1^{p+1})^{1/p} t^{1+1/p}}{\epsilon^{1/p}} \rceil\)

When order is set to 0 (auto), the builder dynamically sweeps orders 2, 3, and 4 to select the order that minimizes the total gate count (product of step reps and step terms) for the target accuracy.

Settings

Setting	Type	Description
order	int	Zassenhaus order (0 for auto, or 1 to 4). Default is 0.
target_accuracy	float	Target approximation error \(\epsilon\). When set to 0.0 (default), automatic step-count estimation is disabled.
num_divisions	int	Explicit number of Zassenhaus steps \(N\). When set to 0 (default), determined from target_accuracy.
error_bound	str	Error bound strategy: "commutator" (default, tighter) or "naive" (simpler).
weight_threshold	float	Coefficient threshold below which Pauli terms are discarded. Default is 1e-12.

Consuming term partitions

When the input QubitOperator carries a populated term_partition, the Trotter builder consumes it directly:

• LayeredPartition (group → layer → index) is used as-is — the outer level controls the Strang/Suzuki splitting and each inner layer becomes one parallelisable sub-step.

• FlatPartition (group → index) is interpreted as a layered partition with one layer per group.

In both cases groups are sorted by ascending layer count so that the smallest groups sit on the outside of the Strang/Suzuki splitting, which maximises merging at recursion boundaries. This typically reduces the number of distinct exponentials per Trotter step and the saving compounds through the recursion at higher orders.

When term_partition is None each Pauli term is exponentiated as its own group. Pre-populate the partition using the term_grouper algorithm or one of the spin model Hamiltonian builders to enable group-aware scheduling.

For lattice models, the edge coloring stored on LatticeGraph feeds directly into the partition: edges of the same color have disjoint qubit supports and form a single parallelisable layer. The model Hamiltonian builders consume this coloring automatically, so no manual geometry handling is required.

Example:

from qdk_chemistry.data import LatticeGraph
from qdk_chemistry.utils.model_hamiltonians import create_ising_hamiltonian
from qdk_chemistry.algorithms import registry

# 4-site open Ising chain: H = Σ Z_i Z_{i+1} + 0.5 Σ X_i
graph = LatticeGraph.chain(4, periodic=False)
hamiltonian = create_ising_hamiltonian(graph, j=1.0, h=0.5)

# The edge coloring partitions ZZ terms into two layers by color:
#   group 0 (field):  1 layer  → [X₀, X₁, X₂, X₃]
#   group 1 (ZZ):     2 layers → [{Z₀Z₁, Z₂Z₃}, {Z₁Z₂}]
print(hamiltonian.term_partition)
# LayeredPartition(strategy='geometry_coloring', num_groups=2)

# Second-order Trotter with 1 division produces the Strang splitting:
#   S₂(t) = e^{fields·t/2}  e^{ZZ_layer0·t}  e^{ZZ_layer1·t}
#           e^{fields·t/2}
# where same-layer ZZ terms (e.g. Z₀Z₁ and Z₂Z₃) have disjoint
# qubit support and are exponentiated independently within one step.
trotter = registry.create("time_evolution_builder", "trotter")
trotter.settings().update({"order": 2, "num_divisions": 1})
evolution = trotter.run(hamiltonian, time=1.0)
container = evolution.get_container()

# The grouped schedule uses 11 exponentiated terms per step,
# vs. 13 for the ungrouped fallback — a 15% reduction that
# compounds at higher Suzuki orders.
print(f"{len(container.step_terms)} terms per Trotter step")
for term in container.step_terms:
    label = ['I'] * 4
    for q, p in term.pauli_term.items():
        label[q] = p
    print(f"  exp(-i * {term.angle:+.4f} * {''.join(reversed(label))})")
# Output:
#   exp(-i * +0.2500 * IIIX)
#   exp(-i * +0.2500 * IIXI)
#   exp(-i * +0.2500 * IXII)
#   exp(-i * +0.2500 * XIII)
#   exp(-i * +1.0000 * IIZZ)
#   exp(-i * +1.0000 * ZZII)
#   exp(-i * +1.0000 * IZZI)
#   exp(-i * +0.2500 * IIIX)
#   exp(-i * +0.2500 * IIXI)
#   exp(-i * +0.2500 * IXII)
#   exp(-i * +0.2500 * XIII)

Block encoding

Block encoding is a technique for embedding a non-unitary operator (such as \(H / \lambda\)) into a larger unitary circuit. Given an \(n\)-qubit Hamiltonian \(H\), a block encoding uses \(a\) ancilla qubits to construct a unitary \(U\) such that:

\[(\langle 0|^{\otimes a} \otimes I)\, U\, (|0\rangle^{\otimes a} \otimes I) = \frac{H}{\lambda}\]

where \(\lambda \geq \|H\|\) is a normalization factor (typically the L1 norm of the coefficients). The Hamiltonian is encoded in the top-left block of the larger unitary — hence the name “block encoding.”

Block encodings are the foundation for quantum algorithms such as qubitization, QSVT, and quantum linear system solvers. QDK/Chemistry currently provides the LCU (Linear Combination of Unitaries) method for constructing block encodings.

LCU (Linear Combination of Unitaries)

Factory name: "lcu"

The LCU builder constructs a block encoding of \(H / \lambda\) using the PREPARE–SELECT pattern [CW12]. For a Hamiltonian \(H = \sum_{j=1}^{L} \alpha_j P_j\) with \(\lambda = \sum_j |\alpha_j|\):

PREPARE oracle

Prepares an ancilla register in the state \(|\psi\rangle = \sum_j \sqrt{|\alpha_j| / \lambda}\, |j\rangle\) using \(\lceil \log_2 L \rceil\) qubits.

SELECT oracle

Applies the corresponding Pauli string controlled on the ancilla index:

\[\text{SELECT} = \sum_{j=0}^{L-1} |j\rangle\langle j| \otimes P_j\]

with sign corrections absorbed into a phase vector.

Block encoding circuit

The composite circuit \(\text{PREP}^\dagger \cdot \text{SEL} \cdot \text{PREP}\) realizes the block encoding:

\[(\langle 0| \otimes I)\, \text{PREP}^\dagger \cdot \text{SEL} \cdot \text{PREP}\, (|0\rangle \otimes I) = \frac{H}{\lambda}\]

Quantum walk operator

When quantum_walk=True, the builder wraps the block encoding with a reflection to form the qubitization walk operator:

\[W = (2|0\rangle\langle 0| - I) \cdot \text{PREP}^\dagger \cdot \text{SEL} \cdot \text{PREP}\]

The walk operator has eigenvalues \(e^{\pm i \arccos(E_k/\lambda)}\) where \(E_k\) are the eigenvalues of \(H\).

Settings

Setting	Type	Description
power	int	The power to which the walk operator is raised. Default is 1.
quantum_walk	bool	If True, wrap the block encoding with a reflection to form the qubitization walk operator. If False, produce the plain block encoding. Default is False.
tolerance	float	Minimum L1 norm below which the decomposition is considered ill-defined. Default is 1e-12.

Example

from qdk_chemistry.algorithms import registry

lcu = registry.create("hamiltonian_unitary_builder", "lcu")
lcu.settings().update({"quantum_walk": True})
unitary = lcu.run(qubit_hamiltonian)

The resulting UnitaryRepresentation wraps an LCUContainer containing the Prepare and Select oracles.

Related classes

• UnitaryRepresentation: Output data class wrapping the exponentiated Pauli terms or LCU container

• QubitOperator: Input qubit Hamiltonian

• PhaseEstimation: Consumer of the hamiltonian unitary

Further reading

• The above examples can be downloaded as a complete Python script.

• PhaseEstimation: Quantum phase estimation algorithms

• QubitMapper: Map fermionic Hamiltonians to qubit operators

• Settings: Configuration settings for algorithms

• Factory Pattern: Understanding algorithm creation