Hamiltonian Simulation

Overview

The unitarylab_algorithms.hamiltonian_simulation package provides five methods for approximating the time-evolution operator $e^{-iHt}$ of a Hamiltonian $H$. Except for the Cartan decomposition, the other four methods all apply a unified input formatting check (_format_system) and accept the same core inputs (H, t, error), while differing in approximation strategy, circuit depth, and accuracy guarantees.

Naming Ambiguity Warning: The QSPHSAlgorithm in this package is located at unitarylab_algorithms.hamiltonian_simulation.qsp.algorithm and is used for approximating Hamiltonian time evolution; the unitarylab_algorithms.linear_algebra package also has a folder with the same name, qsp/, whose QSPAlgorithm is used for solving systems of linear equations (see Linear Algebra Algorithms). The two are completely different in functionality and have different import paths — please be careful to distinguish between them when using them.

Choosing a Method

Unified Note: The status and error Parameters

• For all 5 algorithms in this package, .run() unconditionally sets self.status to "success" upon successful completion of the computation, and the success argument passed to _build_return_dict() is always True — in other words, the top-level status in the returned dictionary is always 'ok', regardless of whether the approximation error actually reaches the expected level of error. If you need to assess accuracy, you should compare the Frobenius norm of error in the returned result (or Final total error for Cartan) against the error threshold you passed in yourself.
• error participates in precision or scale control in QSP-HS, Taylor, and Trotter; qDrift controls precision mainly through steps.

Suzuki-Trotter Product Formula Evolution Algorithm

Background

The Trotter–Suzuki product formula approximates $e^{-i(A+B)t}$ as $(e^{-iAt/r}{e}^{-iBt/r}{)}^r$ (first order or higher). The method is simple to implement and is the most widely used Hamiltonian simulation method in practice.

First-order error: $O(t^2\Vert [A,B]\Vert /r)$ Second order (Suzuki): $O(t^3/r^2)$

Import

from unitarylab_algorithms import TrotterAlgorithm

.run() Parameters

def run(self, H: np.ndarray, t: float, error: float, order: int = 1, steps: int = 1000,
        backend='torch', device='cpu', dtype=np.complex128)

Return Value

circuit_path is a list of two paths: [the full repeated-circuit path, the single time-slice circuit path], unlike the other 4 algorithms in this package (which use a single string).

{
    'status': 'ok',   # always 'ok'; does not reflect the actual approximation accuracy
    'circuit_path': ['/path/.../trotter_full.svg', '/path/.../trotter_slice.svg'],
    'plot': [{'format': 'txt', 'filename': 'trotter_hamiltonian_simulation_algorithm_result.txt'}],
    'circuit': <Circuit>,   # the full repeated circuit
    'Approximate evolution matrix': array([...]),
    'Exact evolution matrix': array([...]),
    'Frobenius norm of error': 1.2e-05,
}

Example

import numpy as np
from unitarylab_algorithms import TrotterAlgorithm
 
H = np.array([[2, 1], [1, 3]])
algo = TrotterAlgorithm()
result = algo.run(H=H, t=1.0, error=1e-8, order=1, steps=1000)
print(result['Frobenius norm of error'])

Quick Demo

from unitarylab_algorithms.hamiltonian_simulation.trotter.algorithm import test
test()

Notes

• steps is actually an upper bound, not a fixed step count — the real step count is automatically computed by a formula based on t, order, error, and the spectral norm of $H$, and the smaller of that value and the passed-in steps is used; therefore passing a very large steps does not necessarily yield that many steps.
• error participates in computing the upper bound on the step count.
• When order is a higher order (e.g., 2, 4), the Suzuki recursive formula (_recurse) is used; in theory only order 1 or even orders are supported.

qDrift Algorithm

Background

qDrift is a randomized product formula that samples Pauli terms according to the magnitude of their coefficients and applies them with appropriately rescaled angles. This produces an unbiased stochastic approximation that converges to the exact evolution as the number of samples $N$ increases.

Gate complexity: $O({\lambda }^2{t}^2/ϵ)$, where $\lambda ={\sum }_j|{\alpha }_j|$

Import

from unitarylab_algorithms import QDriftAlgorithm

.run() Parameters

def run(self, H: np.ndarray, t: float, error: float, steps: int = 5000,
        backend='torch', device='cpu', dtype=np.complex128)

Return Value

{
    'status': 'ok',
    'circuit_path': '/path/to/qdrift_algorithm_circuit.svg',
    'plot': [{'format': 'txt', 'filename': 'qdrift_algorithm_result.txt'}],
    'circuit': <Circuit>,
    'Approximate evolution matrix': array([...]),
    'Exact evolution matrix': array([...]),
    'Frobenius norm of error': 3.4e-03,
}

Example

import numpy as np
from unitarylab_algorithms import QDriftAlgorithm
 
H = np.array([[2, 1], [1, 3]])
algo = QDriftAlgorithm()
result = algo.run(H=H, t=1.0, error=1e-8, steps=5000)
print(result['Frobenius norm of error'])

Quick Demo

from unitarylab_algorithms.hamiltonian_simulation.qdrift.algorithm import test
test()

Notes

• error is used only for validity checking; steps is the key parameter that controls accuracy.
• Because this relies on random sampling (np.random.choice, with no fixed seed set), the specific circuit and the Frobenius norm of error will differ on every run and are not reproducible; to reproduce results, you need to set np.random.seed(...) yourself before calling it.

Taylor-Series Hamiltonian Simulation

Background

The Taylor method expands $e^{-iHt}$ as a truncated Taylor series up to degree $d$, and then implements each term using a linear combination of unitaries (LCU). It can achieve precision $ϵ$ with $d=O(\log (1/ϵ)/\log \log (1/ϵ))$ terms.

Import

from unitarylab_algorithms import TaylorAlgorithm

.run() Parameters

def run(self, H: np.ndarray, t: float, error: float, degree: int = 15,
        backend='torch', device='cpu', dtype=np.complex128)

The default value in the web front end’s parameters.json is d=10, which is inconsistent with .run()’s own default value of degree=15; if the Python API is called without explicitly passing this parameter, 15 is used.

Return Value

{
    'status': 'ok',
    'circuit_path': '/path/to/taylor_hamiltonian_simulation_algorithm_circuit.svg',
    'plot': [{'format': 'txt', 'filename': 'taylor_hamiltonian_simulation_algorithm_result.txt'}],
    'circuit': <Circuit>,   # LCU circuit object (not decomposed)
    'Approximate evolution matrix': array([...]),
    'Exact evolution matrix': array([...]),
    'Frobenius norm of error': 5.6e-06,
}

Example

import numpy as np
from unitarylab_algorithms import TaylorAlgorithm
 
H = np.array([[2, 1], [1, 3]])
algo = TaylorAlgorithm()
result = algo.run(H=H, t=1.0, error=1e-8, degree=15)
print(result['Frobenius norm of error'])

Quick Demo

from unitarylab_algorithms.hamiltonian_simulation.taylor.algorithm import test
test()

Notes

• The effective expansion degree is computed as degree = min(max(degree, ⌈1.5λ + 1.5·ln(1/error)⌉), 15) (where $\lambda =\Vert H{\Vert }_2\cdot t$), with a hard upper bound of 15 — passing degree=100 still results in a computation using no more than 15 degrees. This is an inherent limitation of the current implementation, not a configuration oversight.
• error participates in the effective-degree calculation above.
• The saved circuit diagram corresponds to circuit.decompose() (the decomposed gate sequence), whereas the returned circuit field is the undecomposed LCU circuit object — the two have different structures.

Quantum Signal Processing Hamiltonian Simulation (QSP-HS)

Background

QSP-based Hamiltonian simulation constructs a quantum circuit that approximates the time-evolution operator $e^{-iHt}$ by encoding the eigenvalues of $H$ as a signal and applying a sequence of controlled unitaries and single-qubit rotations. For sparse Hamiltonians, QSP can achieve near-optimal gate complexity. Internally, $H$ is first block-encoded (block_encode(H, method="nagy")), then the two components $\cos (tH)$ and $\sin (tH)$ are constructed via Chebyshev/Bessel coefficients, and finally combined using LCU.

Gate complexity: $O((\Vert H{\Vert }_{\max }t+\log (1/ϵ))\text{ poly}(\log n))$

Import

from unitarylab_algorithms import QSPHSAlgorithm

The class name for this module is QSPHSAlgorithm (distinct from the QSPAlgorithm in the identically named folder under linear_algebra.qsp — see the notice at the top of this page).

.run() Parameters

def run(self, H: np.ndarray, t: float, error: float, degree: int = 15, beta: float = 0.7,
        backend='torch', device='cpu', dtype=np.complex128) -> Dict[str, Any]

Return Value

{
    'status': 'ok',
    'circuit_path': '/path/to/qsp_hamiltonian_simulation_algorithm_circuit.svg',
    'plot': [{'format': 'txt', 'filename': 'qsp_hamiltonian_simulation_algorithm_result.txt'}],
    'circuit': <Circuit>,   # the QSP+LCU circuit for a single time slice
    'Approximate evolution matrix': array([...]),
    'Exact evolution matrix': array([...]),
    'Frobenius norm of error': 8.9e-07,
}

Example

import numpy as np
from unitarylab_algorithms import QSPHSAlgorithm
 
H = np.array([[2, 1], [1, 3]])
algo = QSPHSAlgorithm()
result = algo.run(H=H, t=1.0, error=1e-8, degree=15, beta=0.7)
print(result['Frobenius norm of error'])

Quick Demo

from unitarylab_algorithms.hamiltonian_simulation.qsp.algorithm import test
test()

Notes

• error participates in estimating both the number of time slices and the degree required per slice (_estimate_required_degree).
• In the web front end’s parameter panel, this algorithm’s degree parameter is named d, which is inconsistent with the degree parameter name of the Python .run() method — this is a naming inconsistency; when calling the Python API, use degree.

Cartan Decomposition Algorithm

Background

The Cartan–Lax flow algorithm decomposes the time-evolution operator by decomposing the Lie algebra $\mathfrak{su}(2^n)=\mathfrak{k}\oplus \mathfrak{m}$ into a symmetric subalgebra $\mathfrak{k}$ and an antisymmetric space $\mathfrak{m}$. The resulting circuit has the form $K\cdot e^{-i\eta }\cdot K^{†}$. Internally this delegates directly to unitarylab.library.hamiltonian.hamiltonian_simulation(H, evol_time, method='cartan-lax', ...).

Requirement: $H$ must be a real symmetric matrix (however, unlike the other 4 algorithms, this algorithm does not perform explicit Hermitian / square-matrix / power-of-2 dimension validation at the Python layer — that validation happens inside the underlying hamiltonian_simulation() to which it delegates).

Import

from unitarylab_algorithms import CartanDecompositionAlgorithm

.run() Parameters

def run(self, H: Union[np.ndarray, list], t: float, error: float,
        backend: str = "torch", device: str = "cpu", dtype=np.complex128, **kwargs: Any)

Return Value

The field names in the returned output are completely different from the other 4 algorithms (the other 4 uniformly use Approximate evolution matrix/Exact evolution matrix/Frobenius norm of error, whereas this algorithm uses):

{
    'status': 'ok',
    'circuit_path': '/path/to/cartan_algorithm_circuit.svg',
    'plot': [{'format': 'txt', 'filename': 'cartan_algorithm_result.txt'}],
    'circuit': <Circuit>,
    'Evolution result': <the evolution-result object returned by the underlying simulator>,
    'Final total error': 4.2e-08,
    'Computation time (s)': 0.0345,   # the only algorithm that includes computation time in its output
    'Exact evolution': array([...]),  # the exact evolution matrix exp(-iHt); the key name also differs from the other algorithms (Exact evolution matrix)
}

Example

import numpy as np
from unitarylab_algorithms import CartanDecompositionAlgorithm
 
H = np.array([[2, 1], [1, 2]])
algo = CartanDecompositionAlgorithm()
result = algo.run(H=H, t=1.0, error=1e-3)
print(result['Final total error'])

Quick Demo

from unitarylab_algorithms.hamiltonian_simulation.cartan.algorithm import test
test()

Notes

• The Cartan method currently supports only matrix-form Hamiltonians; Pauli-string input is not yet supported.
• For ill-conditioned Hamiltonians, reducing lr and increasing reps can improve convergence.
• The output field naming (Evolution result/Final total error/Exact evolution) is inconsistent with the other 4 algorithms (Approximate evolution matrix/Exact evolution matrix/Frobenius norm of error). When mixing multiple methods, be careful to distinguish between the key names — do not assume that all hamiltonian_simulation algorithms share a unified output field structure.