Quantum Phase Estimation for Solving Matrix Eigenvalues¶
Quantum Phase Estimation (QPE) is a key algorithm in quantum computing, allowing one to estimate the phase (or eigenvalue) relating to a Hermitian matrix. The algorithm is designed so, that given the inputs of a matrix $M$, and an eigenvalue ${|\psi\rangle}$, the output that will be obtained is $\theta$ where
$ U{|\psi\rangle} = e^{2\pi i\theta}{|\psi\rangle} , U = e^{2\pi iM} $.
By measuring the phase accumulated, the QPE algorithm allows calculating the eigenvalues relating to the chosen inputted vector. To read more about the QPE algorithm and the method it achieves the phase, we recommend [1].
In the general case, where the eigenvectors of the matrix are not known in advance yet the eigenvalues are sought, one can chose a random vector ${|v\rangle}$ for the algorithm’s initial state. For that case, some eigenvalues will be found as the vector can be described in the matrix's basis. Let us define {$\psi_i$} to be a set of eigenvalues of $M$. Generally, any vector can be written as a superposition of any basis set, thus -
${|v\rangle} = \sum_i a_i{|\psi_i\rangle}$
and -
$U{|v\rangle} = \sum_i a_i e^{2\pi i\theta_i}{|\psi_i\rangle}$.
Using execution with enough shots, one will obtain this set of $\theta_i$, i.e., a subset of the matrix's eigenvalues.
During this tutorial, we will follow the QPE algorithms steps as following:¶
a. Define a matrix.
b. Chose either it's eigenstate or a random vector.
c. Chose a resolution for our solution.
d. Find the related eigenvectors by using QPE and analyze the results.
0. Pre-requirments¶
The model is using several Classiq's libraries in addition to IBM's simulating tool.
import matplotlib.pyplot as plt
import numpy as np
import scipy
from numpy import linalg as LA
# for chemistry
from classiq import Model
from classiq.applications.chemistry import Molecule, MoleculeProblem
# for exponentiation
from classiq.builtin_functions import (
UCC,
Exponentiation,
HardwareEfficientAnsatz,
HartreeFock,
PhaseEstimation,
StatePreparation,
UnitaryGate,
)
from classiq.builtin_functions.exponentiation import (
ExponentiationConstraints,
ExponentiationOptimization,
PauliOperator,
)
from classiq.builtin_functions.qpe import (
ExponentiationScaling,
ExponentiationSpecification,
)
from classiq.execution import IBMBackendPreferences
# for state preperation and phase estimation
from classiq.model import Constraints, Preferences
1. Chose the QPE Parameters¶
1.1. Set the Matrix¶
Here, we will define the matrix we would like to submit. This can be any harmitian matrix with size $2^n$ by $2^n$ with $n$ a positive integer. We will use this matrix throughout our code under the name M.
# M = np.array([[0.25, 0], [0, 0.5]])
# M = np.array([[0.3, 0], [0, 0.35]])
# M = np.array([[0,0,0,0], [0, 3,0,0],[0, 0,0.5,0],[0, 0,0,0.75]])
# M = np.array([[0,0,0,0], [0, 3,0,0],[0, 0,-0.5,0],[0, 0,0,-0.75]])
M = np.array([[0, 3, 4, 0], [-0.8, 3, 0, 0], [1, 0, -0.5, 5], [0, 0, 0, -0.75]])
M = np.array(
[
[0.38891555, 0.23315811, 0.21499372, 0.06119557],
[0.23315811, 0.44435328, 0.25197881, -0.13087919],
[0.21499372, 0.25197881, 0.44116509, -0.01961855],
[0.06119557, -0.13087919, -0.01961855, 0.32556608],
]
)
M_t = M.transpose()
M = (M + M_t) / 2
1.2. Set the Vector¶
Here we will choose the vector which will be defined later as the initial condition for the run. There are two options - you can define any vector by using the parameter eigen_vec to be False, and changing the value int_vec.
If you take eigen_vec to be True, then you can choose a number of eigenvalue that will be set as the initial state under ev.
eigen_vec = False
if eigen_vec:
w, v = LA.eig(M)
print("the eigenvalues are", w)
print("the eigenvectors are", v, sep="\n")
ev = 1
int_vec = v[:, ev]
else:
int_vec = np.random.rand(np.shape(M)[0])
print("Your initial state is", int_vec)
2. Preparing the Matrix for QPE¶
2.1 Translating the Matrix into Pauli_ops¶
In order to translate the matrix into quantum circuit language, the matrix should be written in the form of list of strings, of composite pauli operators.
###this is the code to take a matrix and tranform it into pauli string
import itertools # noqa
from itertools import product
from numpy import kron
Paulidict = {
"I": np.array([[1, 0], [0, 1]], dtype=np.complex128),
"Z": np.array([[1, 0], [0, -1]], dtype=np.complex128),
"X": np.array([[0, 1], [1, 0]], dtype=np.complex128),
"Y": np.array([[0, -1j], [1j, 0]], dtype=np.complex128),
}
# generate all combinations of Pauli strings of size n
def generateAllPauliStrings(seq, n):
for s in product(seq, repeat=n):
yield "".join(s)
# convert a Paulistring of size n to 2**n X 2**n matrix
def PauliString2mat(seq):
myPmat = Paulidict[seq[0]]
for p in seq[1:]:
myPmat = kron(myPmat, Paulidict[p])
return myPmat
# Hilbert-Schmidt-Product of two matrices M1, M2
def HilbertSchmidt(M1, M2):
return (np.dot(M1.conjugate().transpose(), M2)).trace()
# Naive decomposition, running over all HS products for all Pauli strings
def LCU_naive(H):
assert H.shape[0] == H.shape[1], "matrix is not square"
assert H.shape[0] != 0, "matrix is of size 0"
assert H.shape[0] & (H.shape[0] - 1) == 0, "matrix size is not 2**n"
n = int(np.log2(H.shape[0]))
myPualiList = list(generateAllPauliStrings("IZXY", n))
mylist = []
for pstr in myPualiList:
co = (1 / 2**n) * HilbertSchmidt(PauliString2mat(pstr), H)
if co != 0:
mylist = mylist + [(pstr, co)]
return mylist
pauli_ops = LCU_naive(M)
print(pauli_ops)
pauli_operator = PauliOperator(pauli_list=pauli_ops)
N = pauli_operator.num_qubits
print(pauli_operator.pauli_list)
2.2 Chose the Algorithm's Precision¶
For QPE algorithms, the precision is set by the number of qubits chosen $n$, such that the resolution is $1/{2^n}$. In case the matrix needs to be normalized, the resolution will be distorted. In the case of normalization, the span of results for the QPE is stretched between the lowest and highest possible phase, thus the resolution will be mapped to $normalization-coefficient/{2^n} ~\sim 1/{((\lambda_{max}-\lambda_{min})*2^n)}$.
Here you can choose the precision in the parameter n_qpe, or by setting your desired resolution. If you choose your desired resolution and set the parameter get_recommended_n to be True, the number of qubits will be calculated for you accordingly.
n_qpe = 8
# recommanded n_qpe:
get_recommended_n = False
import math
desired_resolution = 0.02
def get_nqpe(pauli_operator, desired_resolution):
N = pauli_operator.num_qubits
A = 0
for a, b in pauli_operator.pauli_list:
A = A + abs(b)
nqpe = math.log2(2 * N * A / desired_resolution)
return math.ceil(nqpe)
if get_recommended_n:
n_qpe = get_nqpe(pauli_operator, desired_resolution)
print("number of qubits for QPE is", n_qpe)
2.3 Normalize the Matrix¶
As QPE obtains a phase in the form $e^{2\pi i\theta}$, there is meaning only for $\theta \in [0,2\pi)$. Generally, our matrix M can have any eigenvalue, thus $\theta$ can have any value. In order to fix this discrepancy, the values of the matrix are stretched to be rescaled. We assume $\theta \in [\lambda_{min}, \lambda_{max}]$ and use a normalization function in order to map those values into $[0, 1-1/{2^n}]$, where $n$ is the number of qubits chosen for the QPE process in section 2.2.
We perform the normalization procedure as following:
a. We evaluate $\lambda_{min},\lambda_{max}$ (in the function normalization_params() below). In order to do so we use rough estimation of the absolute max value that can take place by adding together all the pauli coefficients and multiplying by the matrix's dimensions. That will yield us a value $\lambda$ (which is referred in the code as normalization_coeff) and we now assume that the domain is $\theta \in [-\lambda, \lambda]$.
In general, one can build a more accurate assessment, which will decrease the span of solutions and thus achieve a better resolution.
b. We make sure only positive values are available by adding $\lambda*I^n$ to the pauli list. Now our evaluated span is $[0, 2*\lambda]$.
c. We normalize our matrix by multiplying all of the pauli coefficients by $(1-1/2^n)/(2*\lambda)$. Now the span of $\theta$ is $[0, 1-1/2^n]$, as required for proper QPE process.
The values of the matrix's eigenvalues should be now between $0$ to $1-(1/2^n)$. The QPE procedure will be performed on this new normalized matrix. After the phases are obtained, the original phases of the pre-normalized matrix will be gathered by performing opposite steps to this normalization procedure.
- Note that in case your matrix's eigenvalues are naturally between the values $0$ to $1-(1/2^n)$, you may not want to normalize it, as the normalization procedure may enlarge the span, thus lowering the resolution of the algorithm. In that case, you may skip those lines or change the value
normalizeto False.
pauli_operator = PauliOperator(pauli_list=pauli_ops)
N = pauli_operator.num_qubits
print(pauli_operator.pauli_list)
# normalizing the operator
## we need to create a matrix such that its normalized version will have eigenvalues of [0,1/2^k] when k is the resolution of the QPE
normalize = True
def normalization_params(pauli_operator, N):
A = 0
for a, b in pauli_operator.pauli_list:
A = A + abs(b)
return N * A
def normalize_hamiltonian(pauli_operator, normalization_coeff, k):
new_pauli_operator = []
for a, b in pauli_operator.pauli_list:
if a == "I" * N:
new_pauli_operator.append(
(
a,
(b + normalization_coeff)
* (1 - 1 / (2**k))
/ (2 * normalization_coeff),
)
)
else:
new_pauli_operator.append(
(a, b * (1 - 1 / (2**k)) / (2 * normalization_coeff))
)
return new_pauli_operator
if normalize:
normalization_coeff = normalization_params(pauli_operator, N)
new_pauli_operator = normalize_hamiltonian(
pauli_operator, normalization_coeff, n_qpe
)
pauli_ops = new_pauli_operator
print(pauli_ops)
size = math.sqrt(M.size)
I = np.eye(int(size))
Mnew = (
(M + normalization_coeff * I)
* (1 - 1 / (2**n_qpe))
/ (2 * normalization_coeff)
)
else:
Mnew = M
3. Creating the Quantum Circuit¶
We will now create a quantum circuit of the QPE algorithm using the Classiq platform. The user is able to fill in their constraints and perefences as desired.
constraints = Constraints()
preferences = Preferences()
model = Model(preferences=preferences, constraints=constraints)
3.1. Initializing the State¶
StatePreperation function is used to initialize that state, where the amplitudes are determined by the state set by the user in section 1.2.
my_amp = tuple(
int_vec / np.linalg.norm(int_vec)
) # amplitude is given by the eignevector
sp_upper = 0.00 # precision of the State Preparation
sp_params = StatePreparation(
amplitudes=my_amp, error_metric={"L2": {"upper_bound": sp_upper}}
)
sp_out = model.StatePreparation(params=sp_params)
3.2. Calling the Phase Estimation¶
There are generally 2 methods for inseting the matrix into the qft - unitary implementation is exact but long, and exponantiation which is approximated but shorter in depth. We will chose the parameter is exact to indicate what method chosen.
is_exact = False
exp_max_depth = 8000
if is_exact:
my_unitary = scipy.linalg.expm(1j * 2 * np.pi * Mnew)
exp_params = UnitaryGate(data=my_unitary.tolist())
Qreg_name = {"IN": "TARGET", "OUT": "TARGET"}
else:
po = pauli_ops
exp_params = Exponentiation(
pauli_operator=PauliOperator(pauli_list=po),
evolution_coefficient=-2 * np.pi,
# optimization=ExponentiationOptimization.MINIMIZE_ERROR,
# constraints=ExponentiationConstraints(max_depth=exp_max_depth)
)
Qreg_name = {"IN": "IN", "OUT": "OUT"}
if is_exact:
qpe_params = PhaseEstimation(size=n_qpe, unitary_params=exp_params)
else:
qpe_params = PhaseEstimation(
size=n_qpe,
unitary_params=exp_params,
exponentiation_specification=ExponentiationSpecification(
scaling=ExponentiationScaling(max_depth=200, max_depth_scaling_factor=1.8)
),
)
qpe_out = model.PhaseEstimation(
params=qpe_params, in_wires={Qreg_name["IN"]: sp_out["OUT"]}
)
3.3. Defining the Registers¶
The following step is made in order to have a map of the relevant outputted qubits for analysis.
model.set_outputs({"phase_result": qpe_out["PHASE_ESTIMATION"]})
3.4. Synthesizing the Circuit¶
Hereby we syntesize the circuit and show it using the analyzer.
from classiq.execution import ExecutionPreferences
num_shots = 10000
model.execution_preferences = ExecutionPreferences(num_shots=num_shots)
model.sample()
qmod = model.get_model()
with open("qpe_for_matrix.qmod", "w") as f:
f.write(qmod)
from classiq import synthesize
qprog = synthesize(qmod)
from classiq import GeneratedCircuit
circuit = GeneratedCircuit.from_qprog(qprog)
circuit.show()
4. Measurment and GeneratedCircuitalysis¶
4.1. Circuit Execution¶
The circuit is now sent to execution by a chosen backend.
from classiq import execute
results = execute(qprog).result()
from classiq.execution import ExecutionDetails
results = results[0].value
4.2. Presenting the Result's Count¶
4.3. Chosing the Most Probable Solution¶
Hereby the user will choose the number of eigenvalues they wish to extract from the poll of results. The value number_of_solutions will determine how many results out of qpe_results will be analyzed. The algorithm solution_search will save only the number_of_solutions highest count results from the execution.
number_of_solutions = 2 # number of phases sought
solution = [
sampled_state.state["phase_result"]
for sampled_state in results.parsed_counts[:number_of_solutions]
]
print("Your qubit solution is", solution, sep="\n")
4.4. Translating into Eigenvalue (Phase)¶
Here, the value in the vector results are translated from a binary number into a full solution for the eigenvalues.
At first - the binary number is translated into a decimal value.
dec_sol_vec = [sol / (2**n_qpe) for sol in solution]
print("Your decimal solutions are", dec_sol_vec, sep="\n")
Secondy, the decimal value is mapped back into the original values, i.e. renormalized into it's original span.
# renormalize into the "real" solution -
if normalize:
solution = [
(
(value * 2 * normalization_coeff / (1 - (1 / 2**n_qpe)))
- normalization_coeff
)
for value in dec_sol_vec
]
else:
solution = dec_sol_vec
The final results of the phases (matrix's eigenvlues) is therefore:
print(solution)
And the results including the error contributed from the resolution (number of qubits participating in the QPE) are:
if normalize:
energy_resolution = (
(1 / (2**n_qpe)) * 2 * normalization_coeff / (1 - (1 / 2**n_qpe))
)
else:
energy_resolution = 1 / (2**n_qpe)
print("the resolution of results is", energy_resolution)
for sol in solution:
print(
"the solutions are between",
sol - energy_resolution,
"and",
sol + energy_resolution,
)
###if zero or exceed the normalization range need to add conditions
4.5. Compare to Exact Results:¶
w, v = LA.eig(M)
print("the eigenvalues are", w)
print("the eigenvectors are", v, sep="\n")
4.6. Find the Solution's Histogram¶
import matplotlib.pyplot as plt
import numpy as np
energy_vec = []
energy_prob = []
for sampled_state in results.parsed_counts:
temp = sampled_state.state["phase_result"] / (2**n_qpe)
if normalize:
temp2 = (
temp * 2 * normalization_coeff / (1 - (1 / 2**n_qpe))
) - normalization_coeff
else:
temp2 = temp
energy_vec.append(temp2)
energy_prob.append(sampled_state.shots / num_shots)
plt.plot(energy_vec, energy_prob, ".")
plt.show()
5. Appendices¶
This code is converting pauli list of operators to a matrix array. It can be useful for debugging purposes.
# for debug
# make back pauli_ops into M
import numpy as np
from numpy import kron
Paulidict = {
"I": np.array([[1, 0], [0, 1]], dtype=np.complex128),
"Z": np.array([[1, 0], [0, -1]], dtype=np.complex128),
"X": np.array([[0, 1], [1, 0]], dtype=np.complex128),
"Y": np.array([[0, -1j], [1j, 0]], dtype=np.complex128),
}
def PauliString2mat(seq):
myPmat = Paulidict[seq[0]]
for p in seq[1:]:
myPmat = kron(myPmat, Paulidict[p])
return myPmat
pauli_example = pauli_ops
mat_new = 0
for s in pauli_example:
mat_new = mat_new + s[1] * PauliString2mat(s[0])
print(mat_new)
w, v = LA.eig(mat_new)
print("the eigenvalues are", w)
print("the eigenvectors are", v, sep="\n")