Solving the Quantum Linear Systems Problem (QLSP) using AQC
Overview
Adiabatic Quantum Computing (AQC) leverages the adiabatic theorem to solve computational problems by gradually evolving a quantum system from an initial ground state to the ground state of a problem-specific Hamiltonian (see our AQC tutorial).
This tutorial focuses on applying the AQC approach to solve the Quantum Linear Systems Problem (QLSP), a cornerstone problem in quantum computing with significant applications in fields like machine learning, physics, and optimization.
Specifically, we aim to demonstrate how AQC can be utilized to approximate the solution to the QLSP [1]. This problem involves finding a quantum state that corresponds to the solution of a linear system of equations. The tutorial provides a structured overview of the QLSP, its mathematical formulation, and the steps needed to transform it into an eigenvalue problem, laying the foundation for solving it within the AQC framework.
The following demonstration we will follow the paper [1]. The notebook was written in collaboration with Prof. Lin Lin and Dr. Dong An, the authors of the paper.
1. Problem Statement
Given a Hermitian positive-definite matrix
We are given:
-
Matrix
, an invertible Hermitian and positive-definite matrix with condition number and . -
Vector
, a normalized vector. -
Target Error
, specifying the desired accuracy.
The goal is to prepare a quantum state
2. Transformation into AQC
The QLSP is converted into an equivalent eigenvalue problem to leverage quantum computation. This involves the following steps:
2.1 Constructing
Define:
where
Key properties:
-
is Hermitian. -
The null space of
: , where
The null space of a matrix
is the set of all vectors such that:
With regards to eigenstates and eigenvalues: - The null space corresponds to the eigenspace of
associated with the eigenvalue . - Any vector in the null space is an eigenvector of with eigenvalue .
2.2 Constructing
Define:
where
Key properties:
-
If
, then . -
Null space of
: , where
2.3 Adiabatic Interpolation
Construct an interpolation Hamiltonian:
where
2.4 Spectral Gap
-
is a projection operator, and the spectral gap between and the rest of the eigenvalues of is . -
For
, the gap between and the rest of the eigenvalues is bounded from below by . [1]
2.5 Adiabatic Evolution and Null Space
Notice that there is a degeneracy in the number of Null states (unlike the regular adiabatic algorithm usage where we are typically looking at a single grounde state)
We also note that for any
Therefore, there exist an additional statestate
In particular:
-
At
, , the initial state. -
At
, , the solution state.
Thus,
3. AQC Approach to solve QLSP
The adiabatic quantum algorithm prepares the zero-energy state
-
Initialize in the ground state of
, i.e., . -
Slowly evolve the system by varying
from to . -
At the end of the evolution, the system will approximate
, embedding - the solution of the QLSP.
4. Goals
-
Set up a QLSP example: Derive
, and define the interpolation Hamiltonian -
Quantum Circuit Design: Implement Hamiltonian simulation for
. -
Evaluate results: Compare quantum simulation results with the numeric calculation.
Let’s begin with the mathematical setup and proceed to implementation!
Setting up a QLSP example where A is a 4x4 matrix:
Setting an example
For simplicity, we first assume A is Hermitian and positive definite
import numpy as np
# Define matrix A and vector b
A = np.array([[4, 1, 2, 0], [1, 3, 0, 1], [2, 0, 3, 1], [0, 1, 1, 2]])
b = np.array([12, 10, 17, 26])
As a purely mathematical pre-processing step, we will calculate the condition number
In practical scenarios, the condition number is often approximated or known beforehand based on external factors or prior knowledge.
import numpy as np
def compute_condition_number(A):
"""
Computes the condition number (κ) of a matrix A.
Parameters:
A (numpy.ndarray): Input square matrix.
Returns:
float: The condition number of A.
"""
try:
# Compute the norm of A (2-norm, largest singular value)
norm_A = np.linalg.norm(A, 2)
# Compute the inverse of A
A_inv = np.linalg.inv(A)
# Compute the norm of A^-1 (2-norm)
norm_A_inv = np.linalg.norm(A_inv, 2)
# Compute the condition number κ
condition_number = norm_A * norm_A_inv
return condition_number
except np.linalg.LinAlgError:
return float("inf") # Return infinity if the matrix is singular
condition_number = compute_condition_number(A)
print("Condition number:", condition_number)
Condition number: 14.472337634948623
Construct and :
The setup_QLSP
function prepares the necessary Hamiltonians and normalized components to solve the Quantum Linear Systems Problem (QLSP).
The built-in function matrix_to_hamiltonian
, used within setup_QLSP
function, encodes the hamitonian matrix into a sum of Pauli strings that will be used to exponentiate the hamiltonians with a product formula (Suzuki-Trotter) in the next step.
from classiq import *
from classiq.execution import *
def setup_QLSP(A, b):
# Normalize A
norm_A = np.linalg.norm(A, "fro")
A_normalized = A / norm_A
# Normalize vector b
b_normalized = b / np.linalg.norm(b)
# Create the outer product of b
outer_product_b = np.outer(b_normalized, b_normalized)
# Define the identity matrix I with the same size as b
identity_matrix = np.eye(len(b))
# Compute Qb = I - outer_product_b
Qb = identity_matrix - outer_product_b
# Define the Pauli-X (σx) and Pauli-Y (σy) matrices
pauli_x = np.array([[0, 1], [1, 0]])
pauli_y = np.array([[0, -1j], [1j, 0]])
# Define Pauli plus and minus operators
pauli_plus = 0.5 * (pauli_x + 1j * pauli_y)
pauli_minus = 0.5 * (pauli_x - 1j * pauli_y)
# Compute the tensor product of Pauli-X and Qb
H0 = np.kron(pauli_x, Qb)
# Compute A*Qb and Qb*A
A_Qb = np.dot(A, Qb)
Qb_A = np.dot(Qb, A)
# Compute the tensor products
tensor_plus = np.kron(pauli_plus, A_Qb)
tensor_minus = np.kron(pauli_minus, Qb_A)
# Define H1 as the sum of the two tensor products
H1 = tensor_plus + tensor_minus
HO_HAMILTONIAN = matrix_to_hamiltonian(H0)
H1_HAMILTONIAN = matrix_to_hamiltonian(H1)
return H0, H1, HO_HAMILTONIAN, H1_HAMILTONIAN, A_normalized, b_normalized
# Setup
H0, H1, HO_HAMILTONIAN, H1_HAMILTONIAN, A_normalized, b_normalized = setup_QLSP(A, b)
Define the interpolation Hamiltonian:
For the sake of simplicity, we will first use the "vanilla AQC" linear scheduling function
As
# Define the time-dependent Interpolated Hamiltonian, where T is the total evolution time
def hamiltonian_t(H0, H1, t, T):
s = t / T
return ((1 - s) * H0) + (s * H1)
Spectral gap analysis:
From the quantum adiabatic theorem [3, Theorem 3] the formula for the adiabatic error bound
See Appendix A for a detailed explanation of the components.
The formula shows that the adiabatic error is minimized when:
-
The total runtime
is large (slow evolution). -
The spectral gap
is large (well-separated ground and excited states). -
The derivatives
and are small (smooth Hamiltonian changes).
The following function plots the spectral gap
import matplotlib.pyplot as plt
def plot_eigenvalues_evolution(ylim=None):
time_steps = np.linspace(0, 1, 100) # Discrete time steps
# Store eigenvalues at each time step
eigenvalues = []
# Calculate eigenvalues across all time steps
for t in time_steps:
H_t = hamiltonian_t(H0, H1, t, 1)
eigvals = np.linalg.eigvalsh(H_t) # Sorted real eigenvalues
eigenvalues.append(eigvals)
# Convert eigenvalues list to a NumPy array for easier manipulation
eigenvalues = np.array(eigenvalues)
# Add small offsets to separate close eigenvalues visually
offsets = np.linspace(
-0.05, 0.05, eigenvalues.shape[1]
) # Small offsets for each eigenvalue line
# Plot the eigenvalues across time steps
plt.figure(figsize=(10, 6))
for i in range(eigenvalues.shape[1]):
plt.plot(time_steps, eigenvalues[:, i] + offsets[i], label=f"Eigenvalue {i+1}")
# Highlight degenerate eigenvalues (if any)
for step_idx, t in enumerate(time_steps):
unique_vals, counts = np.unique(eigenvalues[step_idx], return_counts=True)
# Apply y-axis limits if provided
if ylim:
plt.ylim(ylim)
# Customize the plot
plt.xlabel("Time (t)", fontsize=12)
plt.ylabel("Eigenvalues", fontsize=12)
plt.title("Eigenvalues Evolution Across $s$", fontsize=14)
plt.grid()
# Move the legend to the side
plt.legend(bbox_to_anchor=(1.05, 1), loc="upper left", borderaxespad=0.0)
plt.tight_layout()
# Show the plot
plt.show()
plot_eigenvalues_evolution()
Should we want to focus on the spectral gap from the Null states, we can zoom-in on our plot:
plot_eigenvalues_evolution(ylim=(-1.5, 1.5))
We can visually observe from the above that although the spectral gap does change throughout TOTAL_EVOLUTION_TIME
). However, if we simply assume
TOTAL_EVOLUTION_TIME = 7
AQC Implementation:
Since
\(
One straightforward approach to achieve this is by using the Trotter splitting method. The lowest order approximation takes the form:
which can further be approximated as:
where
\(
It is proved in [2] that the error of such an approximation is
which indicates that to achieve an
(Note that in our case
Building the Quantum model:
Using the Classiq platform, we implement the adiabatic path with Suzuki-Trotter decomposition for Hamiltonian exponentiation.
In our model TOTAL_EVOLUTION_TIME
and NUM_STEPS
.
In the Trotter implementation
NUM_STEPS = 50
We are now ready to build our quantum model.
The function adiabatic_evolution_qfunc
implements the adiabatic path with Suzuki-Trotter decomposition for Hamiltonian exponentiation:
@qfunc(generative=True)
def adiabatic_evolution_qfunc(
H0: CArray[PauliTerm],
H1: CArray[PauliTerm],
evolution_time: CInt,
num_steps: CInt,
qba: QArray[QBit],
):
# Time step for each increment
delta_t = evolution_time / num_steps
for step in range(num_steps):
t = step * delta_t
suzuki_trotter(
HO_HAMILTONIAN,
evolution_coefficient=delta_t * (1 - t / evolution_time),
order=1,
repetitions=1,
qbv=qba,
)
suzuki_trotter(
H1_HAMILTONIAN,
evolution_coefficient=delta_t * (t / evolution_time),
order=1,
repetitions=10,
qbv=qba,
)
To solve the QLSP we first prepare
def get_model(H0, H1, b, evolution_time, num_steps):
@qfunc
def main(qba: Output[QArray[QBit]]):
prepare_state(
probabilities=(np.abs(np.kron(np.array([1, 0]), b)) ** 2).tolist(),
bound=0,
out=qba,
)
adiabatic_evolution_qfunc(H0, H1, evolution_time, num_steps, qba)
execution_preferences = ExecutionPreferences(
num_shots=1,
backend_preferences=ClassiqBackendPreferences(
backend_name=ClassiqSimulatorBackendNames.SIMULATOR_STATEVECTOR
),
)
return create_model(main, execution_preferences=execution_preferences)
qmod = get_model(
HO_HAMILTONIAN, H1_HAMILTONIAN, b_normalized, TOTAL_EVOLUTION_TIME, NUM_STEPS
)
Synthesize, verify and execute
Synthesize the model into a Quantum Program, verify it and execute it on a state vector simulator:**
qprog = synthesize(qmod)
write_qmod(qmod, "solving_qlsp_with_aqc", decimal_precision=5)
show(qprog)
result_state_vector = execute(qprog).result_value().state_vector
Opening: https://platform.classiq.io/circuit/2rwe34neCtLDeNlc1hOH6HYVXdn?version=0.66.0
Evaluate results:
def plot_state_probabilities(title, x, color="b"):
# Ensure x is a numpy array and normalized
x = np.array(x)
# Calculate probabilities
probabilities = np.abs(x) ** 2
# Create labels for the states
labels = [f"|{i}>" for i in range(len(x))]
# Plot the probabilities
plt.bar(labels, probabilities, color=color, alpha=0.7)
plt.xlabel("States")
plt.ylabel("Probabilities")
plt.title(title)
plt.xticks(rotation=45)
plt.ylim(0, 1)
plt.show()
def compare_states(
state1, state1_label, state2, state1_labe2, color1="gold", color2="b"
):
# Plot a histogram of each state probabilities
plot_state_probabilities(state1_label, state1, color1)
plot_state_probabilities(state1_labe2, state2, color2)
# Check the overlap between states
overlap = np.abs(np.vdot(state1, state2)) ** 2
print(f"Similarity of results: {overlap:.4f}")
# Print the solution vector x
x = np.linalg.solve(A_normalized, b_normalized)
print("Solution vector x:")
normalized_x = x / np.linalg.norm(x)
print(normalized_x)
# Convert dictionary values to complex numbers
print("State vector:")
state_vector = np.array([complex(value) for value in result_state_vector.values()])
print(state_vector)
compare_states(
state_vector,
"Quantum simulator state_vector",
np.kron(np.array([1, 0]), normalized_x),
"Classical solution to normalized_x",
)
Solution vector x:
[ 0.32300564 -0.23491319 -0.22423532 0.88893288]
State vector:
[ 0.01631432-3.08220039e-01j -0.01056133+1.99531093e-01j
-0.01705876+3.22284489e-01j 0.04545515-8.58766326e-01j
0.13241748+7.00895691e-03j -0.02525213-1.33661442e-03j
0.00707045+3.74244064e-04j -0.05315043-2.81329231e-03j]
Similarity of results: 0.9672
By comparing the quantum-computed results with the mathematically expected solution, we observe a good alignment - showcasing the potential of the AQC approach for solving linear problems.
We can observe the runtime of the above implementation by analyzing the depth parameter from the transpiled circuit data of our quantum program
print("Program depth:", QuantumProgram.from_qprog(qprog).transpiled_circuit.depth)
Program depth: 33308
For alternative QLSP configurations, selecting different values for
# Define matrix A and vector b
A = np.array(
[
[3.9270525, 1.06841123, 2.09661281, -0.10400811],
[1.06841123, 2.93584295, -0.0906049, 1.09754032],
[2.09661281, -0.0906049, 2.87204449, 1.13774997],
[-0.10400811, 1.09754032, 1.13774997, 1.85170585],
]
)
b = np.array([12, 10, 17, 26])
condition_number = compute_condition_number(A)
print("Condition number:", condition_number)
H0, H1, HO_HAMILTONIAN, H1_HAMILTONIAN, A_normalized, b_normalized = setup_QLSP(A, b)
Condition number: 5722470278.425136
plot_eigenvalues_evolution(ylim=(-1.5, 1.5))
Although the condition number is higher, the spectral gap remains relatively large. However, using the same values of
qmod = get_model(
HO_HAMILTONIAN, H1_HAMILTONIAN, b_normalized, TOTAL_EVOLUTION_TIME, NUM_STEPS
)
qprog = synthesize(qmod)
show(qprog)
result_state_vector = execute(qprog).result_value().state_vector
Opening: https://platform.classiq.io/circuit/2rweHDjGqIVTawYYsk3Uwic8S9J?version=0.66.0
# Print the solution vector x
x = np.linalg.solve(A_normalized, b_normalized)
print("Solution vector x:")
normalized_x = x / np.linalg.norm(x)
print(normalized_x)
# Convert dictionary values to complex numbers
print("State vector:")
state_vector = np.array([complex(value) for value in result_state_vector.values()])
print(state_vector)
compare_states(
state_vector,
"Quantum simulator state_vector",
np.kron(np.array([1, 0]), normalized_x),
"Classical solution to normalized_x",
)
Solution vector x:
[ 0.42009163 -0.39396804 -0.55637591 0.59896415]
State vector:
[-0.02338635-0.37524301j 0.02140731+0.34348851j 0.04152884+0.6663462j
-0.02503025-0.40161994j 0.28004751-0.01745346j -0.13146114+0.00819308j
-0.18310263+0.01141155j 0.05280268-0.00329083j]
Similarity of results: 0.8209
print("Program depth:", QuantumProgram.from_qprog(qprog).transpiled_circuit.depth)
Program depth: 33308
Choosing
TOTAL_EVOLUTION_TIME = 10
NUM_STEPS = 100
qmod = get_model(
HO_HAMILTONIAN, H1_HAMILTONIAN, b_normalized, TOTAL_EVOLUTION_TIME, NUM_STEPS
)
qprog = synthesize(qmod)
show(qprog)
result_state_vector = execute(qprog).result_value().state_vector
Opening: https://platform.classiq.io/circuit/2rwei3h5uCQJMppqbZQUFzX2v9y?version=0.66.0
# Print the solution vector x
x = np.linalg.solve(A_normalized, b_normalized)
print("Solution vector x:")
normalized_x = x / np.linalg.norm(x)
print(normalized_x)
# Convert dictionary values to complex numbers
print("State vector:")
state_vector = np.array([complex(value) for value in result_state_vector.values()])
print(state_vector)
compare_states(
state_vector,
"Quantum simulator state_vector",
np.kron(np.array([1, 0]), normalized_x),
"Classical solution to normalized_x",
)
Solution vector x:
[ 0.42009163 -0.39396804 -0.55637591 0.59896415]
State vector:
[-0.44075799-0.27456092j 0.2986469 +0.18603581j 0.45432412+0.28301165j
-0.45996295-0.28652425j 0.06242484-0.10021182j -0.02261307+0.0363012j
-0.05203803+0.08353766j 0.01609333-0.02583494j]
Similarity of results: 0.9587
print("Program depth:", QuantumProgram.from_qprog(qprog).transpiled_circuit.depth)
Program depth: 67208
As expected, the overall circuit depth increased to achieve a high degree of similarity in the results.
Important to note: The above implementation example aims to provide an intuitive understanding of the principles behind solving quantum linear solver problems with the adiabatic quantum evolution and the associated error bounds. However, for simplicity and accessibility, several aspects of the implementation are not optimal:
- Suzuki-Trotter Decomposition: We utilized the first order Suzuki-Trotter approximation for time evolution. As such (as mentioned above), a more optimal approach based on the discrete adiabatic theorem [4] could be implemented.
- Schedule Function: In the above example we used the vanilla AQC scheduling function. Using more sophisticated scheduling functions as suggested in [1] will improve runtime.
- Brute-Force Encoding: The encoding of Pauli operators in this tutorial is direct and unoptimized, scaling exponentially with system size. Other encoding techniques (such as suggested in [4]) will be more efficient.
These choices were made to prioritize conceptual clarity over computational efficiency. Next step would be to apply state pf the art techniques and show improvement in gate complexity and runtime.
Appendix A: Explanation of the adiabatic error bound components
The adiabatic error bound
Detailed Explanation of Components:
-
: -
Represents the adiabatic error bound at a specific point
. -
Quantifies the deviation of the quantum state from the ground state during evolution.
-
: -
A proportionality constant that depends on system-specific properties, such as the dimensionality and scaling of norms.
-
: -
The operator norm of the first derivative of the Hamiltonian,
, evaluated at . -
Indicates how quickly the Hamiltonian changes initially.
-
: -
The total runtime of the adiabatic evolution.
-
Larger
values reduce the error, as slower evolution aids adiabaticity. -
: -
The spectral gap at
, defined as the energy difference between the ground state and the first excited state of . -
Larger gaps improve the adiabatic process.
-
: -
The operator norm of the first derivative of
at an intermediate point . -
Reflects how fast the Hamiltonian is changing during evolution.
-
: -
The spectral gap at the point
, mapped via the function . -
: -
The operator norm of the second derivative of the Hamiltonian,
, at a point . -
Captures the curvature or acceleration of the Hamiltonian's evolution.
-
: -
An integral from
to , summing contributions of the Hamiltonian’s derivatives over the path. -
Accounts for cumulative effects of the Hamiltonian's changes during evolution.
-
and :-
Higher powers of the spectral gap at
. -
Larger gaps (in
and ) significantly reduce the adiabatic error.
-
[1]: An, D. and Lin, L. “Quantum Linear System Solver Based on Time-Optimal Adiabatic Quantum Com- puting and Quantum Approximate Optimization Algorithm.” ACM Trans. Quantum Comput. 3 (2022). arXiv:1909.05500..
[2]: Wim van Dam, Michele Mosca, and Umesh Vazirani. 2001. How powerful is adiabatic quantum computation? In Proceedings 42nd IEEE Symposium on Foundations of Computer Science. IEEE, Piscataway, NJ, 279–287
[3]: Sabine Jansen, Mary-Beth Ruskai, and Ruedi Seiler. 2007. Bounds for the adiabatic approximation with applications to quantum computation. J. Math. Phys. 48, 10 (2007), 102111
[4]: The discrete adiabatic quantum linear system solver has lower constant factors than the randomized adiabatic solver Pedro C.S. Costa, Dong An, Ryan Babbush, Dominic Berry