Algorithms and Utility Library

Overview

This chapter introduces the high-level algorithm modules in unitarylab.library, grouped by functionality with descriptions of each module’s purpose, applicable scenarios, and basic usage. After reading this chapter, you will be able to:

• Understand the functional categories of the modules in the algorithm library
• Know which module to use for common quantum computing tasks
• Build quantum circuits using high-level interfaces such as QFT, QPE, LCU, quantum signal processing/QSVT, block encoding, Hamiltonian simulation, and linear system solving

These modules all provide algorithm-level abstractions on top of Circuit, so you generally do not need to write gate sequences by hand.

Top-Level Exports and Import Conventions

The unitarylab.library package’s top level directly exports only the following 10 names (the __all__ in library/__init__.py):

from unitarylab.library import (
    QFT, IQFT,                      # Quantum Fourier Transform / inverse transform
    QPE,                             # Quantum Phase Estimation
    LCU,                             # Linear Combination of Unitaries
    block_encode,                    # Unified entry point for block encoding
    hamiltonian_simulation,          # Unified entry point for Hamiltonian simulation
    solve,                           # Unified entry point for solving linear systems
    QSVT,                            # Quantum Singular Value Transformation (scalar function transform)
    QSP, QSP_hamiltonian_simulation, # Quantum Signal Processing
)

The remaining, lower-level method classes (such as Trotter, HHLSolver, QSVTHermitian, CartanOptimization) are not exported at the top level and must be imported explicitly from their respective submodules; the correct import path is given in each subsequent section of this chapter.

Module Functional Groups

Basic Algorithm Primitives

QFT / IQFT: Quantum Fourier Transform

Purpose: Build a Quantum Fourier Transform or its inverse circuit; a foundational module for QPE, Shor’s algorithm, and others.

from unitarylab import Circuit
from unitarylab.library import QFT, IQFT
 
# Build a 4-qubit QFT circuit
qft_circuit = QFT(n=4)
 
# Build the inverse QFT circuit
iqft_circuit = IQFT(n=4)
 
# Can be embedded into a larger circuit
qc = Circuit(4)
qc.append(qft_circuit, target=[0, 1, 2, 3])
qc.draw()

QPE: Quantum Phase Estimation

Purpose: Estimate the eigenvalue phase $\theta$ of a unitary operator $U$, satisfying $U|\psi ⟩=e^{2\pi i\theta }|\psi ⟩$. Commonly used in quantum chemistry and quantum algorithms to estimate Hamiltonian eigen-energies.

from unitarylab import Circuit
from unitarylab.library import QPE
 
# The unitary operator to be estimated (T gate, phase = 1/8)
U = Circuit(1)
U.t(0)
 
# T|1⟩ = exp(2πi/8)|1⟩, so first prepare the eigenstate |1⟩
prepare_target = Circuit(1)
prepare_target.x(0)
 
# Run QPE, returns (circuit, phase estimate, probability of this result)
qpe_circuit, phase, probability = QPE(
    U=U,
    d=4,                  # number of phase register qubits, precision 1/2^d
    prepare_target=prepare_target,
    device='cpu',
    return_circuit=False
)
print(phase, probability)  # 0.125 1.0 (allowing for numerical rounding error)

LCU: Linear Combination of Unitaries

Purpose: Implement the operator $A={\sum }_j{\alpha }_j{U}_j$, commonly used as a foundational building block for Hamiltonian simulation and block encoding.

from unitarylab import Circuit
from unitarylab.library import LCU
 
# Build two unitary sub-circuits
U1 = Circuit(1)
U1.x(0)
U2 = Circuit(1)
U2.z(0)
 
# LCU: implement 0.6 * U1 + 0.8 * U2
lcu_circuit = LCU([(U1, 0.6), (U2, 0.8)])

Note: All input unitary circuits must have the same number of qubits, otherwise a ValueError is raised. Coefficients must be finite, non-negative floating-point numbers; terms with a coefficient of 0 are automatically skipped.

Hamiltonian Simulation

hamiltonian_simulation() is the unified entry point. Internally it dispatches to different algorithm implementations based on the method parameter, and all return values are instances of HamiltonianSimulationResult (or its subclasses), with a consistent interface.

Unified Entry Point

import numpy as np
from unitarylab.library import hamiltonian_simulation
 
H = np.array([[0, 1], [1, 0]], dtype=float)  # Real symmetric Pauli-X Hamiltonian
 
sim = hamiltonian_simulation(
    H, t=1.0,
    method='trotter',       # 'trotter' (default) / 'qdrift' / 'taylor' / 'qsp' (alias 'qsvt') /
                             # 'cartan-lax' / 'cartan-optimization'
    target_error=1e-3,
    backend='torch', device='cpu', dtype=np.complex128,
    order=1, steps=100,      # method-specific parameters, passed via **kwargs (here corresponding to 'trotter')
)
 
circuit = sim.circuit               # evolution circuit
evolution = sim.evolution_result    # approximate evolution matrix (lazily computed, built on first access)
print(sim.total_error)              # actual approximation error
circuit.draw(title="Hamiltonian Simulation")

Methods and their corresponding dedicated **kwargs:

Passing keyword arguments that are not on the corresponding method’s allowed list raises a ValueError (rather than being silently ignored).

Using the Underlying Method Classes Directly

If you need to bypass the unified entry point and instantiate a method directly (for example, to access additional method-specific attributes), you can import from unitarylab.library.hamiltonian.method:

import numpy as np
from unitarylab.library.hamiltonian.method import Trotter, CartanOptimization
 
H = np.array([[0, 1], [1, 0]], dtype=float)
 
# Equivalent to hamiltonian_simulation(H, t=1.0, method='trotter', order=1, steps=100)
sim = Trotter(H=H, t=1.0, target_error=1e-3, order=1, steps=100, device="cpu")
circuit = sim.circuit
evolution = sim.evolution_result   # Note: the public attribute name is evolution_result, without a leading underscore
print(sim.total_error)

# Cartan decomposition: suitable for exact Hamiltonian evolution on 1-2 qubits; error decreases as optimization converges
H = np.array([[1.0, 0.5], [0.5, -1.0]], dtype=float)  # must be a real symmetric matrix
 
sim = CartanOptimization(
    H=H, t=1.0, target_error=1e-3,
    lr=1e-3, max_steps=10000, optimizer="SD",  # options: 'SD', 'BB', 'BFGS'
    device="cpu",
)
circuit = sim.circuit
evolution = sim.evolution_result  # public attribute, not sim._evolution_result
circuit.draw(title="Cartan Hamiltonian Simulation")

Note: Trotter, QDrift, Taylor, QSP, CartanLax, and CartanOptimization all inherit from HamiltonianSimulationResult. Results should always be accessed through the public attributes .circuit, .evolution_result, and .total_error (internally cached as _circuit/_evolution_result/_total_error, but users should not access the underscore-prefixed internal attributes directly).

Quantum Signal Processing (QSP) and Quantum Singular Value Transformation (QSVT)

These two modules provide quantum algorithm frameworks based on polynomial transformations. They achieve higher precision but produce deeper circuits, making them suitable for scenarios with demanding precision requirements. Both can be imported directly from the unitarylab.library top level.

QSP: Block-Encoding-Based Hamiltonian Simulation

Purpose: Apply a polynomial transformation to an already block-encoded unitary operator. QSP_hamiltonian_simulation is a dedicated wrapper for approximating $e^{-iHt}$.

from unitarylab import Circuit
from unitarylab.library import QSP_hamiltonian_simulation
 
# 1. Construct a minimal block-encoding circuit U_H
# Here we use 1 system qubit + 1 ancilla qubit
# Applying a Z gate on system qubit q0 serves as a simple example of block-encoding a Pauli-Z Hamiltonian
n = 1      # number of system register qubits
m = 1      # number of ancilla qubits for block encoding
 
U_H = Circuit(n + m, name="Block Encoding of Z")
U_H.z(0)
 
# 2. Construct a QSP-based Hamiltonian simulation circuit
circuit, factor, n_ancilla, n_qubits, degree = QSP_hamiltonian_simulation(
    U_H=U_H,
    n=n,
    alpha=1.0,       # block-encoding normalization coefficient
    m=m,
    t=1.0,           # evolution time
    epsilon=1e-3,    # target approximation error
    beta=0.5,        # must satisfy 0 < beta < 1
    flag=True,       # True means approximating exp(-iHt)
)
 
print("factor:", factor)
print("n_ancilla:", n_ancilla)
print("n_qubits:", n_qubits)
print("degree:", degree)
print("circuit qubits:", circuit.get_num_qubits())
 
circuit.draw(title="QSP Hamiltonian Simulation")

QSVT (Quantum Singular Value Transformation): Arbitrary Scalar Function Transformation

Purpose: Apply an arbitrary scalar polynomial function transformation to a Hermitian matrix, applicable to Hamiltonian simulation, matrix inversion, and more. The top-level function QSVT() is the recommended entry point; internally it constructs a QSVTHermitian and immediately completes the fit, so the return value is already a fitted result object from which you can read all result attributes directly.

import numpy as np
from unitarylab.library import QSVT
 
# 1. Define the target Hermitian matrix
hamiltonian_matrix = np.array([
    [0.8, 0.0],
    [0.0, 0.4],
], dtype=complex)
 
# 2. Define the target scalar function, e.g. f(x) = exp(-i x)
target_function = lambda x: np.exp(-1j * x)
 
# 3. Call the unified entry point (the return value has already been fitted and can be used directly)
qsvt = QSVT(
    hamiltonian_matrix,
    function=target_function,
    target_error=1e-6,
    block_encoding_method="nagy",  # options: "nagy" or "fable"
    device="cpu",
)
 
circuit = qsvt.circuit
error = qsvt.total_error
evolution = qsvt.evolution_result
 
print("degree:", qsvt.degree)
print("alpha:", qsvt.alpha)
print("m:", qsvt.m)
print("total_error:", error)
 
circuit.draw(title="QSVT Hermitian Approximation")

If you need finer-grained control (for example, manually fitting in stages), you can use the underlying class directly via from unitarylab.library._qsvt.algorithm import QSVTHermitian. Its constructor parameters are identical to those of QSVT(), but you must additionally call ._run() manually to trigger the fit.

Block Encoding

Purpose: Embed a non-unitary matrix $A$ into the top-left subblock of a unitary matrix ($A/\alpha$). This is a foundational building block for algorithms such as QSVT and LCU.

from unitarylab.library import block_encode
import numpy as np
 
# Block-encode an arbitrary matrix
A = np.array([[0.5, 0.1], [0.2, 0.3]], dtype=complex)
result = block_encode(A, method='nagy', eps=1e-3, verbose=False)
 
circuit = result.circuit          # block-encoding circuit
alpha = result.alpha              # normalization coefficient
total_qubits = result.total_qubits
target_qubits = result.target_qubits

The underlying classes FABLE and Nagy can be imported separately via from unitarylab.library.block_encoding import FABLE, Nagy, for scenarios that require finer-grained control.

Linear System Solving

solve() is the unified entry point for solving linear systems of the form $Ax=b$. Internally it dispatches to different solvers based on the method parameter, and always returns a LinearSolverResult.

Unified Entry Point

from unitarylab.library import solve
import numpy as np
 
A = np.array([[2, 1], [1, 3]], dtype=float)
b = np.array([1, 0], dtype=float)
 
hhl_result = solve(A, b)                       # default method='hhl'
hhl_precise = solve(A, b, method='hhl', d=10) # HHL, 10 phase register qubits
qsvt_result = solve(A, b, method='qsvt')       # QSVT-based solver
vqls_result = solve(A, b, method='vqls')       # variational quantum linear solver
 
print(hhl_result.solution)   # approximate solution vector x
print(hhl_result.circuit)    # the constructed quantum circuit (may be None)

solve() also supports an optional precondition parameter (None / 'diagonal' / 'symmetric' / 'ilu', etc.), which preprocesses $(A,b)$ before calling the quantum solver, and automatically restores the solution to the original problem before returning:

result = solve(A, b, method='hhl', precondition='diagonal')
print(result.precondition_mode)   # preprocessing metadata

Using the Underlying Solvers Directly

from unitarylab.library.linear_solver import HHLSolver
import numpy as np
 
A = np.array([[2, 1], [1, 3]], dtype=float)
b = np.array([1, 0], dtype=float)
 
# Note: HHLSolver returns a tuple (circuit, solution, scaling_factor), not an object
circuit, solution, scale = HHLSolver(A, b, d=6)
print(solution)

QSVTSolver, SCHROSolver, AQCSolver, VQLSSolver, and CKSSolver can likewise be imported from unitarylab.library.linear_solver; their return value format is consistent with HHLSolver (all are (circuit, solution, scaling_factor) triples). Internally, the solve() unified entry point calls these same functions and wraps the result into a LinearSolverResult.

Pauli Utilities

Purpose: Decompose an arbitrary Hermitian matrix into a linear combination of Pauli tensor products. This is a preprocessing step for scenarios such as Hamiltonian simulation and observable estimation.

from unitarylab.library.pauli_operator.pauli_string_decomposition import pauli_string_decomposition
import numpy as np
 
# Decompose a 2x2 Hermitian matrix into Pauli terms
H = np.array([[1, 0.5], [0.5, -1]], dtype=complex)
terms = pauli_string_decomposition(H)
print(terms)
# Should include Z: 1 and X: 0.5; the return order is not guaranteed, and terms with magnitude below 1e-10 are automatically dropped

Return value note: pauli_string_decomposition() returns a list[tuple[str, complex]] (a list of tuples of Pauli strings and coefficients), not a dictionary; the leftmost character of a Pauli string corresponds to qubit 0 (the lowest-order qubit). Optional parameters: sets (a pre-specified set of Pauli strings), partition_commuting (whether to reorder into commuting groups, default True), real_symmetric_hint (if the input is known to be a real symmetric matrix, set to True to skip terms with an odd number of Y’s for speedup).

PDE Solving Framework

unitarylab.library.equation provides a complete quantum partial differential equation (PDE) solving framework based on the Schrödingerization method, including equation parsing, differential operator construction, and built-in equation examples. This submodule is designed for configuration-driven invocation (JSON parameters + an algorithm registry) rather than simple constructor calls, making it suitable for integration with upper-layer services. If you just want to call a specific algorithm from a script, we recommend reading the algorithm.py source in the corresponding example directory directly.

Sub-module Structure

equation/
├── differential_operator/    # Differential operator construction (finite differences, block encoding)
├── equation_parser/          # Equation parser (boundary conditions, discretization, solvers)
├── examples/                 # Built-in equation examples (each subdirectory is an independent algorithm module)
└── schrodingerization/       # Schrödingerization transform (time evolution)

Built-in Equation Examples

There are 20 built-in example directories under equation/examples/. Each directory is an independent algorithm module containing algorithm.py (a BaseAlgorithm subclass), setup.json (default parameter configuration), and __init__.py (declaring ALGORITHM_NAME/ALGORITHM_CLASS):

equation_heat, equation_heat2d, equation_heatVariableCoefficient, equation_backHeat, equation_backHeat2d, equation_advection, equation_burgers, equation_burgers2d, equation_maxwell, equation_Helmholtz, equation_blackScholes, equation_elasticWave, equation_elasticWave2d, equation_OUprocess, equation_SchrABC, equation_general, equation_hamiltonjacobi, equation_multiElliptic, equation_multiTransport, equation_traffic.

The run() method of every algorithm class accepts the common parameters params=None, algo_dir=None: when params is None, it automatically reads the setup.json in that example’s directory (a structured JSON describing boundary conditions, discretization scheme, and equation parameters); you can also pass an equivalent dict or JSON string to override the default configuration. algo_dir is used to store run logs and result files. Note: the signatures are not fully uniform — 11 examples (equation_heat, equation_heat2d, equation_backHeat, equation_backHeat2d, equation_advection, equation_burgers, equation_burgers2d, equation_elasticWave, equation_general, equation_traffic, equation_blackScholes) additionally support the three simulation parameters backend='torch', device='cpu', dtype=np.complex128, while the other 9 (equation_Helmholtz, equation_OUprocess, equation_SchrABC, equation_elasticWave2d, equation_hamiltonjacobi, equation_heatVariableCoefficient, equation_maxwell, equation_multiElliptic, equation_multiTransport) currently accept only the params/algo_dir parameters. Before calling, we recommend using inspect.signature() or checking the corresponding algorithm.py directly to confirm the exact signature for that example.

from unitarylab.library.equation.examples.equation_heat import HeatEquationAlgorithm
 
solver = HeatEquationAlgorithm()
result = solver.run()   # when params=None, automatically reads the default configuration in equation_heat/setup.json

Note: The parameter structure of this framework (fields such as boundary_condition/discrete_format/equation in setup.json) is considerably more complex than a typical Python constructor. We recommend first reading the setup.json in the target example directory to understand the configurable options before deciding whether to override them via params.

Recommended Reading Order

If this is your first time using the algorithm library, we recommend exploring in the following order:

1. QFT / IQFT: understand basic quantum algorithm primitives
2. QPE: use QFT to estimate eigenvalues
3. LCU: combine multiple unitary operators
4. block_encode: encode a classical matrix into a quantum circuit
5. hamiltonian_simulation(): simulate quantum system evolution (default method='trotter')
6. solve(): solve linear systems (default method='hhl')
7. QSVT(): high-precision matrix function transformation