Quantum linear solver based on QSVT
The code here can be integrated as part of a larger CFD solver, e.g., as in qc-cfd repository.
In particular, instead of calling a classical solver, e.g., x = sparse.linalg.spsolve(mat_raw_scr, b_raw), one can call the quantum solver qsvt_solver(mat_raw_scr, b_raw,...).
We implemented two versions for block-encoding, one based on Pauli decomposition of the matrix, and another one based on decomposing the matrix to a finite set of diagonals.
import time
import matplotlib.pyplot as plt
import numpy as np
from cheb_utils import *
from classical_functions_be import *
from quantum_functions_be import *
from scipy import sparse
from classiq import *
from classiq.applications.chemistry.op_utils import qubit_op_to_pauli_terms
np.random.seed(53)
Define functions to get properties of the block-encoding.
def get_pauli_be(mat_raw_scr):
"""
Get relevant block-encoding properties for `lcu_paulis_graycode` block encoding,
Parameters
----------
mat_raw_scr : scipy.sparse.spmatrix
Square sparse matrix of shape (N, N), real or complex, to be block-encoded.
Returns
-------
data_size : int
Size of the data variable.
block_size : int
Size of the block variable.
be_scaling_factor : float
The scaling factor of the block-encoding unitary
BlockEncodedState : QStruct
QSVT-compatible QStruct holding the quantum variables, with fields:
- data : QNum[data_size]
- block : QNum[block_size]
be_qfunc : qfunc
Quantum function that implements the block encoding. Signature:
be_qfunc(be: BlockEncodedState) → None
"""
rval = mat_raw_scr.data
col = mat_raw_scr.indices
rowstt = mat_raw_scr.indptr
nr = mat_raw_scr.shape[0]
raw_size = mat_raw_scr.shape[0]
data_size = max(1, (raw_size - 1).bit_length())
# Set to_symmetrize=False, since we are working with QSVT
paulis_list, transform_matrix = initialize_paulis_from_csr(
rowstt, col, data_size, to_symmetrize=False
)
qubit_op = eval_pauli_op(paulis_list, transform_matrix, rval)
qubit_op.compress(1e-12)
hamiltonian = qubit_op_to_pauli_terms(qubit_op)
be_scaling_factor = sum([np.abs(term.coefficient) for term in hamiltonian.terms])
block_size = size = max(1, (len(hamiltonian.terms) - 1).bit_length())
hamiltonian = hamiltonian * (1 / be_scaling_factor)
class BlockEncodedState(QStruct):
data: QNum[data_size]
block: QNum[block_size]
@qfunc
def be_qfunc(be: BlockEncodedState):
lcu_paulis_graycode(hamiltonian.terms, be.data, be.block)
return data_size, block_size, be_scaling_factor, BlockEncodedState, be_qfunc
def get_banded_diags_be(mat_raw_scr):
"""
Get relevant block-encoding properties for `block_encode_banded` block encoding,
Parameters
----------
mat_raw_scr : scipy.sparse.spmatrix
Square sparse matrix of shape (N, N), real or complex, to be block-encoded.
Returns
-------
data_size : int
Size of the data variable.
block_size : int
Size of the block variable.
be_scaling_factor : float
The scaling factor of the block-encoding unitary
BlockEncodedState : QStruct
QSVT-compatible QStruct holding the quantum variables, with fields:
- data : QNum[data_size]
- block : QNum[block_size]
be_qfunc : qfunc
Quantum function that implements the block encoding. Signature:
be_qfunc(be: BlockEncodedState) → None
"""
raw_size = mat_raw_scr.shape[0]
data_size = max(1, (raw_size - 1).bit_length())
offsets, diags, diags_maxima, prepare_norm = get_be_banded_data(mat_raw_scr)
block_size = int(np.ceil(np.log2(len(offsets)))) + 1
be_scaling_factor = prepare_norm
class BlockEncodedState(QStruct):
data: QNum[data_size]
block: QNum[block_size]
@qfunc
def be_qfunc(be: BlockEncodedState):
block_encode_banded(
offsets=offsets,
diags=diags,
prep_diag=diags_maxima,
block=be.block,
data=be.data,
)
return data_size, block_size, be_scaling_factor, BlockEncodedState, be_qfunc
The solvers were developed in the framework of exploring their performance in hybrid CFD schemes. For simplicity, it is assumed that all the properties of the matrices are known explicitly. In particular, we calculate its sigular values for identyfing the range in which we apply the inversion polnomial.
def get_svd_range(mat_raw_scr):
mat_raw = mat_raw_scr.toarray()
svd = np.linalg.svd(mat_raw)[1]
w_min = min(svd)
w_max = max(svd)
return w_min, w_max
def qsvt_solver(
mat_raw_scr,
b_raw,
poly_degree,
be_method="banded",
cheb_approx_type="optimized",
preferences=Preferences(),
constraints=Constraints(),
qmod_name=None,
):
SCALE = 0.5
b_norm = np.linalg.norm(b_raw) # b normalization
b_normalized = b_raw / b_norm
# Define block encoding
if be_method == "pauli":
data_size, block_size, be_scaling_factor, BlockEncodedState, be_qfunc = (
get_pauli_be(mat_raw_scr)
)
print(
f"Pauli block encoding with block size {block_size} and scaling factor {be_scaling_factor}"
)
elif be_method == "banded":
data_size, block_size, be_scaling_factor, BlockEncodedState, be_qfunc = (
get_banded_diags_be(mat_raw_scr)
)
print(
f"Banded diagonal block encoding with block size {block_size} and scaling factor {be_scaling_factor}"
)
# Get SVD range
w_min, w_max = get_svd_range(mat_raw_scr / be_scaling_factor)
# Get Chebyshev polynomial and the corresponding QSVT angles
pcoefs = get_cheb_coeff(
w_min, poly_degree, w_max, scale=SCALE, method=cheb_approx_type, epsilon=0.01
)
inv_phases = get_qsvt_phases(pcoefs)
# Define QSVT projector
@qfunc
def projector(be: BlockEncodedState, res: QBit):
res ^= be.block == 0
@qfunc
def main(
qsvt_aux: Output[QBit],
data: Output[QNum[data_size]],
block: Output[QNum[block_size]],
):
allocate(qsvt_aux)
allocate(block)
prepare_amplitudes(b_normalized.tolist(), 0, data)
be_state = BlockEncodedState()
within_apply(
lambda: bind([data, block], be_state),
lambda: qsvt_inversion(inv_phases, projector, be_qfunc, be_state, qsvt_aux),
)
if qmod_name is not None:
write_qmod(main, qmod_name, symbolic_only=False)
start_time_syn = time.time()
qprog = synthesize(main, preferences=preferences, constraints=constraints)
print("time to syn:", time.time() - start_time_syn)
execution_preferences = ExecutionPreferences(
num_shots=1,
backend_preferences=ClassiqBackendPreferences(
backend_name=ClassiqSimulatorBackendNames.SIMULATOR_STATEVECTOR
),
)
start_time_exe = time.time()
with ExecutionSession(qprog, execution_preferences) as es:
es.set_measured_state_filter("block", lambda state: state == 0.0)
es.set_measured_state_filter("qsvt_aux", lambda state: state == 0.0)
results = es.sample()
print("time to exe:", time.time() - start_time_exe)
resulting_state = get_projected_state_vector(results, "data")
normalization_factor = (be_scaling_factor * SCALE * w_min) / b_norm
return resulting_state / normalization_factor, qprog
# For faster results use the commented out preferences
# prefs = Preferences(
# transpilation_option="none", optimization_level=0, debug_mode=False, qasm3=True
# )
prefs = Preferences()
poly_degree = 101
We run two examples, for the two different block encoding, but using the same polynomial degree.
import pathlib
path = (
pathlib.Path(__file__).parent.resolve()
if "__file__" in locals()
else pathlib.Path(".")
)
mat_name = "nozzle_small_scr"
matfile = "matrices/" + mat_name + ".npz"
mat_raw_scr = sparse.load_npz(path / matfile)
b_raw = np.load(path / "matrices/b_nozzle_small.npy")
raw_mat_qubits = len(b_raw).bit_length() - 1 # matrix raw size
print(f"Raw matrix size: {raw_mat_qubits} qubits.")
Raw matrix size: 3 qubits.
qsol_banded, qprog_banded = qsvt_solver(
mat_raw_scr,
b_raw,
poly_degree,
be_method="banded",
preferences=prefs,
constraints=Constraints(optimization_parameter="width"),
qmod_name="qsvt_solver_banded_be",
)
show(qprog_banded)
Banded diagonal block encoding with block size 3 and scaling factor 3.8927056451476174
For error 0.01, and given kappa, the needed polynomial degree is: 2*295+1
The Chebyshev LCU requires 9 qubits
Performing convex optimization for the Chebyshev interpolation, with degree 101
Max relative error value: 0.017025127162188403
time to syn: 146.27435111999512
time to exe: 33.56966781616211
Quantum program link: https://platform.classiq.io/circuit/32BPbtdbOD4uIQw80Xgel5Ejvjt
qsol_pauli, qprog_pauli = qsvt_solver(
mat_raw_scr,
b_raw,
poly_degree,
be_method="pauli",
preferences=prefs,
constraints=Constraints(optimization_parameter="width"),
qmod_name="qsvt_solver_pauli_be",
)
show(qprog_pauli)
Pauli block encoding with block size 5 and scaling factor 5.6428068318664355
For error 0.01, and given kappa, the needed polynomial degree is: 2*449+1
The Chebyshev LCU requires 9 qubits
Performing convex optimization for the Chebyshev interpolation, with degree 101
Max relative error value: 0.06140234945729006
time to syn: 180.07041907310486
time to exe: 36.84007692337036
Quantum program link: https://platform.classiq.io/circuit/32BQ5DoazyihYKRCXX8zQCo9htb
mat_raw = mat_raw_scr.toarray()
expected_sol = np.linalg.solve(mat_raw, b_raw)
plt.plot(expected_sol, "o", label="Classical")
ext_idx = np.argmax(np.abs(expected_sol))
correct_sign = np.sign(expected_sol[ext_idx]) / np.sign(qsol_pauli[ext_idx])
qsol_pauli *= correct_sign
plt.plot(qsol_pauli, ".", label=f"QSVT-inv; Pauli BE; degree {poly_degree}")
correct_sign = np.sign(expected_sol[ext_idx]) / np.sign(qsol_banded[ext_idx])
qsol_banded *= correct_sign
plt.plot(qsol_banded, ".", label=f"QSVT-inv; Banded BE; degree {poly_degree}")
plt.legend()
<matplotlib.legend.Legend at 0x312817c10>
We can see that the results is better for the banded diagonal block-encoding. This is because it has a better scaling factor, and thus a better approximation for a given polynomial degree.
assert np.linalg.norm(qsol_banded - expected_sol) < 0.2
assert np.linalg.norm(qsol_pauli - expected_sol) < 0.2