Hi, I have been following the two demos on implementing QSVT in Pennylane: QSVT in practice and How to implement QSVT on hardware using angles calculated using the `laurent` method in PyQSP. As I understand it (and as I will later demonstrate), this should produce a target polynomial that is purely real and thus should not require to layer in LCU with two controlled QSVT circuits to extract the real part as is the case in these two demos. It should only require extracting the upper-left block of the matrix representation of the bare QSVT circuit. However, when I compute this block as in the code example below, I cannot get it to match the exact polynomial transformed operator. If someone could point me to where I am going wrong, it would be greatly appreciated. I am not super familiar with Pennylane, so I may be misunderstanding some code conventions such as the wiring order.
from pyqsp.angle_sequence import QuantumSignalProcessingPhases
from pyqsp import angle_sequence, response
from pyqsp.poly import (polynomial_generators, PolyTaylorSeries, PolyOneOverX, PolyCosineTX, PolySineTX)
import numpy as np
from scipy.linalg import sinm, cosm
import pennylane as qp
import matplotlib.pyplot as plt
#### Define polynomial, calculate phase factors, and plot polynomial as sanity check #######
poly = PolyCosineTX().generate(tau=3, epsilon=10**-4, ensure_bounded=False)
target_func = lambda x: np.cos(3*x)
ang_seq = angle_sequence.QuantumSignalProcessingPhases(
poly,
method='laurent',
chebyshev_basis=False,
signal_operator="Wx")
# Convert for QSVT
angles_pyqsp = qp.transform_angles(ang_seq, "QSP", "QSVT")
# Verify that imaginary part of polynomial is actually zero.
response.PlotQSPResponse(
ang_seq,
pcoefs=poly,
target=target_func,
sym_qsp=False,
simul_error_plot=False)
d = len(ang_seq) - 1
print(f'Number of phase factors: {len(ang_seq)}')
print(f'QSP degree d = {len(ang_seq) - 1}')
print(f'QSP angles: {ang_seq}')
print(f'QSVT angles: {angles_pyqsp}')
reflection_angles = []
reflection_angles.append(ang_seq[0] + (2*d - 1)*np.pi/4)
for _ in range(1,d):
reflection_angles.append(ang_seq[_] - np.pi/2)
reflection_angles.append(ang_seq[d] - np.pi/4)
print(f'Reflection angles calculated explicitly from Wx convention QSP angles:')
print(reflection_angles)
######## Explicitly verify that the <+| * |+> block of the QSP matrix is correct ##########
def Wx(x):
s = np.sqrt(1 - x**2)
return np.array([
[x, 1j * s],
[1j * s, x]
], dtype=complex)
def Rz(phi):
return np.array([
[np.exp(1j * phi), 0],
[0, np.exp(-1j * phi)]
], dtype=complex)
def qsp_matrix(x, phases):
U = Rz(phases[0])
Had = (1/np.sqrt(2)) * np.array([[1,1],
[1,-1]])
for phi in phases[1:]:
#U = U @ Wx(x) @ Rz(phi)
U = Rz(phi) @ Wx(x) @ U
U = Had @ U @ Had
return U
real_vals = []
imag_vals = []
points = np.linspace(start=-1, stop=1, num=200)
for point in points:
mat_elem = qsp_matrix(point, ang_seq)[0,0]
real_vals.append(np.real(mat_elem))
imag_vals.append(np.imag(mat_elem))
plt.plot(points, real_vals, label='Re[<+|U_qsp|+>]')
plt.plot(points, imag_vals, label='Im[<+|U_qsp|+>]')
plt.legend()
plt.show()
########## Generate Hamiltonian to be block-encoded ##############
coeffs = np.array([0.2, -0.7, -0.6])
coeffs /= np.linalg.norm(coeffs, ord=1) # Normalize the coefficients
obs = [qp.X(2), qp.X(2) @ qp.Z(3), qp.Z(2) @ qp.Y(3)]
H = qp.dot(coeffs, obs)
H_mat = qp.matrix(H, wire_order=[2, 3])
############ Generate circuit and check output with exact Cos(3*H) #########
control_wires = [0,1]
block_encode = qp.PrepSelPrep(H, control=control_wires)
projectors = [
qp.PCPhase(angles_pyqsp[i], dim=2 ** len(H.wires), wires=control_wires + H.wires)
for i in range(len(angles_pyqsp))
]
#projectors = [
# qp.PCPhase(reflection_angles[i], dim=2 ** len(H.wires), wires=control_wires + H.wires)
# for i in range(len(reflection_angles))
#]
dev = qp.device("default.qubit")
@qp.qnode(dev)
def circuit():
qp.QSVT(block_encode, projectors)
return qp.state()
matrix = qp.matrix(circuit, wire_order=control_wires + H.wires)()
sub_matrix = np.round(matrix[: 2 ** len(H.wires), : 2 ** len(H.wires)], 4)
# exact cos(3*H)
print(f'Cos(3*H):')
print(np.round(cosm(3*H_mat),4))
print()
# block-encoded P(H)
print('Block-encoded P(H):')
print(sub_matrix)
This produces the output:
Cos(3*H):
[[ 0.3483+0.j 0. +0.j 0. +0.j 0.8915+0.j]
[ 0. +0.j -0.246 +0.j -0.8915+0.j 0. +0.j]
[ 0. +0.j -0.8915+0.j 0.3483+0.j 0. +0.j]
[ 0.8915+0.j 0. +0.j 0. +0.j -0.246 +0.j]]
Block-encoded P(H):
[[ 0.3483-0.1478j 0. -0.j 0. +0.j 0.8914+0.17j ]
[-0. +0.j -0.246 -0.2611j -0.8914-0.17j -0. -0.j ]
[-0. +0.j -0.8914-0.17j 0.3483-0.1478j -0. +0.j ]
[ 0.8914+0.17j -0. +0.j 0. -0.j -0.246 -0.2611j]]
Output of qml.about() :
Name: pennylane
Version: 0.45.0
Summary: PennyLane is a cross-platform Python library for quantum computing, quantum machine learning, and quantum chemistry. Train a quantum computer the same way as a neural network.
Home-page:
Author:
License:
Location: /home/joel/miniconda3/envs/qsp-pennylane/lib/python3.12/site-packages
Platform info: Linux-6.14.5-300.fc42.x86_64-x86_64-with-glibc2.41
Python version: 3.12.0
Numpy version: 2.4.6
Scipy version: 1.17.1
JAX version: None
Catalyst version: None
Installed devices:
- lightning.qubit (pennylane_lightning-0.45.0)
- default.clifford (pennylane-0.45.0)
- default.gaussian (pennylane-0.45.0)
- default.mixed (pennylane-0.45.0)
- default.qubit (pennylane-0.45.0)
- default.qutrit (pennylane-0.45.0)
- default.qutrit.mixed (pennylane-0.45.0)
- default.tensor (pennylane-0.45.0)
- null.qubit (pennylane-0.45.0)
- reference.qubit (pennylane-0.45.0)
Any insight or help would be greatly appreciated.
EDIT: It does looks like ang_seqis symmetric (mod pi) so maybe PyQSP is doing something under the hood that I do not understand. In any case, the above code shows that the imaginary part of the target polynomial is zero, so I’m not quite sure this explains what is going wrong.