Linear Algebra Algorithms

Overview

The unitarylab_algorithms.linear_algebra package provides quantum algorithms for linear algebra tasks, including the Quantum Fourier Transform, linear system solving, non-unitary operator implementation, signal processing, and variational methods.

Return Values and Error Handling

When an algorithm completes normally, it returns status, circuit_path, plot, circuit, and algorithm-specific results. If conditions such as matrix dimension, Hermiticity, or post-selection are not satisfied, a ValueError or RuntimeError is raised, which the caller can catch as needed. Numerical quality should be judged in combination with the error, residual, fidelity, or optimizer results.

Quantum Fourier Transform (QFT)

Background

The Quantum Fourier Transform is the quantum analogue of the discrete Fourier transform, mapping the computational basis state $|j⟩$ to:

$\text{QFT}|j⟩=\frac{1}{\sqrt{N}}{\sum }_{k=0}^{N-1}{e}^{2\pi ijk/N}|k⟩$

QFT is a subroutine of Shor’s algorithm, QPE, HHL, and many other algorithms. This module also supports the inverse QFT (inverse=True).

Import

from unitarylab_algorithms import QFTAlgorithm

.run() Parameters

def run(self, n: int, state: np.ndarray = None, inverse: bool = False,
        backend='torch', device='cpu', dtype=np.complex128) -> Dict[str, Any]

Return Value

The actual output fields (merged into the final return dictionary):

{
    "status": "ok",
    "circuit_path": "....svg",   # save_circuit(qc.decompose()), i.e. the fully decomposed circuit diagram
    "plot": [{"format": "txt", "filename": "quantum_fourier_transform_algorithm_result.txt"}],  # from save_txt(), a result text file rather than a circuit diagram
    "circuit": <Circuit>,
    "Final state": ...,          # the final state vector obtained from simulation
    "Expected state": ...,       # the theoretical expected state computed via numpy's fft/ifft
    "Verification error": ...,   # the L2-norm error between the two (for logging/reference only, does not affect status)
    "Computation time (s)": ...
}

Example

import numpy as np
from unitarylab_algorithms import QFTAlgorithm
 
algo = QFTAlgorithm()
result = algo.run(n=4)
print(result['status'])  # 'ok'
 
# Inverse QFT
result_inv = algo.run(n=4, inverse=True)

Quick Demo

from unitarylab_algorithms.linear_algebra.qft.algorithm import test
test(n=4)

Notes

• The status field is always 'ok'; the Verification error is merely a comparison against the classical numpy.fft/numpy.ifft computation and does not affect status.
• The saved circuit diagram uses the circuit after qc.decompose() (all sub-gates are expanded and displayed), whereas the circuit field in the returned dictionary is the undecomposed original circuit object — the two have different structures.

HHL Algorithm for Solving Linear Systems

Background

The Harrow–Hassidim–Lloyd (HHL) algorithm solves the linear system $Ax=b$ with a time complexity of $O({\kappa }^2\log N/ϵ)$, where $\kappa$ is the condition number of $A$ and $N$ is the matrix dimension. The optimal classical algorithm is $O(N\kappa )$.

The algorithm uses QPE to encode the eigenvalues of $A$ into a phase register, then applies a controlled rotation to compute the reciprocal ${\lambda }^{-1}$ of each eigenvalue, then un-computes the phase register, and finally performs post-selection on the ancilla qubit and the phase register to extract the solution vector.

Requirements:

• $A$ must be a Hermitian square matrix, and its dimension must be a power of 2
• $b$ cannot be an all-zero vector
• The eigenvalues of $A$ cannot all be approximately zero at the same time

Import

from unitarylab_algorithms import HHLAlgorithm

.run() Parameters

def run(self, A: np.ndarray, b: np.ndarray, d: int,
        backend='torch', device='cpu', dtype=np.complex128) -> Dict[str, Any]

The evolution time $t$, the phase-grid starting point k_start, and whether to enable “signed phase mode” (automatically turned on when $A$ contains negative eigenvalues), among others, are all automatically determined by the algorithm based on the eigenspectrum of $A$, and are not parameters the user can pass in.

Return Value

{
    "status": "ok",
    "circuit_path": "....svg",
    "plot": [...],
    "circuit": <Circuit>,
    "Estimated solution (quantum)": ...,   # the solution vector obtained via quantum post-selection
    "Exact solution (classical)": ...,     # the classical solution from np.linalg.solve(A, b)
    "L2 error": ...,                       # the L2-norm difference between the two
    "Post-selection probability": ...,     # the post-selection success probability P(anc=1, phase=|0>)
    "Computation time (s)": ...
}

Example

import numpy as np
from unitarylab_algorithms import HHLAlgorithm
 
A = np.array([[0.8, 0], [0, 0.4]])
b = np.array([1, 2])
 
algo = HHLAlgorithm()
result = algo.run(A=A, b=b, d=11)
print(result['status'])

Quick Demo

from unitarylab_algorithms.linear_algebra.hhl.algorithm import test
test(A=[[0.8, 0], [0, 0.4]], b=[1, 2], d=11)

Notes

• A larger d yields higher precision but significantly increases circuit depth; if d is too small, the algorithm directly raises a ValueError (e.g. “Phase bits d=… insufficient to resolve minimum phase”), rather than returning a failure status.
• If the post-selection success probability is zero (i.e. the branch where ancilla=1 and the phase register is |0...0> never occurs), a ValueError("Post-selection probability is zero...") is raised.
• It is strongly recommended to use matrices with a small condition number (well-conditioned) for experimentation.

Linear Combination of Unitaries (LCU) Algorithm

Background

LCU implements the non-unitary operator $M={\sum }_j{\alpha }_j{U}_j$ by encoding the coefficients $\{{\alpha }_j\}$ into an ancilla “prepare” register (number of qubits $=\lceil {\log }_{2}m\rceil$, where $m$ is the number of terms), and conditionally applying the unitaries $\{U_j\}$ on that register. LCU is a core primitive for Hamiltonian simulation, QSVT, and block encoding.

Import

from unitarylab_algorithms import LCUAlgorithm

.run() Parameters

def run(self, alphas: List[float], unitaries: List[Circuit], n_sys: int,
        initial_state: Circuit = None,
        backend='torch', device='cpu', dtype=np.complex128) -> Dict[str, Any]

Return Value

{
    "status": "ok",
    "circuit_path": "....svg",
    "plot": [...],
    "circuit": <Circuit>,
    "Success probability": ...,       # the probability that all ancilla qubits are measured as |0...0>
    "Computation time (s)": ...,
    "Result state": ...                # the (unnormalized) state slice of the system register in the branch where the ancilla is |0...0>
}

Example

from unitarylab_algorithms import LCUAlgorithm

# Build the unitary circuits (e.g. Pauli operators), then pass them in together with the coefficients
algo = LCUAlgorithm()
result = algo.run(alphas=[0.6, 0.4], unitaries=[U_I, U_X], n_sys=1)
print(result['status'])

Quick Demo

from unitarylab_algorithms.linear_algebra.lcu.algorithm import test
test(n=1, alphas=[0.6, 0.4], paulis=['I', 'X'])

test() automatically constructs the corresponding unitary circuits based on paulis (a list of Pauli strings, e.g. 'IX', 'ZZ'): when a coefficient is negative, its absolute value is taken and a global phase gate U.gp(np.angle(-1)) (representing the sign flip) is appended so that the term is correctly represented, rather than being skipped; only terms whose coefficient is exactly 0 are dropped outright and excluded from the circuit. The paulis parameter also supports being passed as a string separated by the half-width comma ,, the full-width (Chinese) comma ，, the half-width semicolon ;, or the full-width (Chinese) semicolon ； (format_pauli_list automatically converts it into a list); the current implementation does not support the Chinese enumeration comma 、.

Notes

• is_success = success_prob > 1e-6 is the success criterion actually computed internally in the source code; it is used only for log printing (“Valid” / “Very low success rate…”), and is not reflected in the status field of the returned dictionary.
• The number of qubits in the ancilla register is automatically determined by n_anc = ceil(log2(m)) (where m is the length of alphas), and does not need to be specified by the user.

Quantum Signal Processing (QSP, Linear Algebra Version)

Background

linear_algebra.qsp.QSPAlgorithm uses classical optimization (scipy.optimize.minimize, L-BFGS-B) to search for a set of phase sequences $\{{\varphi }_k\}$, and implements a polynomial approximation of $\cos (t\cdot x)$ for a given eigenvalue $x$ using a single-qubit circuit (alternately applying the signal rotation $R_x(2\theta )$ and the phase rotation $R_z(2{\varphi }_k)$); it is a subroutine of QSVT and other advanced algorithms.

Naming Ambiguity Note: This module and hamiltonian_simulation.qsp.QSPHSAlgorithm are two completely different classes — both source directories are named qsp/, but they serve different purposes: this module (linear_algebra.qsp.QSPAlgorithm) only performs a $\cos (t\cdot x)$ approximation for a scalar eigenvalue $x\in [-1,1]$ on a single qubit, and is used to demonstrate/verify the QSP phase-finding procedure itself; hamiltonian_simulation.qsp.QSPHSAlgorithm, on the other hand, performs genuine time-evolution simulation $e^{-iHt}$ for an arbitrary Hermitian matrix $H$. Be sure to double-check the import path.

Import

from unitarylab_algorithms import QSPAlgorithm

.run() Parameters

def run(self, t: float, d: int, x: float = 0.5,
        backend='torch', device='cpu', dtype=np.complex128) -> Dict[str, Any]

Return Value

{
    "status": "ok",
    "circuit_path": "....svg",
    "plot": [...],
    "circuit": <Circuit>,
    "Estimated value": ...,     # the complex amplitude obtained from circuit simulation
    "Ideal value": ...,         # the theoretical value np.cos(t * x)
    "Absolute error": ...,
    "Computation time (s)": ...
}

Example

from unitarylab_algorithms import QSPAlgorithm
 
algo = QSPAlgorithm()
result = algo.run(t=1.0, d=10, x=0.5)
print(result['status'])

Quick Demo

from unitarylab_algorithms.linear_algebra.qsp.algorithm import test
test(t=1.0, d=10, x=0.5)

Notes

• The initial guess for the phase sequence search is np.random.randn(d+1) * 0.1, and no random seed is set, so repeated calls with the same (t, d, x) may converge to different phase solutions and yield different Absolute error values (the results are not fully reproducible); this is the same category of cause as the non-reproducibility issue with qDrift in the Hamiltonian simulation package.
• L-BFGS-B is a local optimizer with no multi-start retry mechanism; when d is small or t is large, it may converge to a poor local optimum, causing Absolute error to be large, though status is still always 'ok'.

QSVT-Based Linear System Solver Algorithm

Background

Quantum Singular Value Transformation (QSVT) generalizes QSP to arbitrary matrices via block encoding, enabling polynomial transformations of singular values. The QSVT linear solver uses this framework to implement the matrix inverse function $x\mapsto 1/x$ on the singular values of $A$, solving $Ax=b$ with near-optimal query complexity.

Import

from unitarylab_algorithms import QSVTLinearSolverAlgorithm

.run() Parameters

def run(self, A, b, epsilon,
        backend='torch', device='cpu', dtype=np.complex128) -> Dict[str, Any]

Return Value

{
    "status": "ok",
    "circuit_path": "....svg",
    "plot": [...],
    "circuit": <Circuit>,
    "Solution vector": ...,          # the solution vector
    "Scaling factor applied": ...,   # the scaling factor applied internally by QSVTSolver
    "Simulation time (s)": ...
}

Example

import numpy as np
from unitarylab_algorithms import QSVTLinearSolverAlgorithm
 
A = np.array([[0.8, 0.1], [0.1, 0.6]])
b = np.array([1.0, 0.5])
 
algo = QSVTLinearSolverAlgorithm()
result = algo.run(A=A, b=b, epsilon=1e-3)
print(result['status'])

Quick Demo

from unitarylab_algorithms.linear_algebra.qsvt_qlsa.algorithm import test
test(A=[[0.8, 0], [0, 0.4]], b=[1., 2.], epsilon=0.01)

Notes

• The core logic of this algorithm (matrix validity checks, block encoding, QSVT phase angle solving) is fully delegated to the library’s internal unitarylab.library.linear_solver.QSVTSolver; QSVTLinearSolverAlgorithm.run() itself is merely a thin wrapper and performs no explicit validation of A/b at the Python layer (e.g. whether they are Hermitian, whether dimensions match, etc.) — any validation failure surfaces as an exception raised internally by QSVTSolver.
• The module-level test() function has no independent docstring.

Variational Quantum Linear Solver (VQLS)

Background

VQLS is a variational algorithm suited for near-term quantum devices (NISQ), used to solve the linear system $Ax=b$: it parameterizes the solution state $|x(\theta )⟩$ with a hardware-efficient ansatz (RY-RZ rotations plus a ring of entangling CNOT layers), and iteratively optimizes the parameters by minimizing a residual-based cost function. The coefficient matrix is automatically Pauli-decomposed internally via unitarylab.library.pauli_operator.pauli_string_decomposition ($A={\sum }_l{c}_l{A}_l$), so the user does not need to supply the decomposition manually.

Import

from unitarylab_algorithms import VQLSAlgorithm

.run() Parameters

def run(self, A: np.ndarray, b: np.ndarray, cost_function: str = "local_ht",
        n_layers: int = 4, maxiter: int = 500, tol: float = 1e-6,
        seed: int = 42, epsilon: Optional[float] = None,
        backend='torch', device='cpu', dtype=np.complex128) -> Dict[str, Any]

Return Value

{
    "status": "ok",
    "circuit_path": "....svg",
    "plot": [...],
    "circuit": <Circuit>,
    "Fidelity": ...,                     # fidelity between the quantum solution and the classical solution; None if A is singular
    "Ax Fidelity": ...,                  # fidelity between the normalized A|x_quantum> and |b>
    "Cost Function": "local_ht" | "local_classical" | "global",
    "Condition Number": ...,             # np.linalg.cond(A)
    "Solution State (Quantum)": ...,
    "Solution State (Classical)": ...,   # may be None (when A is singular)
    "Computation Time (s)": ...,
    "Cost History": [...],               # the history of every cost-function evaluation
    "Early Stopped": True | False
}

Example

import numpy as np
from unitarylab_algorithms import VQLSAlgorithm
 
A = np.array([[1.5, 0.2], [0.2, 1.8]])
b = np.array([1.0, 0.5])
 
algo = VQLSAlgorithm()
result = algo.run(A=A, b=b, cost_function="local_ht", n_layers=4, maxiter=500)
print(result['status'])       # always 'ok'

Quick Demo

from unitarylab_algorithms.linear_algebra.vqls.algorithm import test
test(cost_function="local_ht", n_layers=4, maxiter=500, tol=1e-6, seed=42)
# When A and b are omitted, they default to [[1.5, 0.2], [0.2, 1.8]] and [1.0, 0.5] respectively

Notes

• The optimizer’s convergence status is written to the saved result text; API callers can also combine Fidelity, Cost History, and Early Stopped to judge the quality of the solution.
• The web interface parameter names differ from the Python parameter names: the UI parameter names exposed in parameters.json are max_iterations and tolerance, but the actual keyword arguments of .run() on the Python side are maxiter and tol, respectively.
• The web interface is missing the A/b/seed/epsilon parameters: unlike HHL and QSVT-QLSA, VQLS’s parameters.json exposes only the four items cost_function, n_layers, max_iterations, and tolerance, and does not provide input fields for the matrix A, the vector b, the random seed seed, or the early-stopping threshold epsilon. This means running VQLS via the current web interface uses the default built-in 2x2 example problem in the code, and the user cannot customize the linear system to be solved through the interface.
• The parameters.json default value for n_layers is 2, which does not match the Python .run()/test() default value of 4.
• When A is a singular matrix (np.linalg.solve raises LinAlgError), both Fidelity and Solution State (Classical) return None rather than raising an error, and the algorithm still continues to output the quantum-side results.
• The early-stopping mechanism has a known limitation: the source code comments explicitly note that COBYLA is implemented in Fortran, so a Python exception cannot actually interrupt its internal loop from within the cost-function callback. As a result, setting _early_stop_flag inside _cost_with_stop does not cause COBYLA to exit immediately — it merely records the flag early for a later check, and the actual optimization may still keep running until maxiter (this behavior is acknowledged in the source code comments as a known design limitation).

Discrete Adiabatic Quantum Linear Solver (AQC)

Background

AQCAlgorithm implements the discrete adiabatic quantum linear system solver (QLSP): via Trotterized discrete adiabatic evolution, it solves $Ax=b$ on $n$ system qubits and 5 ancilla qubits. The algorithm performs an SVD block encoding of the coefficient matrix $A$ (normalized so that $\Vert A{\Vert }_2\le 1$, then embedded into a $2N\times 2N$ unitary matrix), constructs the unitary $U_b$ that prepares $|b⟩$ using a Householder reflection, and then evolves step by step for $T$ steps according to the adiabatic schedule parameter $p$, finally performing post-selection on the ancilla register to extract the solution vector.

This algorithm does not accept a user-supplied $A$ or $b$: unlike the other solvers in this package (HHL, QSVT-QLSA, VQLS), AQCAlgorithm.run() does not expose A/b parameters at all — at runtime it first executes np.random.seed(42), and then automatically generates a well-conditioned Hermitian matrix A = randn(N,N) (symmetrized and then with 10·I added, to ensure diagonal dominance and a small condition number), together with a random right-hand-side vector b = randn(N). Because the seed is fixed, the A and b generated for a given n are exactly the same on every run (reproducible), but the user cannot use it to solve a linear system of their own choosing — it can only be used to demonstrate/benchmark the discrete adiabatic solving procedure itself.

Import

from unitarylab_algorithms import AQCAlgorithm

.run() Parameters

def run(self, n: int = 2, T: int = 0, p: float = 1.4,
        backend='torch', device='cpu', dtype=np.complex128) -> Dict[str, Any]

Return Value

General case ($\kappa \not\approx 1$):

{
    "status": "ok",
    "circuit_path": ["....svg", "....svg"],   # [full circuit diagram, single-step slice circuit diagram] (isomorphic to the Trotter algorithm)
    "plot": [...],
    "circuit": <Circuit>,
    "Quantum Solution (x)": ...,
    "Classical Solution": ...,                 # the classical solution from np.linalg.solve(A, b)
    "Residual Norm ||Ax-b||": ...,
    "Error vs Classical (L2)": ...,
    "Internal Scale Factor": ...,
    "Post-selection Amplitude": ...,
    "Simulation Time (s)": ...,
    "Elapsed Time (s)": ...
}

When the internally generated $A$ happens to have condition number $\kappa \approx 1$ (i.e. it is approximately a scalar multiple of the identity matrix), the algorithm takes a trivial fast path: it directly sets x = b (rescaled appropriately), skipping the quantum circuit evolution entirely. In this case, the structure of the returned dictionary differs slightly — it lacks the two entries Post-selection Amplitude and Simulation Time (s), and circuit_path is a single string rather than a list of length 2 (because only the initialization circuit, without the adiabatic evolution steps, is saved).

Example

from unitarylab_algorithms import AQCAlgorithm
 
algo = AQCAlgorithm()
result = algo.run(n=2, T=0, p=1.4)   # T=0 → automatically choose the number of steps based on the condition number
print(result['status'])

Quick Demo

from unitarylab_algorithms.linear_algebra.aqc.algorithm import test
test(n=2, T=0, p=1.4)

Notes

• Post-selection projects onto the |10000> branch of the ancilla register (4 qubits at 0, 1 qubit at 1); if the amplitude of this branch tends to zero, a RuntimeError("Post-selection returned a near-zero state...") is raised, suggesting an increase in T or indicating that the system is ill-conditioned.
• When T is not specified (i.e. T=0), it is automatically set to $T\approx \lceil 10\kappa \rceil$ and rounded to be even; if a T explicitly passed in is odd, it will not be automatically adjusted.
• Circuit size grows linearly with T (11 sub-operations per step); when the full circuit is quite deep, the single-step slice diagram aqc_one_step.svg can be viewed to understand the structure of a single step.
• While saving the circuit diagram, the source code temporarily redirects sys.stdout to an in-memory buffer (to suppress noisy output during plotting), restoring it once complete — this is a normal implementation detail.

Version Reference

To import multiple algorithm classes from this module at once, you can batch-import directly from the top-level package:

from unitarylab_algorithms import (
    HHLAlgorithm, LCUAlgorithm, QFTAlgorithm, QSPAlgorithm,
    QSVTLinearSolverAlgorithm, VQLSAlgorithm, AQCAlgorithm,
)